Hausdorff Distance between images

Started by
11 comments, last by Dmytry 17 years, 8 months ago
Hello, I want to match an image A and a template B. Therefore I compare a window from A with the template (both have same dimensions) and get the Hausdorff distance. The images are edge-images (binary); a 1 represents an edge. The number of "1" between the images can differ due to occlusions or noise. I thought about using the following algorithm (pseudocode), but it's extremely CPU-intensive:

// Dimensions of the two images
widthA = widthB;
heightA = heightB;

// Generate two arrays (P_A, P_B) containing coordinates of points
// where P_A[0][0] is the x-, and P_A[0][1] is the y-coordinate for instance.
for (i=0; i<widthA; i++)
{
	for (j=0; j<heightA; j++)
	{
		if (A[j] != 0)
			AddPoint(*P_A, i,j); // Function that adds a point to the array P_A
		if (B[j] != 0)
			AddPoint(*P_B, i,j); // Function that adds a point to the array P_B
	}
}

// Number of found points in A, B
sizeA = length(P_A);
sizeB = length(P_B);

// Calculate Hausdorff distance d(A,B)
maxDistAB = 0;
for (i=0; i<sizeA; i++)
{
	minB = 1000000; // Helpervariable holding current minimum B
	for (j=0; j<sizeB; j++)
	{
		tmpDist = sqrt((P_A[0]-P_B[j][0])^2 + (P_A[1]-P_B[j][1])^2);

		if (tmpDist < minB)
		{
			minB = tmpDist;
		}
	}

	maxDistAB += minB;
}
maxDistAB *= 1/widthA*heightA;

// Calculate Hausdorff distance d(B,A)
maxDistBA = 0;
for (i=0; i<sizeB; i++)
{
	minA = 100000; // Helpervariable holding current minimum A
	for (j=0; j<sizeA; j++)
	{
		tmpDist = sqrt((P_B[j][0]-P_A[0])^2 + (P_B[j][1]-P_A[1])^2);

		if (tmpDist < minA)
		{
			minA = tmpDist;
		}
	}

	maxDistBA += minA;
}
maxDistBA *= 1/widthA*heightA;

HausdorffDist = max(maxDistAB,maxDistBA);


Formulas: Can this algorithm be made more reliable, or are there errors in it? Thank you, [Edited by - dnaxx on July 28, 2006 2:50:22 AM]
Advertisement
One thing I see:

tmpDist = sqrt((P_A[0]-P_B[j][0])^2 + (P_A[1]-P_B[j][1])^2);

I don't think this does what you want. The ^ operator in C++/Java is exclusive or, not power.
"What are you trying to tell me? That I can write an O(N^2) recursive solution for a 2-dimensional knapsack?" "No, programmer. I'm trying to tell you that when you're ready, you won't have to." -Adapted from "The Matrix"
The code above is only pseudo-code. I implemented it in C and it works (also the Hausdorff distance is correct). But the speed is not the best.
To find out where your program is slowest, put timers around each double for loop as this is where all your CPU will be going to.
maxDistAB *= 1/widthA*heightA
This line means
maxDistAB *= (1/widthA)*heightA
not
maxDistAB *= 1/(widthA*heightA)
which is what I think you want (and the same with maxDistBA).

Unless I missed something from a brief glance - you can move the sqrts outside the loops. Also check that the compiler is inlining/optimising the (P_A[0]-P_B[j][0])^2 type lines - it might be better to write that as a single inline expression that definitely avoids temporaries.
Yep, definately move that sqrt out of the inner loops, just find the min of the squared distances and do the sqrt once you have the min not only does this remove the expensive sqrt but it lets keep everything in the inner loop as integer math.

Also, how does AddPoint() work? Make sure you're not doing anything silly like dynamically resizing the array, just reserve the maximum possible size (ie. widthA*heightA) to begin with.
"Voilà! In view, a humble vaudevillian veteran, cast vicariously as both victim and villain by the vicissitudes of Fate. This visage, no mere veneer of vanity, is a vestige of the vox populi, now vacant, vanished. However, this valorous visitation of a bygone vexation stands vivified, and has vowed to vanquish these venal and virulent vermin vanguarding vice and vouchsafing the violently vicious and voracious violation of volition. The only verdict is vengeance; a vendetta held as a votive, not in vain, for the value and veracity of such shall one day vindicate the vigilant and the virtuous. Verily, this vichyssoise of verbiage veers most verbose, so let me simply add that it's my very good honor to meet you and you may call me V.".....V
Quote:Original post by Veladon Fireheart
To find out where your program is slowest, put timers around each double for loop as this is where all your CPU will be going to.


Umm, it's fairly obvious the slowest part is the inner loops with the sqrt's.

Quote:Original post by Veladon Fireheart
maxDistAB *= 1/widthA*heightA
This line means
maxDistAB *= (1/widthA)*heightA
not
maxDistAB *= 1/(widthA*heightA)
which is what I think you want (and the same with maxDistBA).


Since there's no way it would work either way with integer math he must be using floating point there, in which case I doubt the precision difference between the 2 would matter. And besides, its only pseudocode.
"Voilà! In view, a humble vaudevillian veteran, cast vicariously as both victim and villain by the vicissitudes of Fate. This visage, no mere veneer of vanity, is a vestige of the vox populi, now vacant, vanished. However, this valorous visitation of a bygone vexation stands vivified, and has vowed to vanquish these venal and virulent vermin vanguarding vice and vouchsafing the violently vicious and voracious violation of volition. The only verdict is vengeance; a vendetta held as a votive, not in vain, for the value and veracity of such shall one day vindicate the vigilant and the virtuous. Verily, this vichyssoise of verbiage veers most verbose, so let me simply add that it's my very good honor to meet you and you may call me V.".....V
@joanusdmentia,
The precision difference? Say widthA = heightA = 10, then (1/10)*10 = 1 whereas 1/(10*10) = 1/100. I'd call that quite a big difference.

Quote:Original post by Veladon Fireheart
@joanusdmentia,
The precision difference? Say widthA = heightA = 10, then (1/10)*10 = 1 whereas 1/(10*10) = 1/100. I'd call that quite a big difference.


Ah crap, my bad [smile] No idea where my brain was on that one.....
"Voilà! In view, a humble vaudevillian veteran, cast vicariously as both victim and villain by the vicissitudes of Fate. This visage, no mere veneer of vanity, is a vestige of the vox populi, now vacant, vanished. However, this valorous visitation of a bygone vexation stands vivified, and has vowed to vanquish these venal and virulent vermin vanguarding vice and vouchsafing the violently vicious and voracious violation of volition. The only verdict is vengeance; a vendetta held as a votive, not in vain, for the value and veracity of such shall one day vindicate the vigilant and the virtuous. Verily, this vichyssoise of verbiage veers most verbose, so let me simply add that it's my very good honor to meet you and you may call me V.".....V
thanks for your suggestions.
The most performance is used by this part:
helperX = (P_A[0]-P_B[j][0]);helperY = (P_A[1]-P_B[j][1]);tmpDist = helperX*helperX + helperY*helperY;//tmpDist = pow(helperX,2) + pow(helperY,2);


The pow-function* from "math.h" is even slower than "normal" multiplication.

This topic is closed to new replies.

Advertisement