Triangulate/tripoly.cxx

   1 // tripoly.cxx -- "Triangle" polygon management 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 <Include/fg_constants.h>
  26 #include <Math/point3d.hxx>
  27
  28 #include "tripoly.hxx"
  29
  30
  31 // Constructor
  32 FGTriPoly::FGTriPoly( void ) {
  33 }
  34
  35
  36 // Destructor
  37 FGTriPoly::~FGTriPoly( void ) {
  38 }
  39
  40
  41 // Given a line segment specified by two endpoints p1 and p2, return
  42 // the slope of the line.
  43 static double slope( const Point3D& p0, const Point3D& p1 ) {
  44     if ( fabs(p0.x() - p1.x()) > FG_EPSILON ) {
  45         return (p0.y() - p1.y()) / (p0.x() - p1.x());
  46     } else {
  47         return 1.0e+999; // really big number
  48     }
  49 }
  50
  51
  52 // Given a line segment specified by two endpoints p1 and p2, return
  53 // the y value of a point on the line that intersects with the
  54 // verticle line through x.  Return true if an intersection is found,
  55 // false otherwise.
  56 static bool intersects( const Point3D& p0, const Point3D& p1, double x,
  57                          Point3D *result ) {
  58     // equation of a line through (x0,y0) and (x1,y1):
  59     //
  60     //     y = y1 + (x - x1) * (y0 - y1) / (x0 - x1)
  61
  62     double y;
  63
  64     if ( fabs(p0.x() - p1.x()) > FG_EPSILON ) {
  65         y = p1.y() + (x - p1.x()) * (p0.y() - p1.y()) / (p0.x() - p1.x());
  66     } else {
  67         return false;
  68     }
  69     result->setx(x);
  70     result->sety(y);
  71
  72     if ( p0.y() <= p1.y() ) {
  73         if ( (p0.y() <= y) && (y <= p1.y()) ) {
  74             return true;
  75         }
  76     } else {
  77         if ( (p0.y() >= y) && (y >= p1.y()) ) {
  78             return true;
  79         }
  80     }
  81
  82     return false;
  83 }
  84
  85
  86 // calculate an "arbitrary" point inside this polygon for assigning
  87 // attribute areas
  88 void FGTriPoly::calc_point_inside( const FGTriNodes& trinodes ) {
  89     Point3D tmp, min, ln, p1, p2, p3, m, result;
  90
  91     // 1. find point, min, with smallest y
  92
  93     // min.y() starts greater than the biggest possible lat (degrees)
  94     min.sety( 100.0 );
  95     // int min_index;
  96     int min_node_index = 0;
  97
  98     tripoly_iterator current, last;
  99     current = poly.begin();
 100     last = poly.end();
 101
 102     int counter = 0;
 103     for ( ; current != last; ++current ) {
 104         tmp = trinodes.get_node( *current );
 105         if ( tmp.y() < min.y() ) {
 106             min = tmp;
 107             // min_index = *current;
 108             min_node_index = counter;
 109
 110             // cout << "min index = " << *current
 111             //      << " value = " << min_y << endl;
 112         } else {
 113             // cout << "  index = " << *current << endl;
 114         }
 115         ++counter;
 116     }
 117     cout << "min node index = " << min_node_index << endl;
 118     cout << "min index = " << poly[min_node_index]
 119          << " value = " << trinodes.get_node( poly[min_node_index] )
 120          << " == " << min << endl;
 121
 122     // 2. take midpoint, m, of min with neighbor having lowest
 123     // fabs(slope)
 124
 125     if ( min_node_index == 0 ) {
 126         p1 = trinodes.get_node( poly[1] );
 127         p2 = trinodes.get_node( poly[poly.size() - 1] );
 128     } else if ( min_node_index == (int)(poly.size()) - 1 ) {
 129         p1 = trinodes.get_node( poly[0] );
 130         p2 = trinodes.get_node( poly[poly.size() - 1] );
 131     } else {
 132         p1 = trinodes.get_node( poly[min_node_index - 1] );
 133         p2 = trinodes.get_node( poly[min_node_index + 1] );
 134     }
 135     double s1 = fabs( slope(min, p1) );
 136     double s2 = fabs( slope(min, p2) );
 137     if ( s1 < s2  ) {
 138         ln = p1;
 139     } else {
 140         ln = p2;
 141     }
 142
 143     m.setx( (min.x() + ln.x()) / 2.0 );
 144     m.sety( (min.y() + ln.y()) / 2.0 );
 145     cout << "low mid point = " << m << endl;
 146
 147     // 3. intersect vertical line through m and all other segments.
 148     // save point, p3, with smallest y > m.y
 149
 150     p3.sety(100);
 151     for ( int i = 0; i < (int)(poly.size()) - 1; ++i ) {
 152         p1 = trinodes.get_node( poly[i] );
 153         p2 = trinodes.get_node( poly[i+1] );
 154
 155         if ( intersects(p1, p2, m.x(), &result) ) {
 156             // cout << "intersection = " << result << endl;
 157             if ( ( result.y() < p3.y() ) &&
 158                  ( fabs(result.y() - m.y()) > FG_EPSILON ) ) {
 159                 p3 = result;
 160             }
 161         }
 162     }
 163     p1 = trinodes.get_node( poly[0] );
 164     p2 = trinodes.get_node( poly[poly.size() - 1] );
 165     if ( intersects(p1, p2, m.x(), &result) ) {
 166         // cout << "intersection = " << result << endl;
 167         if ( ( result.y() < p3.y() ) &&
 168              ( fabs(result.y() - m.y()) > FG_EPSILON ) ) {
 169             p3 = result;
 170         }
 171     }
 172     cout << "low intersection of other segment = " << p3 << endl;
 173
 174     // 4. take midpoint of p2 && m as an arbitrary point inside polygon
 175
 176     inside.setx( (m.x() + p3.x()) / 2.0 );
 177     inside.sety( (m.y() + p3.y()) / 2.0 );
 178     cout << "inside point = " << inside << endl;
 179 }
 180
 181
 182 // $Log$
 183 // Revision 1.1  1999/03/20 13:21:36  curt
 184 // Initial revision.
 185 //