ccontour

Go to the documentation of this file.

// Copyright (C) 2010, Guy Barrand. All rights reserved.
// See the file tools.license for terms.
 
#ifndef tools_ccontour
#define tools_ccontour
 
// G.Barrand : inline version of the one found in Lib of OSC 16.11.
//             This code is not mine and I keep it "as it is".
 
// Contour.h: interface for the ccontour class.
//
// ccontour implements Contour plot algorithm descrided in
//      IMPLEMENTATION OF
//      AN IMPROVED CONTOUR
//      PLOTTING ALGORITHM
//      BY
//
//      MICHAEL JOSEPH ARAMINI
//
//      B.S., Stevens Institute of Technology, 1980
// See http://www.ultranet.com/~aramini/thesis.html
//
// Ported to C++ by Jonathan de Halleux.
//
// Using ccontour :
//
// ccontour is not directly usable. The user has to
//  1. derive the function ExportLine that is
//      supposed to draw/store the segment of the contour
//  2. Set the function draw contour of. (using  SetFieldFn
//      The function must be declared as follows
//      double (*myF)(double x , double y);
//
//  History:
//      31-07-2002:
//          - A lot of contribution from Chenggang Zhou (better strip compressions, merging, area, weight),
//          - Got rid of ugly MFC lists for STL.
 
//G.Barrand :
#include <vector>
#include <cstdio>
#include <cstdlib>
#include <cmath>
 
#ifdef TOOLS_MEM
#include "mem"
#endif
 
#include "mnmx"
 
namespace tools {
 
class ccontour
{
#ifdef TOOLS_MEM
    static const std::string& s_class() {
      static const std::string s_v("tools::ccontour");
      return s_v;
    }
#endif
protected:
        // plots a line from (x1,y1) to (x2,y2)
    virtual void ExportLine(int iPlane,int x1,int y1,int x2,int y2) = 0;
 
public:
    ccontour();
    virtual ~ccontour(){
      CleanMemory();
#ifdef TOOLS_MEM
          mem::decrement(s_class().c_str());
#endif
        }
protected: //G.Barrand
        ccontour(const ccontour&){}
private: //G.Barrand
        ccontour& operator=(const ccontour&){return *this;}
public:
protected: //G.Barrand
    // Initialize memory. Called in generate
    virtual void InitMemory();
    // Clean work arrays
    virtual void CleanMemory();
    // generates contour
    // Before calling this functions you must
    //  1. derive the function ExportLine that is
    //      supposed to draw/store the segment of the contour
    //  2. Set the function draw contour of. (using  SetFieldFn
    //      The function must be declared as follows
    //      double (*myF)(double x , double y);
public: //G.Barrand
    virtual void generate();
 
    // Set the dimension of the primary grid
    void set_first_grid(int iCol, int iRow);
    // Set the dimension of the base grid
    void set_secondary_grid(int iCol, int iRow);
    // Sets the region [left, right, bottom,top] to generate contour
    void set_limits(double pLimits[4]);
    // Sets the isocurve values
    void set_planes(const std::vector<double>& vPlanes);
    // Sets the pointer to the F(x,y) funtion
        // G.Barrand : handle a user data pointer.
    void set_field_fcn(double (*_pFieldFcn)(double, double,void*),void*);
 
    size_t get_number_of_planes() const { return m_vPlanes.size();};
    const std::vector<double>& get_planes() const   {   return m_vPlanes;};
    double get_plane(unsigned int i)    const;
 
    // For an indexed point i on the sec. grid, returns x(i)
    double get_xi(int i) const {    return m_pLimits[0] +  i%(m_iColSec+1)*(m_pLimits[1]-m_pLimits[0])/(double)( m_iColSec );};
    // For an indexed point i on the fir. grid, returns y(i)
    double get_yi(int i) const;
 
    void get_limits(double pLimits[4]);
protected: //G.Barrand
 
    // Retrieve dimension of grids, contouring region and isocurve
    int GetColFir() const       {   return m_iColFir;};
    int GetRowFir() const       {   return m_iRowFir;};
    int GetColSec() const       {   return m_iColSec;};
    int GetRowSec() const       {   return m_iRowSec;};
 
private:
  // A structure used internally by ccontour
  /*G.Barrand :
  struct CFnStr {
    double m_dFnVal;
    short m_sLeftLen;
    short m_sRightLen;
    short m_sTopLen;
    short m_sBotLen;
  };
  */
  class CFnStr {
#ifdef TOOLS_MEM
    static const std::string& s_class() {
      static const std::string s_v("tools::ccontour::CFnStr");
      return s_v;
    }
#endif
  public:
    CFnStr():m_dFnVal(0),m_sLeftLen(0),m_sRightLen(0),m_sTopLen(0),m_sBotLen(0){
#ifdef TOOLS_MEM
      mem::increment(s_class().c_str());
#endif
    }
    ~CFnStr(){
#ifdef TOOLS_MEM
      mem::decrement(s_class().c_str());
#endif
    }
  protected:
    CFnStr(const CFnStr&){}
    CFnStr& operator=(const CFnStr&){return *this;}
  public:
    double m_dFnVal;
    short m_sLeftLen;
    short m_sRightLen;
    short m_sTopLen;
    short m_sBotLen;
  };
 
 
protected:
    // Accesibles variables
    std::vector<double> m_vPlanes;          // value of contour planes
    double m_pLimits[4];                        // left, right, bottom, top
    int m_iColFir;                              // primary  grid, number of columns
    int m_iRowFir;                              // primary  grid, number of rows
    int m_iColSec;                              // secondary grid, number of columns
    int m_iRowSec;                              // secondary grid, number of rows
        void* m_pFieldFcnData; // G.Barrand : handle a user data pointer.
        double (*m_pFieldFcn)(double x, double y,void*); // pointer to F(x,y) function
 
    // Protected function
    //virtual void ExportLine(int iPlane, int x1, int y1, int x2, int y2) = 0; // plots a line from (x1,y1) to (x2,y2)
 
    // Work functions and variables
    double m_dDx;
    double m_dDy;
    CFnStr** m_ppFnData;    // pointer to mesh parts
    CFnStr* FnctData(int i,int j)  {    return (m_ppFnData[i]+j);};
 
    double Field(int x, int y);  /* evaluate funct if we must,  */
    void Cntr1(int x1, int x2, int y1, int y2);
    void Pass2(int x1, int x2, int y1, int y2);   /* draws the contour lines */
 
private:
        //G.Barrand : have the below in the class.
        // A simple test function
        static double ContourTestFunction(double x,double y,void*) {
          return 0.5*(::cos(x+3.14/4)+::sin(y+3.14/4));
        }
 
protected:
        static void _ASSERT_(bool a_what,const char* a_where) {
          if(!a_what) {
            ::printf("debug : Contour : assert failure in %s\n",a_where);
            ::exit(0);
          }
        }
 
        static void _ASSERTP_(void* a_what,const char* a_where) {
          if(!a_what) {
            ::printf("debug : Contour : assert failure in %s\n",a_where);
            ::exit(0);
          }
        }
 
};
 
// implementation :
//
//  Code get on the web at :
//    http://www.codeproject.com/cpp/contour.asp
//
//  G.Barrand
//
 
// Construction/Destruction
 
inline ccontour::ccontour()
{
#ifdef TOOLS_MEM
      mem::increment(s_class().c_str());
#endif
    m_iColFir=m_iRowFir=32;
    m_iColSec=m_iRowSec=256;
    m_dDx=m_dDy=0;
        m_pFieldFcnData = NULL;
    m_pFieldFcn=NULL;
    m_pLimits[0]=m_pLimits[2]=0;
    m_pLimits[1]=m_pLimits[3]=5.;
    m_ppFnData=NULL;
 
    // temporary stuff
    m_pFieldFcn=ContourTestFunction;
    m_vPlanes.resize(20);
    for (unsigned int i=0;i<m_vPlanes.size();i++)
    {
        m_vPlanes[i]=(i-m_vPlanes.size()/2.0)*0.1;
    }
}
 
//G.Barrand : inline
 
inline void ccontour::InitMemory()
{
    if (!m_ppFnData)
    {
        m_ppFnData=new CFnStr*[m_iColSec+1];
        for (int i=0;i<m_iColSec+1;i++)
        {
            m_ppFnData[i]=NULL;
        }
    }
}
 
inline void ccontour::CleanMemory()
{
    if (m_ppFnData)
    {
        int i;
        for (i=0;i<m_iColSec+1;i++)
        {
            if (m_ppFnData[i])
                delete[] (m_ppFnData[i]);
        }
        delete[] m_ppFnData;
        m_ppFnData=NULL;
    }
}
 
inline void ccontour::generate()
{
 
    int i, j;
    int x3, x4, y3, y4, x, y, oldx3, xlow;
    const int cols=m_iColSec+1;
    const int rows=m_iRowSec+1;
    //double xoff,yoff;
 
    // Initialize memroy if needed
    InitMemory();
 
    m_dDx = (m_pLimits[1]-m_pLimits[0])/(double)(m_iColSec);
    //xoff = m_pLimits[0];
    m_dDy = (m_pLimits[3]-m_pLimits[2])/(double)(m_iRowSec);
    //yoff = m_pLimits[2];
 
    xlow = 0;
    oldx3 = 0;
    x3 = (cols-1)/m_iRowFir;
    x4 = ( 2*(cols-1) )/m_iRowFir;
    for (x = oldx3; x <= x4; x++)
    {     /* allocate new columns needed
        */
        if (x >= cols)
            break;
        if (m_ppFnData[x]==NULL)
            m_ppFnData[x] = new CFnStr[rows];
 
        for (y = 0; y < rows; y++)
            FnctData(x,y)->m_sTopLen = -1;
    }
 
    y4 = 0;
    for (j = 0; j < m_iColFir; j++)
    {
        y3 = y4;
        y4 = ((j+1)*(rows-1))/m_iColFir;
        Cntr1(oldx3, x3, y3, y4);
    }
 
    for (i = 1; i < m_iRowFir; i++)
    {
        y4 = 0;
        for (j = 0; j < m_iColFir; j++)
        {
            y3 = y4;
            y4 = ((j+1)*(rows-1))/m_iColFir;
            Cntr1(x3, x4, y3, y4);
        }
 
        y4 = 0;
        for (j = 0; j < m_iColFir; j++)
        {
            y3 = y4;
            y4 = ((j+1)*(rows-1))/m_iColFir;
            Pass2(oldx3,x3,y3,y4);
        }
 
        if (i < (m_iRowFir-1))
        {    /* re-use columns no longer needed */
            oldx3 = x3;
            x3 = x4;
            x4 = ((i+2)*(cols-1))/m_iRowFir;
            for (x = x3+1; x <= x4; x++)
            {
                if (xlow < oldx3)
                {
                    if (m_ppFnData[x])
                        delete[] m_ppFnData[x];
                    m_ppFnData[x] = m_ppFnData[xlow];
                    m_ppFnData[ xlow++ ] = NULL;
                }
                else
                    if (m_ppFnData[x]==NULL)
                        m_ppFnData[x] = new CFnStr[rows];
 
                for (y = 0; y < rows; y++)
                    FnctData(x,y)->m_sTopLen = -1;
            }
        }
    }
 
    y4 = 0;
    for (j = 0; j < m_iColFir; j++)
    {
        y3 = y4;
        y4 = ((j+1)*(rows-1))/m_iColFir;
        Pass2(x3,x4,y3,y4);
    }
}
 
inline void ccontour::Cntr1(int x1, int x2, int y1, int y2)
{
    double f11, f12, f21, f22, f33;
    int x3, y3, i, j;
 
    if ((x1 == x2) || (y1 == y2))   /* if not a real cell, punt */
        return;
    f11 = Field(x1, y1);
    f12 = Field(x1, y2);
    f21 = Field(x2, y1);
    f22 = Field(x2, y2);
    if ((x2 > x1+1) || (y2 > y1+1)) {   /* is cell divisible? */
        x3 = (x1+x2)/2;
        y3 = (y1+y2)/2;
        f33 = Field(x3, y3);
        i = j = 0;
        if (f33 < f11) i++; else if (f33 > f11) j++;
        if (f33 < f12) i++; else if (f33 > f12) j++;
        if (f33 < f21) i++; else if (f33 > f21) j++;
        if (f33 < f22) i++; else if (f33 > f22) j++;
        if ((i > 2) || (j > 2)) /* should we divide cell? */
        {
            /* subdivide cell */
            Cntr1(x1, x3, y1, y3);
            Cntr1(x3, x2, y1, y3);
            Cntr1(x1, x3, y3, y2);
            Cntr1(x3, x2, y3, y2);
            return;
        }
    }
    /* install cell in array */
    FnctData(x1,y2)->m_sBotLen = FnctData(x1,y1)->m_sTopLen = x2-x1;
    FnctData(x2,y1)->m_sLeftLen = FnctData(x1,y1)->m_sRightLen = y2-y1;
}
 
inline void ccontour::Pass2(int x1, int x2, int y1, int y2)
{
    int left = 0, right = 0, top = 0, bot = 0,old, iNew, i, j, x3, y3;
    double yy0 = 0, yy1 = 0, xx0 = 0, xx1 = 0, xx3, yy3;
    double v, f11, f12, f21, f22, f33, fold, fnew, f;
    double xoff=m_pLimits[0];
    double yoff=m_pLimits[2];
 
    if ((x1 == x2) || (y1 == y2))   /* if not a real cell, punt */
        return;
    f11 = FnctData(x1,y1)->m_dFnVal;
    f12 = FnctData(x1,y2)->m_dFnVal;
    f21 = FnctData(x2,y1)->m_dFnVal;
    f22 = FnctData(x2,y2)->m_dFnVal;
    if ((x2 > x1+1) || (y2 > y1+1)) /* is cell divisible? */
    {
        x3 = (x1+x2)/2;
        y3 = (y1+y2)/2;
        f33 = FnctData(x3, y3)->m_dFnVal;
        i = j = 0;
        if (f33 < f11) i++; else if (f33 > f11) j++;
        if (f33 < f12) i++; else if (f33 > f12) j++;
        if (f33 < f21) i++; else if (f33 > f21) j++;
        if (f33 < f22) i++; else if (f33 > f22) j++;
        if ((i > 2) || (j > 2)) /* should we divide cell? */
        {
            /* subdivide cell */
            Pass2(x1, x3, y1, y3);
            Pass2(x3, x2, y1, y3);
            Pass2(x1, x3, y3, y2);
            Pass2(x3, x2, y3, y2);
            return;
        }
    }
 
    for (i = 0; i < (int)m_vPlanes.size(); i++)
    {
        v = m_vPlanes[i];
        j = 0;
        if (f21 > v) j++;
        if (f11 > v) j |= 2;
        if (f22 > v) j |= 4;
        if (f12 > v) j |= 010;
        if ((f11 > v) ^ (f12 > v))
        {
            if ((FnctData(x1,y1)->m_sLeftLen != 0) &&
                (FnctData(x1,y1)->m_sLeftLen < FnctData(x1,y1)->m_sRightLen))
            {
                old = y1;
                fold = f11;
                while (1)
                {
                    iNew = old+FnctData(x1,old)->m_sLeftLen;
                    fnew = FnctData(x1,iNew)->m_dFnVal;
                    if ((fnew > v) ^ (fold > v))
                        break;
                    old = iNew;
                    fold = fnew;
                }
                yy0 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
            }
            else
                yy0 = (v-f11)/(f12-f11);
 
            left = (int)(y1+(y2-y1)*yy0+0.5);
        }
        if ((f21 > v) ^ (f22 > v))
        {
            if ((FnctData(x2,y1)->m_sRightLen != 0) &&
                (FnctData(x2,y1)->m_sRightLen < FnctData(x2,y1)->m_sLeftLen))
            {
                old = y1;
                fold = f21;
                while (1)
                {
                    iNew = old+FnctData(x2,old)->m_sRightLen;
                    fnew = FnctData(x2,iNew)->m_dFnVal;
                    if ((fnew > v) ^ (fold > v))
                        break;
                    old = iNew;
                    fold = fnew;
                }
                yy1 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
            }
            else
                yy1 = (v-f21)/(f22-f21);
 
            right = (int)(y1+(y2-y1)*yy1+0.5);
        }
        if ((f21 > v) ^ (f11 > v))
        {
            if ((FnctData(x1,y1)->m_sBotLen != 0) &&
                (FnctData(x1,y1)->m_sBotLen < FnctData(x1,y1)->m_sTopLen)) {
                old = x1;
                fold = f11;
                while (1) {
                    iNew = old+FnctData(old,y1)->m_sBotLen;
                    fnew = FnctData(iNew,y1)->m_dFnVal;
                    if ((fnew > v) ^ (fold > v))
                        break;
                    old = iNew;
                    fold = fnew;
                }
                xx0 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
            }
            else
                xx0 = (v-f11)/(f21-f11);
 
            bot = (int)(x1+(x2-x1)*xx0+0.5);
        }
        if ((f22 > v) ^ (f12 > v))
        {
            if ((FnctData(x1,y2)->m_sTopLen != 0) &&
                (FnctData(x1,y2)->m_sTopLen < FnctData(x1,y2)->m_sBotLen)) {
                old = x1;
                fold = f12;
                while (1) {
                    iNew = old+FnctData(old,y2)->m_sTopLen;
                    fnew = FnctData(iNew,y2)->m_dFnVal;
                    if ((fnew > v) ^ (fold > v))
                        break;
                    old = iNew;
                    fold = fnew;
                }
                xx1 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
            }
            else
                xx1 = (v-f12)/(f22-f12);
 
            top = (int)(x1+(x2-x1)*xx1+0.5);
        }
 
        switch (j)
        {
            case 7:
            case 010:
                ExportLine(i,x1,left,top,y2);
                break;
            case 5:
            case 012:
                ExportLine(i,bot,y1,top,y2);
                break;
            case 2:
            case 015:
                ExportLine(i,x1,left,bot,y1);
            break;
        case 4:
        case 013:
            ExportLine(i,top,y2,x2,right);
            break;
        case 3:
        case 014:
            ExportLine(i,x1,left,x2,right);
            break;
        case 1:
        case 016:
            ExportLine(i,bot,y1,x2,right);
            break;
        case 0:
        case 017:
            break;
        case 6:
        case 011:
            yy3 = (xx0*(yy1-yy0)+yy0)/(1.0-(xx1-xx0)*(yy1-yy0));
            xx3 = yy3*(xx1-xx0)+xx0;
            xx3 = x1+xx3*(x2-x1);
            yy3 = y1+yy3*(y2-y1);
            xx3 = xoff+xx3*m_dDx;
            yy3 = yoff+yy3*m_dDy;
            f = (*m_pFieldFcn)(xx3, yy3,m_pFieldFcnData);
            if (f == v) {
                ExportLine(i,bot,y1,top,y2);
                ExportLine(i,x1,left,x2,right);
            } else
                if (((f > v) && (f22 > v)) || ((f < v) && (f22 < v))) {
                    ExportLine(i,x1,left,top,y2);
                    ExportLine(i,bot,y1,x2,right);
                } else {
                    ExportLine(i,x1,left,bot,y1);
                    ExportLine(i,top,y2,x2,right);
                }
        }
    }
}
 
inline double ccontour::Field(int x, int y)  /* evaluate funct if we must,*/
{
    double x1, y1;
 
    if (FnctData(x,y)->m_sTopLen != -1)  /* is it already in the array */
        return(FnctData(x,y)->m_dFnVal);
 
    /* not in the array, create new array element */
    x1 = m_pLimits[0]+m_dDx*x;
    y1 = m_pLimits[2]+m_dDy*y;
    FnctData(x,y)->m_sTopLen = 0;
    FnctData(x,y)->m_sBotLen = 0;
    FnctData(x,y)->m_sRightLen = 0;
    FnctData(x,y)->m_sLeftLen = 0;
    return (FnctData(x,y)->m_dFnVal = (*m_pFieldFcn)(x1, y1,m_pFieldFcnData));
}
 
inline void ccontour::set_planes(const std::vector<double>& vPlanes)
{
    // cleaning memory
    CleanMemory();
 
    m_vPlanes = vPlanes;
}
 
inline void ccontour::set_field_fcn(double (*_pFieldFcn)(double, double,void*),void* aData)
{
        m_pFieldFcnData = aData;
    m_pFieldFcn=_pFieldFcn;
}
 
inline void ccontour::set_first_grid(int iCol, int iRow)
{
    m_iColFir=mx<int>(iCol,2);
    m_iRowFir=mx<int>(iRow,2);
}
 
inline void ccontour::set_secondary_grid(int iCol, int iRow)
{
    // cleaning work matrices if allocated
    CleanMemory();
 
    m_iColSec=mx<int>(iCol,2);
    m_iRowSec=mx<int>(iRow,2);
}
 
inline void ccontour::set_limits(double pLimits[])
{
    _ASSERT_(pLimits[0]<pLimits[1],"ccontour::set_limits");
    _ASSERT_(pLimits[2]<pLimits[3],"ccontour::set_limits");
    for (int i=0;i<4;i++)
    {
        m_pLimits[i]=pLimits[i];
    }
}
 
inline void ccontour::get_limits(double pLimits[])
{
    for (int i=0;i<4;i++)
    {
        pLimits[i]=m_pLimits[i];
    }
}
 
//G.Barrand : from .h to .cxx to avoid _ASSERT_ in .h
inline double ccontour::get_plane(unsigned int i) const {
  /*_ASSERT_(i>=0);*/
  _ASSERT_(i<m_vPlanes.size(),"ccontour::get_plane");
  return m_vPlanes[i];
}
 
inline double ccontour::get_yi(int i) const {
  if(i<0) ::printf("ccontour::get_yi : %d\n",i);
  _ASSERT_(i>=0,"ccontour::get_yi");
  return m_pLimits[2] +  i/(m_iColSec+1)*(m_pLimits[3]-m_pLimits[2])/(double)( m_iRowSec );
}
 
}
 
#endif