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

IsoPlot.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 IsoPlot.cc  
 \brief Implementation of the Template class IsoPlot.
 
 Magics Team - ECMWF 2004 
 
 Started: Wed 3-Mar-2004
  
 Changes:
 
 
 CONREC
 ======

 http://astronomy.swin.edu.au/~pbourke/projection/conrec/

 Copyright (c) 1996-1997 Nicholas Yue

 This software is copyrighted by Nicholas Yue. This code is base on the work of
 Paul D. Bourke CONREC.F routine

 The authors hereby grant permission to use, copy, and distribute this
 software and its documentation for any purpose, provided that existing
 copyright notices are retained in all copies and that this notice is included
 verbatim in any distributions. Additionally, the authors grant permission to
 modify this software and its documentation for any purpose, provided that
 such modifications are not distributed without the explicit consent of the
 authors and that existing copyright notices are retained in all copies. Some
 of the algorithms implemented by this software are patented, observe all
 applicable patent law.

 IN NO EVENT SHALL THE AUTHORS OR DISTRIBUTORS BE LIABLE TO ANY PARTY FOR
 DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL DAMAGES ARISING OUT
 OF THE USE OF THIS SOFTWARE, ITS DOCUMENTATION, OR ANY DERIVATIVES THEREOF,
 EVEN IF THE AUTHORS HAVE BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

 THE AUTHORS AND DISTRIBUTORS SPECIFICALLY DISCLAIM ANY WARRANTIES, INCLUDING,
 BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
 PARTICULAR PURPOSE, AND NON-INFRINGEMENT.  THIS SOFTWARE IS PROVIDED ON AN
 "AS IS" BASIS, AND THE AUTHORS AND DISTRIBUTORS HAVE NO OBLIGATION TO PROVIDE
 MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.
 
 */
#include "IsoPlot.h"
#include "Factory.h"
#include "MatrixHandler.h"
#include "Timer.h"
#include "GeoPoint.h"
#include "IsoHighlight.h"
#include "Colour.h"

using namespace magics;

template<class P> IsoPlot<P>::IsoPlot() {
      setTag("isoline");
}

template<class P> IsoPlot<P>::~IsoPlot() {
}

/*!
 Class information are given to the output-stream.
 */
template<class P> void IsoPlot<P>::print(ostream& out) const {
      out << "IsoPlot<P>[";
      IsoPlotAttributes<P>::print(out);
      out << "]";
}

#define distance_back(px, py) ( ((*line)->back().x()- px) * ((*line)->back().x()- px) ) + ( ((*line)->back().y()- py) * ((*line)->back().y()- py) )
#define distance_front(px, py) ( ((*line)->front().x()- px) * ((*line)->front().x()- px) ) + ( ((*line)->front().y()- py) * ((*line)->front().y()- py) )

inline void concatenate_back(vector<Polyline* >& lines,Polyline* poly) 
{
      double x = poly->back().x();
      double y = poly->back().y();

      for (vector<Polyline* >::iterator line = lines.begin(); line
                  != lines.end(); ++line) {
            if ( *line == poly || (*line)->empty() )
                  continue;
            
            if (same((*line)->front().x(), x) && same((*line)->front().y(), y) ) {
                  poly->insert(poly->end(), (*line)->begin()+1, (*line)->end());
                  poly->setBoundingBoxCurrent(false); // the insert() operation will invalidate the bounding box
                  Polyline* todelete = *line;
                  lines.erase(line);
                  delete todelete;
                  break;
            }
      }

}

class HelpSort {
public:
      HelpSort() {
      }
      ~HelpSort() {
      }
      bool operator()(Polyline* first, Polyline* second) {
            if (same(first->surface(), second->surface()) ) {
                  if (first->clockwise() )
                        return first->front().value() < second->front().value();
                  else
                        return first->front().value() > second->front().value();

            }
            return first->surface() > second->surface();
      }

};

inline void concatenate_front(vector<Polyline* >& lines,
            Polyline* poly) {
      double x = poly->front().x();
      double y = poly->front().y();

      for (vector<Polyline* >::iterator line = lines.begin(); line
                  != lines.end(); ++line) {
            if ( *line == poly || (*line)->empty() )
                  continue;
            
            if (same((*line)->back().x(), x) && same((*line)->back().y(), y) ) {
                  poly->insert(poly->begin(), (*line)->begin(), (*line)->end()-1);
                  poly->setBoundingBoxCurrent(false); // the insert() operation will invalidate the bounding box
                  Polyline* todelete = *line;
                  lines.erase(line);
                  delete todelete;
                  break;
            }
      }
      
}
#define equal(a,b) (same(a->first, b->first) && same(a->second, b->second) )
#define diff(a,b) !( same(a.first, b.first) && same(a.second, b.second) )
#define EQUAL(a,b)  (same(a.first, b.first) && same(a.second, b.second) )





namespace magics {

template <class P> 
class CellBox : public VectorOfPointers<vector<CellBox<P>* > >
{
public:
      CellBox(CellArray<P>* parent, int row1, int row2, int column1, int column2) :
            parent_(parent), row1_(row1), row2_(row2), column1_(column1),
                        column2_(column2) {
      }

      CellBox(CellArray<P>* parent) :
            parent_(parent), row1_(0), row2_(parent->rows_-1), column1_(0),
                        column2_(parent->columns_-1) {
      }
      CellBox() :
                  parent_(0), row1_(0), row2_(0), column1_(0),
                              column2_(0) {
            }
      virtual ~CellBox() { }
   
      virtual double value () {
            for (int row = row1_; row <= row2_; row++) {
                              for (int column = column1_; column <= column2_; column++) {
                                    if ( ! (*parent_)(row, column)->out() ) 
                                          return (*parent_)(row, column)->value(0);
                              }
                  }
        return 0;
      }
      virtual bool inrange () {
            for (int row = row1_; row <= row2_; row++) {
                                    for (int column = column1_; column <= column2_; column++) {
                                          if ( ! (*parent_)(row, column)->out() ) 
                                                return (*parent_)(row, column)->inrange();
                                    }
                        }
        return false;
      }
      virtual void levels(set<int>& used) {
            for (int row = row1_; row <= row2_; row++) {
                  for (int column = column1_; column <= column2_; column++) {
                        (*parent_)(row, column)->levels(used);
                  }
            }
      }
      void reduce(Shape& shape) {
            static int level = 0;
            level++;
            
            Shape::iterator first = shape.begin();
              Shape::iterator second = shape.begin();
               while(++second != shape.end()) {             
                        if ( *first == *second) {
                        shape.erase(first);                   
                        shape.erase(second);                  
                        reduce(shape);
                        
                        return;
                    }
                    ++first;      
                }
      }
      bool reshape(Shape& shape1, Shape& shape2) {
          if( !shape1.intersect(shape2) ) 
            return false;
          Shape::iterator where;
          Shape::iterator from;
          for ( Shape::iterator line1= shape1.begin(); line1 != shape1.end(); ++line1) {
             for ( Shape::iterator line2= shape2.begin(); line2 != shape2.end(); ++line2) {
                   
                   if ( *line2 == *line1 ) {
                         where = line1;
                         from = line2;
                         ++from;
                         std::rotate(shape2.begin(), from, shape2.end());
                               shape1.splice(where, shape2);
                               
                                  // Remove common lines!
                                reduce(shape1); 
                                if ( shape1.minx_ > shape2.minx_ ) shape1.minx_ = shape2.minx_;
                                if ( shape1.maxx_ < shape2.maxx_ ) shape1.maxx_ = shape2.maxx_;
                                if ( shape1.miny_ > shape2.miny_ ) shape1.miny_ = shape2.miny_;
                                if ( shape1.maxy_ < shape2.maxy_ ) shape1.maxy_ = shape2.maxy_;
                               
                                return true;
                   }
             }
            }
          
          return false;
        
      }
    
    void reduce(std::list<Shape>& list) 
    {
      static int level = 0;
      level++;
      
      for (std::list<Shape>::iterator entry1 = list.begin(); entry1 != list.end(); ++entry1 ) {
            for (std::list<Shape>::iterator entry2 = list.begin(); entry2 != list.end(); ++entry2) {
               if ( entry1 == entry2 ) continue; 
               if ( reshape(*entry1, *entry2) ) {
                    
                        list.erase(entry2);
                        reduce(list);
                        level--;
                        
                        return;
                  }
            }
      }
    }
     

      void reshape(Shape& cell) {
            Colour colour = cell.colour_;
            if ( colour.name() == "none" ) return;
            map<Colour, list<Shape> >::iterator entry = this->helper_.find(colour);
            if (entry == helper_.end() ) {
                  helper_.insert(make_pair(colour, list<Shape>()));
                  entry = helper_.find(colour);
                  entry->second.push_back(cell);
                  return;
            }

            
            
            list<Shape>::iterator last = entry->second.end();
            
            for (list<Shape>::iterator shape = entry->second.begin(); shape != last; ++shape) {
                  if (reshape(*shape, cell) ) {
                        //reduce(entry->second);
                        return;
                  }
            }

            entry->second.push_back(cell);
                  
      }
      
      //virtual void isoline(const vector<int>&, IsoPlot<P>& owner, CellBox<P>* parent = 0);

      virtual void reshape(CellBox<P>* parent) {
            if ( parent == this) return;
            for (map<Colour, list<Shape> >::iterator entry = helper_.begin(); entry != helper_.end(); ++entry) {                      
                              for (list<Shape>::iterator shape = entry->second.begin(); shape != entry->second.end(); ++shape) {                            
                                          parent->reshape(*shape);                              
                        }
            }
            helper_.clear();
      }
      
      
      virtual void split();
      

      virtual Shape& shape(const Colour& colour) {
            if (shape_.empty()) {
                  // bottom 
                  
                  
                  Cell<P>* cell;
        
                  for (int column = column1_; column <= column2_; column++) {
                        cell = (*parent_)(row1_, column);
                      shape_.push_back(cell->column(0), cell->row(0));
                  }

                  // right
                  for (int row = row1_; row <= row2_; row++) {
                        cell = (*parent_)(row, column2_);
                        shape_.push_back(cell->column(1), cell->row(1));
                  }

                  for (int column = column2_; column >= column1_; column--) {
                        cell = (*parent_)(row2_, column);
                        shape_.push_back(cell->column(2), cell->row(2));
                  }

                  for (int row = row2_; row >= row1_; row--) {
                        cell = (*parent_)(row, column1_);   
                        shape_.push_back(cell->column(3), cell->row(3));
                  }
            
            }
      
            shape_.colour_= colour;

            return shape_;

      }
    
      
      void shade(IsoPlot<P>& owner);
      
      void shade(IsoPlot<P>& owner, CellBox* parent) {
            set<int> used;
            this->levels(used);
            int levels = int(used.size());
            vector<int> contours;
            switch (levels) {
                  case 0:
                        break;
                  case 1:
                      if (this->inrange() ) {
                        
                        reshape(this->shape(owner.colour(this->value())));          
                        reshape(parent);
                      }
                        break;
                  
            
                        
                  
                  default:
                        
                        split();          
                        if (this->empty()) { 
                            assert( row1_ == row2_);
                            assert( column1_ == column2_);
                         
                           
                              owner.isoline(*(*parent_)(row1_, column1_), this);
                        }
                        else {
                              
                              for (typename CellBox<P>::iterator cell = this->begin(); cell != this->end(); ++cell) 
                                    (*cell)->shade(owner, this);
                        }
                        reshape(parent);
                        break;
            }
      }
      
      void contour(IsoPlot<P>& owner) {
               for (int row = row1_; row <= row2_; row++) {
                                    for (int column = column1_; column <= column2_; column++) {                                                                         
                                                      owner.isoline(*(*parent_)(row, column));
                                    }
                  }                 
      }
      
      

      void feed(IsoPlot<P>& owner, BasicGraphicsObjectContainer& out) {
            
            for ( map<Colour, list<Shape> >::iterator colour = helper_.begin(); colour!= helper_.end(); ++colour) {
                  reduce(colour->second);
                  
                  for (list<Shape>::const_iterator cells = colour->second.begin(); cells != colour->second.end(); ++cells) {
                            if ( cells->size() < 3 ) continue;
                           
                           
                            vector<PaperPoint> helper1;
                            stack<PaperPoint> helper;
                            double value = cells->value_;

                                               PaperPoint p1, p2, p3;
                                               Shape::const_iterator cell = cells->begin();
                                               helper1.push_back(PaperPoint(cell->p1x_, cell->p1y_, value));
                                               while ( cell != cells->end() ) {
                                                   p3 = PaperPoint(cell->p1x_, cell->p1y_, value);
                                                   if ( !(p3 == helper1.back() ) ) {
                                                         helper1.push_back(p3);
                                                   }
                                                   cell++;
                                               }

                                               vector<PaperPoint>::const_iterator point =  helper1.begin();
                                               for (point = helper1.begin();  point != helper1.end(); ++point) {
                                                   
                                                   
                                                   if ( helper.empty() ) {
                                                        helper.push(*point);
                                      
                                                        continue;
                                                  }
                                                   
                                                   p2 =  helper.top();
                                                   helper.pop();

                                                   if ( helper.empty() ) {    
                                                         helper.push(p2);
                                                         if ( !( p2 == *point) ) {                                                           
                                                               helper.push(*point);
                                                         }
                                                   }
                                                   else if ( !( *point == helper.top()) ) {
                                                      helper.push(p2);
                                                      helper.push(*point);
                                                   }
                                                 

                                               }

                            
                            if ( helper.size() < 3 ) continue;
                            
                            
                          
                            
                            
                        Polyline* poly = new Polyline();
                        poly->setFilled(true); 

                        poly->setThickness(1);
                        poly->setColour(cells->colour_);

                        
                        poly->setFillColour(cells->colour_);
                        
                        
                        
                        
                        
                        
                        while ( !helper.empty() ) {
                        
                              poly->push_back(helper.top());
                              helper.pop();
                        }
                  
                        
                        if ( !poly->empty() ) {
                              owner.getShading()(poly);
                                
                                out.push_back(poly);
                                
                        }
            }
            }
      }
      
      Shape shape_;
      CellArray<P>* parent_;
      int row1_;
      int row2_;
      int column1_;
      int column2_;
      map<Colour, list<Shape> > helper_;

};



}


class SortFunction {
public:
      SortFunction() {
            }
            ~SortFunction() {
            }
            bool operator()(const Shape& first, const Shape& second) {
                  
                  return first.surface() > second.surface();
            }

};
template <class P>
void CellBox<P>::shade(IsoPlot<P>& owner) {
          shade(owner, this);                 
}
template <class P>
void CellBox<P>::split() {
            set<int> used;
            if ( row1_ == row2_ && column1_ ==  column2_ ) 
                  return; 
            int row = (row2_+ row1_) /2;
            int column = (column2_+ column1_)/2;
            if (row2_- row1_ > 1&& column2_- column1_ > 1) {
                  //try first 2 split in columns ...
                  CellBox<P>* cell = new CellBox(parent_, row1_, row2_, column1_, column);
                  cell->levels(used);
                  if (used.size() < 2) {
                        push_back(cell);
                  }

                  else {
                        delete cell;
                        push_back(new CellBox(parent_, row1_, row, column1_, column));
                        push_back(new CellBox(parent_, row+1, row2_, column1_, column));
                  }
                  cell = new CellBox(parent_, row1_, row2_, column+1, column2_);
                  cell->levels(used);
                  if (used.size() < 2) {
                        push_back(cell);
                  } else {
                        delete cell;
                        push_back(new CellBox(parent_, row1_, row, column+1, column2_));
                        push_back(new CellBox(parent_, row+1, row2_, column+1, column2_));
                  }
                  return;
            }
            if (row2_- row1_ > 0) {
                  push_back(new CellBox(parent_, row1_, row, column1_, column2_));
                  push_back(new CellBox(parent_, row+1, row2_, column1_, column2_));
                  return;
            }
            if (column2_- column1_ > 0) {
                  push_back(new CellBox(parent_, row1_, row2_, column1_, column));
                  push_back(new CellBox(parent_, row1_, row2_, column+1, column2_));
                  return;
            }
            
          
            
            
            

      }






template<class P> 
void IsoPlot<P>::isoline(Cell<P>& cell, CellBox<P>* box) {
      static int castab[3][3][3] = { { { 0, 1, 2 }, { 3, 4, 5 }, { 6, 7, 8 } }, {
                  { 9, 10, 11 }, { 12, 13, 14 }, { 15, 16, 17 } }, { { 18, 19, 20 },
                  { 21, 22, 23 }, { 24, 25, 0 } } };
      int m1, m2, m3, case_value;
      double x1=0, x2=0, y1=0, y2=0;
      
      bool missing1, missing2;
      assert(!cell.flagged());
      
      
      
      std::set<int> used;
      cell.levels(used);
      if ( !box && used.size() < 2) return;
      
      if ( used.size() ==  1) {
            if (*used.begin() == -1 ) return;
            double contour=levels_[*used.begin()];
            Shape x;
            for (int i = 0; i < 4; i++) {
                              x.push_back(cell.column(i), cell.row(i));             
            }
            x.colour_ = (cell.value(0) > contour) ? this->shading_->rightColour(contour) : this->shading_->leftColour(contour);
            box->reshape(x);
            return;
      }
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      

      //=================================================================
      //
      // Note: at this stage the relative heights of the corners and the
      // centre are in the h array, and the corresponding coordinates are
      // in the xh and yh arrays. The centre of the box is indexed by 0
      // and the 4 corners by 1 to 4 as shown below.
      // Each triangle is then indexed by the parameter m, and the 3
      // vertices of each triangle are indexed by parameters m1,m2,and m3.
      // It is assumed that the centre of the box is always vertex 2
      // though this is important only when all 3 vertices lie exactly on
      // the same contour level, in which case only the side of the box
      // is drawn.
      //
      //      vertex 4 +-------------------+ vertex 3
      //               | \                 |
      //               |   \               |
      //               |     \         m=1 |
      //               |       \           |
      //               |         \         |      
      //               |           \       |
      //               |             \     |
      //               |    m=0        \   |
      //               |                 \ |
      //      vertex 1 +-------------------+ vertex 2
      //
      //               Scan each triangle in the box
      //=================================================================
      for (int m=0; m<2; m++) {
            m1 = m;
            m2 = m1+1;
            m3 = 3;
            
            
            //int count = 0;
            if (cell.missing(m1)) continue;
            if (cell.missing(m2)) continue;
            if (cell.missing(m3)) continue;
            
            
            
            
      
            
            vector<pair<Shape, Shape> > shapes;
      
            std::set<int>::const_iterator last = used.end();
            if (  used.size() > 1 ) --last;
            
      
            
            // Fisrt, build the list of isolines in this triangle...
            vector<int> rused;            
            vector<int> out;
            
            
            for (std::set<int>::const_iterator l = used.begin(); l != used.end(); ++l) {
                  if (*l == -1) continue; 
                  
                  double contour=levels_[*l];
                  
                  int out  = castab[cell.coef(m1, contour)][cell.coef(m2, contour)][cell.coef(m3, contour)];
                  
                  if ( out != 0 ) {
                        rused.push_back(*l);
                  }           
            }
      
            
            if ( rused.empty() && box ) {
                  Shape shape(this->shading_->colour(cell.value(m1)), cell.value(m1));
                  shape.push_back(cell.column(m1), cell.row(m1));
                  shape.push_back(cell.column(m2), cell.row(m2));
                  shape.push_back(cell.column(m3), cell.row(m3));
                  
                  box->reshape(shape);
            }
            bool complex = ( rused.size() > 1);
            for (vector<int>::const_iterator l = rused.begin(); l != rused.end(); ++l) {                 
                        int level = *l;                     
                        // First make a quich check to see if there is at l
                        
                        double contour=levels_[level];
                        Colour leftcolour = ( box ) ? this->shading_->leftColour(contour) : Colour("none"); 
                        Colour rightcolour = ( box ) ? this->shading_->rightColour(contour) : Colour("none"); 

            case_value = castab[cell.coef(m1, contour)][cell.coef(m2, contour)][cell.coef(m3, contour)];
            
            int add = 2;
            missing1 = false;
            missing2 = false;
            Shape leftcell(leftcolour, contour);
            Shape rightcell(rightcolour, contour);
            
      
      
            
            
            

            switch (case_value) {
            //===========================================================
            //     Case 0 -point out!
            //===========================================================

            case 0:
                  add = 0;
                  break;
                  if ( !box ) break;
                  
                  if (cell.height(m1, contour) < 0) {
                        leftcell.push_back(cell.column(m1), cell.row(m1));
                        leftcell.push_back(cell.column(m2), cell.row(m2));
                        leftcell.push_back(cell.column(m3), cell.row(m3));
                  }
                  else {
                        rightcell.push_back(cell.column(m1), cell.row(m1));
                        rightcell.push_back(cell.column(m2), cell.row(m2));
                        rightcell.push_back(cell.column(m3), cell.row(m3));
                  }
            
                  break;
                  //===========================================================
                  //     Case 1 -Single point 3
                  //===========================================================

            case 1:
                  missing1=cell.missing(m3);
                  x1=cell.column(m3);
                  y1=cell.row(m3);
                  add = 1;
                  if ( !box) break;
                  if (cell.height(m1, contour) < 0) {
                        leftcell.push_back(cell.column(m1), cell.row(m1));
                        leftcell.push_back(cell.column(m2), cell.row(m2));
                        leftcell.push_back(cell.column(m3), cell.row(m3));
                  }
                  else {
                        rightcell.push_back(cell.column(m1), cell.row(m1));
                        rightcell.push_back(cell.column(m2), cell.row(m2));
                        rightcell.push_back(cell.column(m3), cell.row(m3));
                  }
                  
                  break;
                  //===========================================================
                  //     Case 2 - Line between sides 2-3 and 3-1
                  //===========================================================
            case 2:
                  missing1 = cell.missing(m3) || cell.missing(m2);
                  missing2 = cell.missing(m3) || cell.missing(m1);

                  cell.xysect(m2,m3, contour, x1, y1);
                  cell.xysect(m3,m1, contour, x2, y2);
            
                  if ( !box ) break;
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);
                  leftcell.push_back(cell.column(m1), cell.row(m1));
                  leftcell.push_back(cell.column(m2), cell.row(m2));
                  
                  rightcell.push_back(x1, y1);
                  rightcell.push_back(cell.column(m3), cell.row(m3));
                  rightcell.push_back(x2, y2);

                  
                  break;
                  //===========================================================
                  //     Case 3 -Single point 2
                  //===========================================================
            case 3:
                  missing1 = cell.missing(m2);
                  x1=cell.column(m2);
                  y1=cell.row(m2);
                  add = 1;
                  if ( !box) break;
                  
                  if (cell.height(m1, contour) < 0) {
                                          leftcell.push_back(cell.column(m1), cell.row(m1));
                                          leftcell.push_back(cell.column(m2), cell.row(m2));
                                          leftcell.push_back(cell.column(m3), cell.row(m3));
                                    }
                                    else {
                                          rightcell.push_back(cell.column(m1), cell.row(m1));
                                          rightcell.push_back(cell.column(m2), cell.row(m2));
                                          rightcell.push_back(cell.column(m3), cell.row(m3));
                                    }
                              
                  break;
                  //===========================================================
                  //     Case 4 - Line between vertex 2 and vertex 3
                  //===========================================================
            case 4:
                  missing1 = cell.missing(m2);
                  missing2 = cell.missing(m3);
                  x1=cell.column(m2);
                  y1=cell.row(m2);
                  x2=cell.column(m3);
                  y2=cell.row(m3);
                  if ( !box ) break;
                  
                                          leftcell.push_back(cell.column(m1), cell.row(m1));
                                          leftcell.push_back(cell.column(m2), cell.row(m2));
                                          leftcell.push_back(cell.column(m3), cell.row(m3));
                        
            

                  break;
                  //===========================================================
                  //     Case 5 - Line between vertex 2 and side 3-1
                  //===========================================================
            case 5:
                  missing1=cell.missing(m2);
                  x1=cell.column(m2);
                  y1=cell.row(m2);
                  missing2 = cell.missing(m3) || cell.missing(m1);
                  cell.xysect(m3,m1, contour,x2, y2);
                  if ( !box ) break;
                  leftcell.push_back(cell.column(m1), cell.row(m1));
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);
                  
                  rightcell.push_back(x1, y1);
                  rightcell.push_back(cell.column(m3), cell.row(m3));
                  rightcell.push_back(x2, y2);

                  
                  break;
                  //===========================================================
                  //     Case 6 - Line between side 1-2  and side 2.3
                  //===========================================================
            case 6:
                  missing1 = cell.missing(m1) || cell.missing(m2);
                  missing2 = cell.missing(m2) || cell.missing(m3);
                  cell.xysect(m1,m2, contour, x1, y1);
                  cell.xysect(m2,m3, contour, x2, y2);
      
                  if ( !box ) break;
                  leftcell.push_back(cell.column(m1), cell.row(m1));
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);
                  leftcell.push_back(cell.column(m3), cell.row(m3));
                  
                  rightcell.push_back(x1, y1);
                  rightcell.push_back(cell.column(m2), cell.row(m2));
                  rightcell.push_back(x2, y2);
                  
                  break;
                  //===========================================================
                  //     Case 7 - Line between sides 1-2 aanve vertex 3
                  //===========================================================
            case 7:
                  missing2=cell.missing(m3);
                  missing1 = cell.missing(m1) || cell.missing(m2);
                  cell.xysect(m1,m2, contour, x1, y1);
                  x2=cell.column(m3);
                  y2=cell.row(m3);
                  
                  if ( !box ) break;
                  leftcell.push_back(cell.column(m1), cell.row(m1));
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);

                  
                  rightcell.push_back(x1, y1);
                  rightcell.push_back(cell.column(m2), cell.row(m2));
                  rightcell.push_back(x2, y2);
                  
                  break;
                  //===========================================================
                  //     Case 8 - Line between sides 1-2 and 3-1
                  //===========================================================
            case 8:
                  missing1 = cell.missing(m1) || cell.missing(m2);
                  missing2 = cell.missing(m3) || cell.missing(m1);
                  cell.xysect(m1,m2, contour, x1, y1);
                  cell.xysect(m3,m1, contour, x2, y2);
                  
                  if ( !box ) break;
                  leftcell.push_back(cell.column(m1), cell.row(m1));
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);

            
                  rightcell.push_back(x1, y1);
                  rightcell.push_back(cell.column(m2), cell.row(m2));
                  rightcell.push_back(cell.column(m3), cell.row(m3));
                  rightcell.push_back(x2, y2);
                  
                  break;
                  //===========================================================
                  //     Case 9 -single point 1
                  //===========================================================
            case 9:
                  missing1=cell.missing(m1);
                  x1=cell.column(m1);
                  y1=cell.row(m1);
                  add = 1;
                  if ( !box) break;
                  if (cell.height(m2, contour) < 0) {
                                          leftcell.push_back(cell.column(m1), cell.row(m1));
                                          leftcell.push_back(cell.column(m2), cell.row(m2));
                                          leftcell.push_back(cell.column(m3), cell.row(m3));
                                    }
                                    else {
                                          rightcell.push_back(cell.column(m1), cell.row(m1));
                                          rightcell.push_back(cell.column(m2), cell.row(m2));
                                          rightcell.push_back(cell.column(m3), cell.row(m3));
                                    }
                  break;
                  //===========================================================
                  //     Case 10 - Line between vertex 3 and vertex 1
                  //===========================================================
            case 10:
                  missing1=cell.missing(m3);
                  missing2=cell.missing(m1);
                  x1=cell.column(m3);
                  y1=cell.row(m3);
                  x2=cell.column(m1);
                  y2=cell.row(m1);
                  if ( !box) break;

                  
                                                            leftcell.push_back(cell.column(m1), cell.row(m1));
                                                            leftcell.push_back(cell.column(m2), cell.row(m2));
                                                            leftcell.push_back(cell.column(m3), cell.row(m3));
                                                      
            
                  break;
                  //===========================================================
                  //     Case 11 - Line between side 3-2 and vertex 1
                  //===========================================================
            case 11:
                  missing2=cell.missing(m1);
                  missing1 = cell.missing(m3) || cell.missing(m2);
            
                  cell.xysect(m3,m2, contour, x1, y1);
                  x2=cell.column(m1);
                  y2=cell.row(m1);
                  if ( !box ) break;
                  leftcell.push_back(cell.column(m1), cell.row(m1));
                  leftcell.push_back(cell.column(m2), cell.row(m2));
                  leftcell.push_back(x1, y1);
                  
                  rightcell.push_back(x1, y1);
                  rightcell.push_back(cell.column(m3), cell.row(m3));
                  rightcell.push_back(x2, y2);
            
                  break;
                  //===========================================================
                  //     Case 12 - Line between vertex 1 and vertex 2
                  //===========================================================
            case 12:
                  missing1=cell.missing(m1);
                  missing2=cell.missing(m2);
                  x1=cell.column(m1);
                  y1=cell.row(m1);
                  x2=cell.column(m2);
                  y2=cell.row(m2);
                  if ( !box ) break;
            
                                                            leftcell.push_back(cell.column(m1), cell.row(m1));
                                                            leftcell.push_back(cell.column(m2), cell.row(m2));
                                                            leftcell.push_back(cell.column(m3), cell.row(m3));
                                          
            
                  break;
                  //===========================================================
                  //     Case 13 - Flat Area all vertex have the isoline value!
                  //===========================================================
            case 13:
                  add = 0;
                  if ( !box ) break;
                  if (cell.height(m1, contour) < 0) {
                                                            leftcell.push_back(cell.column(m1), cell.row(m1));
                                                            leftcell.push_back(cell.column(m2), cell.row(m2));
                                                            leftcell.push_back(cell.column(m3), cell.row(m3));
                                                      }
                                                      else {
                                                            rightcell.push_back(cell.column(m1), cell.row(m1));
                                                            rightcell.push_back(cell.column(m2), cell.row(m2));
                                                            rightcell.push_back(cell.column(m3), cell.row(m3));
                                                      }
            
                  break;
                  //===========================================================
                  //     Case 14 - Line between vertex 2 and vertex 1
                  //===========================================================
            case 14:
                  missing1=cell.missing(m2);
                  missing2=cell.missing(m1);
                  x1=cell.column(m2);
                  y1=cell.row(m2);
                  x2=cell.column(m1);
                  y2=cell.row(m1);
                  if ( !box ) break;
                  
                                                            rightcell.push_back(cell.column(m1), cell.row(m1));
                                                            rightcell.push_back(cell.column(m2), cell.row(m2));
                                                            rightcell.push_back(cell.column(m3), cell.row(m3));
                              
                  
                  break;
                  //===========================================================
                  //     Case 15 - Line between vertex 1 and side3-2
                  //===========================================================
            case 15:
                  missing1=cell.missing(m1);
                  x1=cell.column(m1);
                  y1=cell.row(m1);
                  
                  cell.xysect(m3,m2, contour, x2, y2);
                  
                  missing2 = cell.missing(m3) || cell.missing(m2);
                  if ( !box ) break;
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);
                  leftcell.push_back(cell.column(m3), cell.row(m3));
                  
                  rightcell.push_back(x1, y1);
                  rightcell.push_back(cell.column(m2), cell.row(m2));
                  rightcell.push_back(x2, y2);
                  
                  break;
                  //===========================================================
                  //     Case 16 - Line between vertex 1 and vertex 3
                  //===========================================================
            case 16:
                  missing1=cell.missing(m1);
                  missing2=cell.missing(m3);
                  x1=cell.column(m1);
                  y1=cell.row(m1);
                  x2=cell.column(m3);
                  y2=cell.row(m3);
                  if ( !box ) break;
                  rightcell.push_back(x1, y1);
                  rightcell.push_back(cell.column(m2), cell.row(m2));
                  rightcell.push_back(x2, y2);
                  
                  
            
                  break;
                  //===========================================================
                  //     Case 17 - single point 1
                  //===========================================================
            case 17:
                  missing1=cell.missing(m1);
                  x1=cell.column(m1);
                  y1=cell.row(m1);
                  add = 1;
                  if ( !box ) break;
                  if (cell.height(m2, contour) < 0) {
                                                            leftcell.push_back(cell.column(m1), cell.row(m1));
                                                            leftcell.push_back(cell.column(m2), cell.row(m2));
                                                            leftcell.push_back(cell.column(m3), cell.row(m3));
                                                      }
                                                      else {
                                                            rightcell.push_back(cell.column(m1), cell.row(m1));
                                                            rightcell.push_back(cell.column(m2), cell.row(m2));
                                                            rightcell.push_back(cell.column(m3), cell.row(m3));
                                                      }
                  
                  break;
                  //===========================================================
                  //     Case 18 - Line between side3-1 and side1-2
                  //===========================================================
            case 18:
            
                  cell.xysect(m3,m1, contour, x1, y1);
            
                  cell.xysect(m1,m2, contour, x2, y2);
                  missing1 = cell.missing(m3) || cell.missing(m1);
                  missing2 = cell.missing(m1) || cell.missing(m2);
                  if ( !box ) break;
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);
                  leftcell.push_back(cell.column(m2), cell.row(m2));
                  leftcell.push_back(cell.column(m3), cell.row(m3));
            
                  rightcell.push_back(x1, y1);

                  rightcell.push_back(cell.column(m1), cell.row(m1));
                  rightcell.push_back(x2, y2);
            
                  break;
                  //===========================================================
                  //     Case 19 - Line between vertex 3 and side1-2
                  //===========================================================
            case 19:
                  missing1=cell.missing(m3);
                  x1=cell.column(m3);
                  y1=cell.row(m3);
            
                  cell.xysect(m1,m2, contour, x2, y2);
                  missing2 = cell.missing(m1) || cell.missing(m2);
                  if ( !box ) break;
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);
                  leftcell.push_back(cell.column(m2), cell.row(m2));
            
                  rightcell.push_back(x1, y1);

                  rightcell.push_back(cell.column(m1), cell.row(m1));
                  rightcell.push_back(x2, y2);
                  
                  break;

                  //===========================================================
                  //     Case 20 - Line between side3-2 and side1-2
                  //===========================================================
            case 20:
            
                  cell.xysect(m3,m2, contour, x1, y1);
            
                  cell.xysect(m1,m2, contour, x2, y2);
                  missing1 = cell.missing(m3) || cell.missing(m2);
                  missing2 = cell.missing(m1) || cell.missing(m2);
                  if ( !box ) break;
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);
                  leftcell.push_back(cell.column(m2), cell.row(m2));
            
                  rightcell.push_back(x1, y1);

                  rightcell.push_back(cell.column(m3), cell.row(m3));
                  rightcell.push_back(cell.column(m1), cell.row(m1));
                  rightcell.push_back(x2, y2);
                  
                  break;
                  //===========================================================
                  //     Case 21 - Line between side 3-1 and vertex 2 
                  //===========================================================
            case 21:
                  missing2=cell.missing(m2);
                  
                  cell.xysect(m3,m1, contour, x1, y1);
                  missing1 = cell.missing(m3) || cell.missing(m2);
                  x2=cell.column(m2);
                  y2=cell.row(m2);
                  if ( !box ) break;
                  
                        leftcell.push_back(x1, y1);
                        leftcell.push_back(x2, y2);
                        leftcell.push_back(cell.column(m3), cell.row(m3));
                        
                        rightcell.push_back(x1, y1);
                        rightcell.push_back(cell.column(m1), cell.row(m1));
                        rightcell.push_back(cell.column(m2), cell.row(m2));
                        
                  
                  
                  break;
                  //===========================================================
                  //     Case 22 - Line between vertex 3 and vertex 2
                  //===========================================================
            case 22:
                  missing1=cell.missing(m3);
                  missing2=cell.missing(m2);
                  x1=cell.column(m3);
                  y1=cell.row(m3);
                  x2=cell.column(m2);
                  y2=cell.row(m2);
                  if ( !box) break;
                  rightcell.push_back(cell.column(m1), cell.row(m1));
                  rightcell.push_back(cell.column(m2), cell.row(m2));
                  rightcell.push_back(cell.column(m3), cell.row(m3));
                  
                  
                  break;
                  //===========================================================
                  //     Case 23 - single point 2
                  //===========================================================
            case 23:
                  missing1=cell.missing(m2);
                  x1=cell.column(m2);
                  y1=cell.row(m2);
                  add = 1;
                  if ( !box) break;
                  if (cell.height(m1, contour) < 0) {
                                                            leftcell.push_back(cell.column(m1), cell.row(m1));
                                                            leftcell.push_back(cell.column(m2), cell.row(m2));
                                                            leftcell.push_back(cell.column(m3), cell.row(m3));
                                                      }
                                                      else {
                                                            rightcell.push_back(cell.column(m1), cell.row(m1));
                                                            rightcell.push_back(cell.column(m2), cell.row(m2));
                                                            rightcell.push_back(cell.column(m3), cell.row(m3));
                                                      }
                  break;
                  //===========================================================
                  //     Case 24 - Line between side1-3  and side3-2
                  //===========================================================
            case 24:
                  missing1 = cell.missing(m1) || cell.missing(m3);
                  missing2 = cell.missing(m3) || cell.missing(m2);
                  cell.xysect(m1,m3, contour, x1, y1);
                  cell.xysect(m3,m2, contour, x2, y2);
                  if ( !box) break;
                  leftcell.push_back(x1, y1);
                  leftcell.push_back(x2, y2);
                  leftcell.push_back(cell.column(m3), cell.row(m3));
                  
                  rightcell.push_back(x1, y1);

                  rightcell.push_back(cell.column(m1), cell.row(m1));
                  rightcell.push_back(cell.column(m2), cell.row(m2));
                  rightcell.push_back(x2, y2);
                  
                  break;
                  //===========================================================
                  //     Case 25 - single point C
                  //===========================================================
            case 25:
                  missing1= cell.missing(m3);
                  x1=cell.column(m3);
                  y1=cell.row(m3);
                  add = 1;
                  if ( !box) break;
                  
                   
                  if (cell.height(m1, contour) < 0) {
                                                                  leftcell.push_back(cell.column(m1), cell.row(m1));
                                                                  leftcell.push_back(cell.column(m2), cell.row(m2));
                                                                  leftcell.push_back(cell.column(m3), cell.row(m3));
                                                            }
                                                            else {
                                                                  rightcell.push_back(cell.column(m1), cell.row(m1));
                                                                  rightcell.push_back(cell.column(m2), cell.row(m2));
                                                                  rightcell.push_back(cell.column(m3), cell.row(m3));
                                                            }
                  
                  
                  break;
                  
            default:
                  break;
            }

            
            if (add == 2 ) {
            
                        bool done = false;
                        for (typename VectorOfPointers<vector<Polyline* > >::iterator
                                    line = helper_[contour].begin(); line
                                    != helper_[contour].end(); ++line) {
                              // try to add to the back ! 
                              if ((*line)->empty() ) {
                                    (*line)->push_back(PaperPoint(x1, y1, contour, missing1));
                                    (*line)->push_back(PaperPoint(x2, y2, contour, missing2));
                                    done = true;
                                    break;
                              }
                              if (same((*line)->back().x(), x1) && same((*line)->back().y(), y1) ) {
                                    
                                    (*line)->push_back(PaperPoint(x2, y2, contour, missing2));
                                    // Check if we can concatenate with another line...
                                    concatenate_back(helper_[contour], *line);
                                    done = true;
                                    break;
                              }
                              if (same((*line)->front().x(), x2) && same((*line)->front().y(), y2) ) {
                                    (*line)->push_front(PaperPoint(x1, y1, contour, missing1));
                                    concatenate_front(helper_[contour], *line);
                                    done = true;
                                    break;
                              }
                        }
      
                              if ( !done ) {
                                    helper_[contour].push_back(new Polyline());
                                    helper_[contour].back()->push_back(PaperPoint(x1, y1, contour, missing1));
                                    helper_[contour].back()->push_back(PaperPoint(x2, y2, contour, missing2));
                              }
                  
                  }
                  
                if ( !box) continue;
                
                if ( !complex ) {
                              box->reshape(leftcell);
                              box->reshape(rightcell);
                        }
                else {
                  //if ( !leftcell.empty() && !rightcell.empty() ) {
                        shapes.push_back(make_pair(leftcell, rightcell)); 
                  //}
                  //else if ( !rightcell.empty() ) {
                  //    box->reshape(rightcell);
                  //}
                  //else if ( !leftcell.empty() ) {
                  //    box->reshape(leftcell);
                  
                }
                            
            } // end of levels...
            
            //
            if ( shapes.empty() ) continue;
            
            // Now we reduce the shape! 
            vector<pair<Shape, Shape> >::iterator current = shapes.begin();   
            vector<pair<Shape, Shape> >::iterator next = shapes.begin();
            
            box->reshape(current->first);
            next++;
            
            while ( next != shapes.end() ) {
                  current->second.intersection(next->second);
                  box->reshape(current->second);
                  current++;
                  next++ ;
            }

            box->reshape(current->second);
            
            
      } // next triangle
                  
}

template<class P> 
void IsoPlot<P>::isoline(MatrixHandler<P>& data, BasicGraphicsObjectContainer& parent)  {

      const Transformation& transformation = parent.transformation();
      levels_.clear();
      // Find the used levels!
      helper_ = map<double, vector<Polyline* > >();
      const vector<double>::const_iterator end = (*this->levelSelection_).end();
      vector<double>::const_iterator previous = (*this->levelSelection_).end();
      vector<double>::const_iterator last = (*this->levelSelection_).end();
      
      if ( (*this->levelSelection_).front()  > data.min() )
            levels_.push_back(data.min());
      
      for (vector<double>::const_iterator l = (*this->levelSelection_).begin(); l != end; ++l) {
            if ( (*l <= data.min()) ) {
                  previous = l;
                  continue;
            }
            if ( (*l >= data.max()) ) {
                         last = l;
                              continue;
                              
            }
            if ( previous != (*this->levelSelection_).end() ) {
                   levels_.push_back(*previous);
                   previous = (*this->levelSelection_).end();
            }
            levels_.push_back(*l);
      }
      
      if ( last != (*this->levelSelection_).end())
                         levels_.push_back(*last);
      else 
             levels_.push_back(data.max());
                        
      
      
      missing_ = data.missing();

      //===========================================================================
      // Note that castab is arranged differently from the FORTRAN code because
      // Fortran and C/C++ arrays are transposed of each other, in this case
      // it is more tricky as castab is in 3 dimension
      //===========================================================================


      IntervalMap<int> range;
      int r= 1;
      for (vector<double>::const_iterator level = levels_.begin(); level
                  != levels_.end(); ++level) {

            helper_.insert(make_pair(*level, vector<Polyline*>()));
            if (level+1!= levels_.end() )
                  range.insert(make_pair(Interval(*level, *(level+1)), r++));
      }
      
      CellArray<P> array(data, range, transformation);

      CellBox<P> view(&array);

       ( this->shading_->shadingMode() ) ? view.shade(*this) : view.contour(*this);

     view.feed(*this,parent);

       
      for (typename map<double, vector<Polyline* > >::iterator lines = helper_.begin(); lines != helper_.end(); ++lines)
            for (typename vector<Polyline* >::iterator line = lines->second.begin(); line != lines->second.end(); ++line)
                  this->push_back(*line);
   
}

template<class P> bool IsoPlot<P>::prepare(MatrixHandler<P>& data) {
      PointsHandler<P> points(data);
      (*this->levelSelection_).clear();
      (*this->levelSelection_).calculate(points.min() , points.max() , this->shading_->shadingMode());
      (*this->label_).prepare(*this->levelSelection_, (*this->colour_).name());
      return (*this->shading_)(*this->levelSelection_);
}

template<class P> void IsoPlot<P>::operator()(MatrixHandler<P>& data, BasicGraphicsObjectContainer& parent)
{
      // Get the triangles list ...
      // Create the isolines...    
      
      
      prepare(data);
      {
            Timer timer("contouring", "Time spent in contouring");
            isoline(data, parent);
      }
      

      (*this->shading_)(data, parent);
      (*this->highlight_).prepare(*this->levelSelection_);

      // Now we feed the task...   
      
      for (typename vector<Polyline* >::const_iterator poly = this->begin(); poly != this->end(); ++poly)
      {
            if ( (*poly)->empty() ) continue;
            (*poly)->setLineStyle(this->style_);
            (*poly)->setThickness(this->thickness_);
            (*poly)->setColour(*this->colour_);

            (*this->highlight_)(*(*poly));
            (*this->label_)(**poly, (*poly)->front().value());

            parent.push_back(*poly);
            
      
      }
      
      


      this->clear();
}

template<class P> void NoIsoPlot<P>::operator()(MatrixHandler<P>& data, BasicGraphicsObjectContainer& parent)
{
      // Create the isolines...        
      if ( !prepare(data) ) {
            (*this->shading_)(data, parent);
            // do not send the isolines...  
            return;
      }
      // The shading needs the isolines..
      // WE will calculate them but will not send them to the printer     
      {
            Timer timer("contouring", "Time spent in contouring");
            isoline(data, parent);
      }

      (*this->shading_)(data, parent);

      typename vector<Polyline* >::const_iterator polyend = this->end();
      for (typename vector<Polyline* >::iterator poly = this->begin(); poly != polyend; ++poly) {
            if ( (*poly)->empty() ) continue;
            (*poly)->setThickness(0);
            (*this->label_)(**poly, (*poly)->front().value());
            parent.push_back(*poly);
      }
      this->clear();
}

template<class P> void IsoPlot<P>::visit(LegendVisitor& legend) {
      (*this->shading_).visit(legend);
}

Generated by  Doxygen 1.6.0   Back to index