geant4/source/analysis/g4tools/include/tools/ccontour

// 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