Warning, /include/Geant4/tools/ccontour is written in an unsupported language. File is not indexed.
0001 // Copyright (C) 2010, Guy Barrand. All rights reserved.
0002 // See the file tools.license for terms.
0003
0004 #ifndef tools_ccontour
0005 #define tools_ccontour
0006
0007 // G.Barrand : inline version of the one found in Lib of OSC 16.11.
0008 // This code is not mine and I keep it "as it is".
0009
0010 // Contour.h: interface for the ccontour class.
0011 //
0012 // ccontour implements Contour plot algorithm descrided in
0013 // IMPLEMENTATION OF
0014 // AN IMPROVED CONTOUR
0015 // PLOTTING ALGORITHM
0016 // BY
0017 //
0018 // MICHAEL JOSEPH ARAMINI
0019 //
0020 // B.S., Stevens Institute of Technology, 1980
0021 // See http://www.ultranet.com/~aramini/thesis.html
0022 //
0023 // Ported to C++ by Jonathan de Halleux.
0024 //
0025 // Using ccontour :
0026 //
0027 // ccontour is not directly usable. The user has to
0028 // 1. derive the function ExportLine that is
0029 // supposed to draw/store the segment of the contour
0030 // 2. Set the function draw contour of. (using SetFieldFn
0031 // The function must be declared as follows
0032 // double (*myF)(double x , double y);
0033 //
0034 // History:
0035 // 31-07-2002:
0036 // - A lot of contribution from Chenggang Zhou (better strip compressions, merging, area, weight),
0037 // - Got rid of ugly MFC lists for STL.
0038 //////////////////////////////////////////////////////////////////////
0039
0040 //G.Barrand :
0041 #include <vector>
0042 #include <cstdio>
0043 #include <cstdlib>
0044 #include <cmath>
0045
0046 #include "mnmx"
0047
0048 namespace tools {
0049
0050 class ccontour
0051 {
0052 protected:
0053 // plots a line from (x1,y1) to (x2,y2)
0054 virtual void ExportLine(int iPlane,int x1,int y1,int x2,int y2) = 0;
0055
0056 public:
0057 ccontour();
0058 virtual ~ccontour(){
0059 CleanMemory();
0060 }
0061 protected: //G.Barrand
0062 ccontour(const ccontour&){}
0063 private: //G.Barrand
0064 ccontour& operator=(const ccontour&){return *this;}
0065 public:
0066 protected: //G.Barrand
0067 // Initialize memory. Called in generate
0068 virtual void InitMemory();
0069 // Clean work arrays
0070 virtual void CleanMemory();
0071 // generates contour
0072 // Before calling this functions you must
0073 // 1. derive the function ExportLine that is
0074 // supposed to draw/store the segment of the contour
0075 // 2. Set the function draw contour of. (using SetFieldFn
0076 // The function must be declared as follows
0077 // double (*myF)(double x , double y);
0078 public: //G.Barrand
0079 virtual void generate();
0080
0081 // Set the dimension of the primary grid
0082 void set_first_grid(int iCol, int iRow);
0083 // Set the dimension of the base grid
0084 void set_secondary_grid(int iCol, int iRow);
0085 // Sets the region [left, right, bottom,top] to generate contour
0086 void set_limits(double pLimits[4]);
0087 // Sets the isocurve values
0088 void set_planes(const std::vector<double>& vPlanes);
0089 // Sets the pointer to the F(x,y) funtion
0090 // G.Barrand : handle a user data pointer.
0091 void set_field_fcn(double (*_pFieldFcn)(double, double,void*),void*);
0092
0093 size_t get_number_of_planes() const { return m_vPlanes.size();};
0094 const std::vector<double>& get_planes() const { return m_vPlanes;};
0095 double get_plane(unsigned int i) const;
0096
0097 // For an indexed point i on the sec. grid, returns x(i)
0098 double get_xi(int i) const { return m_pLimits[0] + i%(m_iColSec+1)*(m_pLimits[1]-m_pLimits[0])/(double)( m_iColSec );};
0099 // For an indexed point i on the fir. grid, returns y(i)
0100 double get_yi(int i) const;
0101
0102 void get_limits(double pLimits[4]);
0103 protected: //G.Barrand
0104
0105 // Retrieve dimension of grids, contouring region and isocurve
0106 int GetColFir() const { return m_iColFir;};
0107 int GetRowFir() const { return m_iRowFir;};
0108 int GetColSec() const { return m_iColSec;};
0109 int GetRowSec() const { return m_iRowSec;};
0110
0111 private:
0112 class CFnStr {
0113 public:
0114 CFnStr():m_dFnVal(0),m_sLeftLen(0),m_sRightLen(0),m_sTopLen(0),m_sBotLen(0){
0115 }
0116 ~CFnStr(){
0117 }
0118 protected:
0119 CFnStr(const CFnStr&){}
0120 CFnStr& operator=(const CFnStr&){return *this;}
0121 public:
0122 double m_dFnVal;
0123 short m_sLeftLen;
0124 short m_sRightLen;
0125 short m_sTopLen;
0126 short m_sBotLen;
0127 };
0128
0129
0130 protected:
0131 // Accesibles variables
0132 std::vector<double> m_vPlanes; // value of contour planes
0133 double m_pLimits[4]; // left, right, bottom, top
0134 int m_iColFir; // primary grid, number of columns
0135 int m_iRowFir; // primary grid, number of rows
0136 unsigned int m_iColSec; // secondary grid, number of columns
0137 unsigned int m_iRowSec; // secondary grid, number of rows
0138 void* m_pFieldFcnData; // G.Barrand : handle a user data pointer.
0139 double (*m_pFieldFcn)(double x, double y,void*); // pointer to F(x,y) function
0140
0141 // Work functions and variables
0142 double m_dDx;
0143 double m_dDy;
0144 CFnStr** m_ppFnData; // pointer to mesh parts
0145 CFnStr* FnctData(int i,int j) { return (m_ppFnData[i]+j);};
0146
0147 double Field(int x, int y); /* evaluate funct if we must, */
0148 void Cntr1(int x1, int x2, int y1, int y2);
0149 void Pass2(int x1, int x2, int y1, int y2); /* draws the contour lines */
0150
0151 private:
0152 //G.Barrand : have the below in the class.
0153 // A simple test function
0154 static double ContourTestFunction(double x,double y,void*) {
0155 return 0.5*(::cos(x+3.14/4)+::sin(y+3.14/4));
0156 }
0157
0158 protected:
0159 static void _ASSERT_(bool a_what,const char* a_where) {
0160 if(!a_what) {
0161 ::printf("debug : Contour : assert failure in %s\n",a_where);
0162 ::exit(0);
0163 }
0164 }
0165
0166 static void _ASSERTP_(void* a_what,const char* a_where) {
0167 if(!a_what) {
0168 ::printf("debug : Contour : assert failure in %s\n",a_where);
0169 ::exit(0);
0170 }
0171 }
0172
0173 };
0174
0175 // implementation :
0176 //
0177 // Code get on the web at :
0178 // http://www.codeproject.com/cpp/contour.asp
0179 //
0180 // G.Barrand
0181 //
0182
0183 //////////////////////////////////////////////////////////////////////
0184 // Construction/Destruction
0185 //////////////////////////////////////////////////////////////////////
0186
0187 inline ccontour::ccontour()
0188 {
0189 m_iColFir=m_iRowFir=32;
0190 m_iColSec=m_iRowSec=256;
0191 m_dDx=m_dDy=0;
0192 m_pFieldFcnData = NULL;
0193 m_pFieldFcn=NULL;
0194 m_pLimits[0]=m_pLimits[2]=0;
0195 m_pLimits[1]=m_pLimits[3]=5.;
0196 m_ppFnData=NULL;
0197
0198 // temporary stuff
0199 m_pFieldFcn=ContourTestFunction;
0200 m_vPlanes.resize(20);
0201 for (unsigned int i=0;i<m_vPlanes.size();i++)
0202 {
0203 m_vPlanes[i]=(i-m_vPlanes.size()/2.0)*0.1;
0204 }
0205 }
0206
0207 //G.Barrand : inline
0208
0209 inline void ccontour::InitMemory()
0210 {
0211 if (!m_ppFnData)
0212 {
0213 m_ppFnData=new CFnStr*[m_iColSec+1];
0214 for (unsigned int i=0;i<m_iColSec+1;i++)
0215 {
0216 m_ppFnData[i]=NULL;
0217 }
0218 }
0219 }
0220
0221 inline void ccontour::CleanMemory()
0222 {
0223 if (m_ppFnData)
0224 {
0225 for (unsigned int i=0;i<m_iColSec+1;i++)
0226 {
0227 if (m_ppFnData[i])
0228 delete[] (m_ppFnData[i]);
0229 }
0230 delete[] m_ppFnData;
0231 m_ppFnData=NULL;
0232 }
0233 }
0234
0235 inline void ccontour::generate()
0236 {
0237
0238 int i, j;
0239 int x3, x4, y3, y4, x, y, oldx3, xlow;
0240 const unsigned int cols=m_iColSec+1;
0241 const unsigned int rows=m_iRowSec+1;
0242 //double xoff,yoff;
0243
0244 // Initialize memroy if needed
0245 InitMemory();
0246
0247 m_dDx = (m_pLimits[1]-m_pLimits[0])/(double)(m_iColSec);
0248 //xoff = m_pLimits[0];
0249 m_dDy = (m_pLimits[3]-m_pLimits[2])/(double)(m_iRowSec);
0250 //yoff = m_pLimits[2];
0251
0252 xlow = 0;
0253 oldx3 = 0;
0254 x3 = (cols-1)/m_iRowFir;
0255 x4 = ( 2*(cols-1) )/m_iRowFir;
0256 for (x = oldx3; x <= x4; x++)
0257 { /* allocate new columns needed
0258 */
0259 if (x >= (int)cols)
0260 break;
0261 if (m_ppFnData[x]==NULL)
0262 m_ppFnData[x] = new CFnStr[rows];
0263
0264 for (y = 0; y < (int)rows; y++)
0265 FnctData(x,y)->m_sTopLen = -1;
0266 }
0267
0268 y4 = 0;
0269 for (j = 0; j < m_iColFir; j++)
0270 {
0271 y3 = y4;
0272 y4 = ((j+1)*(rows-1))/m_iColFir;
0273 Cntr1(oldx3, x3, y3, y4);
0274 }
0275
0276 for (i = 1; i < m_iRowFir; i++)
0277 {
0278 y4 = 0;
0279 for (j = 0; j < m_iColFir; j++)
0280 {
0281 y3 = y4;
0282 y4 = ((j+1)*(rows-1))/m_iColFir;
0283 Cntr1(x3, x4, y3, y4);
0284 }
0285
0286 y4 = 0;
0287 for (j = 0; j < m_iColFir; j++)
0288 {
0289 y3 = y4;
0290 y4 = ((j+1)*(rows-1))/m_iColFir;
0291 Pass2(oldx3,x3,y3,y4);
0292 }
0293
0294 if (i < (m_iRowFir-1))
0295 { /* re-use columns no longer needed */
0296 oldx3 = x3;
0297 x3 = x4;
0298 x4 = ((i+2)*(cols-1))/m_iRowFir;
0299 for (x = x3+1; x <= x4; x++)
0300 {
0301 if (xlow < oldx3)
0302 {
0303 if (m_ppFnData[x])
0304 delete[] m_ppFnData[x];
0305 m_ppFnData[x] = m_ppFnData[xlow];
0306 m_ppFnData[ xlow++ ] = NULL;
0307 }
0308 else
0309 if (m_ppFnData[x]==NULL)
0310 m_ppFnData[x] = new CFnStr[rows];
0311
0312 for (y = 0; y < (int)rows; y++)
0313 FnctData(x,y)->m_sTopLen = -1;
0314 }
0315 }
0316 }
0317
0318 y4 = 0;
0319 for (j = 0; j < m_iColFir; j++)
0320 {
0321 y3 = y4;
0322 y4 = ((j+1)*(rows-1))/m_iColFir;
0323 Pass2(x3,x4,y3,y4);
0324 }
0325 }
0326
0327 inline void ccontour::Cntr1(int x1, int x2, int y1, int y2)
0328 {
0329 double f11, f12, f21, f22, f33;
0330 int x3, y3, i, j;
0331
0332 if ((x1 == x2) || (y1 == y2)) /* if not a real cell, punt */
0333 return;
0334 f11 = Field(x1, y1);
0335 f12 = Field(x1, y2);
0336 f21 = Field(x2, y1);
0337 f22 = Field(x2, y2);
0338 if ((x2 > x1+1) || (y2 > y1+1)) { /* is cell divisible? */
0339 x3 = (x1+x2)/2;
0340 y3 = (y1+y2)/2;
0341 f33 = Field(x3, y3);
0342 i = j = 0;
0343 if (f33 < f11) i++; else if (f33 > f11) j++;
0344 if (f33 < f12) i++; else if (f33 > f12) j++;
0345 if (f33 < f21) i++; else if (f33 > f21) j++;
0346 if (f33 < f22) i++; else if (f33 > f22) j++;
0347 if ((i > 2) || (j > 2)) /* should we divide cell? */
0348 {
0349 /* subdivide cell */
0350 Cntr1(x1, x3, y1, y3);
0351 Cntr1(x3, x2, y1, y3);
0352 Cntr1(x1, x3, y3, y2);
0353 Cntr1(x3, x2, y3, y2);
0354 return;
0355 }
0356 }
0357 /* install cell in array */
0358 FnctData(x1,y2)->m_sBotLen = FnctData(x1,y1)->m_sTopLen = x2-x1;
0359 FnctData(x2,y1)->m_sLeftLen = FnctData(x1,y1)->m_sRightLen = y2-y1;
0360 }
0361
0362 inline void ccontour::Pass2(int x1, int x2, int y1, int y2)
0363 {
0364 int left = 0, right = 0, top = 0, bot = 0,old, iNew, i, j, x3, y3;
0365 double yy0 = 0, yy1 = 0, xx0 = 0, xx1 = 0, xx3, yy3;
0366 double v, f11, f12, f21, f22, f33, fold, fnew, f;
0367 double xoff=m_pLimits[0];
0368 double yoff=m_pLimits[2];
0369
0370 if ((x1 == x2) || (y1 == y2)) /* if not a real cell, punt */
0371 return;
0372 f11 = FnctData(x1,y1)->m_dFnVal;
0373 f12 = FnctData(x1,y2)->m_dFnVal;
0374 f21 = FnctData(x2,y1)->m_dFnVal;
0375 f22 = FnctData(x2,y2)->m_dFnVal;
0376 if ((x2 > x1+1) || (y2 > y1+1)) /* is cell divisible? */
0377 {
0378 x3 = (x1+x2)/2;
0379 y3 = (y1+y2)/2;
0380 f33 = FnctData(x3, y3)->m_dFnVal;
0381 i = j = 0;
0382 if (f33 < f11) i++; else if (f33 > f11) j++;
0383 if (f33 < f12) i++; else if (f33 > f12) j++;
0384 if (f33 < f21) i++; else if (f33 > f21) j++;
0385 if (f33 < f22) i++; else if (f33 > f22) j++;
0386 if ((i > 2) || (j > 2)) /* should we divide cell? */
0387 {
0388 /* subdivide cell */
0389 Pass2(x1, x3, y1, y3);
0390 Pass2(x3, x2, y1, y3);
0391 Pass2(x1, x3, y3, y2);
0392 Pass2(x3, x2, y3, y2);
0393 return;
0394 }
0395 }
0396
0397 for (i = 0; i < (int)m_vPlanes.size(); i++)
0398 {
0399 v = m_vPlanes[i];
0400 j = 0;
0401 if (f21 > v) j++;
0402 if (f11 > v) j |= 2;
0403 if (f22 > v) j |= 4;
0404 if (f12 > v) j |= 010;
0405 if ((f11 > v) ^ (f12 > v))
0406 {
0407 if ((FnctData(x1,y1)->m_sLeftLen != 0) &&
0408 (FnctData(x1,y1)->m_sLeftLen < FnctData(x1,y1)->m_sRightLen))
0409 {
0410 old = y1;
0411 fold = f11;
0412 while (1)
0413 {
0414 iNew = old+FnctData(x1,old)->m_sLeftLen;
0415 fnew = FnctData(x1,iNew)->m_dFnVal;
0416 if ((fnew > v) ^ (fold > v))
0417 break;
0418 old = iNew;
0419 fold = fnew;
0420 }
0421 yy0 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
0422 }
0423 else
0424 yy0 = (v-f11)/(f12-f11);
0425
0426 left = (int)(y1+(y2-y1)*yy0+0.5);
0427 }
0428 if ((f21 > v) ^ (f22 > v))
0429 {
0430 if ((FnctData(x2,y1)->m_sRightLen != 0) &&
0431 (FnctData(x2,y1)->m_sRightLen < FnctData(x2,y1)->m_sLeftLen))
0432 {
0433 old = y1;
0434 fold = f21;
0435 while (1)
0436 {
0437 iNew = old+FnctData(x2,old)->m_sRightLen;
0438 fnew = FnctData(x2,iNew)->m_dFnVal;
0439 if ((fnew > v) ^ (fold > v))
0440 break;
0441 old = iNew;
0442 fold = fnew;
0443 }
0444 yy1 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
0445 }
0446 else
0447 yy1 = (v-f21)/(f22-f21);
0448
0449 right = (int)(y1+(y2-y1)*yy1+0.5);
0450 }
0451 if ((f21 > v) ^ (f11 > v))
0452 {
0453 if ((FnctData(x1,y1)->m_sBotLen != 0) &&
0454 (FnctData(x1,y1)->m_sBotLen < FnctData(x1,y1)->m_sTopLen)) {
0455 old = x1;
0456 fold = f11;
0457 while (1) {
0458 iNew = old+FnctData(old,y1)->m_sBotLen;
0459 fnew = FnctData(iNew,y1)->m_dFnVal;
0460 if ((fnew > v) ^ (fold > v))
0461 break;
0462 old = iNew;
0463 fold = fnew;
0464 }
0465 xx0 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
0466 }
0467 else
0468 xx0 = (v-f11)/(f21-f11);
0469
0470 bot = (int)(x1+(x2-x1)*xx0+0.5);
0471 }
0472 if ((f22 > v) ^ (f12 > v))
0473 {
0474 if ((FnctData(x1,y2)->m_sTopLen != 0) &&
0475 (FnctData(x1,y2)->m_sTopLen < FnctData(x1,y2)->m_sBotLen)) {
0476 old = x1;
0477 fold = f12;
0478 while (1) {
0479 iNew = old+FnctData(old,y2)->m_sTopLen;
0480 fnew = FnctData(iNew,y2)->m_dFnVal;
0481 if ((fnew > v) ^ (fold > v))
0482 break;
0483 old = iNew;
0484 fold = fnew;
0485 }
0486 xx1 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
0487 }
0488 else
0489 xx1 = (v-f12)/(f22-f12);
0490
0491 top = (int)(x1+(x2-x1)*xx1+0.5);
0492 }
0493
0494 switch (j)
0495 {
0496 case 7:
0497 case 010:
0498 ExportLine(i,x1,left,top,y2);
0499 break;
0500 case 5:
0501 case 012:
0502 ExportLine(i,bot,y1,top,y2);
0503 break;
0504 case 2:
0505 case 015:
0506 ExportLine(i,x1,left,bot,y1);
0507 break;
0508 case 4:
0509 case 013:
0510 ExportLine(i,top,y2,x2,right);
0511 break;
0512 case 3:
0513 case 014:
0514 ExportLine(i,x1,left,x2,right);
0515 break;
0516 case 1:
0517 case 016:
0518 ExportLine(i,bot,y1,x2,right);
0519 break;
0520 case 0:
0521 case 017:
0522 break;
0523 case 6:
0524 case 011:
0525 yy3 = (xx0*(yy1-yy0)+yy0)/(1.0-(xx1-xx0)*(yy1-yy0));
0526 xx3 = yy3*(xx1-xx0)+xx0;
0527 xx3 = x1+xx3*(x2-x1);
0528 yy3 = y1+yy3*(y2-y1);
0529 xx3 = xoff+xx3*m_dDx;
0530 yy3 = yoff+yy3*m_dDy;
0531 f = (*m_pFieldFcn)(xx3, yy3,m_pFieldFcnData);
0532 if (f == v) {
0533 ExportLine(i,bot,y1,top,y2);
0534 ExportLine(i,x1,left,x2,right);
0535 } else
0536 if (((f > v) && (f22 > v)) || ((f < v) && (f22 < v))) {
0537 ExportLine(i,x1,left,top,y2);
0538 ExportLine(i,bot,y1,x2,right);
0539 } else {
0540 ExportLine(i,x1,left,bot,y1);
0541 ExportLine(i,top,y2,x2,right);
0542 }
0543 }
0544 }
0545 }
0546
0547 inline double ccontour::Field(int x, int y) /* evaluate funct if we must,*/
0548 {
0549 double x1, y1;
0550
0551 if (FnctData(x,y)->m_sTopLen != -1) /* is it already in the array */
0552 return(FnctData(x,y)->m_dFnVal);
0553
0554 /* not in the array, create new array element */
0555 x1 = m_pLimits[0]+m_dDx*x;
0556 y1 = m_pLimits[2]+m_dDy*y;
0557 FnctData(x,y)->m_sTopLen = 0;
0558 FnctData(x,y)->m_sBotLen = 0;
0559 FnctData(x,y)->m_sRightLen = 0;
0560 FnctData(x,y)->m_sLeftLen = 0;
0561 return (FnctData(x,y)->m_dFnVal = (*m_pFieldFcn)(x1, y1,m_pFieldFcnData));
0562 }
0563
0564 inline void ccontour::set_planes(const std::vector<double>& vPlanes)
0565 {
0566 // cleaning memory
0567 CleanMemory();
0568
0569 m_vPlanes = vPlanes;
0570 }
0571
0572 inline void ccontour::set_field_fcn(double (*_pFieldFcn)(double, double,void*),void* aData)
0573 {
0574 m_pFieldFcnData = aData;
0575 m_pFieldFcn=_pFieldFcn;
0576 }
0577
0578 inline void ccontour::set_first_grid(int iCol, int iRow)
0579 {
0580 m_iColFir=max_of<int>(iCol,2);
0581 m_iRowFir=max_of<int>(iRow,2);
0582 }
0583
0584 inline void ccontour::set_secondary_grid(int iCol, int iRow)
0585 {
0586 // cleaning work matrices if allocated
0587 CleanMemory();
0588
0589 m_iColSec=max_of<int>(iCol,2);
0590 m_iRowSec=max_of<int>(iRow,2);
0591 }
0592
0593 inline void ccontour::set_limits(double pLimits[4])
0594 {
0595 _ASSERT_(pLimits[0]<pLimits[1],"ccontour::set_limits");
0596 _ASSERT_(pLimits[2]<pLimits[3],"ccontour::set_limits");
0597 for (int i=0;i<4;i++)
0598 {
0599 m_pLimits[i]=pLimits[i];
0600 }
0601 }
0602
0603 inline void ccontour::get_limits(double pLimits[4])
0604 {
0605 for (int i=0;i<4;i++)
0606 {
0607 pLimits[i]=m_pLimits[i];
0608 }
0609 }
0610
0611 //G.Barrand : from .h to .cxx to avoid _ASSERT_ in .h
0612 inline double ccontour::get_plane(unsigned int i) const {
0613 /*_ASSERT_(i>=0);*/
0614 _ASSERT_(i<m_vPlanes.size(),"ccontour::get_plane");
0615 return m_vPlanes[i];
0616 }
0617
0618 inline double ccontour::get_yi(int i) const {
0619 if(i<0) ::printf("ccontour::get_yi : %d\n",i);
0620 _ASSERT_(i>=0,"ccontour::get_yi");
0621 return m_pLimits[2] + i/(m_iColSec+1)*(m_pLimits[3]-m_pLimits[2])/(double)( m_iRowSec );
0622 }
0623
0624 }
0625
0626 #endif
