SplitTris/splittris.c

   1 /* splittris.c -- read in a .ele/.node file pair generated by the
   2  *                triangle program and output a simple Wavefront .obj
   3  *                file for the north, south, east, and west edge
   4  *                verticies ... including the normals.
   5  *
   6  * Written by Curtis Olson, started January 1998.
   7  *
   8  * Copyright (C) 1997  Curtis L. Olson  - curt@infoplane.com
   9  *
  10  * This program is free software; you can redistribute it and/or modify
  11  * it under the terms of the GNU General Public License as published by
  12  * the Free Software Foundation; either version 2 of the License, or
  13  * (at your option) any later version.
  14  *
  15  * This program is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  18  * GNU General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU General Public License
  21  * along with this program; if not, write to the Free Software
  22  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  23  *
  24  * $Id$
  25  * (Log is kept at end of this file) */
  26
  27
  28 #include <math.h>
  29 #include <stdio.h>
  30 #include <stdlib.h>   /* for atoi() */
  31 #include <string.h>
  32 #include <sys/stat.h> /* for stat() */
  33 #include <unistd.h>   /* for stat() */
  34
  35 #include "splittris.h"
  36
  37 #include "../../Src/Include/constants.h"
  38 #include "../../Src/Include/types.h"
  39 #include "../../Src/Math/fg_geodesy.h"
  40 #include "../../Src/Math/mat3.h"
  41 #include "../../Src/Math/polar.h"
  42 #include "../../Src/Scenery/tileutils.h"
  43
  44
  45 int nodecount, tricount;
  46 double xmin, xmax, ymin, ymax;
  47
  48 double nodes_orig[MAX_NODES][3];
  49 int tris[MAX_TRIS][3];
  50 int new_tris[MAX_TRIS][3];
  51
  52 struct fgCartesianPoint nodes_cart[MAX_NODES];
  53
  54 long int ne_index, nw_index, sw_index, se_index;
  55 long int north_index, south_index, east_index, west_index;
  56
  57 /* convert a geodetic point lon(arcsec), lat(arcsec), elev(meter) to
  58  * a cartesian point */
  59 struct fgCartesianPoint geod_to_cart(double geod[3]) {
  60     struct fgCartesianPoint p;
  61     double gc_lon, gc_lat, sl_radius;
  62
  63     /* printf("A geodetic point is (%.2f, %.2f, %.2f)\n",
  64            geod[0], geod[1], geod[2]); */
  65
  66     gc_lon = geod[0]*ARCSEC_TO_RAD;
  67     fgGeodToGeoc(geod[1]*ARCSEC_TO_RAD, geod[2], &sl_radius, &gc_lat);
  68
  69     /* printf("A geocentric point is (%.2f, %.2f, %.2f)\n", gc_lon,
  70            gc_lat, sl_radius+geod[2]); */
  71
  72     p = fgPolarToCart(gc_lon, gc_lat, sl_radius+geod[2]);
  73
  74     /* printf("A cart point is (%.8f, %.8f, %.8f)\n", p.x, p.y, p.z); */
  75
  76     return(p);
  77 }
  78
  79
  80 /* given three points defining a triangle, calculate the normal */
  81 void calc_normal(struct fgCartesianPoint p1, struct fgCartesianPoint p2,
  82                  struct fgCartesianPoint p3, double normal[3])
  83 {
  84     double v1[3], v2[3];
  85     float temp;
  86
  87     v1[0] = p2.x - p1.x; v1[1] = p2.y - p1.y; v1[2] = p2.z - p1.z;
  88     v2[0] = p3.x - p1.x; v2[1] = p3.y - p1.y; v2[2] = p3.z - p1.z;
  89
  90     MAT3cross_product(normal, v1, v2);
  91     MAT3_NORMALIZE_VEC(normal,temp);
  92
  93 /*  printf("  Normal = %.2f %.2f %.2f\n", normal[0], normal[1], normal[2]); */
  94 }
  95
  96
  97 /* return the file base name ( foo/bar/file.ext = file.ext ) */
  98 void extract_file(char *in, char *base) {
  99     int len, i;
 100
 101     len = strlen(in);
 102
 103     i = len - 1;
 104     while ( (i >= 0) && (in[i] != '/') ) {
 105         i--;
 106     }
 107
 108     in += (i + 1);
 109     strcpy(base, in);
 110 }
 111
 112
 113 /* return the file path name ( foo/bar/file.ext = foo/bar ) */
 114 void extract_path(char *in, char *base) {
 115     int len, i;
 116
 117     len = strlen(in);
 118     strcpy(base, in);
 119
 120     i = len - 1;
 121     while ( (i >= 0) && (in[i] != '/') ) {
 122         i--;
 123     }
 124
 125     base[i] = '\0';
 126 }
 127
 128
 129 /* return the index of all triangles containing the specified node */
 130 void find_tris(int n, int *t1, int *t2, int *t3, int *t4, int *t5) {
 131     int i;
 132
 133     *t1 = *t2 = *t3 = *t4 = *t5 = 0;
 134
 135     i = 1;
 136     while ( i <= tricount ) {
 137         if ( (n == tris[i][0]) || (n == tris[i][1]) || (n == tris[i][2]) ) {
 138             if ( *t1 == 0 ) {
 139                 *t1 = i;
 140             } else if ( *t2 == 0 ) {
 141                 *t2 = i;
 142             } else if ( *t3 == 0 ) {
 143                 *t3 = i;
 144             } else if ( *t4 == 0 ) {
 145                 *t4 = i;
 146             } else {
 147                 *t5 = i;
 148             }
 149         }
 150         i++;
 151     }
 152 }
 153
 154
 155 /* Initialize a new mesh structure */
 156 void triload(char *basename) {
 157     char nodename[256], elename[256];
 158     FILE *node, *ele;
 159     int dim, junk1, junk2;
 160     int i;
 161
 162     strcpy(nodename, basename);
 163     strcat(nodename, ".node");
 164     strcpy(elename, basename);
 165     strcat(elename, ".ele");
 166
 167     printf("Loading node file:  %s ...\n", nodename);
 168     if ( (node = fopen(nodename, "r")) == NULL ) {
 169         printf("Cannot open file '%s'\n", nodename);
 170         exit(-1);
 171     }
 172
 173     fscanf(node, "%d %d %d %d", &nodecount, &dim, &junk1, &junk2);
 174
 175     if ( nodecount > MAX_NODES - 1 ) {
 176         printf("Error, too many nodes, need to increase array size\n");
 177         exit(-1);
 178     } else {
 179         printf("    Expecting %d nodes\n", nodecount);
 180     }
 181
 182     for ( i = 1; i <= nodecount; i++ ) {
 183         fscanf(node, "%d %lf %lf %lf %d\n", &junk1,
 184                &nodes_orig[i][0], &nodes_orig[i][1], &nodes_orig[i][2], &junk2);
 185         /* printf("%d %.2f %.2f %.2f\n", junk1, n[0], n[1], n[2]); */
 186         nodes_cart[i] = geod_to_cart(nodes_orig[i]);
 187         /* printf("%d %.2f %.2f %.2f\n",
 188                junk1, nodes_cart[i].x, nodes_cart[i].y, nodes_cart[i].z); */
 189
 190         if ( i == 1 ) {
 191             xmin = xmax = nodes_orig[i][0];
 192             ymin = ymax = nodes_orig[i][1];
 193         } else {
 194             if ( nodes_orig[i][0] < xmin ) {
 195                 xmin = nodes_orig[i][0];
 196             }
 197             if ( nodes_orig[i][0] > xmax ) {
 198                 xmax = nodes_orig[i][0];
 199             }
 200             if ( nodes_orig[i][1] < ymin ) {
 201                 ymin = nodes_orig[i][1];
 202             }
 203             if ( nodes_orig[i][1] > ymax ) {
 204                 ymax = nodes_orig[i][1];
 205             }
 206         }
 207     }
 208
 209     fclose(node);
 210
 211     printf("Loading element file:  %s ...\n", elename);
 212     if ( (ele = fopen(elename, "r")) == NULL ) {
 213         printf("Cannot open file '%s'\n", elename);
 214         exit(-1);
 215     }
 216
 217     fscanf(ele, "%d %d %d", &tricount, &junk1, &junk2);
 218
 219     if ( tricount > MAX_TRIS - 1 ) {
 220         printf("Error, too many elements, need to increase array size\n");
 221         exit(-1);
 222     } else {
 223         printf("    Expecting %d elements\n", tricount);
 224     }
 225
 226     for ( i = 1; i <= tricount; i++ ) {
 227         fscanf(ele, "%d %d %d %d\n", &junk1,
 228                &tris[i][0], &tris[i][1], &tris[i][2]);
 229         /* printf("%d %d %d %d\n", junk1, tris[i][0], tris[i][1], tris[i][2]);*/
 230     }
 231
 232     fclose(ele);
 233 }
 234
 235
 236 /* check if a file exists */
 237 int file_exists(char *file) {
 238     struct stat stat_buf;
 239
 240     result = stat(file, &stat_buf);
 241
 242     if ( result != 0 ) {
 243         /* stat failed, no file */
 244         return(0);
 245     } else {
 246         /* stat succeeded, file exists */
 247         return(1);
 248     }
 249 }
 250
 251
 252 /* my custom file opening routine ... don't open if a shared edge or
 253  * vertex alread exists */
 254 FILE *my_open(char *basename, char *basepath, char *ext) {
 255     FILE *fp;
 256     char filename[256];
 257
 258     /* create the output file name */
 259     strcpy(filename, basename);
 260     strcpy(filename, ext);
 261
 262     /* check if a shared object already exist from a different tile */
 263
 264     if ( 0 ) {
 265         /* not an actual file open error, but we've already got the
 266          * shared edge, so we don't want to create another one */
 267         return(NULL);
 268     } else {
 269         /* open the file */
 270         fp = fopen(filename, "w");
 271         return(fp);
 272     }
 273 }
 274
 275
 276 /* dump in WaveFront .obj format */
 277 void dump_obj(char *basename, char *basepath) {
 278     char sw_name[256], se_name[256], ne_name[256], nw_name[256];
 279     char north_name[256], south_name[256], east_name[256], west_name[256];
 280     char body_name[256];
 281     double n1[3], n2[3], n3[3], n4[3], n5[3], norm[3], temp;
 282     FILE *sw, *se, *ne, *nw, *north, *south, *east, *west, *body;
 283     int i, t1, t2, t3, t4, t5, count;
 284
 285     sw = my_open(basename, basepath, ".sw");
 286     se = my_open(basename, basepath, ".se");
 287     ne = my_open(basename, basepath, ".ne");
 288     nw = my_open(basename, basepath, ".nw");
 289
 290     north = my_open(basename, basepath, ".north");
 291     south = my_open(basename, basepath, ".south");
 292     east = my_open(basename, basepath, ".east");
 293     west = my_open(basename, basepath, ".west");
 294
 295     body = my_open(basename, basepath, ".body");
 296
 297     printf("Dumping edges file basename:  %s ...\n", basename);
 298
 299     sw = fopen(sw_name, "w");
 300     se = fopen(se_name, "w");
 301     ne = fopen(ne_name, "w");
 302     nw = fopen(nw_name, "w");
 303
 304     north = fopen(north_name, "w");
 305     south = fopen(south_name, "w");
 306     east = fopen(east_name, "w");
 307     west = fopen(west_name, "w");
 308
 309     body = fopen(body_name, "w");
 310
 311     /* dump vertices */
 312     printf("  writing vertices\n");
 313     for ( i = 1; i <= nodecount; i++ ) {
 314
 315         if ( (fabs(nodes_orig[i][1] - ymin) < FG_EPSILON) &&
 316              (fabs(nodes_orig[i][0] - xmin) < FG_EPSILON) ) {
 317             fprintf(sw, "geodn %.2f %.2f %.2f\n",
 318                     nodes_orig[i][0], nodes_orig[i][1], nodes_orig[i][2]);
 319         } else if ( (fabs(nodes_orig[i][1] - ymin) < FG_EPSILON) &&
 320                     (fabs(nodes_orig[i][0] - xmax) < FG_EPSILON) ) {
 321             fprintf(se, "geodn %.2f %.2f %.2f\n",
 322                     nodes_orig[i][0], nodes_orig[i][1], nodes_orig[i][2]);
 323         } else if ( (fabs(nodes_orig[i][1] - ymax) < FG_EPSILON) &&
 324                     (fabs(nodes_orig[i][0] - xmax) < FG_EPSILON)) {
 325             fprintf(ne, "geodn %.2f %.2f %.2f\n",
 326                     nodes_orig[i][0], nodes_orig[i][1], nodes_orig[i][2]);
 327         } else if ( (fabs(nodes_orig[i][1] - ymax) < FG_EPSILON) &&
 328                     (fabs(nodes_orig[i][0] - xmin) < FG_EPSILON) ) {
 329             fprintf(nw, "geodn %.2f %.2f %.2f\n",
 330                     nodes_orig[i][0], nodes_orig[i][1], nodes_orig[i][2]);
 331         } else if ( fabs(nodes_orig[i][0] - xmin) < FG_EPSILON ) {
 332             fprintf(west, "geodn %.2f %.2f %.2f\n",
 333                     nodes_orig[i][0], nodes_orig[i][1], nodes_orig[i][2]);
 334         } else if ( fabs(nodes_orig[i][0] - xmax) < FG_EPSILON ) {
 335             fprintf(east, "geodn %.2f %.2f %.2f\n",
 336                     nodes_orig[i][0], nodes_orig[i][1], nodes_orig[i][2]);
 337         } else if ( fabs(nodes_orig[i][1] - ymin) < FG_EPSILON ) {
 338             fprintf(south, "geodn %.2f %.2f %.2f\n",
 339                     nodes_orig[i][0], nodes_orig[i][1], nodes_orig[i][2]);
 340         } else if ( fabs(nodes_orig[i][1] - ymax) < FG_EPSILON ) {
 341             fprintf(north, "geodn %.2f %.2f %.2f\n",
 342                     nodes_orig[i][0], nodes_orig[i][1], nodes_orig[i][2]);
 343         } else {
 344             fprintf(body, "geodn %.2f %.2f %.2f\n",
 345                     nodes_orig[i][0], nodes_orig[i][1], nodes_orig[i][2]);
 346         }
 347
 348     }
 349
 350     printf("  calculating and writing normals\n");
 351     /* calculate and generate normals */
 352     for ( i = 1; i <= nodecount; i++ ) {
 353 /*      printf("Finding normal\n"); */
 354
 355         find_tris(i, &t1, &t2, &t3, &t4, &t5);
 356
 357         n1[0] = n1[1] = n1[2] = 0.0;
 358         n2[0] = n2[1] = n2[2] = 0.0;
 359         n3[0] = n3[1] = n3[2] = 0.0;
 360         n4[0] = n4[1] = n4[2] = 0.0;
 361         n5[0] = n5[1] = n5[2] = 0.0;
 362
 363         count = 1;
 364         calc_normal(nodes_cart[tris[t1][0]], nodes_cart[tris[t1][1]],
 365                     nodes_cart[tris[t1][2]], n1);
 366
 367         if ( t2 > 0 ) {
 368             calc_normal(nodes_cart[tris[t2][0]], nodes_cart[tris[t2][1]],
 369                         nodes_cart[tris[t2][2]], n2);
 370             count = 2;
 371         }
 372
 373         if ( t3 > 0 ) {
 374             calc_normal(nodes_cart[tris[t3][0]], nodes_cart[tris[t3][1]],
 375                         nodes_cart[tris[t3][2]], n3);
 376             count = 3;
 377         }
 378
 379         if ( t4 > 0 ) {
 380             calc_normal(nodes_cart[tris[t4][0]], nodes_cart[tris[t4][1]],
 381                         nodes_cart[tris[t4][2]], n4);
 382             count = 4;
 383         }
 384
 385         if ( t5 > 0 ) {
 386             calc_normal(nodes_cart[tris[t5][0]], nodes_cart[tris[t5][1]],
 387                         nodes_cart[tris[t5][2]], n5);
 388             count = 5;
 389         }
 390
 391 /*      printf("  norm[2] = %.2f %.2f %.2f\n", n1[2], n2[2], n3[2]); */
 392
 393         norm[0] = ( n1[0] + n2[0] + n3[0] + n4[0] + n5[0] ) / (double)count;
 394         norm[1] = ( n1[1] + n2[1] + n3[1] + n4[1] + n5[1] ) / (double)count;
 395         norm[2] = ( n1[2] + n2[2] + n3[2] + n4[2] + n5[2] ) / (double)count;
 396
 397 /*      printf("  count = %d\n", count); */
 398 /*      printf("  Ave. normal = %.4f %.4f %.4f\n", norm[0], norm[1], norm[2]);*/
 399         MAT3_NORMALIZE_VEC(norm, temp);
 400 /*      printf("  Normalized ave. normal = %.4f %.4f %.4f\n",  */
 401 /*             norm[0], norm[1], norm[2]); */
 402
 403         if ( (fabs(nodes_orig[i][1] - ymin) < FG_EPSILON) &&
 404              (fabs(nodes_orig[i][0] - xmin) < FG_EPSILON) ) {
 405             fprintf(sw, "vn %.4f %.4f %.4f\n", norm[0], norm[1], norm[2]);
 406         } else if ( (fabs(nodes_orig[i][1] - ymin) < FG_EPSILON) &&
 407                     (fabs(nodes_orig[i][0] - xmax) < FG_EPSILON) ) {
 408             fprintf(se, "vn %.4f %.4f %.4f\n", norm[0], norm[1], norm[2]);
 409         } else if ( (fabs(nodes_orig[i][1] - ymax) < FG_EPSILON) &&
 410                     (fabs(nodes_orig[i][0] - xmax) < FG_EPSILON)) {
 411             fprintf(ne, "vn %.4f %.4f %.4f\n", norm[0], norm[1], norm[2]);
 412         } else if ( (fabs(nodes_orig[i][1] - ymax) < FG_EPSILON) &&
 413                     (fabs(nodes_orig[i][0] - xmin) < FG_EPSILON) ) {
 414             fprintf(nw, "vn %.4f %.4f %.4f\n", norm[0], norm[1], norm[2]);
 415         } else if ( fabs(nodes_orig[i][0] - xmin) < FG_EPSILON ) {
 416             fprintf(west, "vn %.4f %.4f %.4f\n", norm[0], norm[1], norm[2]);
 417         } else if ( fabs(nodes_orig[i][0] - xmax) < FG_EPSILON ) {
 418             fprintf(east, "vn %.4f %.4f %.4f\n", norm[0], norm[1], norm[2]);
 419         } else if ( fabs(nodes_orig[i][1] - ymin) < FG_EPSILON ) {
 420             fprintf(south, "vn %.4f %.4f %.4f\n", norm[0], norm[1], norm[2]);
 421         } else if ( fabs(nodes_orig[i][1] - ymax) < FG_EPSILON ) {
 422             fprintf(north, "vn %.4f %.4f %.4f\n", norm[0], norm[1], norm[2]);
 423         }
 424     }
 425
 426     fclose(sw);
 427     fclose(se);
 428     fclose(ne);
 429     fclose(nw);
 430
 431     fclose(north);
 432     fclose(south);
 433     fclose(east);
 434     fclose(west);
 435 }
 436
 437
 438 int main(int argc, char **argv) {
 439     char basename[256], basepath[256], temp[256];
 440     struct bucket p1, p2;
 441     long int index;
 442     int len;
 443
 444     strcpy(basename, argv[1]);
 445
 446     /* find the base path of the file */
 447     extract_path(basename, basepath);
 448     extract_path(basepath, basepath);
 449     extract_path(basepath, basepath);
 450     printf("%s\n", basepath);
 451
 452     /* find the index of the current file */
 453     extract_file(basename, temp);
 454     len = strlen(temp);
 455     if ( len >= 2 ) {
 456         temp[len-2] = '\0';
 457     }
 458     index = atoi(temp);
 459     printf("%ld\n", index);
 460     parse_index(index, &p1);
 461
 462     /* generate the indexes of the neighbors */
 463     offset_bucket(&p1, &p2,  1,  1); ne_index = gen_index(&p2);
 464     offset_bucket(&p1, &p2,  1, -1); nw_index = gen_index(&p2);
 465     offset_bucket(&p1, &p2, -1,  1); se_index = gen_index(&p2);
 466     offset_bucket(&p1, &p2, -1, -1); sw_index = gen_index(&p2);
 467
 468     offset_bucket(&p1, &p2,  0,  1); north_index = gen_index(&p2);
 469     offset_bucket(&p1, &p2,  0, -1); south_index = gen_index(&p2);
 470     offset_bucket(&p1, &p2,  1,  0); east_index  = gen_index(&p2);
 471     offset_bucket(&p1, &p2, -1,  1); west_index  = gen_index(&p2);
 472
 473     printf("Corner indexes = %ld %ld %ld %ld\n",
 474            ne_index, nw_index, sw_index, se_index);
 475     printf("Edge indexes = %ld %ld %ld %ld\n",
 476            north_index, south_index, east_index, west_index);
 477
 478     /* load the input data files */
 479     triload(basename);
 480
 481     /* dump in WaveFront .obj format */
 482     dump_obj(basename, basepath);
 483
 484     return(0);
 485 }
 486
 487
 488 /* $Log$
 489 /* Revision 1.1  1998/01/14 02:11:31  curt
 490 /* Initial revision.
 491 /*
 492  */