Vlib: add Vect_line_geodesic_distance()

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@61978 15284696-431f-4ddb-bdfa-cd5b030d7da7

lib/vector/Vlib/line.c

@@ -630,7 +630,7 @@ double Vect_line_geodesic_length(const struct line_pnts *Points)
 }
 
 /*!
-  \brief Calculate line distance.
+  \brief Calculate distance of point to line.
   
   Sets (if not null):
    - px, py - point on line,
@@ -750,6 +750,140 @@ int Vect_line_distance(const struct line_pnts *points,
     return (segment);
 }
 
+/*!
+  \brief Calculate geodesic distance of point to line in meters.
+  
+  Sets (if not null):
+   - px, py - point on line,
+   - dist   - distance to line,
+   - spdist - distance to point on line from segment beginning,
+   - sldist - distance to point on line form line beginning along line
+
+  \param points pointer to line_pnts structure
+  \param ux,uy,uz point coordinates
+  \param with_z flag if to use z coordinate (3D calculation)
+  \param[out] px,py,pz point on line
+  \param[out] dist distance to line,
+  \param[out] spdist distance of point from segment beginning
+  \param[out] lpdist distance of point from line
+
+  \return nearest segment (first is 1)
+ */
+int Vect_line_geodesic_distance(const struct line_pnts *points,
+		                double ux, double uy, double uz,
+		                int with_z,
+		                double *px, double *py, double *pz,
+		                double *dist, double *spdist,
+				double *lpdist)
+{
+    int i;
+    double distance;
+    double new_dist;
+    double tpx, tpy, tpz, ttpx, ttpy, ttpz; 
+    double tdist, tspdist, tlpdist = 0, tlpdistseg;
+    double dz;
+    int n_points;
+    int segment;
+
+    G_begin_distance_calculations();
+
+    n_points = points->n_points;
+
+    if (n_points == 1) {
+	distance = G_distance(ux, uy, points->x[0], points->y[0]);
+	if (with_z)
+	    distance = hypot(distance, uz - points->z[0]);
+
+	tpx = points->x[0];
+	tpy = points->y[0];
+	tpz = points->z[0];
+	tdist = distance;
+	tspdist = 0;
+	tlpdist = 0;
+	segment = 0;
+
+    }
+    else {
+
+	distance =
+	    dig_distance2_point_to_line(ux, uy, uz, points->x[0],
+					points->y[0], points->z[0],
+					points->x[1], points->y[1],
+					points->z[1], with_z, 
+					&tpx, &tpy, &tpz,
+					NULL, NULL);
+
+	distance = G_distance(ux, uy, tpx, tpy);
+	if (with_z)
+	    distance = hypot(distance, uz - tpz);
+
+	segment = 1;
+
+	for (i = 1; i < n_points - 1; i++) {
+	    new_dist = dig_distance2_point_to_line(ux, uy, uz,
+						   points->x[i], points->y[i],
+						   points->z[i],
+						   points->x[i + 1],
+						   points->y[i + 1],
+						   points->z[i + 1], with_z,
+						   &ttpx, &ttpy, &ttpz,
+						   NULL, NULL);
+	    new_dist = G_distance(ux, uy, ttpx, ttpy);
+	    if (with_z)
+		new_dist = hypot(new_dist, uz - ttpz);
+
+	    if (new_dist < distance) {
+		distance = new_dist;
+		segment = i + 1;
+		tpx = ttpx;
+		tpy = ttpy;
+		tpz = ttpz;
+	    }
+	}
+
+	/* calculate distance from beginning of segment */
+	tspdist = G_distance(points->x[segment - 1], points->y[segment - 1],
+			     tpx, tpy);
+	if (with_z) {
+	    dz = points->z[segment - 1] - tpz;
+	    tspdist += hypot(tspdist, dz);
+	}
+
+	/* calculate distance from beginning of line */
+	if (lpdist) {
+	    tlpdist = 0;
+	    for (i = 0; i < segment - 1; i++) {
+		tlpdistseg = G_distance(points->x[i], points->y[i],
+		                        points->x[i + 1], points->y[i + 1]);
+
+		if (with_z) {
+		    dz = points->z[i + 1] - points->z[i];
+		    tlpdistseg += hypot(tlpdistseg, dz);
+		}
+
+		tlpdist += tlpdistseg;
+	    }
+	    tlpdist += tspdist;
+	}
+	tdist = distance;
+    }
+
+    if (px)
+	*px = tpx;
+    if (py)
+	*py = tpy;
+    if (pz && with_z)
+	*pz = tpz;
+    if (dist)
+	*dist = tdist;
+    if (spdist)
+	*spdist = tspdist;
+    if (lpdist)
+	*lpdist = tlpdist;
+
+    return (segment);
+}
+
 
 /*!
   \brief Calculate distance of 2 points