Triangulate/triangle.cxx

   1 // triangle.cxx -- "Triangle" interface class
   2 //
   3 // Written by Curtis Olson, started March 1999.
   4 //
   5 // Copyright (C) 1999  Curtis L. Olson  - curt@flightgear.org
   6 //
   7 // This program is free software; you can redistribute it and/or
   8 // modify it under the terms of the GNU General Public License as
   9 // published by the Free Software Foundation; either version 2 of the
  10 // License, or (at your option) any later version.
  11 //
  12 // This program is distributed in the hope that it will be useful, but
  13 // WITHOUT ANY WARRANTY; without even the implied warranty of
  14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  15 // General Public License for more details.
  16 //
  17 // You should have received a copy of the GNU General Public License
  18 // along with this program; if not, write to the Free Software
  19 // Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  20 //
  21 // $Id$
  22 // (Log is kept at end of this file)
  23
  24
  25 #include "triangle.hxx"
  26 #include "tripoly.hxx"
  27
  28 // Constructor
  29 FGTriangle::FGTriangle( void ) {
  30 }
  31
  32
  33 // Destructor
  34 FGTriangle::~FGTriangle( void ) {
  35 }
  36
  37
  38 // populate this class based on the specified gpc_polys list
  39 int
  40 FGTriangle::build( const fitnode_list& fit_list,
  41                    const FGgpcPolyList& gpc_polys )
  42 {
  43     FGTriPoly poly;
  44     int index;
  45
  46     // Point3D junkp;
  47     // int junkc = 0;
  48     // char junkn[256];
  49     // FILE *junkfp;
  50
  51     // traverse the dem fit list and gpc_polys building a unified node
  52     // list and converting the polygons so that they reference the
  53     // node list by index (starting at zero) rather than listing the
  54     // points explicitely
  55
  56     const_fitnode_list_iterator f_current, f_last;
  57     f_current = fit_list.begin();
  58     f_last = fit_list.end();
  59     for ( ; f_current != f_last; ++f_current ) {
  60         index = trinodes.unique_add( *f_current );
  61     }
  62
  63     gpc_polygon *gpc_poly;
  64     const_gpcpoly_iterator current, last;
  65
  66     // process polygons in priority order
  67     cout << "prepairing node list and polygons" << endl;
  68
  69     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
  70         cout << "area type = " << i << endl;
  71         current = gpc_polys.polys[i].begin();
  72         last = gpc_polys.polys[i].end();
  73         for ( ; current != last; ++current ) {
  74             gpc_poly = *current;
  75             cout << "processing a polygon, contours = "
  76                  << gpc_poly->num_contours << endl;
  77
  78             if (gpc_poly->num_contours <= 0 ) {
  79                 cout << "FATAL ERROR! no contours in this polygon" << endl;
  80                 exit(-1);
  81             }
  82
  83             if (gpc_poly->num_contours > 1 ) {
  84                 cout << "FATAL ERROR! no multi-contour support" << endl;
  85                 sleep(2);
  86                 // exit(-1);
  87             }
  88
  89             for ( int j = 0; j < gpc_poly->num_contours; j++ ) {
  90
  91                 poly.erase();
  92
  93                 // sprintf(junkn, "g.%d", junkc++);
  94                 // junkfp = fopen(junkn, "w");
  95
  96                 for ( int k = 0; k < gpc_poly->contour[j].num_vertices; k++ ) {
  97                     Point3D p( gpc_poly->contour[j].vertex[k].x,
  98                                gpc_poly->contour[j].vertex[k].y,
  99                                0 );
 100                     index = trinodes.unique_add( p );
 101                     // junkp = trinodes.get_node( index );
 102                     // fprintf(junkfp, "%.4f %.4f\n", junkp.x(), junkp.y());
 103                     poly.add_node(index);
 104                     // cout << index << endl;
 105                 }
 106                 // fprintf(junkfp, "%.4f %.4f\n",
 107                 //    gpc_poly->contour[j].vertex[0].x,
 108                 //    gpc_poly->contour[j].vertex[0].y);
 109                 // fclose(junkfp);
 110
 111                 poly.calc_point_inside( trinodes );
 112
 113                 polylist[i].push_back(poly);
 114             }
 115         }
 116     }
 117
 118     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 119         if ( polylist[i].size() ) {
 120             cout << get_area_name((AreaType)i) << " = "
 121                  << polylist[i].size() << endl;
 122         }
 123     }
 124
 125     // traverse the polygon lists and build the segment (edge) list
 126     // that is used by the "Triangle" lib.
 127
 128     int i1, i2;
 129     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 130         cout << "area type = " << i << endl;
 131         tripoly_list_iterator tp_current, tp_last;
 132         tp_current = polylist[i].begin();
 133         tp_last = polylist[i].end();
 134
 135         // process each polygon in list
 136         for ( ; tp_current != tp_last; ++tp_current ) {
 137             poly = *tp_current;
 138
 139             for ( int j = 0; j < (int)(poly.size()) - 1; ++j ) {
 140                 i1 = poly.get_pt_index( j );
 141                 i2 = poly.get_pt_index( j + 1 );
 142                 trisegs.unique_add( FGTriSeg(i1, i2) );
 143             }
 144             i1 = poly.get_pt_index( 0 );
 145             i2 = poly.get_pt_index( poly.size() - 1 );
 146             trisegs.unique_add( FGTriSeg(i1, i2) );
 147         }
 148     }
 149
 150     return 0;
 151 }
 152
 153
 154 static void write_out_data(struct triangulateio *out) {
 155     FILE *node = fopen("tile.node", "w");
 156     fprintf(node, "%d 2 %d 0\n",
 157             out->numberofpoints, out->numberofpointattributes);
 158     for (int i = 0; i < out->numberofpoints; i++) {
 159         fprintf(node, "%d %.6f %.6f %.2f\n",
 160                 i, out->pointlist[2*i], out->pointlist[2*i + 1], 0.0);
 161     }
 162     fclose(node);
 163
 164     FILE *ele = fopen("tile.ele", "w");
 165     fprintf(ele, "%d 3 0\n", out->numberoftriangles);
 166     for (int i = 0; i < out->numberoftriangles; i++) {
 167         fprintf(ele, "%d ", i);
 168         for (int j = 0; j < out->numberofcorners; j++) {
 169             fprintf(ele, "%d ", out->trianglelist[i * out->numberofcorners + j]);
 170         }
 171         for (int j = 0; j < out->numberoftriangleattributes; j++) {
 172             fprintf(ele, "%.6f ",
 173                     out->triangleattributelist[i
 174                                               * out->numberoftriangleattributes
 175                                               + j]
 176                     );
 177         }
 178         fprintf(ele, "\n");
 179     }
 180     fclose(ele);
 181
 182     FILE *fp = fopen("tile.poly", "w");
 183     fprintf(fp, "0 2 1 0\n");
 184     fprintf(fp, "%d 0\n", out->numberofsegments);
 185     for (int i = 0; i < out->numberofsegments; ++i) {
 186         fprintf(fp, "%d %d %d\n",
 187                 i, out->segmentlist[2*i], out->segmentlist[2*i + 1]);
 188     }
 189     fprintf(fp, "%d\n", out->numberofholes);
 190     for (int i = 0; i < out->numberofholes; i++) {
 191         fprintf(fp, "%d %.6f %.6f\n",
 192                 i, out->holelist[2*i], out->holelist[2*i + 1]);
 193     }
 194     fprintf(fp, "%d\n", out->numberofregions);
 195     for (int i = 0; i < out->numberofregions; i++) {
 196         fprintf(fp, "%d %.6f %.6f %.6f\n",
 197                 i, out->regionlist[4*i], out->regionlist[4*i + 1],
 198                 out->regionlist[4*i + 2]);
 199     }
 200 }
 201
 202
 203 // triangulate each of the polygon areas
 204 int FGTriangle::run_triangulate() {
 205     FGTriPoly poly;
 206     Point3D p;
 207     struct triangulateio in, out, vorout;
 208     int counter;
 209
 210     // point list
 211     trinode_list node_list = trinodes.get_node_list();
 212     in.numberofpoints = node_list.size();
 213     in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL));
 214
 215     trinode_list_iterator tn_current, tn_last;
 216     tn_current = node_list.begin();
 217     tn_last = node_list.end();
 218     counter = 0;
 219     for ( ; tn_current != tn_last; ++tn_current ) {
 220         in.pointlist[counter++] = tn_current->x();
 221         in.pointlist[counter++] = tn_current->y();
 222     }
 223
 224     in.numberofpointattributes = 1;
 225     in.pointattributelist = (REAL *) malloc(in.numberofpoints *
 226                                             in.numberofpointattributes *
 227                                             sizeof(REAL));
 228     for ( int i = 0; i < in.numberofpoints * in.numberofpointattributes; i++) {
 229         in.pointattributelist[i] = 0.0;
 230     }
 231
 232     in.pointmarkerlist = (int *) malloc(in.numberofpoints * sizeof(int));
 233     for ( int i = 0; i < in.numberofpoints; i++) {
 234         in.pointmarkerlist[i] = 0;
 235     }
 236
 237     // triangle list
 238     in.numberoftriangles = 0;
 239
 240     // segment list
 241     triseg_list seg_list = trisegs.get_seg_list();
 242     in.numberofsegments = seg_list.size();
 243     in.segmentlist = (int *) malloc(in.numberofsegments * 2 * sizeof(int));
 244
 245     triseg_list_iterator s_current, s_last;
 246     s_current = seg_list.begin();
 247     s_last = seg_list.end();
 248     counter = 0;
 249     for ( ; s_current != s_last; ++s_current ) {
 250         in.segmentlist[counter++] = s_current->get_n1();
 251         in.segmentlist[counter++] = s_current->get_n2();
 252     }
 253
 254     // hole list (make holes for airport ignore areas)
 255     in.numberofholes = polylist[(int)AirportIgnoreArea].size();
 256     in.holelist = (REAL *) malloc(in.numberofholes * 2 * sizeof(REAL));
 257
 258     tripoly_list_iterator h_current, h_last;
 259     h_current = polylist[(int)AirportIgnoreArea].begin();
 260     h_last = polylist[(int)AirportIgnoreArea].end();
 261     counter = 0;
 262     for ( ; h_current != h_last; ++h_current ) {
 263         poly = *h_current;
 264         p = poly.get_point_inside();
 265         in.holelist[counter++] = p.x();
 266         in.holelist[counter++] = p.y();
 267     }
 268
 269     // region list
 270     in.numberofregions = 0;
 271     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 272         in.numberofregions += polylist[i].size();
 273     }
 274
 275     in.regionlist = (REAL *) malloc(in.numberofregions * 4 * sizeof(REAL));
 276     counter = 0;
 277     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 278         tripoly_list_iterator h_current, h_last;
 279         h_current = polylist[(int)i].begin();
 280         h_last = polylist[(int)i].end();
 281         for ( ; h_current != h_last; ++h_current ) {
 282             poly = *h_current;
 283             p = poly.get_point_inside();
 284             in.regionlist[counter++] = p.x();  // x coord
 285             in.regionlist[counter++] = p.y();  // y coord
 286             in.regionlist[counter++] = i;      // region attribute
 287             in.regionlist[counter++] = -1.0;   // area constraint (unused)
 288         }
 289     }
 290
 291     // prep the output structures
 292     out.pointlist = (REAL *) NULL;        // Not needed if -N switch used.
 293     // Not needed if -N switch used or number of point attributes is zero:
 294     out.pointattributelist = (REAL *) NULL;
 295     out.pointmarkerlist = (int *) NULL;   // Not needed if -N or -B switch used.
 296     out.trianglelist = (int *) NULL;      // Not needed if -E switch used.
 297     // Not needed if -E switch used or number of triangle attributes is zero:
 298     out.triangleattributelist = (REAL *) NULL;
 299     out.neighborlist = (int *) NULL;      // Needed only if -n switch used.
 300     // Needed only if segments are output (-p or -c) and -P not used:
 301     out.segmentlist = (int *) NULL;
 302     // Needed only if segments are output (-p or -c) and -P and -B not used:
 303     out.segmentmarkerlist = (int *) NULL;
 304     out.edgelist = (int *) NULL;          // Needed only if -e switch used.
 305     out.edgemarkerlist = (int *) NULL;    // Needed if -e used and -B not used.
 306
 307     vorout.pointlist = (REAL *) NULL;     // Needed only if -v switch used.
 308     // Needed only if -v switch used and number of attributes is not zero:
 309     vorout.pointattributelist = (REAL *) NULL;
 310     vorout.edgelist = (int *) NULL;       // Needed only if -v switch used.
 311     vorout.normlist = (REAL *) NULL;      // Needed only if -v switch used.
 312
 313     // TEMPORARY
 314     // write_out_data(&in);
 315
 316     // Triangulate the points.  Switches are chosen to read and write
 317     // a PSLG (p), preserve the convex hull (c), number everything
 318     // from zero (z), assign a regional attribute to each element (A),
 319     // and produce an edge list (e), and a triangle neighbor list (n).
 320
 321     triangulate("pczAen", &in, &out, &vorout);
 322
 323     // TEMPORARY
 324     write_out_data(&out);
 325
 326     // free mem allocated to the "Triangle" structures
 327     free(in.pointlist);
 328     free(in.pointattributelist);
 329     free(in.pointmarkerlist);
 330     free(in.regionlist);
 331     free(out.pointlist);
 332     free(out.pointattributelist);
 333     free(out.pointmarkerlist);
 334     free(out.trianglelist);
 335     free(out.triangleattributelist);
 336     // free(out.trianglearealist);
 337     free(out.neighborlist);
 338     free(out.segmentlist);
 339     free(out.segmentmarkerlist);
 340     free(out.edgelist);
 341     free(out.edgemarkerlist);
 342     free(vorout.pointlist);
 343     free(vorout.pointattributelist);
 344     free(vorout.edgelist);
 345     free(vorout.normlist);
 346
 347     return 0;
 348 }
 349
 350
 351 // $Log$
 352 // Revision 1.8  1999/03/21 14:02:06  curt
 353 // Added a mechanism to dump out the triangle structures for viewing.
 354 // Fixed a couple bugs in first pass at triangulation.
 355 // - needed to explicitely initialize the polygon accumulator in triangle.cxx
 356 //   before each polygon rather than depending on the default behavior.
 357 // - Fixed a problem with region attribute propagation where I wasn't generating
 358 //   the hole points correctly.
 359 //
 360 // Revision 1.7  1999/03/20 20:32:55  curt
 361 // First mostly successful tile triangulation works.  There's plenty of tweaking
 362 // to do, but we are marching in the right direction.
 363 //
 364 // Revision 1.6  1999/03/20 13:22:11  curt
 365 // Added trisegs.[ch]xx tripoly.[ch]xx.
 366 //
 367 // Revision 1.5  1999/03/20 02:21:52  curt
 368 // Continue shaping the code towards triangulation bliss.  Added code to
 369 // calculate some point guaranteed to be inside a polygon.
 370 //
 371 // Revision 1.4  1999/03/19 22:29:04  curt
 372 // Working on preparationsn for triangulation.
 373 //
 374 // Revision 1.3  1999/03/19 00:27:10  curt
 375 // Continued work on triangulation preparation.
 376 //
 377 // Revision 1.2  1999/03/18 04:31:11  curt
 378 // Let's not pass copies of huge structures on the stack ... ye might see a
 379 // segfault ... :-)
 380 //
 381 // Revision 1.1  1999/03/17 23:51:59  curt
 382 // Initial revision.
 383 //