GRASS Programmer's Manual  6.4.2(2012)
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     /* 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 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines