tin.c
Go to the documentation of this file.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
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
00090 *tz = -(d + a * tx + b * ty) / c;
00091 G_debug(3, "TIN: z = %f", *tz);
00092
00093 return 1;
00094 }