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     int index;
  44
  45     // traverse the dem fit list and gpc_polys building a unified node
  46     // list and converting the polygons so that they reference the
  47     // node list by index (starting at zero) rather than listing the
  48     // points explicitely
  49
  50     const_fitnode_list_iterator f_current, f_last;
  51     f_current = fit_list.begin();
  52     f_last = fit_list.end();
  53     for ( ; f_current != f_last; ++f_current ) {
  54         index = trinodes.unique_add( *f_current );
  55     }
  56
  57     gpc_polygon *gpc_poly;
  58     const_gpcpoly_iterator current, last;
  59
  60     // process polygons in priority order
  61     cout << "prepairing node list and polygons" << endl;
  62
  63     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
  64         cout << "area type = " << i << endl;
  65         current = gpc_polys.polys[i].begin();
  66         last = gpc_polys.polys[i].end();
  67         for ( ; current != last; ++current ) {
  68             gpc_poly = *current;
  69             cout << "processing a polygon, contours = "
  70                  << gpc_poly->num_contours << endl;
  71
  72             FGTriPoly poly;
  73
  74             if (gpc_poly->num_contours <= 0 ) {
  75                 cout << "FATAL ERROR! no contours in this polygon" << endl;
  76                 exit(-1);
  77             }
  78
  79             if (gpc_poly->num_contours > 1 ) {
  80                 cout << "FATAL ERROR! no multi-contour support" << endl;
  81                 sleep(2);
  82                 // exit(-1);
  83             }
  84
  85             for ( int j = 0; j < gpc_poly->num_contours; j++ ) {
  86                 for ( int k = 0; k < gpc_poly->contour[j].num_vertices; k++ ) {
  87                     Point3D p( gpc_poly->contour[j].vertex[k].x,
  88                                gpc_poly->contour[j].vertex[k].y,
  89                                0 );
  90                     index = trinodes.unique_add( p );
  91                     poly.add_node(index);
  92                     // cout << index << endl;
  93                 }
  94                 poly.calc_point_inside( trinodes );
  95
  96                 polylist[i].push_back(poly);
  97             }
  98         }
  99     }
 100
 101     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 102         if ( polylist[i].size() ) {
 103             cout << get_area_name((AreaType)i) << " = "
 104                  << polylist[i].size() << endl;
 105         }
 106     }
 107
 108     // traverse the polygon lists and build the segment (edge) list
 109     // that is used by the "Triangle" lib.
 110
 111     FGTriPoly poly;
 112     int i1, i2;
 113     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 114         cout << "area type = " << i << endl;
 115         tripoly_list_iterator tp_current, tp_last;
 116         tp_current = polylist[i].begin();
 117         tp_last = polylist[i].end();
 118
 119         // process each polygon in list
 120         for ( ; tp_current != tp_last; ++tp_current ) {
 121             poly = *tp_current;
 122
 123             for ( int j = 0; j < (int)(poly.size()) - 1; ++j ) {
 124                 i1 = poly.get_pt_index( j );
 125                 i2 = poly.get_pt_index( j + 1 );
 126                 trisegs.unique_add( FGTriSeg(i1, i2) );
 127             }
 128             i1 = poly.get_pt_index( 0 );
 129             i2 = poly.get_pt_index( poly.size() - 1 );
 130             trisegs.unique_add( FGTriSeg(i1, i2) );
 131         }
 132     }
 133
 134     return 0;
 135 }
 136
 137
 138 // triangulate each of the polygon areas
 139 int FGTriangle::run_triangulate() {
 140     FGTriPoly poly;
 141     Point3D p;
 142     struct triangulateio in, out, vorout;
 143     int counter;
 144
 145     // point list
 146     trinode_list node_list = trinodes.get_node_list();
 147     in.numberofpoints = node_list.size();
 148     in.pointlist = (REAL *) malloc(in.numberofpoints * 2 * sizeof(REAL));
 149
 150     trinode_list_iterator tn_current, tn_last;
 151     tn_current = node_list.begin();
 152     tn_last = node_list.end();
 153     counter = 0;
 154     for ( ; tn_current != tn_last; ++tn_current ) {
 155         in.pointlist[counter++] = tn_current->x();
 156         in.pointlist[counter++] = tn_current->y();
 157     }
 158
 159     in.numberofpointattributes = 1;
 160     in.pointattributelist = (REAL *) malloc(in.numberofpoints *
 161                                             in.numberofpointattributes *
 162                                             sizeof(REAL));
 163     for ( int i = 0; i < in.numberofpoints * in.numberofpointattributes; i++) {
 164         in.pointattributelist[i] = 0.0;
 165     }
 166
 167     in.pointmarkerlist = (int *) malloc(in.numberofpoints * sizeof(int));
 168     for ( int i = 0; i < in.numberofpoints; i++) {
 169         in.pointmarkerlist[i] = 0;
 170     }
 171
 172     // segment list
 173     triseg_list seg_list = trisegs.get_seg_list();
 174     in.numberofsegments = seg_list.size();
 175     in.segmentlist = (int *) malloc(in.numberofsegments * 2 * sizeof(int));
 176
 177     triseg_list_iterator s_current, s_last;
 178     s_current = seg_list.begin();
 179     s_last = seg_list.end();
 180     counter = 0;
 181     for ( ; s_current != s_last; ++s_current ) {
 182         in.segmentlist[counter++] = s_current->get_n1();
 183         in.segmentlist[counter++] = s_current->get_n2();
 184     }
 185
 186     // hole list (make holes for airport ignore areas)
 187     in.numberofholes = polylist[(int)AirportIgnoreArea].size();
 188     in.holelist = (REAL *) malloc(in.numberofholes * 2 * sizeof(REAL));
 189
 190     tripoly_list_iterator h_current, h_last;
 191     h_current = polylist[(int)AirportIgnoreArea].begin();
 192     h_last = polylist[(int)AirportIgnoreArea].end();
 193     counter = 0;
 194     for ( ; h_current != h_last; ++h_current ) {
 195         poly = *h_current;
 196         p = poly.get_point_inside();
 197         in.holelist[counter++] = p.x();
 198         in.holelist[counter++] = p.y();
 199     }
 200
 201     // region list
 202     in.numberofregions = 0;
 203     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 204         in.numberofregions += polylist[i].size();
 205     }
 206
 207     in.regionlist = (REAL *) malloc(in.numberofregions * 4 * sizeof(REAL));
 208     for ( int i = 0; i < FG_MAX_AREA_TYPES; ++i ) {
 209         tripoly_list_iterator h_current, h_last;
 210         h_current = polylist[(int)AirportIgnoreArea].begin();
 211         h_last = polylist[(int)AirportIgnoreArea].end();
 212         counter = 0;
 213         for ( ; h_current != h_last; ++h_current ) {
 214             poly = *h_current;
 215             p = poly.get_point_inside();
 216             in.regionlist[counter++] = p.x();  // x coord
 217             in.regionlist[counter++] = p.y();  // y coord
 218             in.regionlist[counter++] = i;      // region attribute
 219             in.regionlist[counter++] = -1.0;   // area constraint (unused)
 220         }
 221     }
 222
 223     // prep the output structures
 224     out.pointlist = (REAL *) NULL;        // Not needed if -N switch used.
 225     // Not needed if -N switch used or number of point attributes is zero:
 226     out.pointattributelist = (REAL *) NULL;
 227     out.pointmarkerlist = (int *) NULL;   // Not needed if -N or -B switch used.
 228     out.trianglelist = (int *) NULL;      // Not needed if -E switch used.
 229     // Not needed if -E switch used or number of triangle attributes is zero:
 230     out.triangleattributelist = (REAL *) NULL;
 231     out.neighborlist = (int *) NULL;      // Needed only if -n switch used.
 232     // Needed only if segments are output (-p or -c) and -P not used:
 233     out.segmentlist = (int *) NULL;
 234     // Needed only if segments are output (-p or -c) and -P and -B not used:
 235     out.segmentmarkerlist = (int *) NULL;
 236     out.edgelist = (int *) NULL;          // Needed only if -e switch used.
 237     out.edgemarkerlist = (int *) NULL;    // Needed if -e used and -B not used.
 238
 239     vorout.pointlist = (REAL *) NULL;     // Needed only if -v switch used.
 240     // Needed only if -v switch used and number of attributes is not zero:
 241     vorout.pointattributelist = (REAL *) NULL;
 242     vorout.edgelist = (int *) NULL;       // Needed only if -v switch used.
 243     vorout.normlist = (REAL *) NULL;      // Needed only if -v switch used.
 244
 245     // Triangulate the points.  Switches are chosen to read and write
 246     // a PSLG (p), preserve the convex hull (c), number everything
 247     // from zero (z), assign a regional attribute to each element (A),
 248     // and produce an edge list (e), and a triangle neighbor list (n).
 249
 250     triangulate("pczAen", &in, &out, &vorout);
 251
 252     // TEMPORARY
 253     //
 254
 255     // Write out the triangulated data to files so we can check
 256     // visually that things seem reasonable
 257
 258     FILE *node = fopen("tile.node", "w");
 259     fprintf(node, "%d 2 %d 0\n",
 260             out.numberofpoints, out.numberofpointattributes);
 261     for (int i = 0; i < out.numberofpoints; i++) {
 262         fprintf(node, "%d %.6f %.6f %.2f\n",
 263                 i, out.pointlist[2*i], out.pointlist[2*i + 1], 0.0);
 264     }
 265     fclose(node);
 266
 267     FILE *ele = fopen("tile.ele", "w");
 268     fprintf(ele, "%d 3 0\n", out.numberoftriangles);
 269     for (int i = 0; i < out.numberoftriangles; i++) {
 270         fprintf(ele, "%d ", i);
 271         for (int j = 0; j < out.numberofcorners; j++) {
 272             fprintf(ele, "%d ", out.trianglelist[i * out.numberofcorners + j]);
 273         }
 274         fprintf(ele, "\n");
 275     }
 276     fclose(ele);
 277
 278     // free mem allocated to the "Triangle" structures
 279     free(in.pointlist);
 280     free(in.pointattributelist);
 281     free(in.pointmarkerlist);
 282     free(in.regionlist);
 283     free(out.pointlist);
 284     free(out.pointattributelist);
 285     free(out.pointmarkerlist);
 286     free(out.trianglelist);
 287     free(out.triangleattributelist);
 288     // free(out.trianglearealist);
 289     free(out.neighborlist);
 290     free(out.segmentlist);
 291     free(out.segmentmarkerlist);
 292     free(out.edgelist);
 293     free(out.edgemarkerlist);
 294     free(vorout.pointlist);
 295     free(vorout.pointattributelist);
 296     free(vorout.edgelist);
 297     free(vorout.normlist);
 298
 299     return 0;
 300 }
 301
 302
 303 // $Log$
 304 // Revision 1.7  1999/03/20 20:32:55  curt
 305 // First mostly successful tile triangulation works.  There's plenty of tweaking
 306 // to do, but we are marching in the right direction.
 307 //
 308 // Revision 1.6  1999/03/20 13:22:11  curt
 309 // Added trisegs.[ch]xx tripoly.[ch]xx.
 310 //
 311 // Revision 1.5  1999/03/20 02:21:52  curt
 312 // Continue shaping the code towards triangulation bliss.  Added code to
 313 // calculate some point guaranteed to be inside a polygon.
 314 //
 315 // Revision 1.4  1999/03/19 22:29:04  curt
 316 // Working on preparationsn for triangulation.
 317 //
 318 // Revision 1.3  1999/03/19 00:27:10  curt
 319 // Continued work on triangulation preparation.
 320 //
 321 // Revision 1.2  1999/03/18 04:31:11  curt
 322 // Let's not pass copies of huge structures on the stack ... ye might see a
 323 // segfault ... :-)
 324 //
 325 // Revision 1.1  1999/03/17 23:51:59  curt
 326 // Initial revision.
 327 //