/Users/barrand/private/dev/softinex/old/inexlib-1.2/inlib/inlib/CContour

Go to the documentation of this file.

// Copyright (C) 2010, Guy Barrand. All rights reserved.
// See the file inlib.license for terms.

#ifndef inlib_CContour
#define inlib_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>

#include "mnmx"

namespace inlib {

class CContour  
{
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();}
protected: //G.Barrand
        inline CContour(const CContour&){}
private: //G.Barrand
        inline CContour& operator=(const CContour&){return *this;}
public:
        // 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);
        virtual void Generate();

        // Set the dimension of the primary grid
        void SetFirstGrid(int iCol, int iRow);
        // Set the dimension of the base grid
        void SetSecondaryGrid(int iCol, int iRow);
        // Sets the region [left, right, bottom,top] to generate contour 
        void SetLimits(double pLimits[4]);
        // Sets the isocurve values
        void SetPlanes(const std::vector<double>& vPlanes);
        // Sets the pointer to the F(x,y) funtion
        // G.Barrand : handle a user data pointer.
        void SetFieldFcn(double (*_pFieldFcn)(double, double,void*),void*);

        // Retrieve dimension of grids, contouring region and isocurve
        inline int GetColFir() const            {       return m_iColFir;};
        inline int GetRowFir() const            {       return m_iRowFir;};
        inline int GetColSec() const            {       return m_iColSec;};
        inline int GetRowSec() const            {       return m_iRowSec;};
        void GetLimits(double pLimits[4]);
        inline unsigned int GetNPlanes() const { return m_vPlanes.size();};
        inline const std::vector<double>& GetPlanes() const     {       return m_vPlanes;};
        double GetPlane(unsigned int i) const;

        // For an indexed point i on the sec. grid, returns x(i)
        inline double GetXi(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 GetYi(int i) const;

private:
  // A structure used internally by CContour
  struct CFnStr {
    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
        inline 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
        inline static double ContourTestFunction(double x,double y,void*) {  
          return 0.5*(::cos(x+3.14/4)+::sin(y+3.14/4)); 
        }

protected:
        inline 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);
          }
        }

        inline 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()
{
        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::SetPlanes(const std::vector<double>& vPlanes)
{       
        // cleaning memory
        CleanMemory();

        m_vPlanes = vPlanes;
}

inline void CContour::SetFieldFcn(double (*_pFieldFcn)(double, double,void*),void* aData) 
{       
        m_pFieldFcnData = aData;
        m_pFieldFcn=_pFieldFcn;
}

inline void CContour::SetFirstGrid(int iCol, int iRow)
{
        m_iColFir=mx<int>(iCol,2);
        m_iRowFir=mx<int>(iRow,2);
}

inline void CContour::SetSecondaryGrid(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::SetLimits(double pLimits[])
{
        ASSERT(pLimits[0]<pLimits[1],"CContour::SetLimits");
        ASSERT(pLimits[2]<pLimits[3],"CContour::SetLimits");
        for (int i=0;i<4;i++)
        {
                m_pLimits[i]=pLimits[i];
        }       
}

inline void CContour::GetLimits(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::GetPlane(unsigned int i) const  {
  /*ASSERT(i>=0);*/
  ASSERT(i<m_vPlanes.size(),"CContour::GetPlane");
  return m_vPlanes[i];
}

inline double CContour::GetYi(int i) const {
  if(i<0) ::printf("CContour::GetYi : %d\n",i);
  ASSERT(i>=0,"CContour::GetYi");
  return m_pLimits[2] +  i/(m_iColSec+1)*(m_pLimits[3]-m_pLimits[2])/(double)( m_iRowSec );
}

}

#endif