[vtkusers] Re: vtkusers Digest, Vol 20, Issue 12

[L.J. van Ruijven]

Hello vtk developers,

When tried to replace a part of my program with a call to 
vtkTriangle::EvaluatePosition, this resulted in bad results. So I tried 
to find the differences between the VTK code and my program. I still 
use VTK 4.2 so I am not sure whether the following still applies.
The function vtkTriangle::EvaluatePosition uses 
vtkLine::DistanceToLine. This function seems strange, because it does 
not always assign a value to the variable t. So, in some cases the 
return value of t is simply the value it had before the function was 
called. Also, when a value is assigned, this value is not always in the 
interval [0,1]. I cannot oversee all the consequences when this 
behavior is changed, but if it is changed the function 
vtkTriangle::EvaluatePosition can be made much shorter.
Also the precision of vtkLine::DistanceToLine can be improved very 
simply for the case where all the points are located far away from the 
origin. My suggestion for vtkLine::DistanceToLine is:

float NewDistanceToLine(float x[3], float p1[3], float p2[3], 
                              float &t, float closestPoint[3])
{
  float p21[3], x1[3], denom, num;
  float x0[3] = {0,0,0};
  float *closest;
  float tolerance;
  //
  //   Determine appropriate vectors
  // 
  x1[0] = x[0]- p1[0];
  x1[1] = x[1]- p1[1];
  x1[2] = x[2]- p1[2];
  p21[0] = p2[0]- p1[0];
  p21[1] = p2[1]- p1[1];
  p21[2] = p2[2]- p1[2];

  //
  //   Get parametric location
  //
  num = vtkMath::Dot(p21,x1);
  denom = vtkMath::Dot(p21,p21);

  // trying to avoid an expensive fabs
  tolerance = VTK_TOL*num;
  if (tolerance < 0.0)
    {
    tolerance = -tolerance;
    }
  if ( -tolerance < denom && denom < tolerance ) //numerically bad!
    {
    closest = x0; //arbitrary, point is (numerically) far away
    t = 0.0;
    }
  //
  // If parametric coordinate is within 0<=p<=1, then the point is 
closest to
  // the line.  Otherwise, it's closest to a point at the end of the 
line.
  //
  else if ( (t=num/denom) < 0.0 )
    {
    closest = x0;
    t = 0.0;
    }
  else if ( t > 1.0 )
    {
    closest = p21;
    t = 1.0;
    }
  else
    {
    closest = p21;
    p21[0] *= t;
    p21[1] *= t;
    p21[2] *= t;
    }

  closestPoint[0] = p1[0] + closest[0];
  closestPoint[1] = p1[1] + closest[1];
  closestPoint[2] = p1[2] + closest[2];
  return vtkMath::Distance2BetweenPoints(closest,x1);
}

The effect on the precision can be seen with the program:

int _tmain(int argc, _TCHAR* argv[])
{
	int i,j;
	float x[3]={0.1972,0.1020,0.3015};	// point coordinates
	float p1[3]={0.1002,0.1999,0.4010};	// endpoints of the line
	float p2[3]={0.2000,0.1000,0.3000};

	const int N = 200;					// call 
DistanceToLine N times
	float Closest_old[3],Closest_new[3];
	float t;
	double d_old[N],d_new[N];			// the 
distances calculated from the points x[i] to the line
	double mean_old=0,mean_new=0;		// and their means

	for (i=0;i<N;i++) {
		float xi[3],p1i[3],p2i[3];
		float d[3];
		for (j=0;j<3;j++) {	// translate the points over a 
large distance
			d[j]= rand()%50;
			xi[j]  = x[j]+d[j];
			p1i[j] = p1[j]+d[j];
			p2i[j] = p2[j]+d[j];
		}

		float dis1 = vtkLine::DistanceToLine
(xi,p1i,p2i,t,Closest_old);
		float dis2 = DistanceToLine(xi,p1i,p2i,t,Closest_new);

		d_old[i] = sqrt(dis1);	// store the results of both 
functions
		d_new[i] = sqrt(dis2);

		mean_old+=d_old[i];
		mean_new+=d_new[i];
	}
	mean_old/=N; mean_new/=N;

// Now calculate the standard deviations of the distances
	double square_old=0, square_new=0;
	for (int i=0;i<N;i++) {
		d_old[i]-=mean_old; square_old+=d_old[i]*d_old[i];
		d_new[i]-=mean_new; square_new+=d_new[i]*d_new[i];
		printf("%4d: olddist - mean = %12.8f, newdist - mean = %
12.8f\n",i,d_old[i],d_new[i]);
	}
	float SD_old = sqrt(square_old/N);
	float SD_new = sqrt(square_new/N);

	printf("\n");
	printf("       Mean Dist   ( SD Dist   ) :\n");
	printf("old: %12.8f (%12.8f)\n",mean_old,SD_old);
	printf("new: %12.8f (%12.8f)\n",mean_new,SD_new);

	return 0;
}

If you decide to change vtkLine::DistanceToLine, please let me know. 
Then I will also propose a new version for 
vtkTriangle::EvaluatePosition.

Leo van Ruijven.

-------------- next part --------------
A non-text attachment was scrubbed...
Name: ljvanruijven.vcf
Type: text/x-vcard
Size: 262 bytes
Desc: Card for "L.J. van Ruijven" <L.J.vanRuijven at amc.uva.nl>
URL: <http://www.vtk.org/pipermail/vtkusers/attachments/20051212/8df09074/attachment.vcf>