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

IsoPlot.h

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.h
    \brief Definition of the Template class IsoPlot.
    
    Magics Team - ECMWF 2004
    
    Started: Wed 3-Mar-2004
    
    Changes:
    
*/

#ifndef IsoPlot_H
#define IsoPlot_H

#include "magics.h"

#include "IsoPlotAttributes.h"
#include "BasicSceneObject.h"
#include "Polyline.h"
#include "VectorOfPointers.h"
#include "XmlNode.h"

namespace magics {


struct ShapeLine 
{
      ShapeLine(double p1x, double p1y, double p2x, double p2y) : p1x_(p1x), p1y_(p1y), p2x_(p2x), p2y_(p2y) {}
      ShapeLine(double p1x, double p1y) : p1x_(p1x), p1y_(p1y), p2x_(0), p2y_(0) {}
      double p1x_;
      double p1y_;
      double p2x_;
      double p2y_;      
      bool operator == (const ShapeLine& other) const 
      {
                  return ( p1x_ == other.p2x_ && p2x_ == other.p1x_  && p1y_ == other.p2y_ && p2y_ == other.p1y_ );                             
      }
      bool same(const ShapeLine& other) const {
            return  ( p1x_ == other.p1x_  && p1y_ == other.p1y_  );
                              
      }
      friend ostream& operator<<(ostream& s,const ShapeLine& p) 
                  { s << "[(" << p.p1x_ << ", " << p.p1y_ << ")->("<< p.p2x_ << ", " << p.p2y_ << ")\n"; return s;  }
};




struct Shape : public list<ShapeLine>
{
      Shape() : colour_("none"), 
      minx_(numeric_limits<double>::max()), 
      maxx_(-numeric_limits<double>::max()), 
      miny_(numeric_limits<double>::max()), 
      maxy_(-numeric_limits<double>::max()) {}
      Shape(const Colour& colour, double val) : 
                        value_(val),
                        colour_(colour), 
                        minx_(numeric_limits<double>::max()), 
                        maxx_(-numeric_limits<double>::max()), 
                        miny_(numeric_limits<double>::max()), 
                        maxy_(-numeric_limits<double>::max()) {}
      ~Shape() {}
      void push_back(double x, double y) {
            points_.push_back(make_pair(x, y));
            if ( !empty () ) {
                  if ( back().p1x_ == x && back().p1y_ ==  y      ) 
                        return;
                  back().p2x_ = x;
                  back().p2y_ = y;                          
                  list<ShapeLine>::push_back(ShapeLine(x, y, front().p1x_, front().p1y_));
            }           
            else {
                  list<ShapeLine>::push_back(ShapeLine(x, y));
            }
            if ( x < minx_ ) minx_ = x;
            if ( x > maxx_ ) maxx_ = x;
            if ( y < miny_ ) miny_ = y;
            if ( y > maxy_ ) maxy_ = y;
      }
      double value_;
      Colour  colour_;
      double minx_;
      double maxx_;
      double miny_;     
      double maxy_;
      vector<pair<double, double> > points_;
      
      bool intersect(Shape& other) {
            if ( minx_ > other.maxx_ ) return false;
            if ( maxx_ < other.minx_ ) return false;
            if ( miny_ > other.maxy_ ) return false;
            if ( maxy_ < other.miny_ ) return false;
            return true;
      }
      
      double surface() const {            
            double area = 0;
            for (const_iterator p = begin(); p != end(); ++p) {
                  area += (p->p1x_ * p->p2y_) - (p->p2x_ * p->p1y_);
            }
            return area;
      }
      
      void clean() {
            points_.clear();
            clear();
            minx_ = numeric_limits<double>::max();
            maxx_ = -numeric_limits<double>::max();
            miny_ = numeric_limits<double>::max();
            maxy_ = -numeric_limits<double>::max();
      }
      
      void reduce() {
             vector<pair<double, double>  >points;                                     
                                           stack<pair<double, double> > helper;
            vector<pair<double, double>  >::const_iterator point = points_.begin();
                                          points.push_back(*point);
                                          while ( point != points_.end() ) {
                                                       if (  *point != points.back()  ) {
                                                                  points.push_back(*point);
                                                       }
                                                       ++point;
                                       }
                                          pair<double, double> top;
                                       for (point = points.begin();  point != points.end(); ++point) {
                                                            if ( helper.empty() ) {
                                                                       helper.push(*point);
                                                                       continue;
                                                            }
                                                                                             
                                                            top =  helper.top();
                                                            helper.pop();
                                                            if ( helper.empty() ) { 
                                                                          helper.push(top);
                                                                          if ( top != *point ) {                                                               
                                                                                    helper.push(*point);
                                                                           }
                                                          }
                                                            else if ( !( *point == helper.top()) ) {
                                                                          helper.push(top);
                                                                          helper.push(*point);
                                                          }
                                          }

                                                      
                                                   clean();
                                                      // rebuild!
                                                   //cout << "reduced----------------->" << endl;
                                                      while ( !helper.empty() ) {
                                                            //cout << "[" << helper.top().first << ", " <<  helper.top().second << "]" << endl;
                                                            push_back(helper.top().first, helper.top().second);
                                                            helper.pop();
                                                      }
                                                      //cout << "<-----------------" << endl;
            
      }
      
      void intersection(Shape& other) {
            //reduce();
            //other.reduce();
            
            vector<pair<double, double > > opoints = other.points_;
            
            std::reverse(opoints.begin(), opoints.end());
            vector<pair<double, double>  >::iterator p1= points_.begin();
            
            
            //colour_ = surface() > other.surface() ? colour_ : other.colour_;

            while (p1!= points_.end() ) {
                        vector<pair<double, double>  >::iterator p2= opoints.begin();
                    while ( p2!= opoints.end()) {
                        if ( *p1 == *p2) {
                                     std::rotate(opoints.begin(), p2, opoints.end());
                                     p1++;
                                     
                                   points_.insert(p1, opoints.begin(), opoints.end());
                                     
                                    //now we need to take off the duplicates! 
                                             
                                   stack<pair<double, double> > helper;
                                   vector<pair<double, double>  >::const_iterator point = points_.begin();
                                       
                                   pair<double, double> top;
                                                                     for (point = points_.begin();  point != points_.end(); ++point) {
                                                                                          if ( helper.empty() ) {
                                                                                                     helper.push(*point);
                                                                                                     continue;
                                                                                          }
                                                                                                                           
                                                                                          top =  helper.top();
                                                                                          helper.pop();
                                                                                          
                                                                                                        if ( top != *point ) {
                                                                                                                    helper.push(top);
                                                                                                                  helper.push(*point);
                                                                                                         }
                                                                                    
                                                                        }

                                   
                                   
                                  

                                                                      clean();//cout << "new shape ----------------->" << endl;
                                                                       // rebuild!
                                                                    
                                                                                                                    
                                                                                                                        while ( !helper.empty() ) {
                                                                                                                              //cout << "[" << setprecision(4) << helper.top().first << ", " <<  helper.top().second << "]" << endl;
                                                                                                                              push_back(helper.top().first, helper.top().second);
                                                                                                                              helper.pop();
                                                                                                                        }
                                                                                                                        //cout << "<-----------------" << endl;
                                                                                                                       
                                                                       
                                                
                                     return;
                        }     
                        p2++;
                    }
                        p1++;                   
            }
      
            
      }
      
      bool operator<(const Shape& other) {
            return this->surface() < other.surface();
      }
      bool operator!=(const Shape& other) {
                  return this->surface() != other.surface();
      }
      
      void  reduce(list<pair<double, double> >& shape)
      {
            if ( shape.size() < 3 ) return; // nothing to reduce!
          list<pair<double, double> >::iterator first = shape.begin();
            list<pair<double, double> >::iterator second = shape.begin();
            ++second;
          while(second != shape.end()) {
             
                  if ( *first == *second) {
                  shape.remove(*first);              
                  reduce(shape);
                  return;
              }
              ++first;      
              ++second;
          }
      }

      
      void rebuild( list<pair<double, double> >& shape)
      {
            this->clear();
                  for ( list<pair<double, double> >::iterator p1 = shape.begin(); p1 != shape.end(); ++p1) {
                        this->push_back(p1->first, p1->second);
                  }     
      }

      friend ostream& operator<<(ostream& s,const Shape& p) 
      {
        s << "Shape[\n";
        for (const_iterator i = p.begin(); i != p.end(); ++i) 
                    s << "\t" << *i << endl;
        s<< "]\n";
        return s; 
      }
};



struct x_function
{
      bool operator()(const PaperPoint& p1, const PaperPoint& p2) {
            return p1.x() < p2.x();
      }
};

struct y_functiong
{
      bool operator()(const PaperPoint& p1, const PaperPoint& p2) {
            return p1.y() < p2.y();
      }
};




template <class P > class Cell;

template<class P> class CellArray : public VectorOfPointers<vector<Cell<P>* > > {
public:
      CellArray(MatrixHandler<P>& data, IntervalMap<int>& range, const Transformation& transformation);
      int rows_;
      int columns_;
      double missing_;
      
      vector<PaperPoint> points_;
      IntervalMap<int> rangeFinder_;

      Cell<P>* operator()(int row, int column) {
            return (*this)[row*columns_ + column];
      }
      
      PaperPoint& point(pair<int, int> pos) {
            return points_[pos.first*(columns_+1) + pos.second];
      }
      
      ~CellArray() { }
};

template <class P>
class Cell 
{
public:     
      Cell(CellArray<P>& parent) : flagged_(false), parent_(parent), missing_(parent.missing_), depth_(0) { } 
      Cell(CellArray<P>& parent, int row, int column):  flagged_(false),  parent_(parent), row_(row), column_(column) ,  depth_(0) {
            minx_ = maxx_ = this->column(0);          
            miny_ = maxy_ = this->row(0);
            
            min_ = max_ =     range(0);
            addOut_ =   (min_ == -1); 
            out_ = true;
            for ( int i = 1; i < 4; i++ ) {
                  
                  if (  this->column(i)  < minx_ ) minx_ = this->column(i); 
                  if (  this->column(i)  > maxx_ ) maxx_ = this->column(i); 
                  if (  this->row(i)  < miny_ ) miny_ = this->row(i); 
                  if (  this->row(i)  > maxy_ ) maxy_ = this->row(i);   
                  if ( range(i) == -1) {
                        addOut_ = true;
                        continue;
                  }
                  out_ = false;
                  if (  range(i)  < min_ ) min_ =   range(i); 
                  if (  range(i)  > max_ ) max_ = range(i);       
            }
        if ( !out_ ) 
            out_ =  ( minx_ == maxx_ || miny_ == maxy_) ? true : false;
        missing_ = parent_.missing_;
      }
      
      virtual ~Cell() {  }
      
      bool flagged_; 
      CellArray<P>& parent_;
      int row_;
      int column_;
      bool out_;
      
      bool flagged() { bool flag = flagged_; flagged_ = true; return flag; }
      bool out () { return out_; }
    
      bool in(int l) { return (min_ <= l && l <= max_); }
      bool in(int l, int m1, int m2, int m3) {
      int min = std::min(range(m1), range(m2));
      min = std::min(min, range(m3));
      int max = std::max(range(m1), range(m2));
      max = std::max(max, range(m3));
      return  (min <= l && l <= max);
    }
      
      virtual void levels(set<int>& levels) { 
      
      
                if (out_) return;
                
                if ( addOut_ ) levels.insert(-1);
            
                for ( int i = min_; i <= max_; i++) {
                  levels.insert(i);
                }
      }
      
      virtual bool inrange() {
            return (this->range(0) != -1);
      }
      
      double missing_;
      
      bool missing(double val) { return same(missing_, val); }
      
      
      int min_;
      int max_;
      double minx_;
      double  miny_;
      double maxx_;
      double maxy_;
      bool addOut_;
      int depth_;
      
      virtual bool equalRange() {
            return min_ == max_;
      }
      
      virtual bool equalRange(int m1, int m2, int m3) {
            return ( range(m1) == range(m2) && range(m1) == range(m3) );
      }
      
      virtual int range() { return min_; }
      
      
      
      pair<int, int> index(int i) {
            
            switch (i) {
                  case 0 : return make_pair(row_, column_);
                  case 1 : return make_pair(row_, column_+1);
                  case 2 : return make_pair(row_+1, column_+1);
                  case 3 : return make_pair(row_+1, column_);
            }
            assert(false);
            return make_pair(-1, -1);
      }
      
      virtual double value(int i) {
            return this->parent_.point(this->index(i)).value();         
            
      }
      
      virtual double height(int i, double val) {
                
                  return missing(i) ? 0 : this->parent_.point(this->index(i)).value() -val;           
                  
      }
      virtual int coef(int i, double val) {
            double height = this->parent_.point(this->index(i)).value() -val;       
            
            return ( height ) ? (height/abs(height))+1 : 1;       
      }

      virtual bool missing(int i) {
            return (same(missing_, this->parent_.point(this->index(i)).value()));         
      }
      
      virtual double column(int i) {
            return parent_.point(index(i)).x();
            //return 0;
                  
            }
      void xysect(int i, int j, double value, double& x, double& y) {
            PaperPoint p1 =  this->parent_.point(this->index(i));
            PaperPoint p2 =  this->parent_.point(this->index(j));
            
            double v1 = p1.value()-value;
            double v2 = p2.value()-value;
            x = (v2*p1.x()-v1*p2.x())/(v2-v1);
            y = (v2*p1.y()-v1*p2.y())/(v2-v1);
      }
      
      virtual double row(int i) {
            return parent_.point(index(i)).y();
            //return 0;
      }           
      virtual int range(int i) {
                  return parent_.point(index(i)).range();
                  //return 0;
            }           
      virtual int findRange(double value) {
            return parent_.rangeFinder_.find(value, -1);
      }
};


template <class P>
CellArray<P>::CellArray(MatrixHandler<P>& data, IntervalMap<int>& range, const Transformation& transformation):
       rangeFinder_(range) {
            //data_(data),, transformation_(transformation) 
            
            points_.reserve(data.rows()* data.columns());         
            this->reserve(data.rows()* data.columns());     
            rows_ = data.rows()-1;
            columns_ = data.columns()-1;
            this->reserve(rows_*columns_);
            int i = 0;
            missing_ = data.missing();
            for (int row = 0; row <= rows_; row++)
                        for (int column = 0; column <= columns_; column++) {
                              
                              double value = data(row, column);
                              int r =  range.find(data(row, column), -1);
                              points_[i] = transformation(P(data.column(row, column), data.row(row, column), data(row, column)));
                              points_[i].range_ = r;
                              if ( transformation.clip(points_[i]) ) {
                                    P point;
                                    transformation.revert(PaperPoint(points_[i].x(), points_[i].y()), point);                                                   
                                    
                                    value =  data.interpolate(point.y(), point.x());                  
                                    if (value < data.min() ) value = data.min();
                                    if (value > data.max() ) value = data.max();
                                  points_[i].value_ = value;
                                  points_[i].range_ = range.find(value, -1);
                              }
                              i++;
                        }
      
            
         
            for (int row = 0; row < rows_; row++) 
                  for (int column = 0; column < columns_; column++) {
                        this->push_back(new Cell<P>(*this, row, column));
                        
                  }
            
}

template <class P>
class CellBox;

template <class P>
class IsoPlot: public IsoPlotAttributes<P>, public  vector<Polyline* >{

public:
      IsoPlot();
      virtual ~IsoPlot();
      
      // Implements the Visdef Interface...
      virtual void operator()(MatrixHandler<P>&, BasicGraphicsObjectContainer&);
      virtual void visit(LegendVisitor&);

      void set(const map<string, string>& map ) {
            IsoPlotAttributes<P>::set(map); 
      }

      void set(const XmlNode& node )  {
            IsoPlotAttributes<P>::set(node); 
      }

      void toxml(ostream&, int)  const {
            //void toxml(ostream& out, int tabs)  const {
            //IsoPlotAttributes<P>::toxml(out, tabs); 
      }

      void setTag(const string&)  {
            //IsoPlotAttributes<P>::setTag(tag); 
      }

      void adjust(double min, double max) 
      { min_ = min; max_ = max; }

      virtual IsoPlot<P>* clone() const { 
            IsoPlot<P>* object = new IsoPlot<P>();
            object->copy(*this);
            return object;
      }

      Colour colour( double value ) {
            return this->shading_->colour(value);
      }

      bool reshape(Shape&,Shape&);
      void reshape(const Colour&, Shape&);
      bool reduce( list<Shape>&, list<Shape>::iterator&);
      void isoline(Cell<P>&, CellBox<P>* = 0);


protected:
      //! Method to print string about this class on to a stream of type ostream (virtual).
      virtual void print(ostream&) const; 

      void isoline(MatrixHandler<P>&, BasicGraphicsObjectContainer&);
      
      bool prepare(MatrixHandler<P>&);
      
      double min_;
      double max_;
      map<double, vector<Polyline* > > helper_;
      double missing_;
      vector<double>  levels_;
      bool shadingMode_; 
      
      

private:
      //! Copy constructor - No copy allowed
      IsoPlot(const IsoPlot&);
      //! Overloaded << operator to copy - No copy allowed
      IsoPlot& operator=(const IsoPlot&);

// -- Friends
      //! Overloaded << operator to call print().
      friend ostream& operator<<(ostream& s,const IsoPlot<P>& p)
            { p.print(s); return s; }
    
};

template <class P>
class NoIsoPlot : public IsoPlot<P>
{
public:
      NoIsoPlot() { this->setTag("noisoline"); };
      ~NoIsoPlot() {};
    
      // Implements the Visualiser Interface...
      void operator()(MatrixHandler<P>&, BasicGraphicsObjectContainer&);
      
      void set(const XmlNode& node)  {
            if ( magCompare(node.name(), "noisoline")  ) {
                  XmlNode iso = node;
                  iso.name("isoline");
                IsoPlotAttributes<P>::set(iso);
            } 
            else 
                  IsoPlotAttributes<P>::set(node);
      }

      virtual IsoPlot<P>* clone() const
      { 
            IsoPlot<P>* object = new NoIsoPlot<P>();
            return object;
      }

protected:
      virtual void print(ostream& out) const 
      { out << "NoIsoPlot" << "\n"; }
};

template<class P>
class Translator<string, IsoPlot<P> > { 
public:
      IsoPlot<P>* operator()(const string& val) {
            return SimpleObjectMaker<IsoPlot<P> >::create(val);
      }     

      IsoPlot<P>* magics(const string& param)
      {
            IsoPlot<P>* object;
            ParameterManager::update(param, object);
            return object;
      }
};

} // namespace magics

#include "IsoPlot.cc"
#endif

Generated by  Doxygen 1.6.0   Back to index