Re: Another vectorisation question

[Jaroslav Hajek]

On Fri, Jul 16, 2010 at 3:24 AM, John Reed <address@hidden> wrote:
> Hello,
>
> I am looking at a nearest neighbor problem.  I have a rectangular grid and 
> random coordinates that fall within the 2d grid.  I want to create an array 
> that stores all the random coordinates that are nearest to each respective 
> grid point.
>
> I've taken two different approaches to solving this problem, one uses the 
> find function, and the other the histc function, but in either case I don't 
> know how to get rid of all the for loops and i would like to speed up the 
> code.
>
> Is there a way to vectorise either of these approaches?  or is there some 
> totally different approach that would give me the answer I want without for 
> loops?
>
> The two different approaches I've tried are below.
>
> Any help would be greatly appreciated,
>
> John
> ..........
>
> incr = 70;   % maximum x & y integer value on square unit grid
> pt_num = 40000;
> coord=incr*rand(pt_num,2);
>
> % 1st Approach - find grid point that coordinate is closest to and assign to 
> array
> % holding all coordinates closest to each respective grid point
> tic;
> for i = 0:incr
>   for j = 0:incr
>      x = [ i - 0.5; i + 0.5];
>      y = [ j - 0.5; j + 0.5];
>
>      indices = find( coord(:,1) >= x(1) & ...
>                   coord(:,1) <  x(2) & ...
>                   coord(:,2) >= y(1) & ...
>                   coord(:,2) <  y(2) );
>
> % build array storing random coordinates indices associated with each grid 
> coordinate
>      num_pt_per_grid(i+1,j+1)=length(indices);
>      if ( length(indices) >= 1 )
>         coords_to_grid(i+1,j+1,1:length(indices)) = indices;
>      end
>   end
> end
> toc
>
> fflush(stdout);
>
> % 2nd approach - using histc group all coordinates closest to grid point and 
> assign
> % to array holding all coordinates closest to each grid point
> tic;
> % create edges for histogram
> histo_edges = [-0.5:1:incr+.5];
>
> [x_count, x_index] = histc(coord(:,1),histo_edges);
> [y_count, y_index] = histc(coord(:,2),histo_edges);
>
> % build array storing random coordinates indices associated with each grid 
> coordinate
> num_pt_per_grid2 = zeros(incr+1,incr+1);
> for i=1:length(x_index)
>   num_pt_per_grid2(x_index(i),y_index(i))++;
>   
> coords_to_grid2(x_index(i),y_index(i),num_pt_per_grid2(x_index(i),y_index(i)))
>  = i;
> end
> toc
>
>
>


% 3rd approach - use a cell array of vectors to accumulate.
tic;

edges = [-0.5:1:incr+.5];
idx = lookup (edges, coord);

coords_to_grid3 = accumarray (idx, 1:pt_num, [incr, incr]+1, @(x) {x});
toc


now, the indices falling into grid (i,j) are given by
coords_to_grid3{i,j} instead of coords_to_grid2(i,j,1:num_pt_per_grid(i,j)).


timings in 3.3.51+ (recommended for best performance):

octave:1> tt
Elapsed time is 2.06051 seconds.
Elapsed time is 0.955433 seconds.
Elapsed time is 0.00869894 seconds.

coords_to_grid3 can be converted to coords_to_grid2 with a little
extra effort (if that is needed at all).

hth

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz