Logo Search packages:      
Sourcecode: magics++ version File versions  Download package

Polyline.cc

Go to the documentation of this file.
/******************************** LICENSE ********************************

 Copyright 2007 European Centre for Medium-Range Weather Forecasts (ECMWF)

 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at 

    http://www.apache.org/licenses/LICENSE-2.0

 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.

 ******************************** LICENSE ********************************/

/*! \file Polyline.cc
    \brief Implementation of polyline graphics class (template).
    
    Magics Team - ECMWF 2004
    
    Started: Jan 2004
    
    Changes:
    
*/
#include "Polyline.h"
#include "Transformation.h"
#include "TeCoord2D.h"
#include "TeGeometryAlgorithms.h"


using namespace magics;


//template <class P>
//inline double squaredistance(const P& p1, const P& p2) 
//{
//    return (p1.y()-p2.y())*(p1.y()-p2.y())+(p1.x()-p2.x())*(p1.x()-p2.x());
//}
//
//template <class P>
//inline double side(const P& p1, const P& p2, const P& p3, bool print = false) 
//{
//    //if (print) //Log::dev()<< "isLeft[" << p1.x() << ", " << p2.x() << ", " << p3.x() << "]"; 
//    //double result = -((p2.y() -p1.y()) * (p1.x()-p3.x()) + (p2.x() - p1.x()) * (p3.y() -p1.y()));
//    //if (print) //Log::dev()<< " = " << result << "\n";
//    return -((p2.y() -p1.y()) * (p1.x()-p3.x()) + (p2.x() - p1.x()) * (p3.y() -p1.y()));
//    
//}
//   
//
//template <class P>
//class NotAcceptable 
//{
//public:
//    NotAcceptable(const P& point, double tolerance) : 
//        point_(point), tolerance_(tolerance) {}
//    ~NotAcceptable() {}
//    bool operator()(pair<P, P> & ab) { 
//        //Log::dev()<< "check (" << point_.x() << ", " << point_.y() << ") against ";        
//        //Log::dev()<< "[(" <<  ab.first.x() << ", " << ab.first.y() << ") ---> ";
//        //Log::dev()<< "(" <<  ab.second.x() << ", " << ab.second.y() << ")]";
//        double d = side(point_, ab.first, ab.second);
//        d = d * d;
//        double e = tolerance_ * squaredistance(ab.first, ab.second);
//        //Log::dev()<< " = " << d << " > " << e << " ??? " << "\n";
//        return (d > e); 
//        }
//protected:
//    const P& point_;
//    double tolerance_;
//};
//
//template <class P>
//class Acceptable 
//{
//public:
//    Acceptable(const P& point, double tolerance) : 
//        point_(point), tolerance_(tolerance) {}
//    ~Acceptable() {}
//    bool operator()(pair<P, P> & ab) { 
//        //Log::dev()<< "check (" << point_.x() << ", " << point_.y() << ") against ";        
//        //Log::dev()<< "[(" <<  ab.first.x() << ", " << ab.first.y() << ") ---> ";
//        //Log::dev()<< "(" <<  ab.second.x() << ", " << ab.second.y() << ")]" << "\n";
//        double d = side(point_, ab.first, ab.second);
//        d = d * d;
//        double e = tolerance_ * squaredistance(ab.first, ab.second);
//        return (d <= e); 
//        }
//protected:
//    const P& point_;
//    double tolerance_;
//};
//
//template <class P>
//class PointReduction 
//{
//    
//public: 
//    PointReduction(vector<PaperPoint>& line, vector<PaperPoint>& result, double tolerance = 1) : 
//        
//        tolerance_(tolerance*tolerance), line_(line), result_(result) {}
//   
//   
//    void reduce() {
//        if (line_.size() < 3) return;
//        current_ = line_.begin();
//        result_.push_back(*current_);
//        working_.push_front(*current_); current_++;
//        working_.push_front(*current_); current_++;
//        working_.push_front(*current_); 
//        
//        while (current_ != line_.end() ) {
//          
//            process();
//            
//        }
//    }
//    
//    void process() {
//        check();    
//        print("state of the hull", working_);
//        deque<PaperPoint>::iterator  point = working_.begin();
//        
//        sides_.push_back(std::make_pair(*(point+1), *(point+2)));
//        sides_.remove_if(NotAcceptable<PaperPoint>(*current_, tolerance_));
//        //Log::dev()<< "reste----> " << sides_.size() << "\n";    
//
//        // Check other side! 
//        {
//            for (int i = 0; i < working_.size() -1 ; i++) { 
//                rotate(point, point+1, working_.end());
//                if ( accept() ) {
//                    sides_.push_back(std::make_pair(*(point+1), *(point+2)));
//                    //Log::dev()<< "Accept--->side" << "\n";
//    
//                }
//            
//            }
//         
//                rotate(point, point+1, working_.end());
//            
//            if (sides_.empty()) {
//                // No acceptable sides ...
//                //print("state of the hull", working_);
//                //Log::dev()<< "" << "\n";
//                //Log::dev()<< "Store point ---> (" <<  (*(current_ -1)).x() << ", " << (*(current_ -1)).y() << ")" << "\n";
//                //Log::dev()<< "" << "\n";
//                result_.push_back(*(current_ -1)); 
//                working_.clear();
//                working_.push_front(result_.back()); 
//                working_.push_front(*current_); 
//               
//                if (++current_ == line_.end()) return;
//            
//                //Log::dev()<< "Add point ---> (" <<  (*current_).x() << ", " << (*current_).y() << ")" << "\n";
//                working_.push_front(*current_); 
//                return;
//            }
//            
//            // add a point to the line!!
//           int pos = -1;
//           while (pos == -1)  {
//         
//           if (++current_ == line_.end()) {
//                result_.push_back(*(--current_));
//                ++current_;
//                return;
//            }
//            
//            //Log::dev()<< "Add point to the hhull ---> (" <<  (*current_).x() << ", " << (*current_).y() << ")" << "\n";
//            //print("state of the hull", working_);
//            
//            int nb = working_.size();
//            
//           
//            for (int i = 0; i < nb - 1; i++) {
//                  
//                  if ( side(*current_, working_.back(), working_.front(), true) >= 0 ) {
//                       
//                            pos = nb-i;
//                        
//                                 
//                  }
//                  rotate(working_.begin(), working_.begin()+1, working_.end());
//            }
//            
//            rotate(working_.begin(), working_.begin()+1, working_.end());
//            
//            //print("new state of the hull", working_);
//            //Log::dev()<< "attach to side --->" << pos << "\n";
//            if (pos == -1) continue;
//            if ( pos == nb &&  side(*current_, working_.back(), working_.front(), true) >=0 ) {                
//                working_.pop_back();
//                working_.push_back(*current_);              
//            }
//            else {
//                
//                print("where are we now? new state of the hull", working_);
//                working_.erase(working_.begin()+(nb-pos-1));
//                working_.insert(working_.begin()+(nb-pos-1), *current_); 
//            }
//            while (working_.front() != *current_) rotate(working_.begin(), working_.begin()+1, working_.end());
//            //print("where qre we now? new state of the hull", working_);
//            
//            
//        }
//           
//           sides_.pop_front();  
////           for (list<pair<P, P> >::const_iterator ab = sides_.begin(); ab != sides_.end(); ++ab)
////           {
////               //Log::dev()<< "Do we keep this one? ---> [" << ab->first.x() << ", " << ab->second.x() << "]" << "\n";
////               
////           }
//           
//          
//            
//            
//            
//            sides_.remove_if(NotAcceptable<PaperPoint>(*current_, tolerance_));
//            //Log::dev()<< "nb of sides ---> " << sides_.size() << "\n";
//            
//                 
//         
//            
//            
//              
//            
//            
//        }
//            
//                 
//        
//    }
//protected:
//   double             tolerance_;
//   vector<PaperPoint>&           line_;
//   vector<PaperPoint>&           result_;
//   list<pair<P, P> >  sides_;
//   deque<PaperPoint>            working_;
//   vector<PaperPoint>::iterator  current_;
// 
//   inline bool accept()   
//   {    
//       deque<PaperPoint>::iterator  point = working_.begin(); 
//       //Log::dev()<< "check (" << (*point).x() << ", " << (*point).y() << ") against ";        
//        //Log::dev()<< "[(" <<  (*(point+1)).x() << ", " << (*(point+1)).y() << ") ---> ";
//        //Log::dev()<< "(" <<  (*(point+2)).x() << ", " << (*(point+2)).y() << ")]" << "\n";
//        double d = side((*point), *(point+1), *(point+2));
//        d = d * d;
//        double e = tolerance_ * squaredistance(*(point+1),*(point+2));
//        return (d < e); 
//   }
//   void check() {
//       deque<PaperPoint>::iterator  point = working_.begin();
//       //print("check before " , working_);
//       if ( side(*point, *(point+1), *(point+2) ) > 0 ) 
//            std::iter_swap((point+1), (point+2));
//        //print("check after " , working_);
//   }
//   
//   void print(string info, deque<PaperPoint>& v) 
//   {
//       //Log::dev()<< info << "--->[";
//       string sep = "";
//       for ( deque<PaperPoint>::const_iterator p = v.begin(); p != v.end(); ++p) {
//            //Log::dev()<< sep << "(" << (*p).x() << ", " << (*p).y() << ")";
//            sep = ", ";
//       }
//       //Log::dev()<< "]" << "\n";
//   }
//   
//  
//};


class SimplePointReduction 
{
public:
    SimplePointReduction(vector<PaperPoint*>& in, vector<PaperPoint>& out, double threshold = 1) : 
        threshold_(threshold), in_(in), out_(out) {}
    ~SimplePointReduction() {}
    
    void reduce() {
      
            
            /* We always store the first point as being significant */

           
            out_.push_back(*(in_.front()));
            vector<PaperPoint*>::const_iterator last = in_.begin();
            vector<PaperPoint*>::const_iterator current = in_.begin();

            
            double gradient;
            double projectedValue;
           
            while ((++current) != in_.end() )
            {
               
                #define PT_X(I) (*I)->x_
                #define PT_Y(I) (*I)->y_
                #define S (out_.end() - 1)           /* last-stored sig point */
                
                double x = (*current)->x_;
                double px = (out_.end() - 1)->x_;
                double y = (*current)->y_;
                double py = (out_.end() - 1)->y_;

                //Log::dev()<< "deal with ---> [" <<  PT_X(current) << ", " <<  PT_Y(current) << "]" << "\n";
                //Log::dev()<< "last stored is ---> [" <<  PT_X(S) << ", " <<  PT_Y(S) << "]" << "\n";
               
                /* Get the gradient from our last significant point to the
                   next point to consider */
                //Log::dev()<< PT_X(current) - PT_X(S) << "\n";
                
                if ( x - px != 0.0)
                {
                    gradient = (y - py) / (x - px);
                    //Log::dev()<< "gradient=" << gradient << "\n";
                }
                else
                {
                    gradient = 1.0;
                }


                /* Now, for each point (i) in-between, project a line from the
                   last significant point to the x-value of the i'th point.
                   If the y-distance is greater than our threshold, then we must
                   store the previous point considered because this one is
                   not good enough. */

                //for (i = lastIndex + 1; i < nextPoint; i++)
                for (vector<PaperPoint*>::const_iterator i = last+1; i != current; ++i)
                {
                    //Log::dev()<< "point in between  ---> [" <<  PT_X(i) << ", " <<  PT_Y(i) << "]" << "\n";
                    projectedValue = y  + (gradient * ((*i)->x_ - px));
                    
                    //Log::dev()<< "---->" << fabs (projectedValue - PT_Y(i)) << "\n";

                    /* If this is too far off, store the previous point */

                    if (fabs (projectedValue - (*i)->y_) > threshold_)
                    {
                        //pnIndices [numSigPts] = nextPoint - 1;
                        out_.push_back(*(*(current-1)));
                        last = current-1;
                        break;
                    }
                }
            }
             out_.push_back(*(in_.back()));
    }
protected:
    double             threshold_;
    vector<PaperPoint*>&           in_;
    vector<PaperPoint>&           out_;
};



/********************************************************************************
    DistanceThresholdPointReduction

    A class containing a function to reduce the number of points in a polyline.
    The basic idea is to attempt to remove each point from the polyline in turn;
    the first and last points in the input set are always retained.
    When a point p is 'deleted' from between points a and b, each point between
    a and b is considered: if the distance from point i to line ab is greater
    than the specified threshold, then we cannot delete the point p and so it is
    added to our output set. If in doubt, trust the code more than this
    explanation.
    Returns 1 if it's ok, but 0 if it got into a state where it could not finish
    the polygon because it was trying to avoid creating an intersection.

*********************************************************************************/


class DistanceThresholdPointReduction 
{
public:
      DistanceThresholdPointReduction(deque<PaperPoint>& in, deque<PaperPoint>& out, double threshold = 1) : 
        threshold_(threshold), in_(in), out_(out) {}
      ~DistanceThresholdPointReduction() {}
    
      bool reduce() {


//        printf ("\nPolygon with %d vertices: ", in_.size());

        // We always store the first point as being significant
        out_.push_back((in_.front()));
        deque<PaperPoint>::const_iterator last = in_.begin();
        deque<PaperPoint>::const_iterator current = in_.begin();

        double magnitude;
        double distance_squared;
        double threshold_squared = threshold_ * threshold_;
        int    previous_intersected = 0;
        int    npoints_processed = 0;

        while ((++current) != in_.end() )
        {
            const int max_intersection_tries = 5;
            double ax = (*current).x_;
            double bx = (out_.end() - 1)->x_;
            double ay = (*current).y_;
            double by = (out_.end() - 1)->y_;
            double u;  // distance along line ab that point i is (see below)
            double ix;
            double iy;
            
            // 'pacifier' so the user knows something is happening
            if (npoints_processed++ % 1000 == 0)
            {
                Log::dev()<< ".";
            }


            // now check whether the new line segment A1-A2 (from previous-stored point to this one)
            // intersects any of the existing line segments B1-B2 in the polygon.
            // if so, then that's bad, so we use the previous point .... unless that would also 
            // cause an interesection (previous_intersected>0)... if this continues for a
            // certain number of points, then we just do it anyway and don't worry about the intersection -
            // if we don't have this 'get-out' clause, then we can end up with a very incomplete coastline.

            if (out_.size() > 3)
            {
                bool linesIntersect = false;
                TeCoord2D A1(ax, ay);
                TeCoord2D A2(bx, by);

                // for each line segment already stored...

                for (deque<PaperPoint>::const_iterator j = out_.begin() + 2; (j < out_.end() - 1) && !linesIntersect; ++j)
                {
                    TeCoord2D B1(j->x_, j->y_);
                    TeCoord2D B2((j-1)->x_, (j-1)->y_);
                    TeCoord2D Intersect;

                    linesIntersect = linesIntersect || TeSegmentsIntersectPoint(A1, A2, B1, B2, Intersect);
                }

                if (linesIntersect)
                {
                    if (previous_intersected)    // 'stuck' behind an existing line?
                    {
                       previous_intersected--;   // keep going, but not forever...
                       continue;
                    }
                    
                    else
                    {
                        out_.push_back((*(current-1)));  // this line intersects, but the previous one wouldn't
                        last = current-1;
                        previous_intersected = max_intersection_tries;
                        continue;
                    }
                }
                else
                {
                    bool skip_point = (previous_intersected > 0);
                    previous_intersected = 0;                      // we might be able to use this point next time
                    if (skip_point) continue;
                }
            }


            /* Now, for each point (i) in-between the last-stored point (b)
               and the one we are considering to reject (a), project a line from
               a to b and find the distance from i to this line.
               If the distance is greater than our threshold, then we must
               store the previous point considered (i). */

            for (deque<PaperPoint>::const_iterator i = last+1; i != current; ++i)
            {         
                // magnitude of line ab

                magnitude = sqrt (((bx - ax) * (bx - ax)) + ((by - ay) * (by - ay)));


                // how far along ab is the closest point to i?
                
                u = ( ( ( (*i).x_ - ax ) * ( bx - ax ) ) + ( ( (*i).y_ - ay ) * (by - ay ) ) ) /
                    ( magnitude * magnitude );


                // find the point on the line that is closest to point i

                ix = ax + u * ( bx - ax );
                iy = ay + u * ( by - ay );


                // now find the distance from i to this point

                distance_squared = ((ix - (*i).x_) * (ix - (*i).x_)) + ((iy - (*i).y_) * (iy - (*i).y_));


                // if this is too far off, store the previous point (don't bother with
                // square roots as we are just comparing distances)

                if (distance_squared > threshold_squared)
                {
                    out_.push_back((*(current-1)));
                    last = current-1;
                    break;
                }
            }
        }
        
        // always store the last point

        out_.push_back((in_.back()));
        
        return !previous_intersected;  // 1 if ok, 0 if problem
    }

protected:
    double        threshold_;
    deque<PaperPoint>&       in_;
    deque<PaperPoint>&       out_;
};


Polyline::Polyline() : 
    maxY_(int_MIN), minY_(int_MAX),
    maxX_(int_MIN), minX_(int_MAX), 
    winding_(int_MAX), area_(double_MIN)
{

      bounding_box_is_current_ = true;
}


Polyline::~Polyline() 
{
}





bool Polyline::reproject(BasicGraphicsObjectContainer& out) const
{
      const Transformation& transformation = out.transformation();
      transformation(*this, out);   
      return false;     
}




Polyline* Polyline::getNew() const
{
      Polyline* poly = new Polyline();

      poly->copy(*this);

      return poly;
}

void Polyline::print(ostream& out) const 
{
      out << "Polyline[";
      //out << " Colour: "<< getColour();
      //out << ", box = " << box_;
      //out << ", label = " << label_;
      out << ", nb_points = " << this->size();
        if ( this->size() < 10 ) {
            out << "Vector[";
            string sep = "";
          const unsigned int si = this->size();
            for (unsigned int i = 0; i < si; i++) {
                out << sep << (*this)[i];
                sep = ", ";
            }
            out << "]";
        }
        else {
            unsigned int nb = this->size();
            out << "Vector[" << (*this)[0] << ", " << (*this)[1] << ", " << (*this)[2];
            out << "...." << (*this)[nb-3] << ", " << (*this)[nb-2]  << ", " << (*this)[nb-1];
            out << "(" << nb << " elements)]";
        }
      out << "]";
}



Polyline* Polyline::simplify(double tol)
{
    Polyline* poly;
    double tolerance = tol;

      poly = getNew();
      DistanceThresholdPointReduction method(*this, *poly, tolerance);
    method.reduce();
    return poly;

/*
    bool ok = false;

    while (!ok)
    {
         poly = getNew();

    //      Log::dev()<< "Simplify " << poly->size() << "---->" << poly->capacity()  << "\n";
    //      for (vector<PaperPoint*>::iterator point = working_.begin(); point != working_.end(); ++point) 
    //            poly->push_back(*(*point));
    //      SimplePointReduction<PaperPoint> method(working_, *poly, tol);


          DistanceThresholdPointReduction method(*this, *poly, tolerance);

          if (method.reduce() || tolerance < 0.0004)
        {
            ok = true;
            return poly;
        }
        
        else
        {
            printf ("Polyline::simplify: going from tolerance %f to %f\n", tolerance, tolerance * 0.8);
            delete poly;
            tolerance *= 0.8;
        }
    }
*/


      //Log::dev()<< "Simplify " << poly->size() << "---->" << poly->capacity()  << "\n";
      
      //Log::dev()<< "Simplify " << working_.size() << "---->" << poly->size() <<  "---->" << poly->capacity() << "\n";
}


struct WindingFunction
{
    WindingFunction(const Polyline& ref) : reference_(ref), factor_(0) {}
    ~WindingFunction() {}
    int operator()(const PaperPoint& p1, const PaperPoint& p2) 
    {
        if ( p1.y() <= reference_.y() ) {         // start y <= P.y
            if ( p2.y() >  reference_.y() )      // an upward crossing
                if ( reference_.isLeft( p1, p2) > 0)  // P left of edge
                    factor_++;            // have a valid up intersect
        }
        else {                       // start y > P.y (no test needed)
            if ( p2.y() <= reference_.y() )     // a downward crossing
                if ( reference_.isLeft( p1, p2) < 0)  // P right of edge
                    factor_--;            // have a valid down intersect
        }
        return factor_;   
    }
    const Polyline& reference_;
    int factor_;  
private:
    WindingFunction(const WindingFunction&);
};


int Polyline::windingFactor() const
{
      if (winding_ != int_MAX ) return winding_;
      if (this->empty()) return winding_;
      // Work on the  vector of points.. not on the working array of PaperPoint* 
      WindingFunction tool(*this);
      vector<PaperPoint> result;

      const unsigned int isize = this->size() - 1;
      for (unsigned int i = 0; i < isize; i++)
      {
            tool((*this)[i], (*this)[i+1]);
      }

      winding_ = tool.factor_;
      return winding_;
}


void Polyline::makeCounterClockwise()
{
      area();
      if ( clockwise() ) 
            std::reverse(this->begin(), this->end()); 
}

void Polyline::makeClockwise()
{
      area();
      if ( !clockwise() ) 
            std::reverse(this->begin(), this->end()); 
}
/*
bool Polyline::inside(const Polyline& other) const
{
      unsigned int i, j=0;
      bool  odd = false;
      PaperPoint point = (*this)[this->size()/2];

      for (i=0; i < other.size(); i++)
      {
            j++; 
            if ( j == other.size() ) j=0;
            if ( other[i].y() < point.y() && other[j].y() >= point.y()
              || other[j].y() < point.y() && other[i].y() >= point.y() )
            {
                  if (other[i].x()+(point.y()-other[i].y())/(this->poly[j].y() - other[i].y() )*(other[j].x()-other[i].x())<point.x())
                  {
                        odd =! odd;
                  }
            }
      }
      return odd; 
}
*/


double Polyline::area()
{
      if ( this->area_ != double_MIN ) return this->area_;
      this->area_ = 0 ;       
      if (this->size()) this->push_back((*this)[0]);
      for(unsigned int j = 0; j < this->size(); j++)
      {
            if(j) this->area_ += ((*this)[j-1].x() * (*this)[j].y() - (*this)[j].x() * (*this)[j-1].y() );
      }
      if(this->size()) this->pop_back();
      return this->area_;
}



void Polyline::redisplay(const BaseDriver& driver) const 
{
    driver.redisplay(*this);
}



Generated by  Doxygen 1.6.0   Back to index