I need to select the points with in the regions that closest match the template, with the knowlede that the region may not contain any valid points.

A Voronoi diagram of your template will give you a set of "nearest-neighbour" regions. That is, each template point will have a Voronoi cell associated with it; every input point in that cell is closer to that cell's template point than to any other template point. (That would be so much easier to say in LaTeX math-mode.)

If your points are on a grid, you can build a bitmap from the Voronoi diagram, with a unique colour for each template point's Voronoi cell. Checking inputs is simply a matter of looking up their position in the bitmap (finding the colour, and thus which template point they're closest to), and calculating the distance to that point. Simple and O(n) once you have the Voronoi diagram.

In the continuous domain, things are probably a bit trickier. Find a good computational geometry text in your local university library; it should have plenty of stuff on this matter.

Here is a nice Voronoi applet that graphically shows the cells.

cheers

tachyon

Given there are probably going to be rotation and x-y registration issues it seems to me that reducing the problem from 2 dimensions to a single dimension vector may work very swell in practice. It is easier to explain in code but the guts of the idea is simple. No matter how rotated or displaced, any planar figure (and the points define such a figure) will have a constant geometric centroid. Each point will have a constant distance from this centriod in one dimensional space regardless of rotation or displacement. We can use this to generate a very useful vector. We can also make the matching as precise of fuzzy as desired.

use constant X => 0;
use constant Y => 1;

my @points1   = ( [1,1], [2,2], [3,3] );     # a line like /
my @points2   = ( [2,1], [3,2], [4,3] );     # a line like / displaced
+ R
my @points3   = ( [1,1], [2.41,1], [3.82,1]);# a line like -
my @template  = ( [1,3], [2,2], [3,1] );     # a line like \

my $vp1 = vector(\@points1,2);
my $vp2 = vector(\@points2,2);
my $vp3 = vector(\@points3,2);
my $vt  = vector(\@template,2);
print "@$vp1\n@$vp2\n@$vp3\ntemplate\n@$vt\n";

sub vector {
    my ($points, $precision) = @_;
    $precision ||= 10;
    # calculate the geometric centroid
    my ( $cx, $cy ) = (0,0);
    for (@$points) {
        $cx += $_->[X];
        $cy += $_->[Y];
    }
    $cx /= @$points;
    $cy /= @$points;
    
    # map the 1 dimensional distance of all points from centroid
    # using simple pythagorean geometry
    my @vector = map{ sqrt(($cx-$_->[X])**2 + ($cy-$_->[Y])**2) } @$po
+ints;
    
    # fuzzify to desired precision and sort to create finished vector
    my $template = "%.${precision}f";
    @vector = map{ sprintf $template, $_ } sort { $a <=> $b } @vector;
  return \@vector;
}
__DATA__
0.00 1.41 1.41
0.00 1.41 1.41
0.00 1.41 1.41
template
0.00 1.41 1.41
[download]

So in one dimensional space our vector tells us that the 3 lines all match the template. The effects of rotation and displacement have been removed. You can also deal with scaling simply by using normalisation.

Matching the single dimension vector against a template set can be achieved in many ways. Stringifying the vector at a reliable precision and using a hash lookup against the template vectors is an obvious method.

This vector process has some caveats. For any given 1 dimensional vector there are an infinite number of 2 dimensional shapes that map to it. For example a vector of 1,1,1,1 in 1 dimensional space can represent any number of rectangles or a square in 2 dimensional space.

C = centroid

.       .           .         .          :      C      :
                         C
    C               .         .    

.       .
[download]

I like this solution because it is simple and clever, but the three problems that I see (and I may be wrong) are as follows:

1. I have seven vertices making up a 2d shape, this shape may have suffered rotational and translational displacement. I have a rough idea of the maximum displayment, so I look at a given size ROI around where the points would be with no displacement. After discriminating the noise for size of points I am usually left with 1-2 possible points (not bad), but I have seen as many as twelve. So looking at a worse case (but very possible), if I am left with 10 points in each ROI, then starting with one ROI and working through each set of points I have 10^7 possible geometric objects that roughly look like my template. lots of checks....

2. The distance between the points in the ROI << then the distance to the calculated centroid, so I am concerned about the ability to consistantly choose the seven correct points and not seven almost true points (similar to what you were saying about infinite number of shapes containg the same centroid). In your example it only looks at rotational displacement, there will be translation displacement, but that is easily delt with looking at the relative distances from the points to the centroid....so still valid.

3. What if the ROI contains points, but does not contain the True point?


This does give me some thoughts that I will work through. I am going to take all of the suggestions, try to hash something out in the next couple of days, and repost what I have done.

Cameron

The method copes fine with 360 degrees of rotation or any ammount x-y translation (there is an even an example of that). This is simple a function of mapping from 2D to 1D space. In one dimensional space the concepts or rotation and x-y tranlation simply have no meaning.

However the method shown relies on having *all* the points, no extra points and no missing points. This is because it is using the centriod to create the vector. A single missing or additional point will result in a different centroid, and vastly different vector.

It is not clear what the real problem you are dealing with is. Curious? It seems now that missing and additional points are an issue. If so you can modify the method as follows.

Without an accurate centroid we must choose one (or more) arbitrary points to use for our 2D->1D mapping. You could choose x_min, or y_min or any other arbitrary point but given templates with ~ 10 points I would just choose all of them in sequence, ie arbitratily declare each point the 'centroid' in turn and calculate the 1D distance mapping vector as before. Insert this data into some sort of easy lookup hash structure. You probably won't need the sort anymore which is O(logn) anyway. You now have all 10 possible 1D mappings of your template object from every possible point, thus rotation, translation, and added/missing points are not an issue. For any given point set you can rapidly check that sets correlation against all the templates in O(n) time.

I will be interested to hear how you get on.

cheers

tachyon

The following code and output attempts to demonstrate that it is possible to generate a single numerical "signature" or digest for your sets of points that will cause similar sets of points to sort close together. By generating these digests for each set in your database and then generating one for your new film as you receive it you can use a binary search/tree to locate the insertion point into the database, and so find those that have a similar signature. You then do a more thorough spacial comparison of a few either side to locate the best matches.

This is not an exact science as the following output shows.

The first set of data (above the ----- line) are the signatures and related datasets read from the DATA section of the code.

The second set of data includes all the first set again, plus the same datasets with modifications made, sorted by sugnature. In some I have deleted one or more pairs, other I've adjusted some values at random by +/-10.

The output show that the algorithm I'm using (the third I tried) does a fairly good job of grouping the modified dataset close or adjacent to the original it is derived from.

It's not perfect. With more tweaking and testing it should be possible to adjust it to do a better (more sensitive) job, but that requires a detailed knowledge (and sample datasets) to do properly.

The basic approach (see calcSig()) is :

1. Produce normalised (0-1) value for each X & Y.
2. Sum all the normalised Xs and divide by the number of pairs to produce a "normalised average X".
3. Do the same for the Ys.
4. Multiply the X average by a factor and int to produce an integer large enough to have enough significant digits to distribute the X signatures nicely.

   This should probably be increased slightly.

5. Add the normalised average Y value to produce the floating point signature.

   These last two step are prime candidates for improvement as this causes the signature to favour the Xs. The idea is to find an algorithm that distributes the datasets over a wide hash space fairly evenly without wrapping. I have achieved that here yet.

Here's the code. It's not documented. If you are interested in pursuing the idea further and have questions, fell free to ask.

#! perl -slw
use strict;
use Data::Dumper;
use List::Util qw[ reduce sum ];
use constant {
    MAXX => 100,
    MAXY => 100,
};

sub calcSig {
    my $sum = reduce { 
        $a->[ 0 ] += $b->[ 0 ] / MAXX;
        $a->[ 1 ] += $b->[ 1 ] / MAXY;
        $a;
    } [], @_;
    $sum->[ 0 ] /= @_; 
    $sum->[ 1 ] /= @_;
    return int( $sum->[ 0 ] * MAXX * MAXY ) + $sum->[ 1 ]    
}

my( @sigs, %sigs );
while( <DATA> ) {
    chomp;
    my @points;
    push @points, [ $1, $2 ] while m[\((\d+?),(\d+?)\)(?:,\s|$)]g;
#    print Dumper \@points;
    push @sigs, my $sig = calcSig @points;
    $sigs{ $sig } = \@points;
}

@sigs = sort{ $a <=> $b } @sigs;
printf "%9.5f => %s\n", $_, join ' ', map{ @$_ } @{ $sigs{ $_ } } for 
+@sigs;
print "\n ------ \n";

my @adjusted = map{ [ @$_ ] } values %sigs;
for my $aref ( @adjusted ) {
    if( rand > .3 ) { ## delete a pair
        splice @$aref, rand( @$aref ), 1;
    }
    else { ## Adjust a few values;
        @$aref = map{
            rand > .5 and $_->[ 0 ] += -10+int( rand 20 );
            rand > .5 and $_->[ 1 ] += -10+int( rand 20 );
            $_;
        } @$aref;
    }
    push @sigs, my $sig = calcSig @$aref;
    $sigs{ $sig } = $aref;
}
@sigs = sort{ $a <=> $b } @sigs;
printf "%9.5f => %s\n", $_, join ' ', map{ @$_ } @{ $sigs{ $_ } } for 
+@sigs;

__DATA__
(67,66), (7,45), (77,19), (4,84), (74,53), (0,65), (37,12), (91,65), (
+84,36)
(92,44), (27,30), (94,43), (51,96)
(46,57), (45,2), (64,46), (46,1), (40,54)
(84,26), (93,43), (11,80), (16,50), (66,80), (46,76)
(92,32), (91,28), (80,76), (55,2)
(94,40), (90,13), (57,36), (83,3)
(13,69), (53,24), (61,90), (56,28), (61,82), (17,85), (12,59), (14,55)
(33,44), (1,71), (35,71), (97,26), (39,78), (31,3), (32,65), (87,61)
(70,26), (38,71), (78,36), (84,58), (84,5), (81,75)
(70,98), (20,25), (79,19), (18,49), (48,78), (50,77), (40,53)
(27,32), (54,19), (46,20), (38,93), (61,19), (9,53), (25,1)
(43,83), (41,83), (16,69), (99,49), (22,60)
(98,42), (79,0), (72,26), (95,11), (89,38)
(67,0), (57,24), (46,50), (44,38)
(76,20), (41,22), (93,64), (74,92)
(82,24), (88,16), (49,9), (63,19), (84,77), (25,90), (48,11), (69,45),
+ (30,88)
(53,23), (72,95), (27,29), (94,59), (18,22), (89,33), (57,97), (63,95)
(20,43), (25,36), (33,84), (63,11), (15,2), (70,69), (47,77), (67,91),
+ (96,32)
(59,55), (4,88), (74,63), (20,33), (28,34), (3,93), (9,19), (10,22), (
+42,44)
(8,33), (21,55), (51,48), (6,20), (52,63)
(72,27), (30,52), (37,35), (7,76), (23,45), (71,25), (75,17), (86,29),
+ (95,30)
(51,38), (31,65), (62,4), (37,35), (73,48)
(98,33), (3,22), (75,40), (64,18), (91,67), (8,65), (29,68), (11,96)
(45,32), (65,80), (58,22), (95,44), (51,36), (93,61), (76,94)
(55,90), (88,6), (9,1), (56,45), (26,56)
(26,21), (25,46), (46,14), (79,11), (85,39), (30,9)
(37,18), (90,23), (38,13), (97,63), (17,26), (41,56), (22,85)
(13,42), (79,34), (5,47), (25,20), (78,35), (59,75), (89,42), (0,26), 
+(48,36)
(60,73), (51,45), (95,63), (98,28)
(91,84), (9,40), (13,98), (88,60), (17,83)
[download]

Sorry to take so long to reply, I was working through the code and thinking about implementation. Nifty idea, but not I do not believe practical for my situation. I have 7 ROIs, each could contain as many as 10 points, and of those ten the true point may or may not exist. So coming up with all of the possible combinations/permutations could make this solution prohibitive.

I have been thinking about filtering by point instead of collection (at least initially). For good registration, I need at least 4 points. Knowing the distance between the points I can work through each of the points and check the distance between each of the other points. If I set some tolerance, then I can discriminate the points that are not within the tolereance for at least 4 of the 6 distances. This should be fairly quick { 6 distances per point max ~ 10 points per ROI => 420 checks, and if I splice the array for the ROI then this number should decrease as bad points are thrown out. }


If I can run a few tests then I will see what we typically have left. If it is only 1-3 points per ROI then your solution is much more viable.


Quick question (I have not looked for this yet, thought of it while writting), do you or does anyone know of a good module for doing the permutations I need? i.e. N arrays with M elements, creating new unique arrays containing only 1 element from each of the N arrays? I am sure that it is out there, I will look for it tonight.

Thanks for you help,


Cameron

...does anyone know of a good module for doing the permutations...

Algorithm::Permute and NextPermute() from tye's Algorithm::Loops module are both very good implementations.

I think that the idea of generating a signature could be made to work well for your application. I saw the extra information that you posted in reply to tachyon and I think I see the problems your describing, and also some of a possible solution. The key to taking the thoughts further is having some test data to work with. If you felt like placing some on your doowah2004's scratchpad, and /msg'ing me, I enjoy playing with it.

That said, I'm still not entirely sure that I understand the full scope of the problem, so maybe I should just accept your more informed judgement on the matter :)

---

Examine what is said, not who speaks.

"Efficiency is intelligent laziness." -David Dunham
"Think for yourself!" - Abigail
"Memory, processor, disk in that order on the hardware side. Algorithm, algorithm, algorithm on the code side." - tachyon

How about Algorithm::Points::MinimumDistance

The points are puntures through the film by a physical template, so the distance between the true points should be very close to the measured distance for the template. But when the holes are made the film might be slightly rotated. What this means is that the scanned film x,y values are meaningless when comparing to the template x,y values. BUT the delta values (relative distances) between the points should closely correlate.

You might have more luck with a rigid point correspondence algorithm, then.

Edit: Oops, posted too soon. Provided that the deltas are pretty small, the Voronoi diagram solution should still work. If it does, it will probably be easier to implement and much quicker than a point-correspondence algorithm (most of the PC approaches I'm aware of start by building an affinity matrix between all of the feature points, then do an eigenvalue decomposition on that, so they're going to be pretty slow compared to a Voronoi-cell lookup).

If you're looking for hella speed and have the resources to do it, you could probably do the Voronoi lookup fairly easily with programmable graphics hardware (read: consumer graphics cards no more than, say, two years old). Iterating through, say, a thousand possible rotations of the film onto the template should be a piece of cake; that's exactly the sort of thing that modern graphics cards are optimized for. If this sounds like something you have the resources to do, send me a /msg -- we're getting a bit off-topic here, but we're also getting closer to my area of research.

--
F o x t r o t U n i f o r m
Found a typo in this node? /msg me
% man 3 strfry

use Algorithm::Loops qw/ NestedLoops /;

$|++;


$i = 0;
$toggle = -1;
@p = qw|ROI1 ROI2 ROI3 ROI4 ROI5 ROI6 ROI7|;




$precision = shift;
$match = shift;
$minpoints = shift;
system('cls');

if ((!$precision) || (!$match) || (!$minpoints)) {
    print "Usage: pointfilter.exe x y z\n\n\tx,\tPrecision in pixels f
+or the distances.",
    "\n\t\tDecrease this value to lower total points.",
    "\n\t\tBe aware, lowering it too much will return bad solutions.",
    "\n\t\tTypical values are 10 - 20.",
    "\n\n\ty,\tNumber of distances that a point must be within toelere
+nce to be considered valid.",
    "\n\t\tThis number will typically be 4 or 5, 6 is maximum.",
    "\n\n\tz,\tThe minimum amount of points contained in the best fit.
+",
    "\n\t\tTypical values would be 4 or 5, minimum is 4, maximum is 7.
+";
    exit (0);
}





##### Get Template Distances #####

open (TEMPLATE, "<avedist.txt");
@avedist = <TEMPLATE>;
close (TEMPLATE);
for $j(0..6) {
    for $k(0..6) {
        if ($j == $k) {}
        else {
            chomp $avedist[$i];
            $data->{template}->{$p[$j]}->{$p[$k]}->{AVE} = $avedist[$i
+];
            $i++;
        }
    }
}





##### Get the Film Point Data #####

$i = 0;
$alpha = "A";
$m = 0;
open (POINTS, "<points.txt");
while (<POINTS>) {
    chomp $_;
    ($a,$b,$c) = split(/\t/,$_);
    if ($a == 1) {
        $toggle *= -1;
        if ($toggle < 0) {
            $i++;
            $alpha = "A";
            $toggle = 1;
        }
    }
    if ($b == -1) {
        $i++;
        $alpha = "A";
        $toggle = -1;
    } else {
    $data->{$p[$i]}->{$p[$i]."point".$alpha}->{X} = $b;
    $data->{$p[$i]}->{$p[$i]."point".$alpha++}->{Y} = $c;
    $m++;
    }
}
close(POINTS);

print "Total number of Points = $m\n\n";
print "Running first filtering pass...";





##### First Pass at Filtering the Data #####


$passed = 0;
for $i(0..6) {
    foreach $testpoints (sort keys %{$data->{$p[$i]}}) {
        $j = 0;
        $testX = $data->{$p[$i]}->{$testpoints}->{X};
        $testY = $data->{$p[$i]}->{$testpoints}->{Y};
        for $l(0..6) {
            if ($l == $i) {}
            else {
                $truedist = $data->{template}->{$p[$i]}->{$p[$l]}->{AV
+E};
                foreach $points(sort keys %{$data->{$p[$l]}}) {
                    $X = $data->{$p[$l]}->{$points}->{X};
                    $Y = $data->{$p[$l]}->{$points}->{Y};
                    $distance = sqrt(($testX - $X)**2 + ($testY - $Y)*
+*2);
                    $data->{$p[$i]}->{$testpoints}->{$p[$l]}->{$points
+} = abs($truedist - $distance);
                    if (abs($truedist - $distance) < $precision) {
                        $passed = 1;
                    }
                }
                if ($passed == 1) {
                    $j++;
                    $passed = 0;
                }
            }
        }
        if ($j < $match) {
            delete $data->{$p[$i]}->{$testpoints};
        }
    }
}


$i = 0;
foreach (@p) {
    foreach $key(sort keys %{$data->{$_}}){
        $i++;
    }
}


print "done!\n", $m - $i ," points were filtered.  $i points remaining
+.\n\n";
print "Calculating the total number of combinations....";





##### Calculate the Combinations #####

@arrays = ([keys %{$data->{ROI1}}], [keys %{$data->{ROI2}}], [keys %{$
+data->{ROI3}}],
                 [keys %{$data->{ROI4}}], [keys %{$data->{ROI5}}], [ke
+ys %{$data->{ROI6}}],
                       [keys %{$data->{ROI7}}]);

NestedLoops(
  [ map [ undef, @$_ ], @arrays ],  
  sub {
    $count = 0;
    $string = join '', grep { defined && ++$count } @_;
    push @combinations, [grep { defined && ++$count } @_] if $count >=
+ $minpoints;
  },
);


print "done!\nThere are $#combinations total combinations\n\n";
print "Running final filter...";





##### Final Filter #####   
    
$lowscore = 1000000;       
foreach (@combinations){
    @fit = @$_;
    $i = 0;
    for $j(0..$#fit) {
        ($region1,$trash) = split (/point/,$fit[$j]);
        for $k($j+1..$#fit) {
            ($region2,$trash) = split (/point/,$fit[$k]);
            $i++;
            $tempsum += $data->{$region1}->{$fit[$j]}->{$region2}->{$f
+it[$k]};
            if ($data->{$region1}->{$fit[$j]}->{$region2}->{$fit[$k]} 
+> ($precision * 1.25)) {
                $failed = 1;
            }
        }
    }
    do {
        @bestfit = @fit;
        $lowscore = $tempsum / $i;
    } unless (($lowscore < ($tempsum / $i)) || ($lowscore == 0) || ($f
+ailed == 1));
    undef $tempsum;    
    $failed = 0;

}


print "done!\n\n\n";        


print "Bestfit = \n";
foreach $bf(@bestfit) {
    ($region,$trash) = split (/point/,$bf);
    print "$data->{$region}->{$bf}->{X} | $data->{$region}->{$bf}->{Y}
+\t";
}
print "\n";
[download]