GRASS Programmer's Manual 6.4.1(2011)
|
00001 00020 #include <grass/Vect.h> 00021 00035 int 00036 Vect_tin_get_z(struct Map_info *Map, 00037 double tx, double ty, double *tz, double *angle, double *slope) 00038 { 00039 int i, area, n_points; 00040 struct Plus_head *Plus; 00041 P_AREA *Area; 00042 static struct line_pnts *Points; 00043 static int first_time = 1; 00044 double *x, *y, *z; 00045 double vx1, vx2, vy1, vy2, vz1, vz2; 00046 double a, b, c, d; 00047 00048 /* TODO angle, slope */ 00049 00050 Plus = &(Map->plus); 00051 if (first_time == 1) { 00052 Points = Vect_new_line_struct(); 00053 first_time = 0; 00054 } 00055 00056 area = Vect_find_area(Map, tx, ty); 00057 G_debug(3, "TIN: area = %d", area); 00058 if (area == 0) 00059 return 0; 00060 00061 Area = Plus->Area[area]; 00062 if (Area->n_isles > 0) 00063 return -1; 00064 00065 Vect_get_area_points(Map, area, Points); 00066 n_points = Points->n_points; 00067 if (n_points != 4) 00068 return -1; 00069 00070 x = Points->x; 00071 y = Points->y; 00072 z = Points->z; 00073 for (i = 0; i < 3; i++) { 00074 G_debug(3, "TIN: %d %f %f %f", i, x[i], y[i], z[i]); 00075 } 00076 00077 vx1 = x[1] - x[0]; 00078 vy1 = y[1] - y[0]; 00079 vz1 = z[1] - z[0]; 00080 vx2 = x[2] - x[0]; 00081 vy2 = y[2] - y[0]; 00082 vz2 = z[2] - z[0]; 00083 00084 a = vy1 * vz2 - vy2 * vz1; 00085 b = vz1 * vx2 - vz2 * vx1; 00086 c = vx1 * vy2 - vx2 * vy1; 00087 d = -a * x[0] - b * y[0] - c * z[0]; 00088 00089 /* OK ? */ 00090 *tz = -(d + a * tx + b * ty) / c; 00091 G_debug(3, "TIN: z = %f", *tz); 00092 00093 return 1; 00094 }