added all relevant files

LICENSE

@@ -0,0 +1,165 @@
+GNU LESSER GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+
+  This version of the GNU Lesser General Public License incorporates
+the terms and conditions of version 3 of the GNU General Public
+License, supplemented by the additional permissions listed below.
+
+  0. Additional Definitions.
+
+  As used herein, "this License" refers to version 3 of the GNU Lesser
+General Public License, and the "GNU GPL" refers to version 3 of the GNU
+General Public License.
+
+  "The Library" refers to a covered work governed by this License,
+other than an Application or a Combined Work as defined below.
+
+  An "Application" is any work that makes use of an interface provided
+by the Library, but which is not otherwise based on the Library.
+Defining a subclass of a class defined by the Library is deemed a mode
+of using an interface provided by the Library.
+
+  A "Combined Work" is a work produced by combining or linking an
+Application with the Library.  The particular version of the Library
+with which the Combined Work was made is also called the "Linked
+Version".
+
+  The "Minimal Corresponding Source" for a Combined Work means the
+Corresponding Source for the Combined Work, excluding any source code
+for portions of the Combined Work that, considered in isolation, are
+based on the Application, and not on the Linked Version.
+
+  The "Corresponding Application Code" for a Combined Work means the
+object code and/or source code for the Application, including any data
+and utility programs needed for reproducing the Combined Work from the
+Application, but excluding the System Libraries of the Combined Work.
+
+  1. Exception to Section 3 of the GNU GPL.
+
+  You may convey a covered work under sections 3 and 4 of this License
+without being bound by section 3 of the GNU GPL.
+
+  2. Conveying Modified Versions.
+
+  If you modify a copy of the Library, and, in your modifications, a
+facility refers to a function or data to be supplied by an Application
+that uses the facility (other than as an argument passed when the
+facility is invoked), then you may convey a copy of the modified
+version:
+
+   a) under this License, provided that you make a good faith effort to
+   ensure that, in the event an Application does not supply the
+   function or data, the facility still operates, and performs
+   whatever part of its purpose remains meaningful, or
+
+   b) under the GNU GPL, with none of the additional permissions of
+   this License applicable to that copy.
+
+  3. Object Code Incorporating Material from Library Header Files.
+
+  The object code form of an Application may incorporate material from
+a header file that is part of the Library.  You may convey such object
+code under terms of your choice, provided that, if the incorporated
+material is not limited to numerical parameters, data structure
+layouts and accessors, or small macros, inline functions and templates
+(ten or fewer lines in length), you do both of the following:
+
+   a) Give prominent notice with each copy of the object code that the
+   Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the object code with a copy of the GNU GPL and this license
+   document.
+
+  4. Combined Works.
+
+  You may convey a Combined Work under terms of your choice that,
+taken together, effectively do not restrict modification of the
+portions of the Library contained in the Combined Work and reverse
+engineering for debugging such modifications, if you also do each of
+the following:
+
+   a) Give prominent notice with each copy of the Combined Work that
+   the Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the Combined Work with a copy of the GNU GPL and this license
+   document.
+
+   c) For a Combined Work that displays copyright notices during
+   execution, include the copyright notice for the Library among
+   these notices, as well as a reference directing the user to the
+   copies of the GNU GPL and this license document.
+
+   d) Do one of the following:
+
+       0) Convey the Minimal Corresponding Source under the terms of this
+       License, and the Corresponding Application Code in a form
+       suitable for, and under terms that permit, the user to
+       recombine or relink the Application with a modified version of
+       the Linked Version to produce a modified Combined Work, in the
+       manner specified by section 6 of the GNU GPL for conveying
+       Corresponding Source.
+
+       1) Use a suitable shared library mechanism for linking with the
+       Library.  A suitable mechanism is one that (a) uses at run time
+       a copy of the Library already present on the user's computer
+       system, and (b) will operate properly with a modified version
+       of the Library that is interface-compatible with the Linked
+       Version.
+
+   e) Provide Installation Information, but only if you would otherwise
+   be required to provide such information under section 6 of the
+   GNU GPL, and only to the extent that such information is
+   necessary to install and execute a modified version of the
+   Combined Work produced by recombining or relinking the
+   Application with a modified version of the Linked Version. (If
+   you use option 4d0, the Installation Information must accompany
+   the Minimal Corresponding Source and Corresponding Application
+   Code. If you use option 4d1, you must provide the Installation
+   Information in the manner specified by section 6 of the GNU GPL
+   for conveying Corresponding Source.)
+
+  5. Combined Libraries.
+
+  You may place library facilities that are a work based on the
+Library side by side in a single library together with other library
+facilities that are not Applications and are not covered by this
+License, and convey such a combined library under terms of your
+choice, if you do both of the following:
+
+   a) Accompany the combined library with a copy of the same work based
+   on the Library, uncombined with any other library facilities,
+   conveyed under the terms of this License.
+
+   b) Give prominent notice with the combined library that part of it
+   is a work based on the Library, and explaining where to find the
+   accompanying uncombined form of the same work.
+
+  6. Revised Versions of the GNU Lesser General Public License.
+
+  The Free Software Foundation may publish revised and/or new versions
+of the GNU Lesser General Public License from time to time. Such new
+versions will be similar in spirit to the present version, but may
+differ in detail to address new problems or concerns.
+
+  Each version is given a distinguishing version number. If the
+Library as you received it specifies that a certain numbered version
+of the GNU Lesser General Public License "or any later version"
+applies to it, you have the option of following the terms and
+conditions either of that published version or of any later version
+published by the Free Software Foundation. If the Library as you
+received it does not specify a version number of the GNU Lesser
+General Public License, you may choose any version of the GNU Lesser
+General Public License ever published by the Free Software Foundation.
+
+  If the Library as you received it specifies that a proxy can decide
+whether future versions of the GNU Lesser General Public License shall
+apply, that proxy's public statement of acceptance of any version is
+permanent authorization for you to choose that version for the
+Library.

README

@@ -0,0 +1,29 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+                             COPYRIGHT
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+This package contains Matlab source code for object detection with LDA models.
+
+This repo is based on the original version of who from 
+(Bharath Hariharan, Jitendra Malik and Deva Ramanan. 
+Discriminative decorrelation for clustering and classification. In ECCV 2012 ).
+We significantly altered it, making it more flexible, easier to access, 
+and less focused on hog features in general.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+                           START / SETUP
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Compile mex-files via
+compileWHOEfficient
+
+Setup workspace
+initWorkspaceWHOGeneric
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+                           RUN DETECTION DEMO
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Run a demo which learns an LDA model for detecting bycicles and 
+using HOG features as underlying representation
+demoWhitenedHOG

compileWHOEfficient.m

@@ -0,0 +1,18 @@
+% Compile mex functions
+%
+% feature extraction
+mex -O features/featuresHOGGrayScale.cc -o features/featuresHOGGrayScale
+mex -O features/featuresHOGColor.cc     -o features/featuresHOGColor
+mex -O features/featuresHOGorig.cc      -o features/featuresHOGorig
+%
+% image resizing ( by scalar factors )
+mex -O features/resizeGrayScale.cc      -o features/resizeGrayScale
+mex -O features/resizeColor.cc          -o features/resizeColor
+%
+% image shrinking ( always by factor of 2 )
+mex -O features/reduceGrayScale.cc      -o features/reduceGrayScale
+mex -O features/reduceColor.cc          -o features/reduceColor
+%
+% discrete convolutions
+mex -O features/fconv2D.cc              -o features/fconv2D
+mex -O features/fconv3D.cc              -o features/fconv3D

demos/demoWhitenedHOG.m

@@ -0,0 +1,77 @@
+% author: Alexander Freytag
+% date  : 27-02-2014 (dd-mm-yyyy)
+
+
+%% Training
+% We will train a model from a single instance.
+% So let's create a tiny list of positive examples
+pos(1).im = 'train.jpg';
+pos(1).x1 = 70;
+pos(1).y1 = 202;
+pos(1).x2 = 255;
+pos(1).y2 = 500;
+
+
+
+%read image ... 
+im       = readImage(pos(1).im);
+figTrain = figure;
+set ( figTrain, 'name', 'Training Image');
+% ... and show initial bounding box
+showboxes(im ,[pos(1).x1 pos(1).y1 pos(1).x2 pos(1).y2]);
+
+%settings for feature extraction
+settings.i_binSize = 8;
+settings.interval  = 10; % same as on who demo file
+settings.order     = 20;
+% note:
+% change the representation as you desire. This repo comes along with HOG
+% features as default, however, any diffferent feature type can be plugged
+% in as well given the proper wrapper funtion.
+% Examples can be found in our repository about patch discovery
+settings.fh_featureExtractor = ...
+  struct('name','Compute HOG features using WHO code', ...
+         'mfunction',@computeHOGs_WHOorig, ...
+         'b_leaveBoundary' , true );
+
+
+% try locate previously trained bg-struct
+try
+    fileToBG = fullfile(pwd, 'data/bg11.mat');
+    load( fileToBG );
+    % compatibility to older versions
+    if ( isfield(bg,'sbin') && ~isfield(bg, 'i_binSize') );
+        bg.i_binSize = bg.sbin;
+    end
+catch
+    % if not possible, leave bg empty and compute it from scratch lateron
+    bg=[];
+end
+
+% no negative examples, use 'universal' negative model
+neg   = [];
+model = learn_dataset( pos, neg, bg, settings );
+
+%show learned model, i.e., visualized HOG feature for positive and negative weights
+b_closeImg = false;
+showWeightVectorHOG( model.w, b_closeImg )
+
+%% Testing
+test(1).im='test.jpg';
+
+%perform detection
+boxes=test_dataset( test, model, settings );
+
+%convert output from 1x1-cell to double array
+boxes = boxes{1};
+%only take highest response
+%[~,bestIdx] = max ( boxes(:,5) );
+%boxes = boxes(bestIdx,:);
+
+%show detections
+im      = readImage(test(1).im);
+figTest = figure;
+set ( figTest, 'name', 'Test Image');
+showboxes(im, boxes);
+
+% that's it :)

detect/detectWithGivenFeatures.m

@@ -0,0 +1,66 @@
+function boxes = detectWithGivenFeatures(pyraFeats,model,thresh)
+%function boxes = detectWithGivenFeatures(pyraFeats,model,thresh)
+% 
+% BRIEF:
+%     Detect object in images by searching for large model responses in
+%     terms of convolution scores
+% 
+%     Pretty much the same as the original code of who, only differs in
+%     handing over precomputed features.
+%     Additionally, no hallucinated levels are deleted.
+% 
+% author: Alexander Freytag
+% date:   26-02-2014 (dd-mm-yyyy)
+
+    levels   = 1:length(pyraFeats.feat);
+    %pre-allocate some memory (assume not more then 10k boxes in advance)
+    boxes    = zeros(10000,5);
+    myDetectorFilter  = {model.w};
+
+    % padding of image borders
+    padx = pyraFeats.padx;
+    pady = pyraFeats.pady;
+    
+    sizx = size(model.w,2);
+    sizy = size(model.w,1);
+    cnt  = 0;
+
+    for l = levels,
+        % get features of current scale
+        scale = pyraFeats.scale(l);
+      
+        % do the actual convolution
+        if ( ndims ( myDetectorFilter{1} ) > 2 ) 
+            resp  = fconv3D( pyraFeats.feat{l}, myDetectorFilter, 1,1);
+        else
+            resp  = fconv2D( pyraFeats.feat{l}, myDetectorFilter, 1,1);
+        end        
+        
+        % correction of data structure
+        resp  = resp{1};
+        
+        % only accept scores over thresh
+        [y,x] = find(resp >= thresh);
+        
+        if ( ~isempty(x) )
+            I  = (x-1)*size(resp,1)+y;
+
+            % convert responses to original image coordinates
+            x1 = (x-1-padx)*scale + 1;
+            y1 = (y-1-pady)*scale + 1;
+            x2 = x1 + sizx*scale  - 1;
+            y2 = y1 + sizy*scale  - 1;
+
+            i = cnt+1:cnt+length(I);
+
+            % store corners of box + score
+            boxes(i,:) = [x1 y1 x2 y2 resp(I)];
+
+            % increase number of responses found so far
+            cnt = cnt+length(I);
+        end
+    end
+
+    %that's it, return the found boxes with corresponding scores
+    boxes = boxes(1:cnt,:);
+end

detect/test_dataset.m

@@ -0,0 +1,21 @@
+function boxes=test_dataset(test, model, settings )
+% boxes=test_dataset(test, model, setting )
+% test is struct array with fields:
+%	.im:full path to image
+
+%TODO
+
+for i = 1:length(test),
+  fprintf('testing: %d/%d\n', i, length(test));
+  im = readImage(test(i).im);
+  %tic;
+  pyraFeat = featPyramidGeneric( im, model, settings ) ;
+  boxes{i} = detectWithGivenFeatures( pyraFeat, model, model.thresh );
+  %toc; tic;
+  %%%
+  %TODO replace by own method
+  boxes{i} = nonMaxSupp(boxes{i},.25);
+  %toc;
+ % showboxes(im,boxes{i});
+end
+

features/computeFeaturesForBlocks.m

@@ -0,0 +1,27 @@
+function features = computeFeaturesForBlocks ( blocks, settings )
+%%TODO docu
+
+    %% (1) check input
+    if ( nargin < 2 )
+        settings = [];
+    end
+    
+    %check for feature extractor, if not existing, set to default
+    fh_featureExtractor = struct('name','Compute mean patches', 'mfunction',@computeMeanPatches);
+    settings = addDefaultVariableSetting( settings, 'fh_featureExtractor', fh_featureExtractor, settings );    
+    
+    %check for feature extractor, if not existing, set to []
+    settingsFeatureExtractorDefault = [];
+    settings = addDefaultVariableSetting( settings, 'settingsFeatureExtractor', settingsFeatureExtractorDefault, settings );
+    
+    %% (2) compute features
+    
+    n = length(blocks);
+       
+    for i=n:-1:1 
+        feature = settings.fh_featureExtractor.mfunction ( blocks{i}, settings.settingsFeatureExtractor );
+        features(i).feature = feature;
+    end
+
+
+end

features/computeHOGs_WHO.m

@@ -0,0 +1,29 @@
+function hogFeature = computeHOGs_WHO ( img, settings )
+%%TODO docu
+
+    %% (1) check input
+    if ( ( nargin >= 2 ) && ...
+         ( ~isempty (settings) ) && ...
+         ( isstruct ( settings ) ) && ...
+         ( isfield(settings, 'sbin')  )...
+       )
+        sbin = settings.sbin;
+    %just for backward compatibility we check both options
+    elseif ( ( nargin >= 2 ) && ...
+             ( ~isempty (settings) ) && ...
+             ( isstruct ( settings ) ) && ...
+             ( isfield(settings, 'i_binSize')  )...
+           )
+        sbin = settings.i_binSize;
+    else
+        sbin = 8;
+    end
+    
+    %% (2) compute features
+    if ( ndims(img) == 3 )
+        hogFeature = featuresHOGColor( double(img), sbin );
+    else
+        hogFeature = featuresHOGGrayScale( double(img), sbin );
+    end   
+    
+end

features/computeHOGs_WHOorig.m

@@ -0,0 +1,29 @@
+function hogFeature = computeHOGs_WHOorig ( img, settings )
+%%TODO docu
+
+    %% (1) check input
+    if ( ( nargin >= 2 ) && ...
+         ( ~isempty (settings) ) && ...
+         ( isstruct ( settings ) ) && ...
+         ( isfield(settings, 'sbin')  )...
+       )
+        sbin = settings.sbin;
+    %just for backward compatibility we check both options
+    elseif ( ( nargin >= 2 ) && ...
+             ( ~isempty (settings) ) && ...
+             ( isstruct ( settings ) ) && ...
+             ( isfield(settings, 'i_binSize')  )...
+           )
+        sbin = settings.i_binSize;
+    else
+        sbin = 8;
+    end
+    
+    %% (2) compute features
+    if ( ndims(img) == 3 )
+        hogFeature = featuresHOGorig( double(img), sbin );
+    else
+        hogFeature = featuresHOGorig( double(repmat(img,[1,1,3])), sbin );
+    end   
+    
+end

features/computeOptimalCellNumber.m

@@ -0,0 +1,44 @@
+function i_numCells = computeOptimalCellNumber ( blocks, i_binSize )
+% function i_numCells = computeOptimalCellNumber ( blocks, i_binSize )
+% 
+% author: Alexander Freytag
+% date:   02-03-2014 ( dd-mm-yyyy )
+% 
+
+    % pick mode of aspect ratios
+    h = [blocks(:).y2]' - [blocks(:).y1]' + 1;
+    w = [blocks(:).x2]' - [blocks(:).x1]' + 1;
+
+    xx = -2:.02:2;
+    filter = exp(-[-100:100].^2/400);
+    aspects = hist(log(double(h)./double(w)), xx);
+    aspects = convn(aspects, filter, 'same');
+    [peak, I] = max(aspects);
+    aspect = exp(xx(I));
+
+    % pick 20 percentile area
+    mean(h);
+    mean(w);
+    areas = sort(h.*w);
+    % nasty hack to constrain HoG feature support areas to be not too big or 
+    % too small, even if our input data would tell us otherwise
+    area = areas(max(floor(length(areas) * 0.2),1));
+    area = max(min(area, 7000), 5000);
+
+
+    % how many pixels shall a cell cover in x and y direction?
+    sbin = i_binSize;
+
+    % pick dimensions according to data suggestions
+    w  = sqrt(double(area)/aspect);
+    h  = w*aspect;
+
+
+    % resulting number of cells in x direction
+    sizeX = round(h/sbin);
+    % resulting number of cells in y direction
+    sizeY = round(w/sbin);
+    
+    i_numCells = [sizeX sizeY];
+    i_numCells = max(i_numCells,1);
+end

features/fconv2D.cc

@@ -0,0 +1,179 @@
+#include "mex.h"
+#include <math.h>
+#include <string.h>
+
+/*
+ * This code is used for computing filter responses.  It computes the
+ * response of a set of filters with a feature map.  
+ *
+ * Basic version, relatively slow but very compatible.
+ */
+
+struct thread_data {
+  double *A;
+  double *B;
+  double *C;
+  mxArray *mxC;
+  const mwSize *A_dims;
+  const mwSize *B_dims;
+  mwSize C_dims[2];
+};
+
+// convolve A and B
+void *process(void *thread_arg)
+{
+  thread_data *args = (thread_data *)thread_arg;
+  
+  //input A - e.g., an image
+  double *A = args->A;
+  
+  //input B - e.g., a filter
+  double *B = args->B;
+  
+  //output C - the convolution response map
+  double *C = args->C;
+  
+  const mwSize *A_dims = args->A_dims;
+  const mwSize *B_dims = args->B_dims;
+  const mwSize *C_dims = args->C_dims;
+
+  //new pointers just for comparability with original fconv.cc
+  double *dst = C;
+  double *A_src = A;
+  double *B_src = B;
+  
+  //run over possible locations on x dimension
+  for (int x = 0; x < C_dims[1]; x++)
+  {
+    //run over possible locations on y dimension
+    for (int y = 0; y < C_dims[0]; y++)
+    {
+      // compute convolution score for current position
+      // loop over filter size
+      double val = 0;
+      for (int xp = 0; xp < B_dims[1]; xp++)
+      {
+        double *A_off = A_src + (x+xp)*A_dims[0] + y;
+        double *B_off = B_src + xp*B_dims[0];
+        
+        //hope-to-be efficient version for up to 20 dimensions
+        switch(B_dims[0])
+        {
+          case 20: val += A_off[19] * B_off[19];
+          case 19: val += A_off[18] * B_off[18];
+          case 18: val += A_off[17] * B_off[17];
+          case 17: val += A_off[16] * B_off[16];
+          case 16: val += A_off[15] * B_off[15];
+          case 15: val += A_off[14] * B_off[14];
+          case 14: val += A_off[13] * B_off[13];
+          case 13: val += A_off[12] * B_off[12];
+          case 12: val += A_off[11] * B_off[11];
+          case 11: val += A_off[10] * B_off[10];
+          case 10: val += A_off[9]  * B_off[9];
+          case 9:  val += A_off[8]  * B_off[8];
+          case 8:  val += A_off[7]  * B_off[7];
+          case 7:  val += A_off[6]  * B_off[6];
+          case 6:  val += A_off[5]  * B_off[5];
+          case 5:  val += A_off[4]  * B_off[4];
+          case 4:  val += A_off[3]  * B_off[3];
+          case 3:  val += A_off[2]  * B_off[2];
+          case 2:  val += A_off[1]  * B_off[1];
+          case 1:  val += A_off[0]  * B_off[0];
+            break;
+          // less efficient version with more than 20 dimensions
+          default:              
+            for (int yp = 0; yp < B_dims[0]; yp++)
+            {
+              val += *(A_off++) * *(B_off++);
+            }
+        } // nasty switch-statement
+      } //for-loop over filtersize
+      
+      //write value to current position in output data and increment output pointer
+      *(dst++) += val;
+    } // for-loop over y-positions
+  } // for-loop over x-positions
+}
+
+// matlab entry point
+// C = fconv(A, cell of B, start, end);
+// for convolving a single filter with an image, call fconv2D( A, B, 1, 1)
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
+  if (nrhs != 4)
+    mexErrMsgTxt("Wrong number of inputs"); 
+  if (nlhs != 1)
+    mexErrMsgTxt("Wrong number of outputs");
+
+  // get input A, e.g., an image
+  const mxArray *mxA = prhs[0];
+  // safety check for input A
+  if ( mxGetNumberOfDimensions(mxA) != 2 /* only single dimension feature matrices supported here*/ || 
+       mxGetClassID(mxA) != mxDOUBLE_CLASS
+     )
+  {
+    mexErrMsgTxt("Invalid input: A");
+  }
+
+  // get input B, e.g., a cell array of filters,  and start/end indices for cells to use
+  const mxArray *cellB = prhs[1];
+  mwSize num_bs = mxGetNumberOfElements(cellB);  
+  
+  // with which filter to start
+  int start = (int)mxGetScalar(prhs[2]) - 1;
+  // with which filter to end
+  int end = (int)mxGetScalar(prhs[3]) - 1;
+  
+  // safety check for start and end wrt input cell array B
+  if ( (start < 0) || (end >= num_bs) || (start > end) )
+  {
+    mexErrMsgTxt("Invalid input: start/end");
+  }
+  
+  int i_numOfFilters = end-start+1;
+  
+
+  // output cell
+  plhs[0] = mxCreateCellMatrix(1, i_numOfFilters);
+
+  // do convolutions
+  thread_data td;
+  const mwSize *A_dims = mxGetDimensions(mxA);
+  double *A = (double *)mxGetPr(mxA);
+  
+  for ( int i = 0; i < i_numOfFilters; i++ )
+  {
+    const mxArray *mxB = mxGetCell(cellB, i+start);
+    
+    td.A_dims = A_dims;
+    td.A = A;
+    
+    td.B_dims = mxGetDimensions(mxB);
+    td.B = (double *)mxGetPr(mxB);
+    
+    // safety check for input B
+    if ( mxGetNumberOfDimensions(mxB) != 2 /* only single dimension feature matrices supported here*/ ||
+         mxGetClassID(mxB) != mxDOUBLE_CLASS
+       )
+    {
+      mexErrMsgTxt("Invalid input: B");
+    }
+
+    // compute size of output
+    int height = td.A_dims[0] - td.B_dims[0] + 1;
+    int width  = td.A_dims[1] - td.B_dims[1] + 1;
+    
+    if ( (height < 1) || (width < 1) )
+      mexErrMsgTxt("Invalid input: B should be smaller than A");
+    
+    td.C_dims[0] = height;
+    td.C_dims[1] = width;
+    td.mxC = mxCreateNumericArray(2, td.C_dims, mxDOUBLE_CLASS, mxREAL);
+    td.C = (double *)mxGetPr(td.mxC);
+    
+    // call the main function doing the hacky convolution
+    process( (void *)&td );
+    mxSetCell(plhs[0], i, td.mxC);
+  }
+}
+
+

features/fconv3D.cc

@@ -0,0 +1,183 @@
+#include "mex.h"
+#include <math.h>
+#include <string.h>
+
+/*
+ * This code is used for computing filter responses.  It computes the
+ * response of a set of filters with a feature map.  
+ *
+ * Basic version, relatively slow but very compatible.
+ */
+
+struct thread_data {
+  double *A;
+  double *B;
+  double *C;
+  mxArray *mxC;
+  const mwSize *A_dims;
+  const mwSize *B_dims;
+  mwSize C_dims[2];
+};
+
+// convolve A and B
+void *process(void *thread_arg)
+{
+  thread_data *args = (thread_data *)thread_arg;
+  
+  //input A - e.g., an image
+  double *A = args->A;
+  
+  //input B - e.g., a filter
+  double *B = args->B;
+  
+  //output C - the convolution response map
+  double *C = args->C;
+  
+  const mwSize *A_dims = args->A_dims;
+  const mwSize *B_dims = args->B_dims;
+  const mwSize *C_dims = args->C_dims;
+
+  int num_features = args->A_dims[2];
+
+  // run over all dimensions (feature planes)
+  for (int f = 0; f < num_features; f++)
+  {
+    //new pointers just for comparability with original fconv.cc
+    double *dst = C;
+    double *A_src = A + f*A_dims[0]*A_dims[1];      
+    double *B_src = B + f*B_dims[0]*B_dims[1];
+    
+    //run over possible locations on x dimension
+    for (int x = 0; x < C_dims[1]; x++)
+    {
+      //run over possible locations on y dimension
+      for (int y = 0; y < C_dims[0]; y++)
+      {
+        // compute convolution score for current position
+        // loop over filter size        
+        double val = 0;
+        for (int xp = 0; xp < B_dims[1]; xp++)
+        {
+          double *A_off = A_src + (x+xp)*A_dims[0] + y;
+          double *B_off = B_src + xp*B_dims[0];
+          
+          //hope-to-be efficient version for up to 20 dimensions
+          switch(B_dims[0])
+          {
+            case 20: val += A_off[19] * B_off[19];
+            case 19: val += A_off[18] * B_off[18];
+            case 18: val += A_off[17] * B_off[17];
+            case 17: val += A_off[16] * B_off[16];
+            case 16: val += A_off[15] * B_off[15];
+            case 15: val += A_off[14] * B_off[14];
+            case 14: val += A_off[13] * B_off[13];
+            case 13: val += A_off[12] * B_off[12];
+            case 12: val += A_off[11] * B_off[11];
+            case 11: val += A_off[10] * B_off[10];
+            case 10: val += A_off[9]  * B_off[9];
+            case 9:  val += A_off[8]  * B_off[8];
+            case 8:  val += A_off[7]  * B_off[7];
+            case 7:  val += A_off[6]  * B_off[6];
+            case 6:  val += A_off[5]  * B_off[5];
+            case 5:  val += A_off[4]  * B_off[4];
+            case 4:  val += A_off[3]  * B_off[3];
+            case 3:  val += A_off[2]  * B_off[2];
+            case 2:  val += A_off[1]  * B_off[1];
+            case 1:  val += A_off[0]  * B_off[0];
+            break;
+             // less efficient version with more than 20 dimensions
+            default:
+              for (int yp = 0; yp < B_dims[0]; yp++)
+              {
+                val += *(A_off++) * *(B_off++);
+              }
+          } // nasty switch-statement
+        } //for-loop over filtersize
+        //write value to current position in output data and increment output pointer
+        *(dst++) += val;
+      }// for-loop over y-positions
+    }// for-loop over x-positions
+  }// for-loop over dimensions (feature planes)
+}
+
+// matlab entry point
+// C = fconv(A, cell of B, start, end);
+// for convolving a multi-dim filter with a multi-dim image, call fconv3D( A, B, 1, 1)
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
+  if (nrhs != 4)
+    mexErrMsgTxt("Wrong number of inputs"); 
+  if (nlhs != 1)
+    mexErrMsgTxt("Wrong number of outputs");
+
+ // get input A, e.g., a RGB image
+  const mxArray *mxA = prhs[0];
+  // safety check for input A
+  if (mxGetNumberOfDimensions(mxA) != 3 /* only feature matrices with more than 1 dimension supported here*/|| 
+      mxGetClassID(mxA) != mxDOUBLE_CLASS)
+    mexErrMsgTxt("Invalid input: A");
+
+  // get input B, e.g., a cell array of filters,  and start/end indices for cells to use
+  const mxArray *cellB = prhs[1];
+  mwSize num_bs = mxGetNumberOfElements(cellB);  
+  
+  // with which filter to start
+  int start = (int)mxGetScalar(prhs[2]) - 1;
+  // with which filter to end
+  int end = (int)mxGetScalar(prhs[3]) - 1;
+  
+  // safety check for start and end wrt input cell array B
+  if ( (start < 0) || (end >= num_bs) || (start > end) )
+  {
+    mexErrMsgTxt("Invalid input: start/end");
+  }
+  
+  int i_numOfFilters = end-start+1;
+
+  // output cell
+  plhs[0] = mxCreateCellMatrix(1, i_numOfFilters);
+
+  // do convolutions
+  thread_data td;
+  
+  const mwSize *A_dims = mxGetDimensions(mxA);
+  double *A = (double *)mxGetPr(mxA);
+  
+  for ( int i = 0; i < i_numOfFilters; i++ )
+  {
+    const mxArray *mxB = mxGetCell(cellB, i+start);
+    
+    td.A_dims = A_dims;
+    td.A = A;
+    
+    td.B_dims = mxGetDimensions(mxB);
+    td.B = (double *)mxGetPr(mxB);
+    
+    // safety check for input B
+    if (mxGetNumberOfDimensions(mxB) != 3 /* only feature matrices with more than 1 dimension supported here*/||
+        mxGetClassID(mxB) != mxDOUBLE_CLASS ||
+        td.A_dims[2] != td.B_dims[2])
+    {
+      mexErrMsgTxt("Invalid input: B");
+    }
+
+    // compute size of output
+    int height = td.A_dims[0] - td.B_dims[0] + 1;
+    int width  = td.A_dims[1] - td.B_dims[1] + 1;
+    
+    if ( (height < 1) || (width < 1) )
+      mexErrMsgTxt("Invalid input: B should be smaller than A");
+    
+    td.C_dims[0] = height;
+    td.C_dims[1] = width;
+    td.mxC = mxCreateNumericArray(2, td.C_dims, mxDOUBLE_CLASS, mxREAL);
+    td.C = (double *)mxGetPr(td.mxC);
+    
+    // call the main function doing the hacky convolution
+    process( (void *)&td );
+  
+    /* Assign the new value to the ith cell. */
+    mxSetCell(plhs[0], i, td.mxC);
+  }
+}
+
+

features/featPyramidGeneric.m

@@ -0,0 +1,180 @@
+function pyra = featPyramidGeneric(im, model, settings )
+% function pyra = featPyramidGeneric(im, model, settings )
+% 
+%  author: Alexander Freytag
+%  date:   26-02-2014 ( dd-mm-yyyy )
+% 
+%  BRIEF: 
+%     Computes a feature pyramid from an image. Features are padded such
+%     that at least a single cell is visible. For every scale factor, the
+%     resized image is as often downsampled by factor 2 until less than 
+%     i_numCells *i_binSize pixel are available.
+%
+%  INPUT:
+%     im       -- image to extract feature pyramid from, gray scale and 
+%                 RGB images are supported right now
+% 
+%     model    -- struct, with at least the following fields 
+%         .interval   -- scalar, number of octaves for rescaling the img
+%         .i_numCells -- 2x1, number of cells in y and x direction
+%         .i_binSize  -- scalar, number of pixel of every cell in x and y direction 
+% 
+%     settings -- (optional) struct, with possible fields
+%         .fh_featureExtractor -- (struct), contains name and mfunction for
+%                                 feature extraction
+%         .d_maxRelBoxExtent   -- what shall be the max extent of a box
+%                                 features are extracted from? (1.0 allows
+%                                 for responses covering the whole image)
+% 
+%  OUTPUT: 
+%     pyra.feat{i} is the i-th level of the feature pyramid.
+%     pyra.scales{i} is the scaling factor used for the i-th level.
+%     pyra.feat{i+interval} is computed at exactly half the resolution of feat{i}.
+% 
+
+
+    %% (1) check input
+    if ( nargin < 3 )
+        settings = [];
+    end
+    
+    %check for feature extractor, if not existing, set to default
+    fh_featureExtractor = struct('name','Compute patch means', 'mfunction',@computePatchMeans);
+    settings = addDefaultVariableSetting( settings, 'fh_featureExtractor', fh_featureExtractor, settings );    
+    
+    %check for feature extractor, if not existing, set to []
+    settingsFeatureExtractorDefault = [];
+    settings = addDefaultVariableSetting( settings, 'settingsFeatureExtractor', settingsFeatureExtractorDefault, settings );
+
+    %% (2) compute features
+    interval    = model.interval;
+    i_numCells  = model.i_numCells;
+    i_binSize   = model.i_binSize;
+
+    % Select padding, allowing for at least one cell in model to be visible
+    if (getFieldWithDefault ( settings, 'b_padFeatPyramid', true ) )
+        padx = max(model.i_numCells(2)-1,0);
+        pady = max(model.i_numCells(1)-1,0);
+    else
+        padx = 0;
+        pady = 0;
+    end
+
+    % Even padding allows for consistent spatial relations across 2X scales
+    padx = floor(padx/2)*2;
+    pady = floor(pady/2)*2;
+
+        
+
+    sc = 2 ^(1/interval);
+    imsize = [size(im, 1) size(im, 2)];
+    
+    % a value of 1.0 allows responses covering the whole image, whereas 0.5
+    % results on boxes covering 25% at most
+    d_maxRelBoxExtent = getFieldWithDefault ( settings, 'd_maxRelBoxExtent', 1.0 );
+    imsizeMaxAcceptable = d_maxRelBoxExtent*imsize;
+    
+    % determine maximum scale such that at least i_numCells *i_binSize px
+    % are visible in the scaled image
+    maxScalesProLevel = zeros( interval, 1);
+    for i = 1:interval
+        maxScalesProLevel(i) = ...
+                1 + floor(   log2(    1/sc^(i-1).*min( double(imsizeMaxAcceptable) ...
+                                             ./ ...
+                                         (i_numCells *i_binSize) ...
+                                        ) ...
+                                ) ...
+                         );
+    end
+
+    %pre-allocate memory                     
+    pyra.feat  = cell(  sum( maxScalesProLevel), 1);
+    pyra.scale = zeros( sum( maxScalesProLevel), 1);
+    % our resize function wants floating point values
+    im = double(im);
+    %%TODO check whether double<-> uint8 is still needed for meanPatch features
+
+    if ( length(size(im)) == 2)
+        resizeFct = @resizeGrayScale;
+        reduceFct = @reduceGrayScale;
+    else
+        resizeFct = @resizeColor;
+        reduceFct = @reduceColor;
+    end
+
+    
+    settings.settingsFeatureExtractor.i_binSize = i_binSize;    
+    
+    % every iteration rescales the image with exponentially decreasing
+    % factors, corresponds to the scale space of Lowes SIFT paper
+    for i = 1:interval
+        %rescale image according to current scaling factor
+        scaledImg = resizeFct(im, 1/sc^(i-1));        
+        
+        %initial scaling factor for this iteration
+        pyra.scale( i ) = 1/sc^(i-1);
+        
+        % always doubles support of cells by shrinking the image by a factor of
+        % 2 in every iteration, corresponds to octaves of Lowes SIFT paper        
+        for j = 1:maxScalesProLevel(i)
+            
+            %extracat features for current level
+            pyra.feat{  i+(j-1)*interval } = settings.fh_featureExtractor.mfunction ( scaledImg, settings.settingsFeatureExtractor );
+            
+            % the scale is exactly have the scale of the previous level
+            pyra.scale( i+(j-1)*interval ) = 0.5^(j-1) * pyra.scale(i);
+            
+            % shrink image by factor of 2 for next iteration in inner round
+            scaledImg = reduceFct(scaledImg);            
+        end
+    end
+
+    %% write results to output object
+    % write boundary occlusion feature
+    
+    %Ricks code: td = model.features.truncation_dim;
+    if ( strcmp ( settings.fh_featureExtractor.name, 'Compute HOG features using WHO code' ) || ...
+         strcmp ( settings.fh_featureExtractor.name, 'HOG and Patch Means concatenated' ) ...
+       )  
+        td = 32;
+    else
+        % no truncation for all other known feature types
+        td = 0;
+    end
+    if ( strcmp ( settings.fh_featureExtractor.name, 'Compute HOG features using WHO code' ) )  
+        % add 1 to padding because feature generation deletes a 1-cell
+        % wide border around the feature map        
+        padAdd = 1;
+    else
+        % nothing to change for padding
+        padAdd = 0;
+    end    
+    
+    for i = 1:length(pyra.feat)
+        
+        if ( ndims ( pyra.feat{i} ) > 2 )
+            pyra.feat{i} = myPadArray(pyra.feat{i}, [pady+padAdd padx+padAdd 0], 0);
+        else
+            pyra.feat{i} = myPadArray(pyra.feat{i}, [pady+padAdd padx+padAdd], 0);
+        end
+
+
+        % possibly correct the truncation feature
+        if ( td > 0 )
+            pyra.feat{i}(1:pady,            :,       td) = 1;
+            pyra.feat{i}(end-pady:end,      :,       td) = 1;
+            pyra.feat{i}(:,            1:padx,       td) = 1;
+            pyra.feat{i}(:,            end-padx:end, td) = 1;
+        end
+    end
+      
+
+    %% add further settings to output object
+    pyra.scale    = model.i_binSize./pyra.scale;
+    pyra.interval = interval;
+    pyra.imy      = imsize(1);
+    pyra.imx      = imsize(2);
+    pyra.pady     = pady;
+    pyra.padx     = padx;
+
+end

features/featuresHOGColor.cc

@@ -0,0 +1,407 @@
+#include <math.h>
+#include "mex.h"
+
+// small value, used to avoid division by zero
+#define eps 0.0001
+
+// unit vectors used to compute gradient orientation
+double uu[9] = {1.0000, 
+                0.9397, 
+                0.7660, 
+                0.500, 
+                0.1736, 
+                -0.1736, 
+                -0.5000, 
+                -0.7660, 
+                -0.9397
+                };
+                
+double vv[9] = {0.0000, 
+                0.3420, 
+                0.6428, 
+                0.8660, 
+                0.9848, 
+                0.9848, 
+                0.8660, 
+                0.6428, 
+                0.3420
+               };
+
+static inline double min(double x, double y) { return (x <= y ? x : y); }
+static inline double max(double x, double y) { return (x <= y ? y : x); }
+
+static inline int min(int x, int y) { return (x <= y ? x : y); }
+static inline int max(int x, int y) { return (x <= y ? y : x); }
+
+// main function:
+// takes a double color image and a bin size 
+// returns HOG features
+mxArray *process(const mxArray *mximage, const mxArray *mxsbin)
+{
+  double *im = (double *)mxGetPr(mximage);
+  const int *dims = mxGetDimensions(mximage);
+  if (mxGetNumberOfDimensions(mximage) != 3 ||
+      dims[2] != 3 ||
+      mxGetClassID(mximage) != mxDOUBLE_CLASS)
+    mexErrMsgTxt("Invalid input");
+
+  // size of a cell in pixel, used for both dimensions  
+  int sbin = (int)mxGetScalar(mxsbin);
+
+  // memory for caching orientation histograms & their norms
+  int blocks[2];
+  // compute number of cells fitting into current image given speficied cell size in pixel
+  // use floor to prevent running over image borders
+  blocks[0] = (int)floor( (double)dims[0] / (double)sbin );
+  blocks[1] = (int)floor( (double)dims[1] / (double)sbin ); 
+  
+  // pre-allocate memory  
+  double *hist = (double *)mxCalloc(blocks[0]*blocks[1]*18, sizeof(double));
+  double *norm = (double *)mxCalloc(blocks[0]*blocks[1], sizeof(double));
+
+  // memory for HOG features
+  int out[3];
+  // we compute as many blocks as possible given the specified cell size and the current image
+  out[0] = max(blocks[0], 0);
+  out[1] = max(blocks[1], 0);  
+  // note: previously, a subtraction of 2 was done to guarantee identical bilinear interpolation behaviour in all cells
+  // however, are more interested in obtaining a correct number of returned cells, rather than in slightly more consistent interpolation results...
+  // apart from that, the only cells affected are the ones on the top and left border of the cell array
+  
+  out[2] = 27+4+1;
+  mxArray *mxfeat = mxCreateNumericArray(3, out, mxDOUBLE_CLASS, mxREAL);
+  double *feat = (double *)mxGetPr(mxfeat);
+  
+  int visible[2];
+  visible[0] = blocks[0]*sbin;
+  visible[1] = blocks[1]*sbin;
+  
+  // If you want to be 100% save, init output values here.
+  // Note that this should not be needed, since every value will be set lateron explicitely.
+  //
+  //double *dst = feat;
+  //for ( int tmp = 0; tmp < out[0]*out[1]*out[2]; tmp++, dst++ )
+  //{
+  // *dst = 0.0;
+  //}
+  
+  // start by 1 and end by -1 to ensure safe gradient calculations on boundaries  
+  for (int x = 1; x < visible[1]-1; x++)
+  {
+    for (int y = 1; y < visible[0]-1; y++)
+    {
+      // first color channel
+      // NOTE: why minimum check? boundary safety is given by for loop ends      
+      //               col offset                  current row position      
+      double *s = im + min(x, dims[1]-1)*dims[0] + min(y, dims[0]-1);
+      //gradient in y direction
+      double dy = *(s+1) - *(s-1);
+      // gradient in x direction
+      double dx = *(s+dims[0]) - *(s-dims[0]);
+      // squared gradient strength on current pixel
+      double v = dx*dx + dy*dy;
+
+      // second color channel
+      s += dims[0]*dims[1];
+      //gradient in y direction
+      double dy2 = *(s+1) - *(s-1);
+      // gradient in x direction
+      double dx2 = *(s+dims[0]) - *(s-dims[0]);
+      // squared gradient strength on current pixel
+      double v2 = dx2*dx2 + dy2*dy2;
+
+      // third color channel
+      s += dims[0]*dims[1];
+      //gradient in y direction
+      double dy3 = *(s+1) - *(s-1);
+      // gradient in x direction
+      double dx3 = *(s+dims[0]) - *(s-dims[0]);
+      // squared gradient strength on current pixel
+      double v3 = dx3*dx3 + dy3*dy3;
+
+      // pick channel with strongest gradient
+      if (v2 > v)
+      {
+        v = v2;
+        dx = dx2;
+        dy = dy2;
+      } 
+      if (v3 > v)
+      {
+        v = v3;
+        dx = dx3;
+        dy = dy3;
+      }
+
+      //
+      // now, discretize gradient orientation into one of 18 possible (oriented) bins 
+      //
+      
+      // strength of strongest orientation in this pixel
+      double best_dot = 0;
+      // index of strongest orientation in this pixel
+      int best_o = 0;
+      
+      for (int o = 0; o < 9; o++)
+      {
+        double dot = uu[o]*dx + vv[o]*dy;
+        if (dot > best_dot)
+        {
+          best_dot = dot;
+          best_o = o;
+        }
+        else if (-dot > best_dot)
+        {
+          best_dot = -dot;
+          best_o = o+9;
+        }
+      }
+      
+      // add to 4 histograms around pixel using linear interpolation
+      
+      // current position normalized to cell scale, e.g. xp = 1.4 -> located in the left part of second cell
+      // subtraction of 0.5 to move relative to cell center
+      double xp  = ((double)x+0.5)/(double)sbin - 0.5;      
+      double yp  = ((double)y+0.5)/(double)sbin - 0.5;   
+      
+      // that's the index of the cell, whose center is directly left of current position in x direction
+      int ixp    = (int)floor(xp);
+      // that's the index of the cell, whose center is directly on top of current position in y direction      
+      int iyp    = (int)floor(yp);
+      
+      // distance to left, used for weighting the gradient strength by 1-distance, guaranteed to be in [0,1]      
+      double vx0 = xp-ixp;
+      // distance to top
+      double vy0 = yp-iyp;
+      // distance to right
+      double vx1 = 1.0-vx0;
+      // distance to bottom
+      double vy1 = 1.0-vy0;
+      
+      // normalized gradient strength on current pixel
+      v = sqrt(v);
+
+      // if left upper cell exists
+      if ( (ixp >= 0) && (iyp >= 0) )
+      {
+        *(hist + ixp*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
+        vx1*vy1*v; // i.e., (1-distX0)*(1-distY0)*v
+      }
+
+      // if right upper cell exists
+      if ( (ixp+2 < blocks[1]) && (iyp >= 0) )
+      {
+        *(hist + (ixp+1)*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
+        vx0*vy1*v;
+      }
+
+      // if left lower cell exists
+      if ( (ixp >= 0) && (iyp+2 < blocks[0]) )
+      {
+        *(hist + ixp*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
+        vx1*vy0*v;
+      }
+
+      // if right lower cell exists
+      if ( (ixp+2 < blocks[1]) && (iyp+2 < blocks[0]) )
+      {
+        *(hist + (ixp+1)*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
+        vx0*vy0*v;
+      }
+    }
+  }
+
+  // compute energy in each block by summing over orientations
+  for (int o = 0; o < 9; o++)
+  {
+    double *src1 = hist + o*blocks[0]*blocks[1];
+    double *src2 = hist + (o+9)*blocks[0]*blocks[1];
+    double *dst  = norm;
+    double *end  = norm + blocks[1]*blocks[0];
+    
+    while (dst < end)
+    {
+      *(dst++) += (*src1 + *src2) * (*src1 + *src2);
+      src1++;
+      src2++;
+    }
+  }
+
+  // compute features
+  for (int x = 0; x < out[1]; x++)
+  {
+    for (int y = 0; y < out[0]; y++)
+    {
+      double *dst = feat + x*out[0] + y;      
+      double *src, *p, n1, n2, n3, n4;
+
+      // compute normalization factors for all 4 possible blocks of 2x2 cells
+      
+      
+      // block with current, right, lower, and lower right cell
+      if ( ( (x+1) < out[1] ) && ( (y+1) < out[0] ) )
+      {
+        p = norm + x*blocks[0] + y;
+        n1 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      }
+      else 
+        n1 = 1.0 / sqrt(eps);      
+      
+      // block with current, upper, right, and upper right cell
+      if ( ( (x+1) < out[1] ) && ( (y-1) >= 0 ) )
+      {      
+        p = norm + x*blocks[0] + (y-1);
+        n2 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      }    
+      else 
+        n2 = 1.0 / sqrt(eps);      
+      
+      // block with current, lower, left, and lower left cell
+      if ( ( (x-1) >= 0 ) && ( (y+1) < out[0] ) )
+      {   
+        p = norm + (x-1)*blocks[0] + y;
+        n3 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      }      
+      else 
+        n3 = 1.0 / sqrt(eps);      
+    
+      // block with current, upper, left, and upper left cell
+      if ( ( (x-1) >= 0 ) && ( (y-1) >= 0 ) )
+      {
+        p = norm + (x-1)*blocks[0] + (y-1);
+        n4 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      }
+      else 
+        n4 = 1.0 / sqrt(eps);
+              
+      
+      
+      // copy normalization factors for blocks on boundaries
+      //   -----------------
+      //   | n4 |     | n2 |
+      //   -----------------
+      //   |    | x,y |    |
+      //   -----------------
+      //   | n3 |     | n1 |
+      //   -----------------      
+      if ( (x) == 0 ) // left boundary
+      {
+        if ( (y) == 0 ) // left top corner
+        {
+          n4 = n1; n3 = n1; n2 = n1;
+        }
+        else if ( (y) == (out[0]-1) ) // left lower corner
+        {
+          n4 = n2; n3 = n2; n1 = n2;
+        }
+        else
+        {
+          n4 = n2; n3 = n1;
+        }
+      }
+      else if ( (x) == (out[1]-1) ) // right boundary
+      {
+        if ( (y) == 0 ) // right top corner
+        {
+          n4 = n3; n2 = n3; n1 = n3;
+        }
+        else if ( (y) == (out[0]-1) ) // right lower corner
+        {
+          n3 = n4; n2 = n4; n1 = n4;
+        }
+        else
+        {
+          n2 = n4; n1 = n3; 
+        }
+      }
+      if ( (y) == 0 ) // upper boundary ( corners already tackled)
+      {
+        if ( ( x > 0 ) && ( x < (out[1]-1) ) )
+        {
+          n4 = n3; n2 = n1;
+        }
+      }
+      else if ( (y) == (out[0]-1) ) // lower boundary ( corners already tackled)
+      {
+        if ( ( x > 0 ) && ( x < (out[1]-1) ) )
+        {
+          n3 = n4; n1 = n2;
+        }
+      }     
+      
+      double t1 = 0;
+      double t2 = 0;
+      double t3 = 0;
+      double t4 = 0;
+      
+      // contrast-sensitive features
+      src = hist + (x)*blocks[0] + (y);
+      for (int o = 0; o < 18; o++)
+      {
+        double h1 = min(*src * n1, 0.2);
+        double h2 = min(*src * n2, 0.2);
+        double h3 = min(*src * n3, 0.2);
+        double h4 = min(*src * n4, 0.2);
+        *dst = 0.5 * (h1 + h2 + h3 + h4);
+        t1 += h1;
+        t2 += h2;
+        t3 += h3;
+        t4 += h4;
+        // move pointers to next position
+        dst += out[0]*out[1];
+        src += blocks[0]*blocks[1];
+      }
+      
+      // contrast-insensitive features
+      src = hist + (x)*blocks[0] + (y);
+      for (int o = 0; o < 9; o++)
+      {
+        double sum = *src + *(src + 9*blocks[0]*blocks[1]);
+        double h1 = min(sum * n1, 0.2);
+        double h2 = min(sum * n2, 0.2);
+        double h3 = min(sum * n3, 0.2);
+        double h4 = min(sum * n4, 0.2);
+        *dst = 0.5 * (h1 + h2 + h3 + h4);
+        dst += out[0]*out[1];
+        src += blocks[0]*blocks[1];
+      }
+
+      // texture features
+      //to be complicable to FFLD code
+      *dst = 0.2357 * t4;
+      //*dst = 0.2357 * t1;
+      dst += out[0]*out[1];
+      *dst = 0.2357 * t2;
+      dst += out[0]*out[1];
+      *dst = 0.2357 * t3;
+      dst += out[0]*out[1];
+      //to be complicable to FFLD code
+      *dst = 0.2357 * t1;
+      //*dst = 0.2357 * t4;
+
+      // truncation feature
+      dst += out[0]*out[1];
+      *dst = 0;
+    }
+  }
+
+
+  mxFree(hist);
+  mxFree(norm);
+  return mxfeat;
+}
+
+// matlab entry point
+// F = features(image, bin)
+// image should be color with double values
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{ 
+  if (nrhs != 2)
+    mexErrMsgTxt("Wrong number of inputs"); 
+  if (nlhs != 1)
+    mexErrMsgTxt("Wrong number of outputs");
+  plhs[0] = process(prhs[0], prhs[1]);
+}
+
+
+

features/featuresHOGColor.m

@@ -0,0 +1,25 @@
+% function hogFeature = featuresHOGColor( img, sbin );
+% 
+% BRIEF:
+%    Compute an array of extended HOG features for a given color image. 
+%    Differs from the original implementation (featuresHOGorig.cc) by not 
+%    leaving a boundary of sbin pixels at the border of the image "unused".
+% 
+%    Advantage:    k*sbin pixel result in k cells (orig: k-2)
+%    Disadvantage: block normalization ill posed on boundary
+%
+% INPUT:
+%    img   -- (x,y,3) double array, input image
+%    sbin  -- double scalar, number of pixels each cell covers in x and y
+%             direction
+%
+% OUTPUT:
+%    hogFeature   -- (x/sbin, y/sbin, 32) double array,
+%                    extracted hog array, last dim equals 0
+% 
+% NOTE:
+%    Don't miss to mex (compile) the .cc-file!
+%   
+% author: Alexander Freytag
+% last update: 11-03-2014 ( dd-mm-yyyy )
+% 

features/featuresHOGGrayScale.cc

@@ -0,0 +1,382 @@
+#include <math.h>
+#include "mex.h"
+
+// small value, used to avoid division by zero
+#define eps 0.0001
+
+// unit vectors used to compute gradient orientation
+double uu[9] = {1.0000, 
+                0.9397, 
+                0.7660, 
+                0.500, 
+                0.1736, 
+                -0.1736, 
+                -0.5000, 
+                -0.7660, 
+                -0.9397
+                };
+                
+double vv[9] = {0.0000, 
+                0.3420, 
+                0.6428, 
+                0.8660, 
+                0.9848, 
+                0.9848, 
+                0.8660, 
+                0.6428, 
+                0.3420
+               };
+
+static inline double min(double x, double y) { return (x <= y ? x : y); }
+static inline double max(double x, double y) { return (x <= y ? y : x); }
+
+static inline int min(int x, int y) { return (x <= y ? x : y); }
+static inline int max(int x, int y) { return (x <= y ? y : x); }
+
+// main function:
+// takes a double gray scale image and a bin size 
+// returns HOG features
+mxArray *process(const mxArray *mximage, const mxArray *mxsbin)
+{
+  double *im = (double *)mxGetPr(mximage);
+  const int *dims = mxGetDimensions(mximage);
+  
+    
+  if (mxGetNumberOfDimensions(mximage) != 2 ||
+      mxGetClassID(mximage) != mxDOUBLE_CLASS)
+        mexErrMsgTxt("Invalid input, either not a gray scale img (!= 2 dim), or input not given as double matrix.");
+
+  // size of a cell in pixel, used for both dimensions
+  int sbin = (int)mxGetScalar(mxsbin);
+
+  // memory for caching orientation histograms & their norms
+  int blocks[2];
+  // compute number of cells fitting into current image given speficied cell size in pixel
+  // use floor to prevent running over image borders
+  blocks[0] = (int)floor( (double)dims[0] / (double)sbin );
+  blocks[1] = (int)floor( (double)dims[1] / (double)sbin );  
+  
+  // pre-allocate memory
+  double *hist = (double *)mxCalloc( blocks[0]*blocks[1]*18, sizeof(double) );
+  double *norm = (double *)mxCalloc( blocks[0]*blocks[1],    sizeof(double) );
+  
+  // memory for HOG features
+  int out[3];
+  // we compute as many blocks as possible given the specified cell size and the current image
+  out[0] = max(blocks[0], 0);
+  out[1] = max(blocks[1], 0);  
+  // note: previously, a subtraction of 2 was done to guarantee identical bilinear interpolation behaviour in all cells
+  // however, are more interested in obtaining a correct number of returned cells, rather than in slightly more consistent interpolation results...
+  // apart from that, the only cells affected are the ones on the top and left border of the cell array
+  
+  out[2] = 27+4+1;
+  mxArray *mxfeat = mxCreateNumericArray(3, out, mxDOUBLE_CLASS, mxREAL);
+  double *feat = (double *)mxGetPr(mxfeat);
+  
+   
+  int visible[2];
+  visible[0] = blocks[0]*sbin;
+  visible[1] = blocks[1]*sbin;
+  
+  // If you want to be 100% save, init output values here.
+  // Note that this should not be needed, since every value will be set lateron explicitely.
+  //
+  //double *dst = feat;
+  //for ( int tmp = 0; tmp < out[0]*out[1]*out[2]; tmp++, dst++ )
+  //{
+  // *dst = 0.0;
+  //}  
+  
+  // start by 1 and end by -1 to ensure safe gradient calculations on boundaries
+  for (int x = 1; x < visible[1]-1; x++)
+  {
+    for (int y = 1; y < visible[0]-1; y++)
+    {
+      // simple gray scale channel
+      // NOTE: why minimum check? boundary safety is given by for loop ends      
+      //               col offset                  current row position
+      double *s = im + min(x, dims[1]-1)*dims[0] + min(y, dims[0]-1);
+      //gradient in y direction
+      double dy = *(s+1) - *(s-1);
+      // gradient in x direction
+      double dx = *(s+dims[0]) - *(s-dims[0]);
+      // squared gradient strength on current pixel
+      double v = dx*dx + dy*dy;
+
+
+      //
+      // now, discretize gradient orientation into one of 18 possible (oriented) bins 
+      //
+      
+      // strength of strongest orientation in this pixel
+      double best_dot = 0;
+      // index of strongest orientation in this pixel
+      int best_o = 0;
+      
+      for (int o = 0; o < 9; o++)
+      {
+        double dot = uu[o]*dx + vv[o]*dy;
+        
+        if (dot > best_dot)
+        {
+          best_dot = dot;
+          best_o   = o;
+        }
+        else if (-dot > best_dot)
+        {
+          best_dot = -dot;
+          best_o   = o+9;
+        }
+      }
+      
+      // add to 4 histograms around pixel using linear interpolation
+      
+      // current position normalized to cell scale, e.g. xp = 1.4 -> located in the left part of second cell
+      // subtraction of 0.5 to move relative to cell center
+      double xp  = ((double)x+0.5)/(double)sbin - 0.5;      
+      double yp  = ((double)y+0.5)/(double)sbin - 0.5;      
+      
+      // that's the index of the cell, whose center is directly left of current position in x direction
+      int ixp    = (int)floor(xp);
+      // that's the index of the cell, whose center is directly on top of current position in y direction      
+      int iyp    = (int)floor(yp);
+      
+      // check whether current pixel is nearer to left or right boundary of cell
+      
+      // distance to left, used for weighting the gradient strength by 1-distance, guaranteed to be in [0,1]
+      double vx0 = xp-ixp;
+      // distance to top
+      double vy0 = yp-iyp;
+      // distance to right
+      double vx1 = 1.0-vx0;
+      // distance to bottom
+      double vy1 = 1.0-vy0;
+      
+      // normalized gradient strength on current pixel
+      v = sqrt(v);
+
+      // if left upper cell exists
+      if ( (ixp >= 0) && (iyp >= 0) )
+      {
+        *(hist + ixp*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
+        vx1*vy1*v; // i.e., (1-distX0)*(1-distY0)*v
+      }
+
+      // if right upper cell exists
+      if ( (ixp+2 < blocks[1]) && (iyp >= 0) )
+      {
+        *(hist + (ixp+1)*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
+        vx0*vy1*v;
+      }
+
+      // if left lower cell exists
+      if ( (ixp >= 0) && (iyp+2 < blocks[0]) )
+      {
+        *(hist + ixp*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
+        vx1*vy0*v;
+      }
+
+      // if right lower cell exists
+      if ( (ixp+2 < blocks[1]) && (iyp+2 < blocks[0]) )
+      {
+        *(hist + (ixp+1)*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
+        vx0*vy0*v;
+      }
+    }
+  }
+
+  // compute energy in each block by summing over orientations
+  for (int o = 0; o < 9; o++)
+  {
+    double *src1 = hist + o*blocks[0]*blocks[1];
+    double *src2 = hist + (o+9)*blocks[0]*blocks[1];
+    double *dst  = norm;
+    double *end  = norm + blocks[1]*blocks[0];
+    
+    while (dst < end)
+    {
+      *(dst++) += (*src1 + *src2) * (*src1 + *src2);
+      src1++;
+      src2++;
+    }
+  }
+
+  // compute features
+  for (int x = 0; x < out[1]; x++)
+  {
+    for (int y = 0; y < out[0]; y++)
+    {
+      double *dst = feat + x*out[0] + y;      
+      double *src, *p, n1, n2, n3, n4;
+
+      // compute normalization factors for all 4 possible blocks of 2x2 cells
+      
+      
+       
+      // block with current, right, lower, and lower right cell
+      if ( ( (x+1) < out[1] ) && ( (y+1) < out[0] ) )
+      {
+        p = norm + x*blocks[0] + y;
+        n1 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      }
+      else 
+        n1 = 1.0 / sqrt(eps);      
+      
+      // block with current, upper, right, and upper right cell
+      if ( ( (x+1) < out[1] ) && ( (y-1) >= 0 ) )
+      {      
+        p = norm + x*blocks[0] + (y-1);
+        n2 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      }    
+      else 
+        n2 = 1.0 / sqrt(eps);      
+      
+      // block with current, lower, left, and lower left cell
+      if ( ( (x-1) >= 0 ) && ( (y+1) < out[0] ) )
+      {   
+        p = norm + (x-1)*blocks[0] + y;
+        n3 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      }      
+      else 
+        n3 = 1.0 / sqrt(eps);      
+    
+      // block with current, upper, left, and upper left cell
+      if ( ( (x-1) >= 0 ) && ( (y-1) >= 0 ) )
+      {
+        p = norm + (x-1)*blocks[0] + (y-1);
+        n4 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      }
+      else 
+        n4 = 1.0 / sqrt(eps);
+      
+      
+      // copy normalization factors for blocks on boundaries
+      //   -----------------
+      //   | n4 |     | n2 |
+      //   -----------------
+      //   |    | x,y |    |
+      //   -----------------
+      //   | n3 |     | n1 |
+      //   -----------------      
+      if ( (x) == 0 ) // left boundary
+      {
+        if ( (y) == 0 ) // left top corner
+        {
+          n4 = n1; n3 = n1; n2 = n1;
+        }
+        else if ( (y) == (out[0]-1) ) // left lower corner
+        {
+          n4 = n2; n3 = n2; n1 = n2;
+        }
+        else
+        {
+          n4 = n2; n3 = n1;
+        }
+      }
+      else if ( (x) == (out[1]-1) ) // right boundary
+      {
+        if ( (y) == 0 ) // right top corner
+        {
+          n4 = n3; n2 = n3; n1 = n3;
+        }
+        else if ( (y) == (out[0]-1) ) // right lower corner
+        {
+          n3 = n4; n2 = n4; n1 = n4;
+        }
+        else
+        {
+          n2 = n4; n1 = n3; 
+        }
+      }
+      if ( (y) == 0 ) // upper boundary ( corners already tackled)
+      {
+        if ( ( x > 0 ) && ( x < (out[1]-1) ) )
+        {
+          n4 = n3; n2 = n1;
+        }
+      }
+      else if ( (y) == (out[0]-1) ) // lower boundary ( corners already tackled)
+      {
+        if ( ( x > 0 ) && ( x < (out[1]-1) ) )
+        {
+          n3 = n4; n1 = n2;
+        }
+      }       
+      
+    
+      
+      double t1 = 0;
+      double t2 = 0;
+      double t3 = 0;
+      double t4 = 0;
+      
+      // contrast-sensitive features
+      src = hist + (x)*blocks[0] + (y);
+      for (int o = 0; o < 18; o++)
+      {
+        double h1 = min(*src * n1, 0.2);
+        double h2 = min(*src * n2, 0.2);
+        double h3 = min(*src * n3, 0.2);
+        double h4 = min(*src * n4, 0.2);
+        *dst = 0.5 * (h1 + h2 + h3 + h4);
+        t1 += h1;
+        t2 += h2;
+        t3 += h3;
+        t4 += h4;
+        // move pointers to next position
+        dst += out[0]*out[1];
+        src += blocks[0]*blocks[1];
+      }
+      
+      // contrast-insensitive features
+      src = hist + (x)*blocks[0] + (y);
+      for (int o = 0; o < 9; o++)
+      {
+        double sum = *src + *(src + 9*blocks[0]*blocks[1]);
+        double h1 = min(sum * n1, 0.2);
+        double h2 = min(sum * n2, 0.2);
+        double h3 = min(sum * n3, 0.2);
+        double h4 = min(sum * n4, 0.2);
+        *dst = 0.5 * (h1 + h2 + h3 + h4);
+        dst += out[0]*out[1];
+        src += blocks[0]*blocks[1];
+      }
+
+      // texture features
+      //to be complicable to FFLD code
+      *dst = 0.2357 * t4;
+      //*dst = 0.2357 * t1;
+      dst += out[0]*out[1];
+      *dst = 0.2357 * t2;
+      dst += out[0]*out[1];
+      *dst = 0.2357 * t3;
+      dst += out[0]*out[1];
+      //to be complicable to FFLD code
+      *dst = 0.2357 * t1;
+      //*dst = 0.2357 * t4;
+
+      // truncation feature
+      dst += out[0]*out[1];
+      *dst = 0;
+    }
+  }
+
+  mxFree(hist);
+  mxFree(norm);
+  return mxfeat;
+}
+
+// matlab entry point
+// F = featuresGrayScale(image, bin)
+// image should be gray scale with double values
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
+{ 
+  if (nrhs != 2)
+    mexErrMsgTxt("Wrong number of inputs"); 
+  if (nlhs != 1)
+    mexErrMsgTxt("Wrong number of outputs");
+  plhs[0] = process(prhs[0], prhs[1]);
+}
+
+
+

features/featuresHOGGrayScale.m

@@ -0,0 +1,26 @@
+% function hogFeature = featuresHOGGrayScale( img, sbin );
+% 
+% BRIEF:
+%    Compute an array of extended HOG features for a given gray scale image. 
+%    Original implementation did not supported simply gray scale images.
+%    Additionally, differs from the original implementation (featuresHOGorig.cc)
+%    by not leaving a boundary of sbin pixels at the border of the image "unused".
+% 
+%    Advantage:    k*sbin pixel result in k cells (orig: k-2)
+%    Disadvantage: block normalization ill posed on boundary
+%
+% INPUT:
+%    img   -- (x,y) double array, input image
+%    sbin  -- double scalar, number of pixels each cell covers in x and y
+%             direction
+%
+% OUTPUT:
+%    hogFeature   -- (x/sbin, y/sbin, 32) double array,
+%                    extracted hog array, last dim equals 0
+% 
+% NOTE:
+%    Don't miss to mex (compile) the .cc-file!
+%   
+% author: Alexander Freytag
+% last update: 11-03-2014 ( dd-mm-yyyy )
+% 

features/featuresHOGorig.cc

@@ -0,0 +1,238 @@
+#include <math.h>
+#include "mex.h"
+
+// small value, used to avoid division by zero
+#define eps 0.0001
+
+// unit vectors used to compute gradient orientation
+double uu[9] = {1.0000, 
+		0.9397, 
+		0.7660, 
+		0.500, 
+		0.1736, 
+		-0.1736, 
+		-0.5000, 
+		-0.7660, 
+		-0.9397};
+double vv[9] = {0.0000, 
+		0.3420, 
+		0.6428, 
+		0.8660, 
+		0.9848, 
+		0.9848, 
+		0.8660, 
+		0.6428, 
+		0.3420};
+
+static inline double min(double x, double y) { return (x <= y ? x : y); }
+static inline double max(double x, double y) { return (x <= y ? y : x); }
+
+static inline int min(int x, int y) { return (x <= y ? x : y); }
+static inline int max(int x, int y) { return (x <= y ? y : x); }
+
+// main function:
+// takes a double color image and a bin size 
+// returns HOG features
+mxArray *process(const mxArray *mximage, const mxArray *mxsbin) {
+  double *im = (double *)mxGetPr(mximage);
+  const int *dims = mxGetDimensions(mximage);
+  if (mxGetNumberOfDimensions(mximage) != 3 ||
+      dims[2] != 3 ||
+      mxGetClassID(mximage) != mxDOUBLE_CLASS)
+    mexErrMsgTxt("Invalid input");
+
+  int sbin = (int)mxGetScalar(mxsbin);
+
+  // memory for caching orientation histograms & their norms
+  int blocks[2];
+  blocks[0] = (int)round((double)dims[0]/(double)sbin);
+  blocks[1] = (int)round((double)dims[1]/(double)sbin);
+  double *hist = (double *)mxCalloc(blocks[0]*blocks[1]*18, sizeof(double));
+  double *norm = (double *)mxCalloc(blocks[0]*blocks[1], sizeof(double));
+
+  // memory for HOG features
+  int out[3];
+  out[0] = max(blocks[0]-2, 0);
+  out[1] = max(blocks[1]-2, 0);
+  out[2] = 27+4+1;
+  mxArray *mxfeat = mxCreateNumericArray(3, out, mxDOUBLE_CLASS, mxREAL);
+  double *feat = (double *)mxGetPr(mxfeat);
+  
+  int visible[2];
+  visible[0] = blocks[0]*sbin;
+  visible[1] = blocks[1]*sbin;
+  
+  for (int x = 1; x < visible[1]-1; x++) {
+    for (int y = 1; y < visible[0]-1; y++) {
+      // first color channel
+      double *s = im + min(x, dims[1]-2)*dims[0] + min(y, dims[0]-2);
+      double dy = *(s+1) - *(s-1);
+      double dx = *(s+dims[0]) - *(s-dims[0]);
+      double v = dx*dx + dy*dy;
+
+      // second color channel
+      s += dims[0]*dims[1];
+      double dy2 = *(s+1) - *(s-1);
+      double dx2 = *(s+dims[0]) - *(s-dims[0]);
+      double v2 = dx2*dx2 + dy2*dy2;
+
+      // third color channel
+      s += dims[0]*dims[1];
+      double dy3 = *(s+1) - *(s-1);
+      double dx3 = *(s+dims[0]) - *(s-dims[0]);
+      double v3 = dx3*dx3 + dy3*dy3;
+
+      // pick channel with strongest gradient
+      if (v2 > v) {
+	v = v2;
+	dx = dx2;
+	dy = dy2;
+      } 
+      if (v3 > v) {
+	v = v3;
+	dx = dx3;
+	dy = dy3;
+      }
+
+      // snap to one of 18 orientations
+      double best_dot = 0;
+      int best_o = 0;
+      for (int o = 0; o < 9; o++) {
+	double dot = uu[o]*dx + vv[o]*dy;
+	if (dot > best_dot) {
+	  best_dot = dot;
+	  best_o = o;
+	} else if (-dot > best_dot) {
+	  best_dot = -dot;
+	  best_o = o+9;
+	}
+      }
+      
+      // add to 4 histograms around pixel using linear interpolation
+      double xp = ((double)x+0.5)/(double)sbin - 0.5;
+      double yp = ((double)y+0.5)/(double)sbin - 0.5;
+      int ixp = (int)floor(xp);
+      int iyp = (int)floor(yp);
+      double vx0 = xp-ixp;
+      double vy0 = yp-iyp;
+      double vx1 = 1.0-vx0;
+      double vy1 = 1.0-vy0;
+      v = sqrt(v);
+
+      if (ixp >= 0 && iyp >= 0) {
+	*(hist + ixp*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
+	  vx1*vy1*v;
+      }
+
+      if (ixp+1 < blocks[1] && iyp >= 0) {
+	*(hist + (ixp+1)*blocks[0] + iyp + best_o*blocks[0]*blocks[1]) += 
+	  vx0*vy1*v;
+      }
+
+      if (ixp >= 0 && iyp+1 < blocks[0]) {
+	*(hist + ixp*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
+	  vx1*vy0*v;
+      }
+
+      if (ixp+1 < blocks[1] && iyp+1 < blocks[0]) {
+	*(hist + (ixp+1)*blocks[0] + (iyp+1) + best_o*blocks[0]*blocks[1]) += 
+	  vx0*vy0*v;
+      }
+    }
+  }
+
+  // compute energy in each block by summing over orientations
+  for (int o = 0; o < 9; o++) {
+    double *src1 = hist + o*blocks[0]*blocks[1];
+    double *src2 = hist + (o+9)*blocks[0]*blocks[1];
+    double *dst = norm;
+    double *end = norm + blocks[1]*blocks[0];
+    while (dst < end) {
+      *(dst++) += (*src1 + *src2) * (*src1 + *src2);
+      src1++;
+      src2++;
+    }
+  }
+
+  // compute features
+  for (int x = 0; x < out[1]; x++) {
+    for (int y = 0; y < out[0]; y++) {
+      double *dst = feat + x*out[0] + y;      
+      double *src, *p, n1, n2, n3, n4;
+
+      p = norm + (x+1)*blocks[0] + y+1;
+      n1 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      p = norm + (x+1)*blocks[0] + y;
+      n2 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      p = norm + x*blocks[0] + y+1;
+      n3 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      p = norm + x*blocks[0] + y;      
+      n4 = 1.0 / sqrt(*p + *(p+1) + *(p+blocks[0]) + *(p+blocks[0]+1) + eps);
+      double t1 = 0;
+      double t2 = 0;
+      double t3 = 0;
+      double t4 = 0;
+	  double bh_test=0;
+      // contrast-sensitive features
+      src = hist + (x+1)*blocks[0] + (y+1);
+      for (int o = 0; o < 18; o++) {
+	bh_test=bh_test+(*src * n1)+(*src * n2)+(*src * n3)+(*src * n4);
+	double h1 = min(*src * n1, 0.2);
+	double h2 = min(*src * n2, 0.2);
+	double h3 = min(*src * n3, 0.2);
+	double h4 = min(*src * n4, 0.2);
+	*dst = 0.5 * (h1 + h2 + h3 + h4);
+	t1 += h1;
+	t2 += h2;
+	t3 += h3;
+	t4 += h4;
+	dst += out[0]*out[1];
+	src += blocks[0]*blocks[1];
+      }
+	//printf("%f\n", bh_test);
+      // contrast-insensitive features
+      src = hist + (x+1)*blocks[0] + (y+1);
+      for (int o = 0; o < 9; o++) {
+        double sum = *src + *(src + 9*blocks[0]*blocks[1]);
+        double h1 = min(sum * n1, 0.2);
+        double h2 = min(sum * n2, 0.2);
+        double h3 = min(sum * n3, 0.2);
+        double h4 = min(sum * n4, 0.2);
+        *dst = 0.5 * (h1 + h2 + h3 + h4);
+        dst += out[0]*out[1];
+        src += blocks[0]*blocks[1];
+      }
+
+      // texture features
+      *dst = 0.2357 * t1;
+      dst += out[0]*out[1];
+      *dst = 0.2357 * t2;
+      dst += out[0]*out[1];
+      *dst = 0.2357 * t3;
+      dst += out[0]*out[1];
+      *dst = 0.2357 * t4;
+
+      // truncation feature
+      dst += out[0]*out[1];
+      *dst = 0;
+    }
+  }
+
+  mxFree(hist);
+  mxFree(norm);
+  return mxfeat;
+}
+
+// matlab entry point
+// F = features(image, bin)
+// image should be color with double values
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
+  if (nrhs != 2)
+    mexErrMsgTxt("Wrong number of inputs"); 
+  if (nlhs != 1)
+    mexErrMsgTxt("Wrong number of outputs");
+  plhs[0] = process(prhs[0], prhs[1]);
+}
+
+
+

features/featuresHOGorig.m

@@ -0,0 +1,25 @@
+% function hogFeature = featuresHOGorig( img, sbin );
+% 
+% BRIEF:
+%    Compute an array of extended HOG features for a given color image. 
+%    No support for color images. Leaves a boundary of sbin pixels at 
+%    the border of the image "unused".
+% 
+%    Advantage:    block normalization well posed on boundary
+%    Disadvantage: k*sbin pixel result in k-2 cells
+%
+% INPUT:
+%    img   -- (x,y,3) double array, input image
+%    sbin  -- double scalar, number of pixels each cell covers in x and y
+%             direction
+%
+% OUTPUT:
+%    hogFeature   -- (x/sbin -2, y/sbin -2, 32) double array,
+%                    extracted hog array, last dim equals 0
+%   
+% NOTE:
+%    Don't miss to mex (compile) the .cc-file!
+% 
+% author: Alexander Freytag
+% last update: 11-03-2014 ( dd-mm-yyyy )
+% 

features/reduceColor.cc

@@ -0,0 +1,92 @@
+#include <math.h>
+#include <assert.h>
+#include <string.h>
+#include "mex.h"
+
+// reduce(im) resizes im to half its size, using a 5-tap binomial filter for anti-aliasing
+// (see Burt & Adelson's Laplacian Pyramid paper)
+
+// reduce each column
+// result is transposed, so we can apply it twice for a complete reduction
+void reduce1dtran(double *src, int sheight, double *dst, int dheight, 
+		  int width, int chan) {
+  // resize each column of each color channel
+  bzero(dst, chan*width*dheight*sizeof(double));
+  int y;
+  double *s, *d;
+
+  for (int c = 0; c < chan; c++) {
+    for (int x = 0; x < width; x++) {
+      s  = src + c*width*sheight + x*sheight;
+      d  = dst + c*dheight*width + x;
+
+      // First row
+      *d = s[0]*.6875 + s[1]*.2500 + s[2]*.0625;      
+
+      for (y = 1; y < dheight-2; y++) {	
+	s += 2;
+	d += width;
+	*d = s[-2]*0.0625 + s[-1]*.25 + s[0]*.375 + s[1]*.25 + s[2]*.0625;
+      }
+
+      // Last two rows
+      s += 2;
+      d += width;
+      if (dheight*2 <= sheight) {
+	*d = s[-2]*0.0625 + s[-1]*.25 + s[0]*.375 + s[1]*.25 + s[2]*.0625;
+      } else {
+	*d = s[1]*.3125 + s[0]*.3750 + s[-1]*.2500 + s[-2]*.0625;
+      }
+      s += 2;
+      d += width;
+      *d = s[0]*.6875 + s[-1]*.2500 + s[-2]*.0625;
+    }
+  }
+}
+
+// main function
+// takes a double color image and a scaling factor
+// returns resized image
+mxArray *reduce(const mxArray *mxsrc) {
+  double *src = (double *)mxGetPr(mxsrc);
+
+  const int *sdims = mxGetDimensions(mxsrc);
+  if (mxGetNumberOfDimensions(mxsrc) != 3 || 
+      mxGetClassID(mxsrc) != mxDOUBLE_CLASS)
+    mexErrMsgTxt("Invalid input");  
+
+  int ddims[3];
+  ddims[0] = (int)round(sdims[0]*.5);
+  ddims[1] = (int)round(sdims[1]*.5);
+  ddims[2] = sdims[2];
+
+  if (sdims[0] < 5|| sdims[1] < 5)
+    mexErrMsgTxt("Minimum size of image is 5x5");
+
+  mxArray *mxdst = mxCreateNumericArray(3, ddims, mxDOUBLE_CLASS, mxREAL);
+  double *dst = (double *)mxGetPr(mxdst);
+
+  double *tmp = (double *)mxCalloc(ddims[0]*sdims[1]*sdims[2], sizeof(double));
+  reduce1dtran(src, sdims[0], tmp, ddims[0], sdims[1], sdims[2]);
+  reduce1dtran(tmp, sdims[1], dst, ddims[1], ddims[0], sdims[2]);
+  mxFree(tmp);
+
+  return mxdst;
+}
+
+// matlab entry point
+// dst = resize(src, scale)
+// image should be color with double values
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
+  if (nrhs != 1)
+    mexErrMsgTxt("Wrong number of inputs"); 
+  if (nlhs != 1)
+    mexErrMsgTxt("Wrong number of outputs");
+  plhs[0] = reduce(prhs[0]);
+}
+
+/**********
+d = repmat([1:6]',[1 5 3]);
+a = reduce(d);
+ *********/
+

features/reduceGrayScale.cc

@@ -0,0 +1,89 @@
+#include <math.h>
+#include <assert.h>
+#include <string.h>
+#include "mex.h"
+
+// reduce(im) resizes im to half its size, using a 5-tap binomial filter for anti-aliasing
+// (see Burt & Adelson's Laplacian Pyramid paper)
+
+// reduce each column
+// result is transposed, so we can apply it twice for a complete reduction
+void reduce1dtran(double *src, int sheight, double *dst, int dheight, 
+		  int width) {
+  // resize each column 
+  bzero(dst, width*dheight*sizeof(double));
+  int y;
+  double *s, *d;
+
+    for (int x = 0; x < width; x++) {
+      s  = src + x*sheight;
+      d  = dst + x;
+
+      // First row
+      *d = s[0]*.6875 + s[1]*.2500 + s[2]*.0625;      
+
+      for (y = 1; y < dheight-2; y++) {	
+	s += 2;
+	d += width;
+	*d = s[-2]*0.0625 + s[-1]*.25 + s[0]*.375 + s[1]*.25 + s[2]*.0625;
+      }
+
+      // Last two rows
+      s += 2;
+      d += width;
+      if (dheight*2 <= sheight) {
+	*d = s[-2]*0.0625 + s[-1]*.25 + s[0]*.375 + s[1]*.25 + s[2]*.0625;
+      } else {
+	*d = s[1]*.3125 + s[0]*.3750 + s[-1]*.2500 + s[-2]*.0625;
+      }
+      s += 2;
+      d += width;
+      *d = s[0]*.6875 + s[-1]*.2500 + s[-2]*.0625;
+    }
+}
+
+// main function
+// takes a double color image and a scaling factor
+// returns resized image
+mxArray *reduceGrayScale(const mxArray *mxsrc) {
+  double *src = (double *)mxGetPr(mxsrc);
+
+  const int *sdims = mxGetDimensions(mxsrc);
+  if (mxGetNumberOfDimensions(mxsrc) != 2 || 
+      mxGetClassID(mxsrc) != mxDOUBLE_CLASS)
+    mexErrMsgTxt("Invalid input - no gray scale double image given!");  
+
+  int ddims[2];
+  ddims[0] = (int)round(sdims[0]*.5);
+  ddims[1] = (int)round(sdims[1]*.5);
+
+  if (sdims[0] < 5|| sdims[1] < 5)
+    mexErrMsgTxt("Minimum size of image is 5x5");
+
+  mxArray *mxdst = mxCreateNumericArray(2, ddims, mxDOUBLE_CLASS, mxREAL);
+  double *dst = (double *)mxGetPr(mxdst);
+
+  double *tmp = (double *)mxCalloc(ddims[0]*sdims[1], sizeof(double));
+  reduce1dtran(src, sdims[0], tmp, ddims[0], sdims[1]);
+  reduce1dtran(tmp, sdims[1], dst, ddims[1], ddims[0]);
+  mxFree(tmp);
+
+  return mxdst;
+}
+
+// matlab entry point
+// dst = resize(src, scale)
+// image should be gray scale with double values
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
+  if (nrhs != 1)
+    mexErrMsgTxt("Wrong number of inputs"); 
+  if (nlhs != 1)
+    mexErrMsgTxt("Wrong number of outputs");
+  plhs[0] = reduceGrayScale(prhs[0]);
+}
+
+/**********
+d = repmat([1:6]',[1 5 3]);
+a = reduce(d);
+ *********/
+

features/resizeColor.cc

@@ -0,0 +1,121 @@
+#include <math.h>
+#include <assert.h>
+#include <string.h>
+#include "mex.h"
+
+/*
+ * Fast image subsampling.
+ * This is used to construct the feature pyramid.
+ */
+
+// struct used for caching interpolation values
+struct alphainfo {
+  int si, di;
+  double alpha;
+};
+
+// copy src into dst using pre-computed interpolation values
+void alphacopy(double *src, double *dst, struct alphainfo *ofs, int n) {
+  struct alphainfo *end = ofs + n;
+  while (ofs != end) {
+    dst[ofs->di] += ofs->alpha * src[ofs->si];
+    ofs++;
+  }
+}
+
+// resize along each column
+// result is transposed, so we can apply it twice for a complete resize
+void resize1dtran(double *src, int sheight, double *dst, int dheight, 
+		  int width, int chan) {
+  double scale = (double)dheight/(double)sheight;
+  double invscale = (double)sheight/(double)dheight;
+  
+  // we cache the interpolation values since they can be 
+  // shared among different columns
+  int len = (int)ceil(dheight*invscale) + 2*dheight;
+  alphainfo ofs[len];
+  int k = 0;
+  for (int dy = 0; dy < dheight; dy++) {
+    double fsy1 = dy * invscale;
+    double fsy2 = fsy1 + invscale;
+    int sy1 = (int)ceil(fsy1);
+    int sy2 = (int)floor(fsy2);       
+
+    if (sy1 - fsy1 > 1e-3) {
+      assert(k < len);
+      assert(sy-1 >= 0);
+      ofs[k].di = dy*width;
+      ofs[k].si = sy1-1;
+      ofs[k++].alpha = (sy1 - fsy1) * scale;
+    }
+
+    for (int sy = sy1; sy < sy2; sy++) {
+      assert(k < len);
+      assert(sy < sheight);
+      ofs[k].di = dy*width;
+      ofs[k].si = sy;
+      ofs[k++].alpha = scale;
+    }
+
+    if (fsy2 - sy2 > 1e-3) {
+      assert(k < len);
+      assert(sy2 < sheight);
+      ofs[k].di = dy*width;
+      ofs[k].si = sy2;
+      ofs[k++].alpha = (fsy2 - sy2) * scale;
+    }
+  }
+
+  // resize each column of each color channel
+  bzero(dst, chan*width*dheight*sizeof(double));
+  for (int c = 0; c < chan; c++) {
+    for (int x = 0; x < width; x++) {
+      double *s = src + c*width*sheight + x*sheight;
+      double *d = dst + c*width*dheight + x;
+      alphacopy(s, d, ofs, k);
+    }
+  }
+}
+
+// main function
+// takes a double color image and a scaling factor
+// returns resized image
+mxArray *resize(const mxArray *mxsrc, const mxArray *mxscale) {
+  double *src = (double *)mxGetPr(mxsrc);
+  const int *sdims = mxGetDimensions(mxsrc);
+  if (mxGetNumberOfDimensions(mxsrc) != 3 || 
+      mxGetClassID(mxsrc) != mxDOUBLE_CLASS)
+    mexErrMsgTxt("Invalid input");  
+
+  double scale = mxGetScalar(mxscale);
+  if (scale > 1)
+    mexErrMsgTxt("Invalid scaling factor");   
+
+  int ddims[3];
+  ddims[0] = (int)round(sdims[0]*scale);
+  ddims[1] = (int)round(sdims[1]*scale);
+  ddims[2] = sdims[2];
+  mxArray *mxdst = mxCreateNumericArray(3, ddims, mxDOUBLE_CLASS, mxREAL);
+  double *dst = (double *)mxGetPr(mxdst);
+
+  double *tmp = (double *)mxCalloc(ddims[0]*sdims[1]*sdims[2], sizeof(double));
+  resize1dtran(src, sdims[0], tmp, ddims[0], sdims[1], sdims[2]);
+  resize1dtran(tmp, sdims[1], dst, ddims[1], ddims[0], sdims[2]);
+  mxFree(tmp);
+
+  return mxdst;
+}
+
+// matlab entry point
+// dst = resize(src, scale)
+// image should be color with double values
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
+  if (nrhs != 2)
+    mexErrMsgTxt("Wrong number of inputs"); 
+  if (nlhs != 1)
+    mexErrMsgTxt("Wrong number of outputs");
+  plhs[0] = resize(prhs[0], prhs[1]);
+}
+
+
+

features/resizeGrayScale.cc

@@ -0,0 +1,117 @@
+#include <math.h>
+#include <assert.h>
+#include <string.h>
+#include "mex.h"
+
+/*
+ * Fast image subsampling.
+ * This is used to construct the feature pyramid.
+ */
+
+// struct used for caching interpolation values
+struct alphainfo {
+  int si, di;
+  double alpha;
+};
+
+// copy src into dst using pre-computed interpolation values
+void alphacopy(double *src, double *dst, struct alphainfo *ofs, int n) {
+  struct alphainfo *end = ofs + n;
+  while (ofs != end) {
+    dst[ofs->di] += ofs->alpha * src[ofs->si];
+    ofs++;
+  }
+}
+
+// resize along each column
+// result is transposed, so we can apply it twice for a complete resize
+void resize1dtran(double *src, int sheight, double *dst, int dheight, 
+		  int width) {
+  double scale = (double)dheight/(double)sheight;
+  double invscale = (double)sheight/(double)dheight;
+  
+  // we cache the interpolation values since they can be 
+  // shared among different columns
+  int len = (int)ceil(dheight*invscale) + 2*dheight;
+  alphainfo ofs[len];
+  int k = 0;
+  for (int dy = 0; dy < dheight; dy++) {
+    double fsy1 = dy * invscale;
+    double fsy2 = fsy1 + invscale;
+    int sy1 = (int)ceil(fsy1);
+    int sy2 = (int)floor(fsy2);       
+
+    if (sy1 - fsy1 > 1e-3) {
+      assert(k < len);
+      assert(sy-1 >= 0);
+      ofs[k].di = dy*width;
+      ofs[k].si = sy1-1;
+      ofs[k++].alpha = (sy1 - fsy1) * scale;
+    }
+
+    for (int sy = sy1; sy < sy2; sy++) {
+      assert(k < len);
+      assert(sy < sheight);
+      ofs[k].di = dy*width;
+      ofs[k].si = sy;
+      ofs[k++].alpha = scale;
+    }
+
+    if (fsy2 - sy2 > 1e-3) {
+      assert(k < len);
+      assert(sy2 < sheight);
+      ofs[k].di = dy*width;
+      ofs[k].si = sy2;
+      ofs[k++].alpha = (fsy2 - sy2) * scale;
+    }
+  }
+
+  // resize each column 
+  bzero(dst, width*dheight*sizeof(double));
+    for (int x = 0; x < width; x++) {
+      double *s = src + x*sheight;
+      double *d = dst + x;
+      alphacopy(s, d, ofs, k);
+    }
+
+}
+
+// main function
+// takes a double color image and a scaling factor
+// returns resized image
+mxArray *resize(const mxArray *mxsrc, const mxArray *mxscale) {
+  double *src = (double *)mxGetPr(mxsrc);
+  const int *sdims = mxGetDimensions(mxsrc);
+  if (mxGetNumberOfDimensions(mxsrc) != 2 || 
+      mxGetClassID(mxsrc) != mxDOUBLE_CLASS)
+    mexErrMsgTxt("Invalid input - double gray scale image expected!");  
+
+  double scale = mxGetScalar(mxscale);
+  if (scale > 1)
+    mexErrMsgTxt("Invalid scaling factor");   
+
+  int ddims[2];
+  ddims[0] = (int)round(sdims[0]*scale);
+  ddims[1] = (int)round(sdims[1]*scale);
+  mxArray *mxdst = mxCreateNumericArray(2, ddims, mxDOUBLE_CLASS, mxREAL);
+  double *dst = (double *)mxGetPr(mxdst);
+
+  double *tmp = (double *)mxCalloc(ddims[0]*sdims[1], sizeof(double));
+  resize1dtran(src, sdims[0], tmp, ddims[0], sdims[1]);
+  resize1dtran(tmp, sdims[1], dst, ddims[1], ddims[0]);
+  mxFree(tmp);
+
+  return mxdst;
+}
+
+// matlab entry point
+// dst = resize(src, scale)
+// image should be color with double values
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { 
+  if (nrhs != 2)
+    mexErrMsgTxt("Wrong number of inputs"); 
+  if (nlhs != 1)
+    mexErrMsgTxt("Wrong number of outputs");
+  plhs[0] = resize(prhs[0], prhs[1]);
+}
+

features/visualization/compareHogsVisually.m

@@ -0,0 +1,37 @@
+function compareHogsVisually ( hog1, hog2, idxDimsToShow )
+% function compareHogsVisually ( hog1, hog2, idxStart, idxEnd )
+% 
+% brief:  compute dim-wise differences between two hog features and plot
+%         the results into an image
+% author: Alexander Freytag
+% date:   05-02-2014 (dd-mm-yyyy)
+%
+% INPUT:  hog1 and hog2 two hog features of corresponding size
+%         idxDimsToShow (optionally) vector with dim indices
+
+    if ( nargin < 3 )
+        idxDimsToShow = 1:size(hog1,3);
+    end
+    
+    
+    for i=1:length(idxDimsToShow)
+        
+        % create new figure
+        fig1=figure; 
+        % set title indicating current dimension
+        s_titleHoG = sprintf('Dim %d', idxDimsToShow(i) );
+        set ( fig1, 'name', s_titleHoG);        
+        
+        % plot actual differences
+        imagesc( (hog1(:,:,idxDimsToShow(i) ) - hog2(:,:, idxDimsToShow(i) )));
+        
+        % don't miss a colorbar indicating the actual values
+        colorbar; 
+
+        % wait for user input 
+        pause; 
+
+        % close current figure and switch to next dimension
+        close(fig1);
+    end
+end

features/visualization/computeHoGSizeNormalized.m

@@ -0,0 +1,23 @@
+function hogFeature = computeHoGSizeNormalized ( imgfn )
+
+    img = readImage( imgfn );
+    [xsize,ysize] = size(img);
+
+    load('bg15Scenes.mat');
+
+    pos.x1=1;
+    pos.y1=1;
+    pos.x2=xsize;
+    pos.y2=ysize;
+    pos.im = imgfn;
+    
+    modeltemplate = initmodel('name', pos, bg15ScenesGray); 
+
+
+    warpedImgRegular = warppos('doesntMatter', modeltemplate, pos);
+%     imshow(warpedImgRegular{1})
+
+    % compute HoG feature from the given image (size normalized)
+    hogFeature = computeHOGs_WHO( img );
+
+end

features/visualization/myHOGpicture.m

@@ -0,0 +1,28 @@
+function im = myHOGpicture(w, bsx, bsy )
+    % myHOGpicture(w, bsx, bxy )
+    % Make picture of positive HOG weights.
+
+    % construct a "glyph" for each orientaion
+    bim1 = zeros(bsy, bsx);
+    bim1(:,round(bsx/2):round(bsx/2)+1) = 1;
+    no   = 9;
+    bim  = zeros([size(bim1) no]);
+    bim(:,:,1) = bim1;
+    for i = 2:no,
+      bim(:,:,i) = imrotate(bim1, -(i-1)*(180/no), 'crop');
+    end
+
+    % make pictures of positive weights bs adding up weighted glyphs
+    s = size(w);    
+    w(w < 0) = 0;    
+    im = zeros(bsy*s(1), bsx*s(2));
+    for i = 1:s(1),
+      iis = (i-1)*bsy+1:i*bsy;
+      for j = 1:s(2),
+        jjs = (j-1)*bsx+1:j*bsx;          
+        for k = 1:no,
+          im(iis,jjs) = im(iis,jjs) + bim(:,:,k) * w(i,j,k);
+        end
+      end
+    end
+end

features/visualization/showWeightVectorHOG.m

@@ -0,0 +1,108 @@
+function out = showWeightVectorHOG( w, settings )
+% function out = showWeightVectorHOG( w, settings )
+% 
+% author: Alexander Freytag
+% date  : 27-02-2014 (dd-mm-yyyy)
+% 
+% BRIEF :
+%   Given a weight vector w obtained by training a model with DPM HOG features, 
+%   positive and negative components are displayed separately
+% 
+% INPUT :
+%    w       --  weight vector of model
+%    settings
+%            --  (optional), struct with possible fields, e.g.,
+%                'b_closeImg', ...
+% 
+% OUTPUT :
+%    out     -- (optional), the resulting image of visualized model
+
+    %% ( 0 ) check input
+       
+    if ( nargin  < 2 )
+        settings = [];
+    end
+    
+    b_closeImg    = getFieldWithDefault ( settings, 'b_closeImg',   false);
+    b_removeAxis  = getFieldWithDefault ( settings, 'b_removeAxis', true);
+    s_destination = getFieldWithDefault ( settings, 's_destination', '');
+        
+    if ( isempty (s_destination) ) 
+        b_saveAsEPS = false;
+    else
+        b_saveAsEPS = true;
+    end
+    
+    
+    widthOfCell  = getFieldWithDefault( settings, 'widthOfCell',  20 ); 
+    heightOfCell = getFieldWithDefault( settings, 'heightOfCell', 20 );
+ 
+
+    %% ( 1 ) Make pictures of positive and negative weights        
+
+    
+    % compute the visualization
+    imPos = myHOGpicture( w,  widthOfCell, heightOfCell );
+    imNeg = myHOGpicture( -w, widthOfCell, heightOfCell );
+    
+    scale = max( [w(:); -w(:)] );
+    imPos = imPos ./ scale;
+    imNeg = imNeg ./ scale;    
+    
+    
+    
+    %% ( 2 ) Put pictures together
+    % 
+    % a bit of padding for nice visualization
+    buff = 10;
+    imPos = myPadArray( imPos, [buff buff], 0.5);
+    imNeg = myPadArray( imNeg, [buff buff], 0.5 );
+    
+    im = [imPos imNeg];
+       
+
+    %% ( 3 ) saturate image information out of [0,1]
+    
+    im(im < 0) = 0;
+    im(im > 1) = 1;
+    
+    % scale to [0,255]
+    im = im*255;    
+    
+    % convert to uint8 as usually done for images
+    im = uint8(im);
+    
+  
+    
+    %% ( 4 ) draw figure or output result
+    if ( nargout == 0 )    
+    
+        % create new figure
+        figHandle = figure;
+
+        imagesc(im);
+        colormap gray;
+
+        % make images beeing displayed correctly, i.e., not skewed
+        axis image;
+        %don't show axis ticks
+        set(gca,'Visible','off');    
+
+        if ( b_removeAxis )
+            set(gca,'Visible','off')
+        end
+
+        if ( b_closeImg )
+            pause
+            close ( figHandle );
+        end    
+    else    
+      out = im;
+    end        
+
+    
+    if ( b_saveAsEPS )
+        print( '-depsc', '-r300' , s_destination );
+    end
+    
+end

features/visualization/showWeightVectorHOGPosOnly.m

@@ -0,0 +1,100 @@
+function out = showWeightVectorHOGPosOnly( w, settings )
+% function out = showWeightVectorHOGPosOnly( w, settings )
+% 
+% author: Alexander Freytag
+% date  : 14-02-2014 (dd-mm-yyyy)
+% 
+% BRIEF :
+%   Given a weight vector w obtained by training a model with DPM HOG features, 
+%   positive components are displayed
+% 
+% INPUT :
+%    w       --  weight vector of model
+%    settings
+%            --  (optional), struct with possible fields, e.g.,
+%                'b_closeImg', ...
+% 
+% OUTPUT :
+%    out     -- (optional), the resulting image of visualized model
+
+    %% ( 0 ) check input
+
+    if ( nargin  < 2 )
+        settings = [];
+    end
+    
+    b_closeImg    = getFieldWithDefault ( settings, 'b_closeImg',   false);
+    b_removeAxis  = getFieldWithDefault ( settings, 'b_removeAxis', true);
+    s_destination = getFieldWithDefault ( settings, 's_destination', '');
+        
+    if ( isempty (s_destination) ) 
+        b_saveAsEPS = false;
+    else
+        b_saveAsEPS = true;
+    end
+    
+    
+    widthOfCell  = getFieldWithDefault( settings, 'widthOfCell',  20 ); 
+    heightOfCell = getFieldWithDefault( settings, 'heightOfCell', 20 );
+ 
+
+    %% ( 1 ) Make pictures of positive weights        
+    
+    wwp = foldHOG(w);
+    scale = max(wwp(:));
+    scale = double(255)/scale;
+    
+
+    % compute the visualization
+    im = myHOGpicture(wwp, widthOfCell, heightOfCell );
+    % scale to [0,255]
+    im = im*scale;
+
+    % 
+    % a bit of padding for nice visualization
+    buff = 10;
+    im = myPadArray( im, [buff buff], 200 );
+
+    % convert to uint8 as usually done for images
+    im = uint8(im);    
+    
+    %% ( 4 ) draw figure or output result
+    if ( nargout == 0 )    
+    
+        % create new figure
+        figHandle = figure;
+
+        imagesc(im);
+        colormap gray;
+
+        % make images beeing displayed correctly, i.e., not skewed
+        axis image;
+        %don't show axis ticks
+        set(gca,'Visible','off');    
+
+        if ( b_removeAxis )
+            set(gca,'Visible','off')
+        end
+
+        if ( b_closeImg )
+            pause
+            close ( figHandle );
+        end    
+    else
+      out = im;
+    end        
+
+    
+    if ( b_saveAsEPS )
+        print( '-depsc', '-r300' , s_destination );
+    end
+    
+end
+
+function f = foldHOG(w)
+    % f = foldHOG(w)
+    % Condense HOG features into one orientation histogram.
+    % Used for displaying a feature.
+
+    f=max(w(:,:,1:9),0)+max(w(:,:,10:18),0)+max(w(:,:,19:27),0);
+end

initWorkspaceWHOGeneric.m

@@ -0,0 +1,39 @@
+function initWorkspaceWHOGeneric
+
+
+    %% add paths
+    
+    % pre-computed data
+    addpath( genpath( fullfile(pwd, 'data') ) );
+    
+    % things for detecting objects
+    addpath( genpath( fullfile(pwd, 'detect') ) );
+    
+    % mainly HOG and scripts for image resizing
+    addpath( genpath( fullfile(pwd, 'features') ) );
+    
+    % everything regarding learning the actual LDA detectors
+    addpath( genpath( fullfile(pwd, 'learn') ) );
+    
+    % path to introducing demo scripts
+    addpath( genpath( fullfile(pwd, 'demos') ) );    
+    
+    % minor things, e.g., warping of blocks or padding of arrays
+    addpath( 'misc' );    
+    
+    
+    % read images from given filename
+    if (exist('readImage') ~= 2 )
+        addpath( genpath( fullfile(pwd, 'misc/imageIO') ) );    
+    end
+    
+    % scripts for configuration of setting-objects
+    if (exist('getFieldWithDefault') ~= 2 )
+        addpath( genpath( fullfile(pwd, 'misc/settings') ) );    
+    end    
+    
+    %% 3rd party, untouched
+    
+    % nothing needed here...
+    
+end

learn/initLDAStuff.m

@@ -0,0 +1,49 @@
+function ldaStuff = initLDAStuff( ny, nx, modeltemplate)
+% function ldaStuff = initLDAStuff( ny, nx, modeltemplate)
+%
+% BRIEF:
+%    At the beginning, every detector is built only from a single
+%    block (see the 1-SVM paper for motivation details)
+%
+%    Lateron, further blocks are added to this struct increasing the
+%    number of training samples of every detector
+%
+% INPUT:
+%    ny             -- ysize of our HoG features, assumed to be constant
+%    nx             -- xsize of our HoG features, assumed to be constant
+%    modeltemplate  -- template for best LDA model (size) determined by who code
+%
+% OUTPUT:
+%    ldaStuff   -- struct containing following fields:
+%                   ldaStuff.R             -- covariance matrix of all data
+%                   ldaStuff.neg           -- mean of universal negative data
+%                   ldaStuff.modelTemplate -- model template for LDA models
+
+
+    %this needs to be done only ones since our parameters are always the same
+
+    model = modeltemplate;
+    
+    if ( isfield(modeltemplate.lda,'b_noiseDropOut') && ~isempty(modeltemplate.lda.b_noiseDropOut) )
+        model.lda.b_noiseDropOut = modeltemplate.lda.b_noiseDropOut;
+    end
+  
+    if ( isfield(modeltemplate.lda,'d_dropOutProb') && ~isempty(modeltemplate.lda.d_dropOutProb) )
+        model.lda.d_dropOutProb = modeltemplate.lda.d_dropOutProb;
+    end    
+  
+    if ( isfield(modeltemplate.lda,'lambda') && ~isempty(modeltemplate.lda.lambda) )
+        model.lda.lambda = modeltemplate.lda.lambda;
+    end      
+    
+    [ldaStuff.R,ldaStuff.neg] = whitenWithDropout(model.bg, model.lda, nx,ny);
+    
+    ldaStuff.modelTemplate = modeltemplate;
+    
+    %per default, we work on all dimensions
+    %however, DPM HOG features have as last feature dimension truncation
+    %features constant to 0. In thoses cases, we set the flag to true and
+    %ignore the last dimension in further processing stages
+    ldaStuff.b_ignoreLastDim = false;
+        
+end

learn/initmodel_static.m

@@ -0,0 +1,60 @@
+function model = initmodel_static( settings, i_numDim )
+% function model = initmodel_static( settings, i_numDim )
+% 
+% BRIEF
+%    Initialize model structure.
+% 
+% OUTPUT
+%    model = initmodel( settings )
+%        model.maxsize = [y,x] size of root filter in HOG cells
+% ...TODO
+% 
+% author: Alexander Freytag
+% date:   13-02-2014 (dd-mm-yyyy) (last updated)
+
+    
+%     %% how many dimensions does our feature has?
+    
+    % for DPM-HOG features, this would result in the following layout
+    % how many dimensions does our resulting 'augmented HoG feature' has?
+    % see DPM Paper from 2010 for more details
+    % 2*numberBins for keeping gradient sign + 
+    % 1*numberBins without the sign + 
+    % 4 strange texture feature dimensions + 
+    % 1 dimension constant to zero
+    % -> default: 32 dimensions    
+    
+    sizeOfModel = [ settings.lda.bg.i_numCells, ...
+                    i_numDim...
+                  ];
+
+    %% initialize the rest of the model structure
+    
+    %empty model of according size
+    model.w          = zeros(sizeOfModel);
+    % size of root filter in HOG cells
+    model.i_numCells = sizeOfModel(1:2);
+    % size of each cell in pixel
+    model.i_binSize  = settings.lda.bg.i_binSize;
+    % strange interval 
+    model.interval   = settings.lda.bg.interval;    
+    %threshold to reject detection with score lwoer than that
+    model.d_detectionThreshold ...
+                     = settings.lda.d_detectionThreshold;
+
+    %negative mean, cov matrix, and stuff like that
+    model.bg         = settings.lda.bg;
+    
+    %======== ======== ======== ======== 
+    %      add here noise model for
+    %     modeling de-noising effect
+    %======== ======== ======== ========  
+    
+    % this adds noise on the main diagonal of the covariance matrix
+    model.lda.lambda         = settings.lda.lambda;
+    
+    %%% this additionally adds a drop-out noise model 
+    model.lda.b_noiseDropOut = settings.lda.b_noiseDropOut;
+    model.lda.d_dropOutProb  = settings.lda.d_dropOutProb;
+    
+end

learn/learnWithGivenWhitening.m

@@ -0,0 +1,89 @@
+function modelNew = learnWithGivenWhitening(model,R, neg, pyraFeats, i_truncDim)
+% function modelNew = learnWithGivenWhitening(model,R, neg, pyraFeats, i_truncDim)
+% BRIEF:
+%    Learn model by linear discriminant analysis with already given
+%    negative mean and covariance matrix
+% 
+% INPUT:
+%    model           -- struct, previously untrained model, contains at
+%                       least the following fields:
+%        .w            -- previous (most likely empty) weight vector of
+%                         model, relevant for determining the correct size
+%        .i_numCells   -- number of cells, specified per dimension (copied
+%                         only)
+%        .i_binSize    -- number of pixel of every cell in both directions
+%                         (copied only)
+%        .interval     -- how many octaves are used for feature extraction 
+%                         (copied only)
+%        .d_detectionThreshold -- minimum scores for accepting a response
+% 
+%    R               -- covariance matrix learned previously
+%    neg             -- negative mean learned previously
+%    pyraFeats       -- features of positive examples
+%    i_truncDim      -- int, indicating which dimension, if any,  serves as 
+%                       truncation feature by being constant to zero ( as
+%                       done for DPM HOG features )
+% 
+% OUTPUT:
+%   modelNew         -- struct, learned model, with fields 'w',
+%                       'i_numCells', 'i_binSize', 'interval', 'd_detectionThreshold
+% 
+% author:               Alexander Freytag
+% last time modified:   27-02-2014 (dd-mm-yyyy)
+
+
+
+    [ny,nx,nf] = size(model.w);
+    
+    if ( i_truncDim > 0 )
+        nf = nf - 1;
+    end
+
+
+    %num is the number of blocks we have for this cache
+    num  = length(pyraFeats);
+    feats = zeros(ny*nx*nf,num);
+
+    % flatten features into single vectors
+    for i = 1:num
+        % get current feature
+        feat = pyraFeats(i).feature;
+      
+        % possibly remove unneeded truncation feature
+        if ( i_truncDim > 0 )
+            feat = feat(:, :, 1:end~=i_truncDim );
+        end
+
+        % flatten vector and store in feats-struct
+        feats(:,i) = feat(:);
+    end
+
+    % average all features together
+    pos = mean(feats,2);
+
+    % perform actual LDA training
+    w=R\(R'\(pos-neg));
+    
+    % bring weight vector into correct layout
+    w = reshape(w,[ny nx nf]);
+
+    if ( i_truncDim > 0 )
+        % Add in occlusion feature
+        %note: might only be troublesome if very first dim is the td...
+        w = cat( 3, w(:,:,1:(i_truncDim-1) ), ...
+                    zeros ( size(w,1),size(w,2),1,class(w) ), ...
+                    w(:, :, i_truncDim:end) ); 
+    end
+    
+    modelNew.w                    = w;
+    % size of root filter in HOG cells
+    modelNew.i_numCells           = model.i_numCells;
+    % size of each cell in pixel
+    modelNew.i_binSize            = model.i_binSize;    
+    % strange interval
+    modelNew.interval             = model.interval;
+    %threshold to reject detection with score lower than that
+    modelNew.d_detectionThreshold = model.d_detectionThreshold;   
+    
+end
+

learn/learn_dataset.m

@@ -0,0 +1,98 @@
+function model = learn_dataset( pos, neg, bg, settings )
+% model = learn_dataset(pos, neg, name, settings )
+% 
+% author: Alexander Freytag
+% date:   13-02-2014 (dd-mm-yyyy) (last updated)
+% 
+% BRIEF
+%    Initialize model structure. (note: If you just have image patches, instead of bounding boxes, consider using learn.m directly)
+% 
+% INPUT
+%     pos      -- is a struct array for the positive patches, with fields:
+%	              .im (full path to the image), .x1 (xmin), .y1 (ymin), .x2 (xmax), .y2 (ymax)
+%     neg      -- struct array for the negative patches with field:
+%	              im: full path to the image 
+%                 Used only when the background statistics cannot be found.
+%     bg       -- (optional) pre-computed whitening info (cov mat, negative mean, ...)
+%     settings -- struct with config settings
+%      
+% OUTPUT
+%    model     -- learnt model
+% 
+
+
+    if ( isempty(neg) )
+        allImages = {pos.im};
+    else
+        allImages = [ {pos.im}; {neg.im} ];
+    end
+
+    
+    %use given background statistics if they exist; else build them
+    if ( isempty(bg) )
+      bg  = trainBGwithArbitraryFeatures( allImages, settings );
+    end
+
+
+    % Define model structure
+    % how many dimensions does our feature has?
+    % compute a feature for a small,  empty image to fetch the number of
+    % dimensions every cell-feature has
+    i_numImgChannels = size ( readImage( allImages{1}),3);
+    i_numDim = size( settings.fh_featureExtractor.mfunction ( zeros([3 3 i_numImgChannels]) ),3 );  
+    
+    clear ( 'allImages' );
+    
+    %computed lda variables
+    settings.lda.bg = bg;
+    %threshold to reject detection with score lwoer than th
+    settings.lda.d_detectionThreshold = 0;
+    
+    % how many cells are useful in both directions?
+    i_numCells = computeOptimalCellNumber ( pos, settings.i_binSize );
+    settings.lda.bg.i_numCells = i_numCells;
+    
+    %======== ======== ======== ======== 
+    %      add here noise model for
+    %     modeling de-noising effect
+    %======== ======== ======== ========      
+    % this adds noise on the main diagonal of the covariance matrix
+    settings.lda.lambda = bg.lambda;
+    %%% this additionally adds a drop-out noise model 
+    settings.lda.b_noiseDropOut = false;
+    settings.lda.d_dropOutProb = 0.0;
+    
+    settings.lda.bg.interval = settings.interval;
+    
+    model = initmodel_static(settings, i_numDim);
+    
+    %skip models if the HOG window is too skewed
+    if( max(model.i_numCells)<4*min(model.i_numCells) )
+        
+        %get image patches
+        warpedTrainBlocks = warpBlocksToStandardSize( model, pos, settings.fh_featureExtractor );
+        
+        % pre-compute features from the seeding blocks (size normalized)
+        feats = computeFeaturesForBlocks( warpedTrainBlocks, settings);  
+
+        
+        [ldaStuff.R,ldaStuff.neg] = whitenWithDropout(model.bg, model.lda, model.i_numCells(2),model.i_numCells(1));
+        
+        % for the HOG features computed here, the 32nd dim is constant to
+        % zero serving as truncation dim.
+        i_truncDim = 32;
+        
+        model = learnWithGivenWhitening(model, ...
+            ldaStuff.R, ldaStuff.neg, ...
+            feats, i_truncDim );
+
+    else
+        model.maxsize
+        error('HOG window is too skewed!');
+    end
+    
+    model.w=model.w./(norm(model.w(:))+eps);
+    model.thresh = 0.5;
+    model.bg=[];
+end
+

learn/learn_dataset_fixedmodel.m

@@ -0,0 +1,44 @@
+function model = learn_dataset_fixedmodel(pos, neg, name, model)
+
+%TODO
+
+% Load background statistics if they exist; else build them
+file = bg_file_name;
+try
+  load(file);
+catch
+  all = rmfield(pos,{'x1','y1','x2','y2'});
+  all = [all neg];
+  bg  = trainBG(all,20,5,8);
+  save(file,'bg');
+end
+bg
+
+% Define model structure
+%model = initmodel(name,pos,bg);
+model.bg=bg;
+%skip models if the HOG window is too skewed
+if(max(model.maxsize)<4*min(model.maxsize))
+
+
+warped=warppos(name, model, pos);
+
+%flip if necessary
+if(isfield(pos, 'flipped'))
+    fprintf('Warning: contains flipped images. Flipping\n');
+    for k=1:numel(warped)
+        if(pos(k).flipped)
+            warped{k}=warped{k}(:,end:-1:1,:);
+        end
+    end
+end
+
+
+% Learn by linear discriminant analysis
+model = learn(name,model,warped);
+end
+model.w=model.w./norm(model.w(:));
+model.bg=[];
+model.thresh = 0.5;
+model.name=name;
+

learn/trainBGwithArbitraryFeatures.m

@@ -0,0 +1,246 @@
+function BG = trainBGwithArbitraryFeatures( allImages, settings )
+% function BG = trainBGwithArbitraryFeatures( allImages, settings )
+% 
+% BRIEF
+%    Trains a spatial autocorrelation function from a list of images
+% 
+% OUTPUT
+% ...TODO
+% 
+% author: Alexander Freytag
+% date:   13-02-2014 (dd-mm-yyyy)
+
+    %% (1) check input
+    
+    order               = settings.order;
+    
+    interval            = settings.interval;
+    
+    % number of cells in every dimension
+    i_binSize          = settings.i_binSize;
+    
+    % those three assignments need to be done here already, since the BG
+    % object is passed to pyramid construction methods
+    % (featPyramidGeneric), e.g., in line 84 (currently)
+    BG.i_binSize  = i_binSize;
+    BG.interval   = interval;  
+    % size of root filter in cells
+    % we want to extract at least two cells per dim to ensure save
+    % calcaluations (strange things happened if set to [1,1] only)
+    BG.i_numCells = [2,2];      
+    
+    i_truncDim = getFieldWithDefault ( settings, 'i_truncDim', -1 );
+    
+    
+    
+    %NOTE perhaps we should explicitely specify to NOT use any padding at
+    %all... 
+    
+
+    %% (2) set necessary variables
+    % Ignoring an empty feature dimension? For DPM HOG features, this is 
+    % an empty (===0) truncation feature
+    
+    % if no dim was specified, we compare against the feature types we
+    % know...
+    if ( i_truncDim < 0 )
+        if ( strcmp ( settings.fh_featureExtractor.name, 'Compute HOG features using WHO code' ) || ...
+             strcmp ( settings.fh_featureExtractor.name, 'HOG and Patch Means concatenated' )    || ...
+             strcmp ( settings.fh_featureExtractor.name, 'HOG and Color Names' )...
+           ) 
+            i_truncDim = 32;
+        else
+            i_truncDim = 0;
+        end
+    end
+
+    % compute a feature for a small, empty image to fetch the number of
+    % dimensions every cell-feature has
+    i_numImgChannels = size ( readImage(allImages{1}),3);
+    i_numDim = size( settings.fh_featureExtractor.mfunction ( zeros([3 3 i_numImgChannels]) ),3 );    
+    
+    if ( i_truncDim > 0 )
+      display('Ignoring last truncation feature');
+      i_numDim  = i_numDim-1;
+    end
+    
+    neg = zeros(i_numDim,1);
+    
+    % will be the total number of cells extracted from all images on all
+    % scales
+    n   = 0;
+    
+    
+
+
+    %% (3) start learning negative mean
+    fprintf('\nLearning negative mean\n');
+
+    % average features over all images, and all possible locations and
+    % scales
+    for i = 1:length(allImages)
+        
+      % progressbar-like output  
+      if(rem(i,100)==0)
+          fprintf('%d / %d\n', i,length(allImages) );
+      end
+   
+      im = readImage( allImages{i} );
+      if ( i_numImgChannels ~= size ( im,3) )
+            % if by chance there are some images not fitting to the other
+            % ones...
+            continue;
+      end      
+
+      % Extract feature pryamid, removing hallucinated octave
+      pyra = featPyramidGeneric(im, BG, settings );
+            
+      % run over all scales
+      for s = 1:length(pyra.feat)
+        featIm = pyra.feat{s};
+        
+        % possibly remove last dimension 
+        if ( i_truncDim > 0 )
+            featIm = featIm(:, :, 1:end~=i_truncDim );
+        end        
+   
+        [imy,imx,imz] = size(featIm);
+        
+        % total number of extractable from the image at the current scale
+        t = imy*imx;
+        
+        % compress all cells in each channel into a single long feature
+        % vector
+        feat = reshape(featIm,t,imz);
+        % increase number of totally inspected features
+        n    = n + t;
+        % add features to previous universal mean
+        try
+            neg  = neg + sum(feat)';        
+        catch err
+            err
+        end
+      end  
+    end
+
+    % normalize mean accordingly
+    neg = neg'/n;
+
+    w    = order;
+    h    = order;
+    dxy = [];
+    for x = 0:w-1,
+      for y = 0:h-1,
+        dxy = [dxy; [x y]];
+        if x > 0 & y > 0,
+          dxy = [dxy; [x -y]];
+        end
+      end
+    end
+
+
+
+    %% (4) start learning covariance matrix
+    k    = size(dxy,1);
+    ns   = zeros(k,1);
+    cov  = zeros(i_numDim,i_numDim,k);    
+    
+    fprintf('\nLearning stationairy negative covariance\n');
+
+    for i = 1:length(allImages)
+                
+        % progressbar-like output
+        if( rem(i,100)==0 )
+            fprintf('%d / %d\n', i,length(allImages) );
+        end
+        im        = readImage( allImages{i} );
+        if ( i_numImgChannels ~= size ( im,3) )
+            % if by chance there are some images not fitting to the other
+            % ones...
+            continue;
+        end
+
+        % Extract feature pryamid
+        pyra      = featPyramidGeneric(im, BG, settings );
+      
+        % Subtract mean from all features extracted from current image
+        for s = 1:length(pyra.feat),
+            featIm = pyra.feat{s};
+
+            if ( i_truncDim > 0 )
+                featIm = featIm(:, :, 1:end~=i_truncDim );
+            end
+        
+            [imy,imx,imz] = size(featIm);
+            featIm = reshape(featIm,imy*imx,imz);
+            featIm = bsxfun(@minus,featIm,neg);
+            if ( imz > 1 )
+                featIm = reshape(featIm,[imy imx imz]);
+            else
+                featIm = reshape(featIm,[imy imx]);
+            end
+            pyra.feat{s} = featIm;
+        end
+      
+        for s = 1:length(pyra.feat),
+          for i = 1:k,
+            dx = dxy(i,1);
+            dy = dxy(i,2);
+            [imy,imx,~] = size(pyra.feat{s});
+            
+            if ( dy > 0 )
+              y11 = 1;
+              y12 = imy - dy;
+            else
+              y11 = -dy + 1;
+              y12 = imy;
+            end
+            
+            if ( dx > 0 )
+              x11 = 1;
+              x12 = imx - dx;
+            else
+              x11 = -dx + 1;
+              x12 = imx;
+            end
+            
+            if ( ( y12 < y11 ) || ( x12 < x11 ) )
+              continue;
+            end
+            
+            y21 = y11 + dy;
+            y22 = y12 + dy;
+            x21 = x11 + dx;
+            x22 = x12 + dx;        
+            assert(y11 >= 1 && y12 <= imy && ...
+                   x21 >= 1 && x22 <= imx);
+               
+            t    = (y12 - y11 + 1)*(x12 - x11 + 1);
+            feat1 = reshape(pyra.feat{s}(y11:y12,x11:x12,:),t,i_numDim);        
+            feat2 = reshape(pyra.feat{s}(y21:y22,x21:x22,:),t,i_numDim);
+            
+            cov(:,:,i) = cov(:,:,i) + feat1'*feat2;
+            ns(i) = ns(i) + t;
+          end
+        end
+      
+    end %for i = 1:length(allImages) 
+
+    fprintf('\n');
+    neg = neg';
+
+    for i = 1:k
+      cov(:,:,i) = cov(:,:,i) / ns(i);
+    end
+
+    
+    %% (5) write results to output object
+    
+    BG.neg        = neg;
+    BG.cov        = cov;
+    BG.dxy        = dxy;
+    BG.ns         = ns;
+    BG.lambda     = .01;
+  
+    
+end

learn/whitenWithDropout.m

@@ -0,0 +1,93 @@
+% function [R,neg] = whitenWithDropout(bg,lda,nx,ny)
+% Obtain whitenixng matrix and mean from a general HOG model   
+% by a cholesky decompoition on a stationairy covariance matrix 
+% feat' = R\(feat - neg) has zero mean and unit covariance
+%
+% bg.neg: negative mean (nf by 1)
+% bg.cov: covariance for k spatial offsets (nf by nf by k)
+% bg.dxy: k spatial offsets (k by 2)
+% lda.lambda: regularizer
+function [R,neg] = whitenWithDropout(bg,lda,nx,ny)
+  neg = repmat(bg.neg',ny*nx,1);
+  neg = neg(:);
+  p=1;
+  
+  
+%   if ( ~isfield(lda,'b_noiseDropOut') || isempty(lda.b_noiseDropOut) )
+%       b_noiseDropOut = true;
+%   else
+%       b_noiseDropOut = lda.b_noiseDropOut;
+%   end
+%   
+%   if ( ~isfield(lda,'d_dropOutProb') || isempty(lda.d_dropOutProb) )
+%       d_dropOutProb = 0.1;
+%   else
+%       d_dropOutProb = lda.d_dropOutProb;
+%   end
+  
+
+  
+  while(p~=0)
+      sig = reconstructSig(nx,ny,bg.cov,bg.dxy);
+
+      % drop-out like noise model, as described by Chen et al. (Marginalized Denoising Autoencoders for Domain Adaptation)
+      if ( lda.b_noiseDropOut )
+          d=size(sig,1);
+          q=ones(d,1) .* (1-lda.d_dropOutProb);
+          sig=sig.*(q*q');
+      end
+
+      % Gaussian noise model for every dimension
+      sig = sig + lda.lambda*eye(size(sig));  
+
+      [R,p] = chol(sig);
+      if p ~= 0,
+        disp('Increasing lambda');
+        lda.lambda = lda.lambda+0.01;
+        %display('Sig is not positive definite, add a larger regularizer');
+        %keyboard;
+      end
+  end
+end
+  
+function w = reconstructSig(nx,ny,ww,dxy)
+% W = reconstructSig(nx,ny,ww,dxy)
+% W = n x n 
+% n = ny * nx * nf
+
+  k  = size(dxy,1);
+  nf = size(ww,1);
+  n  = ny*nx;  
+  w  = zeros(nf,nf,n,n);
+  
+  for x1 = 1:nx,
+    for y1 = 1:ny,
+      i1 = (x1-1)*ny + y1;
+      for i = 1:k,
+        x = dxy(i,1);
+        y = dxy(i,2);
+        x2 = x1 + x;        
+        y2 = y1 + y;
+        if x2 >= 1 && x2 <= nx && y2 >= 1 && y2 <= ny,
+          i2 = (x2-1)*ny + y2;
+          w(:,:,i1,i2) = ww(:,:,i); 
+        end
+        x2 = x1 - x;        
+        y2 = y1 - y;
+        if x2 >= 1 && x2 <= nx && y2 >= 1 && y2 <= ny,
+          i2 = (x2-1)*ny + y2; 
+          w(:,:,i1,i2) = ww(:,:,i)'; 
+        end
+      end
+    end
+  end
+  
+  % Permute [nf nf n n] to [n nf n nf]
+  w = permute(w,[3 1 4 2]);
+  w = reshape(w,n*nf,n*nf);
+  
+  % Make sure returned matrix is close to symmetric
+  assert(sum(sum(abs(w - w'))) < 1e-5);
+  
+  w = (w+w')/2;
+end

misc/imageIO/readImage.m

@@ -0,0 +1,44 @@
+function im = readImage( s_imageName, settings )
+% function im = readImage( s_imageName, settings )
+% 
+% author: Alexander Freytag
+% date:   31-03-2014 ( dd-mm-yyyy )
+% 
+% BRIEF:  read an image from filename, possible repeate a gray value image
+%         to result in 3 dimensions (pseudo-RGB)
+% 
+% INPUT: 
+%         s_imageName -- char array ( filename to image )
+%         settings    -- struct (optional), with useable fields
+%                        'b_resizeImageToStandardSize',
+%                        'i_standardImageSize'
+% 
+% OUTPUT:
+%         im          -- h x w x 3 uint8 image
+% 
+    if( nargin < 2 ) 
+        settings = [];
+    end
+
+    %% (1) READ IMAGE FROM FILENAME
+    im = imread( s_imageName );
+    
+    %% (2) RESIZING TO STANDARD SIZE IF DESIRED
+    if ( getFieldWithDefault (settings, 'b_resizeImageToStandardSize', false ) )
+        
+        i_standardImageSize = getFieldWithDefault ( settings, 'i_standardImageSize', [128,128]);
+        
+        if ( ndims(i_standardImageSize) == 1)
+            i_standardImageSize = repmat(i_standardImageSize, [1,2]);
+        end
+        
+       im = imresize(im, i_standardImageSize);  
+    end
+
+    %% (3) MAKE IMAGE ALWAYS TO BE COLOR, OR PSEUDO COLOR
+    if( ndims( im ) == 2 )
+        im = repmat ( im, [1,1,3] );
+    end
+    
+end
+

misc/myPadArray.m

@@ -0,0 +1,53 @@
+function paddedArray = myPadArray( myArray, padSize, valForPadding )
+% function newf = myPadArray( myArray, padSize, valForPadding )
+% 
+% BRIEF:
+%   Quite the same thing as matlabs padarray, but significantly faster ( about 
+%    ~10 times)
+% 
+% author: Alexander Freytag
+% date  : 28-02-2014 ( dd-mm-yyyy )
+
+    if ( numel ( padSize ) < 2 )
+        padSize = [padSize, padSize];
+    end
+    
+    if ( (ndims ( myArray ) > 2 ) && (numel ( valForPadding ) ~= size(myArray,3) ) )
+        valForPadding = repmat ( valForPadding, size(myArray,3),1 );
+    end    
+
+    classOfMyArray = class ( myArray );
+    
+    newsize = size( myArray );
+    newsize(1) = newsize(1)+2*padSize(1);
+    newsize(2) = newsize(2)+2*padSize(2);
+    
+    startpos = [padSize(1)+1, padSize(2)+1];
+    endpos   = startpos+[size( myArray,1 ), size( myArray,2 )]-1;
+    
+    % be safe about not changing the classtype of the padded array
+    valForPadding = cast ( valForPadding, classOfMyArray);    
+    
+
+    
+    
+    if ( ndims ( myArray ) == 3 )
+        % create 'plain' array
+        paddedArray         = ones( newsize, classOfMyArray );
+        for i=1:size(myArray,3)
+            paddedArray(:,:,i)  = valForPadding(i)*paddedArray(:,:,i);
+        end
+        
+        % set original content to according position    
+        paddedArray (startpos(1):endpos(1), startpos(2):endpos(2), :) = myArray;        
+        
+    elseif ( ndims ( myArray ) == 2 )
+        % create 'plain' array
+        paddedArray         = valForPadding*ones( newsize, classOfMyArray );
+        
+        % set original content to according position
+        paddedArray (startpos(1):endpos(1), startpos(2):endpos(2) ) = myArray;        
+    else
+        disp('Number of feature dimensions does not match!')
+    end
+end

misc/nonMaxSupp.m

@@ -0,0 +1,59 @@
+function [top, pick] = nonMaxSupp(boxes,overlap,numpart)
+% Non-maximum suppression.
+% Greedily select high-scoring detections and skip detections
+% that are significantly covered by a previously selected detection.
+pick=[];
+if nargin < 2
+    overlap = 0.5;
+end
+if nargin < 3
+    numpart = floor(size(boxes,2)/4);
+end
+
+if isempty(boxes)
+  top = [];
+else
+  x1 = zeros(size(boxes,1),numpart);
+  y1 = zeros(size(boxes,1),numpart);
+  x2 = zeros(size(boxes,1),numpart);
+  y2 = zeros(size(boxes,1),numpart);
+  area = zeros(size(boxes,1),numpart);
+  for p = 1:numpart
+    x1(:,p) = boxes(:,1+(p-1)*4);
+    y1(:,p) = boxes(:,2+(p-1)*4);
+    x2(:,p) = boxes(:,3+(p-1)*4);
+    y2(:,p) = boxes(:,4+(p-1)*4);
+    area(:,p) = (x2(:,p)-x1(:,p)+1) .* (y2(:,p)-y1(:,p)+1);
+  end
+  
+  s = boxes(:,5);
+  [vals, I] = sort(s);
+  pick = [];
+  while ~isempty(I)
+	
+    last = length(I);
+    i = I(last);
+	%similar=find(abs(s(I)-s(i))<1e-5);
+	%[m,i]=max(area(I(similar)));
+	%i=I(similar(i));
+	%disp('here');
+	
+	%pause;
+    pick = [pick; i];
+
+    xx1 = bsxfun(@max,x1(i,:), x1(I,:));
+    yy1 = bsxfun(@max,y1(i,:), y1(I,:));
+    xx2 = bsxfun(@min,x2(i,:), x2(I,:));
+    yy2 = bsxfun(@min,y2(i,:), y2(I,:));
+    
+    w = xx2-xx1+1;
+    w(w<0) = 0;
+    h = yy2-yy1+1;
+    h(h<0) = 0;    
+    inter  = sum(w.*h,2);
+    o = double(inter) ./double(sum(area(I,:),2));
+    I(o > overlap) = [];
+	
+  end  
+  top = boxes(pick,:);
+end

misc/settings/addDefaultVariableSetting.m

@@ -0,0 +1,13 @@
+function newSetting = addDefaultVariableSetting( setting, strSettingName, value, newSetting )
+% function newSetting = addDefaultVariableSetting( setting, strSettingName, value, newSetting )
+% 
+% author: Alexander Freytag
+% data:   04-02-2014 (dd-mm-yyyy)
+
+  if ( ( ~isfield(setting,strSettingName))  || isempty(setting.(strSettingName) ) )
+        newSetting.(strSettingName) = value;
+  else
+        newSetting.(strSettingName) = setting.(strSettingName);
+  end  
+  
+end

misc/settings/getFieldWithDefault.m

@@ -0,0 +1,26 @@
+function myOut = getFieldWithDefault ( myStruct, myField, myDefault )
+% function myOut = getFieldWithDefault ( myStruct, myField, myDefault )
+% 
+%  BRIEF:
+%    Get the content of a named field of a struct if existing, or return a
+%    specified default value instead
+%    ...inspired by NICE::Config.gI ('section','name', default )
+% 
+%  INPUT:
+%    myStruct  -- a struct
+%    myField   -- string with desired field name
+%    myDefault -- default value to use if field is non-existing or empty
+% 
+%  OUTPUT:
+%    myOut     -- content of field or default value
+% 
+% 
+% author: Alexander Freytag
+% date  : 04-03-2014 ( dd-mm-yyyy )
+
+    if ( ~isempty(myStruct) && isfield(myStruct, myField) && ~isempty( getfield( myStruct, myField )))
+        myOut = getfield( myStruct, myField );
+    else
+        myOut = myDefault;
+    end
+end

misc/showboxes.m

@@ -0,0 +1,38 @@
+function figHandle = showboxes(im, boxes, partcolor, b_WnHgiven)
+% showboxes(im, boxes)
+% Draw boxes on top of image.
+
+if ( nargin < 4 ) 
+    b_WnHgiven = false;
+end
+
+
+if ( (nargin < 3) || isempty ( partcolor) )
+  partcolor(1)    = {'r'};
+  partcolor(2:20) = {'b'};
+end
+
+%imagesc(im); axis image; axis off;
+handle = imshow(im);
+hold on;
+if ~isempty(boxes)
+  numparts = floor(size(boxes, 2)/4);
+  for i = 1:numparts
+    x1 = boxes(:,1+(i-1)*4);
+    y1 = boxes(:,2+(i-1)*4);
+    if ( b_WnHgiven )
+        x2 = x1+boxes(:,3+(i-1)*4);
+        y2 = y1+boxes(:,4+(i-1)*4);
+    else
+        x2 = boxes(:,3+(i-1)*4);
+        y2 = boxes(:,4+(i-1)*4);        
+    end
+    line([x1 x1 x2 x2 x1]',[y1 y2 y2 y1 y1]','color',partcolor{i},'linewidth',5);
+  end
+end
+drawnow;
+hold off;
+
+if ( nargout > 0 )
+    figHandle = handle;
+end

misc/warpBlocksToStandardSize.m

@@ -0,0 +1,88 @@
+function warped = warpBlocksToStandardSize( model, pos, fh_featureExtractor )
+% function warped = warpBlocksToStandardSize( model, pos, fh_featureExtractor )
+%
+% BRIEF: 
+%     Warp positive examples to fit model dimensions.
+%     Used for training root filters from positive bounding boxes.
+%
+% author: Alexander Freytag
+% date:   13-03-2014 ( last updated)
+
+    if ( nargin < 3 ) 
+        fh_featureExtractor = [];
+    end
+
+    i_modelSize = size(model.w);
+    i_modelSize = i_modelSize(1:2);
+    
+    numpos = length(pos);
+    
+    b_leaveBoundary = getFieldWithDefault ( fh_featureExtractor, 'b_leaveBoundary', false );
+    
+    if ( b_leaveBoundary )
+        cropsize = (i_modelSize+2) .* model.i_binSize;
+        %only needed in this case
+        pixels   = double(i_modelSize * model.i_binSize);
+        heights  = double([pos(:).y2]' - [pos(:).y1]' + 1);
+        widths   = double([pos(:).x2]' - [pos(:).x1]' + 1);
+    else
+        cropsize = (i_modelSize) .* model.i_binSize;
+    end
+    
+    warped  = [];
+    lastreadimg='';
+    
+    for i = 1:numpos
+    %  fprintf('%s: warp: %d/%d\n', name, i, numpos);
+      if(~strcmp(pos(i).im, lastreadimg))	    
+        im          = readImage(pos(i).im);
+        lastreadimg = pos(i).im;
+      end
+      
+      if ( b_leaveBoundary )
+        padx = model.i_binSize * widths(i) / pixels(2);
+        pady = model.i_binSize * heights(i) / pixels(1);
+      else
+        padx=0;
+        pady=0;
+      end
+      
+      x1 = round(double(pos(i).x1)-padx);
+      x2 = round(double(pos(i).x2)+padx);
+      y1 = round(double(pos(i).y1)-pady);
+      y2 = round(double(pos(i).y2)+pady);
+      
+      window = subarray(im, y1, y2, x1, x2, 1);%note: 0 as last option is currently not supported
+      warped{end+1} = imresize(window, cropsize, 'bilinear');%, 'Antialiasing', false);
+    end
+
+    if numpos == 1,
+      assert(~isempty(warped));
+    end
+end
+
+
+
+
+function B = subarray(A, i1, i2, j1, j2, pad)
+% B = subarray(A, i1, i2, j1, j2, pad)
+% Extract subarray from array
+% pad with boundary values if pad = 1
+% pad with zeros if pad = 0
+
+    dim = size(A);
+    %i1
+    %i2
+    is = i1:i2;
+    js = j1:j2;
+
+    if pad
+      is = max(is,1);
+      js = max(js,1);
+      is = min(is,dim(1));
+      js = min(js,dim(2));
+      B  = A(is,js,:);
+    else
+      % todo
+    end
+end