% THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION IS RELEASED "AS IS." THE U.S. GOVERNMENT MAKES NO WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, CONCERNING THIS SOFTWARE AND ANY ACCOMPANYING DOCUMENTATION, INCLUDING, WITHOUT LIMITATION, ANY WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE. IN NO EVENT WILL THE U.S. GOVERNMENT BE LIABLE FOR ANY DAMAGES, INCLUDING ANY LOST PROFITS, LOST SAVINGS OR OTHER INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE, OR INABILITY TO USE, THIS SOFTWARE OR ANY ACCOMPANYING DOCUMENTATION, EVEN IF INFORMED IN ADVANCE OF THE POSSIBILITY OF SUCH DAMAGES.
%
% file: group_detection_analyzer.m
%
% Inputs:
% object_locations - a column of (x,y) location pairs
% group_min_object_count - the fewest number of objects that constitute a "group"
% group_max_separation - the largest distance between objects within a group
% dist_measure - the type of distance measure to use
% 1 for Euclidean
% 2 for taxicab
% max_group_size - the largest number of objects ever expected to occur in a group (used for "memory allocation")
%
% Outputs:
% group_report - a row vector containing:
% group_count - the number of groups detected
% object_count_post_group - number of objects that are in the detected groups
% object_count_pre_group - number of input objects (whether in a group or not)
% group_freq - a vector of length max_group_size where the value of group_freq(i) is the number of groups with i members
%
% group_detector's run time (PC, c.2002):
% objects seconds
% 100 0.02
% 300 0.44
% 500 1.6
% 1000 9.0
%
% for group detection background, see T. D. Ross, M. Bryant, R. Dilsavor, J. Jackson, "Performance Assessment of Group Detection Algorithms", Proc. SPIE, Vol. 5427, Algorithms for Synthetic Aperture Radar Imagery XI, April 2004
% --------- change log ---------------------
% 030306 tdr built from group detection code of Dilsavor, 11/02, and Bryant 2/03
% 030314 tdr accepts object_locations rather than a filename
% 030325 tdr added distance_measure parameter
% 040130 tdr updates for delivery to Erik
% added comments
% included subfunctions in this file (rather than leaving them as separate files)
% 040309 cleared for Public Release by ASC Public Affairs. Disposition Date: 3/9/2004. Document Number ASC 04-0679.
function group_report = group_detection_analyzer(object_locations, group_min_object_count, ...
group_max_separation, dist_measure, max_group_size)
%tic;
object_count_pre_group = size(object_locations,1);
switch dist_measure
case 1 % Euclidean
if object_count_pre_group > 3000
D = sq_dist2(object_locations, object_locations);
else
D = sq_dist1(object_locations);
end;
case 2 % taxicab
D = sq_dist3(object_locations);
end;
D = D < (group_max_separation^2);
group_count = 0;
object_count_post_group = 0;
[biggest_group_size, biggest_group_index] = max(sum(D,1));
% sum(D,1) gets the sum of columns, which is the number of objects within
% group_max_separation of the object associated with a particular column,
% including that object itself.
group_freq = zeros(max_group_size,1);
if biggest_group_size > max_group_size error('group size exceeded max_group_size'); end;
while biggest_group_size >= group_min_object_count,
group_count = group_count + 1;
object_count_post_group = object_count_post_group + biggest_group_size;
group_freq(biggest_group_size) = group_freq(biggest_group_size) + 1;
unused_rows = (D(:,biggest_group_index)==0);
D = D(unused_rows,:);
[biggest_group_size, biggest_group_index] = max(sum(D,1));
end
% report group count, size, ...
group_size_array = 1:max_group_size;
object_group_count = group_size_array' .* group_freq;
group_freq(1) = object_count_pre_group - sum(object_group_count);
% finish filling in the group_freq for smaller groups
while biggest_group_size >= 2,
group_freq(biggest_group_size) = group_freq(biggest_group_size) + 1;
unused_rows = (D(:,biggest_group_index)==0);
D = D(unused_rows,:);
[biggest_group_size, biggest_group_index] = max(sum(D,1));
end
group_report = [[group_count object_count_post_group object_count_pre_group] group_freq'];
%toc
return
%==============================================================
% distance function
% input is a list of x,y locations
% output is an inner distance matrix
% This is only good for up-to a few thousand locations (1.8GHz Pentium 4, 0.5GB RAM):
% - 1000 locations - takes a few seconds
% - 2500 locations - takes a few tens of seconds
% - 4000 locations - takes ten minutes or so
% - 5300 locations - runs out of memory
% so we compute the distance matrix in 1000x1000 blocks (roughly
% speaking)
function dist_squared_array = sq_dist(locations_rows, locations_columns)
dist_squared_array = 0;
object_count_rows = size(locations_rows,1);
object_count_columns = size(locations_columns,1);
if (object_count_rows == 0) | (object_count_columns == 0) return; end;
if max([object_count_rows object_count_columns]) > 1000
mid_point_rows = ceil(object_count_rows/2);
mid_point_columns = ceil(object_count_columns/2);
locations_rows_a = locations_rows(1:mid_point_rows,:);
locations_rows_b = locations_rows(mid_point_rows+1:object_count_rows,:);
locations_columns_a = locations_columns(1:mid_point_columns,:);
locations_columns_b = locations_columns(mid_point_columns+1:object_count_columns,:);
dist_squared_array_aa = sq_dist(locations_rows_a, locations_columns_a);
dist_squared_array_ab = sq_dist(locations_rows_a, locations_columns_b);
dist_squared_array_ba = sq_dist(locations_rows_b, locations_columns_a);
dist_squared_array_bb = sq_dist(locations_rows_b, locations_columns_b);
dist_squared_array = [dist_squared_array_aa dist_squared_array_ab; dist_squared_array_ba dist_squared_array_bb];
clear dist_squared_array_*;
else
cmplx_locations_rows = locations_rows(:,1) + locations_rows(:,2)*i;
cmplx_locations_columns = locations_columns(:,1) + locations_columns(:,2)*i;
clear locations_*;
cmplx_dist_matrix = (cmplx_locations_rows * ones(1,object_count_columns)) - (ones(object_count_rows,1) * cmplx_locations_columns.');
clear cmplx_locations_rows cmplx_locations_columns;
dist_squared_array = cmplx_dist_matrix .* conj(cmplx_dist_matrix);
clear cmplx_dist_matrix;
end
return
%==============================================================
% distance function 1
% input is a list of x,y locations
% output is an inner distance matrix
% This is only good for up-to a few thousand locations (1.8GHz Pentium 4, 0.5GB RAM):
% - 1000 locations - takes a few seconds
% - 2500 locations - takes a few tens of seconds
% - 4000 locations - takes ten minutes or so
% - 5300 locations - runs out of memory
% use sq_dist2.m for larger lists (takes a long time, but doesn't run out
% of memory on the 5300 long list.
function dist_squared_array = sq_dist1(locations)
dist_squared_array = 0;
object_count = size(locations,1); % total number of detections
if object_count == 0 return; end;
cmplx_locations = locations(:,1) + locations(:,2)*i;
clear locations;
cmplx_dist_matrix = (cmplx_locations * ones(1,object_count)) - (ones(object_count,1) * cmplx_locations.');
clear cmplx_locations;
dist_squared_array = cmplx_dist_matrix .* conj(cmplx_dist_matrix);
clear cmplx_dist_matrix;
return
%==============================================================
% distance function 2
% input is a list of x,y locations
% output is an inner distance matrix
%
% This seems to run for a list of 5300 (takes 15 minutes though).
% If the list of locations is small (less than a couple thousand) then sq_dist1.m is better.
% The approach here is to break up the list and recursively call itself
% until the list is < 1000 then it runs a sq_dist1 type of routine.
% There's some extra generality required here though, so this isn't as fast
% as sq_dist1.m even for small lists.
function dist_squared_array = sq_dist2(locations_rows, locations_columns)
dist_squared_array = 0;
object_count_rows = size(locations_rows,1);
object_count_columns = size(locations_columns,1);
if (object_count_rows == 0) | (object_count_columns == 0) return; end;
if max([object_count_rows object_count_columns]) > 1000
mid_point_rows = ceil(object_count_rows/2);
mid_point_columns = ceil(object_count_columns/2);
locations_rows_a = locations_rows(1:mid_point_rows,:);
locations_rows_b = locations_rows(mid_point_rows+1:object_count_rows,:);
locations_columns_a = locations_columns(1:mid_point_columns,:);
locations_columns_b = locations_columns(mid_point_columns+1:object_count_columns,:);
dist_squared_array_aa = sq_dist(locations_rows_a, locations_columns_a);
dist_squared_array_ab = sq_dist(locations_rows_a, locations_columns_b);
dist_squared_array_ba = sq_dist(locations_rows_b, locations_columns_a);
dist_squared_array_bb = sq_dist(locations_rows_b, locations_columns_b);
dist_squared_array = [dist_squared_array_aa dist_squared_array_ab; dist_squared_array_ba dist_squared_array_bb];
clear dist_squared_array_*;
else
cmplx_locations_rows = locations_rows(:,1) + locations_rows(:,2)*i;
cmplx_locations_columns = locations_columns(:,1) + locations_columns(:,2)*i;
clear locations_*;
cmplx_dist_matrix = (cmplx_locations_rows * ones(1,object_count_columns)) - (ones(object_count_rows,1) * cmplx_locations_columns.');
clear cmplx_locations_rows cmplx_locations_columns;
dist_squared_array = cmplx_dist_matrix .* conj(cmplx_dist_matrix);
clear cmplx_dist_matrix;
end
return
%==============================================================
% distance function 3 - taxicab distances
% input is a list of x,y locations
% output is an inner distance matrix
function dist_squared_array = sq_dist3(locations)
dist_squared_array = 0;
object_count = size(locations,1); % total number of detections
if object_count == 0 return; end;
cmplx_locations = locations(:,1) + locations(:,2)*i;
clear locations;
cmplx_dist_matrix = (cmplx_locations * ones(1,object_count)) - (ones(object_count,1) * cmplx_locations.');
clear cmplx_locations;
taxicab_dist_matrix = abs(real(cmplx_dist_matrix)) + abs(imag(cmplx_dist_matrix));
dist_squared_array = taxicab_dist_matrix .^ 2;
clear taxicab_dist_matrix;
return
