title: Efficient histogram computation from image ROIs
body: <p>I am looking for some already existing functions/tools to compute standard Bag of Visual Words histograms from multiple ROI (Regions of Interest) in an image. Let me explain:</p>
  
  <p>(1) Suppose you have an image where each "pixel" carries an integer: 1 ... K
  Each such "pixel" has following information</p>
  
  <ol>
  <li>x,y coordinate</li>
  <li>a value from 1 to K</li>
  </ol>
  
  <p>(2) Suppose a LARGE amount of fixed size regions are sample from all the image in format:</p>
  
  <ol>
  <li>(x1,y1) - top,left coordinate</li>
  <li>(x2,y2) - bottom,right coordinate</li>
  </ol>
  
  <p>(3) For every region: Compute a K bin histogram that counts number of occurrences of the "pixel" values that fall in that region</p>
  
  <p>I have implemented a following function in MATLAB but due to multiple for loops in the code, it is very slow</p>
  
  <pre><code>function [H words] = sph_roi( wind, tree, desc, feat, bins )
  % FUNCTION computes an SPH histogram for a collection of windows. Spatial
  % information is captured by splitting the window in bins horizontally.
  %
  %   [H words] = sph_roi( obj_wind, tree, desc, feat, [ bins ] );
  %
  % INPUT :
  %   wind        - sampled ROI windows
  %                   [left_x, top_y, right_x, bottom_y] - see sample_roi()
  %   tree        - vocabulary tree
  %   desc        - descriptors matrix
  %   feat        - features matrix
  %   bins        - number of horizontal cells (1=BOVW, 2... SPH)
  %                   by default set to the multiples of window height.
  %
  % OUTPUT :
  %   H           - SPH histograms
  %   words       - word IDs found for every descriptor
  %
  
  verbose = 0;
  
  % input argument number check
  if nargin &lt; 4
  error( 'At least 4 input arguments required.' );
  end
  
  % default number of horizontal cells
  if nargin &lt; 5
  bins = -1;  % will be set in multiples of each window height corresp.
  end
  
  % number of windows
  num_wind = size( wind, 1 );
  
  % number of visual words
  num_words = tree.K;
  
  % pre-compute all visual words
  words = vl_hikmeanspush( tree, desc );
  
  % initialize SPH histograms matrix
  H = zeros( num_words * bins, num_wind );
  
  % compute BOVW for each ROI
  for i = 1 : num_wind
  
  if verbose == 1
      fprintf( 'sph_roi(): processing %d / %d\n', i, num_wind );
  end
  
  % pick a window
  wind_i = wind( i, : );
  
  % get the dimensions of the window
  [w h] = wind_size( wind_i );
  
  % if was not set - the number of horizontal bins
  if bins == -1
      bins = round( w / h );
  end
  
  % return a list of subcell windows
  scw = create_sph_wind( wind_i, bins );
  
  for j = 1 : bins
  
      % pick a cell
      wind_tmp = scw( j, : );
  
      % get the descriptor ids falling in that cell
      ids = roi_feat_ids( wind_tmp, feat );
  
      % compute the BOVW histogram for the current cell
      h = vl_hikmeanshist( tree, words(ids) );
  
      % assemble the SPH histogram in the output matrix directly
      H( 1+(j-1)*num_words : j*num_words, i ) = h( 2:end );
  
  end
  
  end
  
  function ids = roi_feat_ids( w, f )
  % FUNCTION returns those feature ids that fall in the window.
  %
  %   ids = roi_feat_ids( w, f );
  %
  % INPUT :
  %   w   - window
  %   f   - all feature points
  %
  % OUTPUT :
  %   ids - feature ids
  %
  
  % input argument number check
  if nargin ~= 2
  error( 'Two input arguments required.' );
  end
  
  left_x = 1;
  top_y = 2;
  right_x = 3;
  bottom_y = 4;
  
  % extract and round the interest point coordinates
  x = round( f(1,:) );
  y = round( f(2,:) );
  
  % bound successively the interest points
  s1 = ( x &gt; w(left_x) ); % larger than left_x
  s2 = ( x &lt; w(right_x) ); % smaller than right_x
  s3 = ( y &gt; w(top_y) ); % larger than top_y
  s4 = ( y &lt; w(bottom_y) ); % smaller than bottom_y
  
  % intersection of these 4 sets are the ROI enclosed interest points
  ids = s1 &amp; s2 &amp; s3 &amp; s4;
  
  % convert ids to real
  ids = find( ids );
  </code></pre>
  
  <p>I've looked at routines proposed by OpenCV and even in Intel's MKL but found nothing appropriate. Using the Matlab's profiler, I found that considerable time is spent in the roi_feat_ids() and the outer loop over each region in the function sph_roi() is slow too. Before trying to implement a MEX file, I would like to see if I could recycle some existing code. </p>