[postgis-users] distance function only supports 2d?

[David Blasby]

We could probably construct a distance3d() function that would work on 
point and lines, but it will be extreamly difficult to do correctly for 
polygons.

We'll need the following functions:

a) 3d point to 3d point distance (easy)
b) 3d point to 3d line segement distance (not too difficult)
c) 3d line to 3d line distance (hopefuly we can find this on the internet)

Here's the code for (b) and (c) in the 2d case.  The 3d case will be 
much much more difficult:



//distance from p to line A->B
double distance_pt_seg(POINT3D *p, POINT3D *A, POINT3D *B)
{
	double	r,s;


	//if start==end, then use pt distance
	if (  ( A->x == B->x) && (A->y == B->y) )
		return distance_pt_pt(p,A);

	//otherwise, we use comp.graphics.algorithms Frequently Asked Questions 
method

	/*(1)     	      AC dot AB
                    r = ---------
                          ||AB||^2
		r has the following meaning:
		r=0 P = A
		r=1 P = B
		r<0 P is on the backward extension of AB
		r>1 P is on the forward extension of AB
		0<r<1 P is interior to AB
	*/

	r = ( (p->x-A->x) * (B->x-A->x) + (p->y-A->y) * (B->y-A->y) )/( 
(B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) );

	if (r<0)
		return (distance_pt_pt(p,A));
	if (r>1)
		return(distance_pt_pt(p,B));


	/*(2)
		     (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
		s = -----------------------------
		             	L^2

		Then the distance from C to P = |s|*L.

	*/

	s = ((A->y-p->y)*(B->x-A->x)- (A->x-p->x)*(B->y-A->y) )/ 
((B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) );

	return abs(s) * sqrt(((B->x-A->x)*(B->x-A->x) +(B->y-A->y)*(B->y-A->y) ));
}



// find the distance from AB to CD
double distance_seg_seg(POINT3D *A, POINT3D *B, POINT3D *C, POINT3D *D)
{

	double	s_top, s_bot,s;
	double	r_top, r_bot,r;


//printf("seg_seg [%g,%g]->[%g,%g]  by [%g,%g]->[%g,%g] 
\n",A->x,A->y,B->x,B->y, C->x,C->y, D->x, D->y);
		//A and B are the same point

	if (  ( A->x == B->x) && (A->y == B->y) )
		return distance_pt_seg(A,C,D);

		//U and V are the same point

	if (  ( C->x == D->x) && (C->y == D->y) )
		return distance_pt_seg(D,A,B);

	// AB and CD are line segments
	/* from comp.graphics.algo

	Solving the above for r and s yields
				(Ay-Cy)(Dx-Cx)-(Ax-Cx)(Dy-Cy)
	           r = ----------------------------- (eqn 1)
				(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)

		 	(Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
		s = ----------------------------- (eqn 2)
			(Bx-Ax)(Dy-Cy)-(By-Ay)(Dx-Cx)
	Let P be the position vector of the intersection point, then
		P=A+r(B-A) or
		Px=Ax+r(Bx-Ax)
		Py=Ay+r(By-Ay)
	By examining the values of r & s, you can also determine some other 
limiting conditions:
		If 0<=r<=1 & 0<=s<=1, intersection exists
		r<0 or r>1 or s<0 or s>1 line segments do not intersect
		If the denominator in eqn 1 is zero, AB & CD are parallel
		If the numerator in eqn 1 is also zero, AB & CD are collinear.

	*/
	r_top = (A->y-C->y)*(D->x-C->x) - (A->x-C->x)*(D->y-C->y) ;
	r_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x) ;

	s_top = (A->y-C->y)*(B->x-A->x) - (A->x-C->x)*(B->y-A->y);
	s_bot = (B->x-A->x)*(D->y-C->y) - (B->y-A->y)*(D->x-C->x);

	if  ( (r_bot==0) || (s_bot == 0) )
	{
		return (
			min(distance_pt_seg(A,C,D),
				min(distance_pt_seg(B,C,D),
					min(distance_pt_seg(C,A,B),
						distance_pt_seg(D,A,B)    ) ) )
			 );
	}
	s = s_top/s_bot;
	r=  r_top/r_bot;

	if ((r<0) || (r>1) || (s<0) || (s>1) )
	{
		//no intersection
		return (
			min(distance_pt_seg(A,C,D),
				min(distance_pt_seg(B,C,D),
					min(distance_pt_seg(C,A,B),
						distance_pt_seg(D,A,B)    ) ) )
			 );

	}
	else
		return -0; //intersection exists

}