/pgl/algotools/Contour.h
C Header | 606 lines | 438 code | 57 blank | 111 comment | 96 complexity | fcd06c81f57f8c0510355749eac93f72 MD5 | raw file
- /* ****************************************************************************
- *
- * Copyright (c) Microsoft Corporation.
- *
- * This source code is subject to terms and conditions of the Microsoft Public License. A
- * copy of the license can be found in the License.html file at the root of this distribution. If
- * you cannot locate the Microsoft Public License, please send an email to
- * dlr@microsoft.com. By using this source code in any fashion, you are agreeing to be bound
- * by the terms of the Microsoft Public License.
- *
- * You must not remove this notice, or any other, from this software.
- *
- * THIS CODE AND INFORMATION ARE PROVIDED "AS IS" WITHOUT WARRANTY OF ANY
- * KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR FITNESS FOR A
- * PARTICULAR PURPOSE.
- * ***************************************************************************/
-
-
- #if !defined(AFX_CONTOUR_H__45A84754_43CD_414A_80A0_9211F954B0DA__INCLUDED_)
- #define AFX_CONTOUR_H__45A84754_43CD_414A_80A0_9211F954B0DA__INCLUDED_
-
- #if _MSC_VER > 1000
- #pragma once
- #endif // _MSC_VER > 1000
-
- #include "Field.h"
-
-
- //! Contouring algorithm
- namespace contour
- {
-
- /*! \brief A contouring class
-
- TContour 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 TContour :
-
- TContour 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. Attach a field object to it
-
- \ingroup ContourLibGroup
-
- */
- template<class T>
- class TContour
- {
- public:
- typedef std::vector<T> PlaneContainer;
- struct SFnStr
- {
- T m_dFnVal;
- short m_sLeftLen;
- short m_sRightLen;
- short m_sTopLen;
- short m_sBotLen;
- };
-
- //! \name Constructors
- //@{
- //! Default constructor
- TContour();
- virtual ~TContour();
- //@}
-
- //! \name Main Methods
- //@{
- //! Initialize memory. Called in Generate
- virtual void InitMemory();
- //! Clean work arrays
- virtual void CleanMemory();
-
- /*! \brief 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();
- //@}
-
- //! \name Grids
- //@{
- //! Sets the pointer to the F(x,y) funtion
- void SetField(TField<T>* _pField);
- //! Returns pointer to the current field function
- TField<T>* GetField() { return m_pField;};
- //@}
-
- //! \name Grid Methods
- //@{
- //! Set the dimension of the primary grid
- void SetFirstGrid(int iRow, int iCol);
- //! Set the dimension of the base grid
- void SetSecondaryGrid(int iRow, int iCol);
- /*! \brief Sets the number of isocurve and their values
- pPlanes must be allocated on the heap. It will be deleted by TContour.
- */
- //! Retrieve number of columns of primary grid
- int GetColFir() const { return m_iColFir;};
- //! Retrieve number of rows of primary grid
- int GetRowFir() const { return m_iRowFir;};
- //! Retrieve number of columns of secondary grid
- int GetColSec() const { return m_iColSec;};
- //! Retrieve number of rows of secondary grid
- int GetRowSec() const { return m_iRowSec;};
- //@}
-
- //! \name Plane methods
- //@{
- void SetPlanes(const PlaneContainer vPlanes){ CleanMemory();m_vPlanes=vPlanes;};
- //! Return current number of planes
- size_t GetPlanesSize() const { return m_vPlanes.size();};
- //! Returns planes
- const PlaneContainer& GetPlanes() const { return m_vPlanes;};
- //! Returns i-th plane
- double GetPlane(UINT i) const { ASSERT(i>=0); ASSERT(i<m_vPlanes.size()); return m_vPlanes[i];};
- //@}
-
- //! \name Approximated point coordinates
- //@{
- //! For an indexed point i on the sec. grid, returns x(i)
- T GetXi(int i) const
- {
- const TField<T>::SLimits& limits = m_pField->GetLimits();
- return (T)(limits.left+(i%(m_iColSec+1))*(limits.GetWidth())/(T)( m_iColSec ));
- };
- //! For an indexed point i on the fir. grid, returns y(i)
- T GetYi(int i) const
- {
- const TField<T>::SLimits& limits = m_pField->GetLimits();
- UINT index=i/(m_iColSec+1);
- return (T)(limits.bottom+index*limits.GetHeight()/(T)(m_iRowSec));
- };
- //@}
-
- protected:
- //! \name Attributes
- //@{
- //! value of contour planes
- PlaneContainer m_vPlanes;
- //! primary grid, number of columns
- int m_iColFir;
- //! primary grid, number of rows
- int m_iRowFir;
- //! secondary grid, number of columns
- int m_iColSec;
- //! secondary grid, number of rows
- int m_iRowSec;
- T m_dDx;
- T m_dDy;
- //! pointer to F(x,y) function
- TField<T>* m_pField;
- //@}
-
- //! Protected function
- virtual void ExportLine(UINT iPlane,int x1, int y1, int x2, int y2) = 0; // plots a line from (x1,y1) to (x2,y2)
-
- // Work functions and variables
- //! pointer to mesh parts
- SFnStr** m_ppFnData;
- SFnStr* FnctData(int i,int j) { return (m_ppFnData[i]+j);};
- T 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 */
- };
-
- //////////////////////////////////////////////////////////////////////
- // Construction/Destruction
- //////////////////////////////////////////////////////////////////////
-
- template <class T>
- TContour<T>::TContour()
- : m_vPlanes(1,0)
- {
- m_iColFir=m_iRowFir=32;
- m_iColSec=m_iRowSec=256;
- m_dDx=m_dDy=0;
- m_pField=NULL;
- m_ppFnData=NULL;
- }
-
- template <class T>
- TContour<T>::~TContour()
- {
- CleanMemory();
- }
-
- template <class T>
- void TContour<T>::InitMemory()
- {
- if (!m_ppFnData)
- {
- m_ppFnData=new SFnStr*[m_iColSec+1];
- for (int i=0;i<m_iColSec+1;i++)
- {
- m_ppFnData[i]=NULL;
- }
- }
- }
-
- template <class T>
- void TContour<T>::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;
- }
- }
-
- template <class T>
- void TContour<T>::Generate()
- {
- if (!m_pField)
- return;
-
- int i, j;
- int x3, x4, y3, y4, x, y, oldx3, xlow;
- const int cols=m_iColSec+1;
- const int rows=m_iRowSec+1;
- T xoff,yoff;
-
- // Initialize memroy if needed
- InitMemory();
-
- const TField<T>::SLimits& limits = m_pField->GetLimits();
-
- m_dDx = limits.GetWidth()/(T)(m_iColSec);
- xoff = limits.left;
- m_dDy = limits.GetHeight()/(T)(m_iRowSec);
- yoff = limits.bottom;
-
- 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 SFnStr[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 SFnStr[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);
- }
- }
-
- template <class T>
- void TContour<T>::Cntr1(int x1, int x2, int y1, int y2)
- {
- T 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;
- }
-
- template <class T>
- void TContour<T>::Pass2(int x1, int x2, int y1, int y2)
- {
- int left, right, top, bot,old, iNew, x3, y3;
- UINT i,j;
- T yy0, yy1, xx0, xx1, xx3, yy3;
- T v, f11, f12, f21, f22, f33, fold, fnew, f;
- _ASSERT(m_pField);
- const TField<T>::SLimits& limits = m_pField->GetLimits();
-
- T xoff=limits.left;
- T yoff=limits.bottom;
-
- 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 < 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 = (T)((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;
-
- _ASSERT(m_pField);
- f = m_pField->Eval(xx3, yy3);
- 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);
- }
- }
- }
- }
-
- template <class T>
- T TContour<T>::Field(int x, int y) /* evaluate funct if we must, */
- {
- T x1, y1;
- _ASSERT(m_pField);
- const TField<T>::SLimits& limits = m_pField->GetLimits();
-
- 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 = limits.left+m_dDx*x;
- y1 = limits.bottom+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_pField->Eval(x1, y1));
- }
-
- template <class T>
- void TContour<T>::SetFirstGrid(int iRow, int iCol)
- {
- m_iColFir=__max(iCol,2);
- m_iRowFir=__max(iRow,2);
- }
-
- template <class T>
- void TContour<T>::SetSecondaryGrid(int iRow, int iCol)
- {
- // cleaning work matrices if allocated
- CleanMemory();
-
- m_iColSec=__max(iCol,2);
- m_iRowSec=__max(iRow,2);
- }
-
- template <class T>
- void TContour<T>::SetField(TField<T>* _pField)
- {
- m_pField=_pField;
- };
-
- };
-
- #endif // !defined(AFX_CONTOUR_H__45A84754_43CD_414A_80A0_9211F954B0DA__INCLUDED_)