| Geant4 Cross Reference |
1 // Copyright (C) 2010, Guy Barrand. All rights reserved.
2 // See the file tools.license for terms.
3
4 #ifndef tools_ccontour
5 #define tools_ccontour
6
7 // G.Barrand : inline version of the one found in Lib of OSC 16.11.
8 // This code is not mine and I keep it "as it is".
9
10 // Contour.h: interface for the ccontour class.
11 //
12 // ccontour implements Contour plot algorithm descrided in
13 // IMPLEMENTATION OF
14 // AN IMPROVED CONTOUR
15 // PLOTTING ALGORITHM
16 // BY
17 //
18 // MICHAEL JOSEPH ARAMINI
19 //
20 // B.S., Stevens Institute of Technology, 1980
21 // See http://www.ultranet.com/~aramini/thesis.html
22 //
23 // Ported to C++ by Jonathan de Halleux.
24 //
25 // Using ccontour :
26 //
27 // ccontour is not directly usable. The user has to
28 // 1. derive the function ExportLine that is
29 // supposed to draw/store the segment of the contour
30 // 2. Set the function draw contour of. (using SetFieldFn
31 // The function must be declared as follows
32 // double (*myF)(double x , double y);
33 //
34 // History:
35 // 31-07-2002:
36 // - A lot of contribution from Chenggang Zhou (better strip compressions, merging, area, weight),
37 // - Got rid of ugly MFC lists for STL.
38 //////////////////////////////////////////////////////////////////////
39
40 //G.Barrand :
41 #include <vector>
42 #include <cstdio>
43 #include <cstdlib>
44 #include <cmath>
45
46 #ifdef TOOLS_MEM
47 #include "mem"
48 #endif
49
50 #include "mnmx"
51
52 namespace tools {
53
54 class ccontour
55 {
56 #ifdef TOOLS_MEM
57 static const std::string& s_class() {
58 static const std::string s_v("tools::ccontour");
59 return s_v;
60 }
61 #endif
62 protected:
63 // plots a line from (x1,y1) to (x2,y2)
64 virtual void ExportLine(int iPlane,int x1,int y1,int x2,int y2) = 0;
65
66 public:
67 ccontour();
68 virtual ~ccontour(){
69 CleanMemory();
70 #ifdef TOOLS_MEM
71 mem::decrement(s_class().c_str());
72 #endif
73 }
74 protected: //G.Barrand
75 ccontour(const ccontour&){}
76 private: //G.Barrand
77 ccontour& operator=(const ccontour&){return *this;}
78 public:
79 protected: //G.Barrand
80 // Initialize memory. Called in generate
81 virtual void InitMemory();
82 // Clean work arrays
83 virtual void CleanMemory();
84 // generates contour
85 // Before calling this functions you must
86 // 1. derive the function ExportLine that is
87 // supposed to draw/store the segment of the contour
88 // 2. Set the function draw contour of. (using SetFieldFn
89 // The function must be declared as follows
90 // double (*myF)(double x , double y);
91 public: //G.Barrand
92 virtual void generate();
93
94 // Set the dimension of the primary grid
95 void set_first_grid(int iCol, int iRow);
96 // Set the dimension of the base grid
97 void set_secondary_grid(int iCol, int iRow);
98 // Sets the region [left, right, bottom,top] to generate contour
99 void set_limits(double pLimits[4]);
100 // Sets the isocurve values
101 void set_planes(const std::vector<double>& vPlanes);
102 // Sets the pointer to the F(x,y) funtion
103 // G.Barrand : handle a user data pointer.
104 void set_field_fcn(double (*_pFieldFcn)(double, double,void*),void*);
105
106 size_t get_number_of_planes() const { return m_vPlanes.size();};
107 const std::vector<double>& get_planes() const { return m_vPlanes;};
108 double get_plane(unsigned int i) const;
109
110 // For an indexed point i on the sec. grid, returns x(i)
111 double get_xi(int i) const { return m_pLimits[0] + i%(m_iColSec+1)*(m_pLimits[1]-m_pLimits[0])/(double)( m_iColSec );};
112 // For an indexed point i on the fir. grid, returns y(i)
113 double get_yi(int i) const;
114
115 void get_limits(double pLimits[4]);
116 protected: //G.Barrand
117
118 // Retrieve dimension of grids, contouring region and isocurve
119 int GetColFir() const { return m_iColFir;};
120 int GetRowFir() const { return m_iRowFir;};
121 int GetColSec() const { return m_iColSec;};
122 int GetRowSec() const { return m_iRowSec;};
123
124 private:
125 // A structure used internally by ccontour
126 /*G.Barrand :
127 struct CFnStr {
128 double m_dFnVal;
129 short m_sLeftLen;
130 short m_sRightLen;
131 short m_sTopLen;
132 short m_sBotLen;
133 };
134 */
135 class CFnStr {
136 #ifdef TOOLS_MEM
137 static const std::string& s_class() {
138 static const std::string s_v("tools::ccontour::CFnStr");
139 return s_v;
140 }
141 #endif
142 public:
143 CFnStr():m_dFnVal(0),m_sLeftLen(0),m_sRightLen(0),m_sTopLen(0),m_sBotLen(0){
144 #ifdef TOOLS_MEM
145 mem::increment(s_class().c_str());
146 #endif
147 }
148 ~CFnStr(){
149 #ifdef TOOLS_MEM
150 mem::decrement(s_class().c_str());
151 #endif
152 }
153 protected:
154 CFnStr(const CFnStr&){}
155 CFnStr& operator=(const CFnStr&){return *this;}
156 public:
157 double m_dFnVal;
158 short m_sLeftLen;
159 short m_sRightLen;
160 short m_sTopLen;
161 short m_sBotLen;
162 };
163
164
165 protected:
166 // Accesibles variables
167 std::vector<double> m_vPlanes; // value of contour planes
168 double m_pLimits[4]; // left, right, bottom, top
169 int m_iColFir; // primary grid, number of columns
170 int m_iRowFir; // primary grid, number of rows
171 unsigned int m_iColSec; // secondary grid, number of columns
172 unsigned int m_iRowSec; // secondary grid, number of rows
173 void* m_pFieldFcnData; // G.Barrand : handle a user data pointer.
174 double (*m_pFieldFcn)(double x, double y,void*); // pointer to F(x,y) function
175
176 // Protected function
177 //virtual void ExportLine(int iPlane, int x1, int y1, int x2, int y2) = 0; // plots a line from (x1,y1) to (x2,y2)
178
179 // Work functions and variables
180 double m_dDx;
181 double m_dDy;
182 CFnStr** m_ppFnData; // pointer to mesh parts
183 CFnStr* FnctData(int i,int j) { return (m_ppFnData[i]+j);};
184
185 double Field(int x, int y); /* evaluate funct if we must, */
186 void Cntr1(int x1, int x2, int y1, int y2);
187 void Pass2(int x1, int x2, int y1, int y2); /* draws the contour lines */
188
189 private:
190 //G.Barrand : have the below in the class.
191 // A simple test function
192 static double ContourTestFunction(double x,double y,void*) {
193 return 0.5*(::cos(x+3.14/4)+::sin(y+3.14/4));
194 }
195
196 protected:
197 static void _ASSERT_(bool a_what,const char* a_where) {
198 if(!a_what) {
199 ::printf("debug : Contour : assert failure in %s\n",a_where);
200 ::exit(0);
201 }
202 }
203
204 static void _ASSERTP_(void* a_what,const char* a_where) {
205 if(!a_what) {
206 ::printf("debug : Contour : assert failure in %s\n",a_where);
207 ::exit(0);
208 }
209 }
210
211 };
212
213 // implementation :
214 //
215 // Code get on the web at :
216 // http://www.codeproject.com/cpp/contour.asp
217 //
218 // G.Barrand
219 //
220
221 //////////////////////////////////////////////////////////////////////
222 // Construction/Destruction
223 //////////////////////////////////////////////////////////////////////
224
225 inline ccontour::ccontour()
226 {
227 #ifdef TOOLS_MEM
228 mem::increment(s_class().c_str());
229 #endif
230 m_iColFir=m_iRowFir=32;
231 m_iColSec=m_iRowSec=256;
232 m_dDx=m_dDy=0;
233 m_pFieldFcnData = NULL;
234 m_pFieldFcn=NULL;
235 m_pLimits[0]=m_pLimits[2]=0;
236 m_pLimits[1]=m_pLimits[3]=5.;
237 m_ppFnData=NULL;
238
239 // temporary stuff
240 m_pFieldFcn=ContourTestFunction;
241 m_vPlanes.resize(20);
242 for (unsigned int i=0;i<m_vPlanes.size();i++)
243 {
244 m_vPlanes[i]=(i-m_vPlanes.size()/2.0)*0.1;
245 }
246 }
247
248 //G.Barrand : inline
249
250 inline void ccontour::InitMemory()
251 {
252 if (!m_ppFnData)
253 {
254 m_ppFnData=new CFnStr*[m_iColSec+1];
255 for (unsigned int i=0;i<m_iColSec+1;i++)
256 {
257 m_ppFnData[i]=NULL;
258 }
259 }
260 }
261
262 inline void ccontour::CleanMemory()
263 {
264 if (m_ppFnData)
265 {
266 for (unsigned int i=0;i<m_iColSec+1;i++)
267 {
268 if (m_ppFnData[i])
269 delete[] (m_ppFnData[i]);
270 }
271 delete[] m_ppFnData;
272 m_ppFnData=NULL;
273 }
274 }
275
276 inline void ccontour::generate()
277 {
278
279 int i, j;
280 int x3, x4, y3, y4, x, y, oldx3, xlow;
281 const unsigned int cols=m_iColSec+1;
282 const unsigned int rows=m_iRowSec+1;
283 //double xoff,yoff;
284
285 // Initialize memroy if needed
286 InitMemory();
287
288 m_dDx = (m_pLimits[1]-m_pLimits[0])/(double)(m_iColSec);
289 //xoff = m_pLimits[0];
290 m_dDy = (m_pLimits[3]-m_pLimits[2])/(double)(m_iRowSec);
291 //yoff = m_pLimits[2];
292
293 xlow = 0;
294 oldx3 = 0;
295 x3 = (cols-1)/m_iRowFir;
296 x4 = ( 2*(cols-1) )/m_iRowFir;
297 for (x = oldx3; x <= x4; x++)
298 { /* allocate new columns needed
299 */
300 if (x >= (int)cols)
301 break;
302 if (m_ppFnData[x]==NULL)
303 m_ppFnData[x] = new CFnStr[rows];
304
305 for (y = 0; y < (int)rows; y++)
306 FnctData(x,y)->m_sTopLen = -1;
307 }
308
309 y4 = 0;
310 for (j = 0; j < m_iColFir; j++)
311 {
312 y3 = y4;
313 y4 = ((j+1)*(rows-1))/m_iColFir;
314 Cntr1(oldx3, x3, y3, y4);
315 }
316
317 for (i = 1; i < m_iRowFir; i++)
318 {
319 y4 = 0;
320 for (j = 0; j < m_iColFir; j++)
321 {
322 y3 = y4;
323 y4 = ((j+1)*(rows-1))/m_iColFir;
324 Cntr1(x3, x4, y3, y4);
325 }
326
327 y4 = 0;
328 for (j = 0; j < m_iColFir; j++)
329 {
330 y3 = y4;
331 y4 = ((j+1)*(rows-1))/m_iColFir;
332 Pass2(oldx3,x3,y3,y4);
333 }
334
335 if (i < (m_iRowFir-1))
336 { /* re-use columns no longer needed */
337 oldx3 = x3;
338 x3 = x4;
339 x4 = ((i+2)*(cols-1))/m_iRowFir;
340 for (x = x3+1; x <= x4; x++)
341 {
342 if (xlow < oldx3)
343 {
344 if (m_ppFnData[x])
345 delete[] m_ppFnData[x];
346 m_ppFnData[x] = m_ppFnData[xlow];
347 m_ppFnData[ xlow++ ] = NULL;
348 }
349 else
350 if (m_ppFnData[x]==NULL)
351 m_ppFnData[x] = new CFnStr[rows];
352
353 for (y = 0; y < (int)rows; y++)
354 FnctData(x,y)->m_sTopLen = -1;
355 }
356 }
357 }
358
359 y4 = 0;
360 for (j = 0; j < m_iColFir; j++)
361 {
362 y3 = y4;
363 y4 = ((j+1)*(rows-1))/m_iColFir;
364 Pass2(x3,x4,y3,y4);
365 }
366 }
367
368 inline void ccontour::Cntr1(int x1, int x2, int y1, int y2)
369 {
370 double f11, f12, f21, f22, f33;
371 int x3, y3, i, j;
372
373 if ((x1 == x2) || (y1 == y2)) /* if not a real cell, punt */
374 return;
375 f11 = Field(x1, y1);
376 f12 = Field(x1, y2);
377 f21 = Field(x2, y1);
378 f22 = Field(x2, y2);
379 if ((x2 > x1+1) || (y2 > y1+1)) { /* is cell divisible? */
380 x3 = (x1+x2)/2;
381 y3 = (y1+y2)/2;
382 f33 = Field(x3, y3);
383 i = j = 0;
384 if (f33 < f11) i++; else if (f33 > f11) j++;
385 if (f33 < f12) i++; else if (f33 > f12) j++;
386 if (f33 < f21) i++; else if (f33 > f21) j++;
387 if (f33 < f22) i++; else if (f33 > f22) j++;
388 if ((i > 2) || (j > 2)) /* should we divide cell? */
389 {
390 /* subdivide cell */
391 Cntr1(x1, x3, y1, y3);
392 Cntr1(x3, x2, y1, y3);
393 Cntr1(x1, x3, y3, y2);
394 Cntr1(x3, x2, y3, y2);
395 return;
396 }
397 }
398 /* install cell in array */
399 FnctData(x1,y2)->m_sBotLen = FnctData(x1,y1)->m_sTopLen = x2-x1;
400 FnctData(x2,y1)->m_sLeftLen = FnctData(x1,y1)->m_sRightLen = y2-y1;
401 }
402
403 inline void ccontour::Pass2(int x1, int x2, int y1, int y2)
404 {
405 int left = 0, right = 0, top = 0, bot = 0,old, iNew, i, j, x3, y3;
406 double yy0 = 0, yy1 = 0, xx0 = 0, xx1 = 0, xx3, yy3;
407 double v, f11, f12, f21, f22, f33, fold, fnew, f;
408 double xoff=m_pLimits[0];
409 double yoff=m_pLimits[2];
410
411 if ((x1 == x2) || (y1 == y2)) /* if not a real cell, punt */
412 return;
413 f11 = FnctData(x1,y1)->m_dFnVal;
414 f12 = FnctData(x1,y2)->m_dFnVal;
415 f21 = FnctData(x2,y1)->m_dFnVal;
416 f22 = FnctData(x2,y2)->m_dFnVal;
417 if ((x2 > x1+1) || (y2 > y1+1)) /* is cell divisible? */
418 {
419 x3 = (x1+x2)/2;
420 y3 = (y1+y2)/2;
421 f33 = FnctData(x3, y3)->m_dFnVal;
422 i = j = 0;
423 if (f33 < f11) i++; else if (f33 > f11) j++;
424 if (f33 < f12) i++; else if (f33 > f12) j++;
425 if (f33 < f21) i++; else if (f33 > f21) j++;
426 if (f33 < f22) i++; else if (f33 > f22) j++;
427 if ((i > 2) || (j > 2)) /* should we divide cell? */
428 {
429 /* subdivide cell */
430 Pass2(x1, x3, y1, y3);
431 Pass2(x3, x2, y1, y3);
432 Pass2(x1, x3, y3, y2);
433 Pass2(x3, x2, y3, y2);
434 return;
435 }
436 }
437
438 for (i = 0; i < (int)m_vPlanes.size(); i++)
439 {
440 v = m_vPlanes[i];
441 j = 0;
442 if (f21 > v) j++;
443 if (f11 > v) j |= 2;
444 if (f22 > v) j |= 4;
445 if (f12 > v) j |= 010;
446 if ((f11 > v) ^ (f12 > v))
447 {
448 if ((FnctData(x1,y1)->m_sLeftLen != 0) &&
449 (FnctData(x1,y1)->m_sLeftLen < FnctData(x1,y1)->m_sRightLen))
450 {
451 old = y1;
452 fold = f11;
453 while (1)
454 {
455 iNew = old+FnctData(x1,old)->m_sLeftLen;
456 fnew = FnctData(x1,iNew)->m_dFnVal;
457 if ((fnew > v) ^ (fold > v))
458 break;
459 old = iNew;
460 fold = fnew;
461 }
462 yy0 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
463 }
464 else
465 yy0 = (v-f11)/(f12-f11);
466
467 left = (int)(y1+(y2-y1)*yy0+0.5);
468 }
469 if ((f21 > v) ^ (f22 > v))
470 {
471 if ((FnctData(x2,y1)->m_sRightLen != 0) &&
472 (FnctData(x2,y1)->m_sRightLen < FnctData(x2,y1)->m_sLeftLen))
473 {
474 old = y1;
475 fold = f21;
476 while (1)
477 {
478 iNew = old+FnctData(x2,old)->m_sRightLen;
479 fnew = FnctData(x2,iNew)->m_dFnVal;
480 if ((fnew > v) ^ (fold > v))
481 break;
482 old = iNew;
483 fold = fnew;
484 }
485 yy1 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
486 }
487 else
488 yy1 = (v-f21)/(f22-f21);
489
490 right = (int)(y1+(y2-y1)*yy1+0.5);
491 }
492 if ((f21 > v) ^ (f11 > v))
493 {
494 if ((FnctData(x1,y1)->m_sBotLen != 0) &&
495 (FnctData(x1,y1)->m_sBotLen < FnctData(x1,y1)->m_sTopLen)) {
496 old = x1;
497 fold = f11;
498 while (1) {
499 iNew = old+FnctData(old,y1)->m_sBotLen;
500 fnew = FnctData(iNew,y1)->m_dFnVal;
501 if ((fnew > v) ^ (fold > v))
502 break;
503 old = iNew;
504 fold = fnew;
505 }
506 xx0 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
507 }
508 else
509 xx0 = (v-f11)/(f21-f11);
510
511 bot = (int)(x1+(x2-x1)*xx0+0.5);
512 }
513 if ((f22 > v) ^ (f12 > v))
514 {
515 if ((FnctData(x1,y2)->m_sTopLen != 0) &&
516 (FnctData(x1,y2)->m_sTopLen < FnctData(x1,y2)->m_sBotLen)) {
517 old = x1;
518 fold = f12;
519 while (1) {
520 iNew = old+FnctData(old,y2)->m_sTopLen;
521 fnew = FnctData(iNew,y2)->m_dFnVal;
522 if ((fnew > v) ^ (fold > v))
523 break;
524 old = iNew;
525 fold = fnew;
526 }
527 xx1 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
528 }
529 else
530 xx1 = (v-f12)/(f22-f12);
531
532 top = (int)(x1+(x2-x1)*xx1+0.5);
533 }
534
535 switch (j)
536 {
537 case 7:
538 case 010:
539 ExportLine(i,x1,left,top,y2);
540 break;
541 case 5:
542 case 012:
543 ExportLine(i,bot,y1,top,y2);
544 break;
545 case 2:
546 case 015:
547 ExportLine(i,x1,left,bot,y1);
548 break;
549 case 4:
550 case 013:
551 ExportLine(i,top,y2,x2,right);
552 break;
553 case 3:
554 case 014:
555 ExportLine(i,x1,left,x2,right);
556 break;
557 case 1:
558 case 016:
559 ExportLine(i,bot,y1,x2,right);
560 break;
561 case 0:
562 case 017:
563 break;
564 case 6:
565 case 011:
566 yy3 = (xx0*(yy1-yy0)+yy0)/(1.0-(xx1-xx0)*(yy1-yy0));
567 xx3 = yy3*(xx1-xx0)+xx0;
568 xx3 = x1+xx3*(x2-x1);
569 yy3 = y1+yy3*(y2-y1);
570 xx3 = xoff+xx3*m_dDx;
571 yy3 = yoff+yy3*m_dDy;
572 f = (*m_pFieldFcn)(xx3, yy3,m_pFieldFcnData);
573 if (f == v) {
574 ExportLine(i,bot,y1,top,y2);
575 ExportLine(i,x1,left,x2,right);
576 } else
577 if (((f > v) && (f22 > v)) || ((f < v) && (f22 < v))) {
578 ExportLine(i,x1,left,top,y2);
579 ExportLine(i,bot,y1,x2,right);
580 } else {
581 ExportLine(i,x1,left,bot,y1);
582 ExportLine(i,top,y2,x2,right);
583 }
584 }
585 }
586 }
587
588 inline double ccontour::Field(int x, int y) /* evaluate funct if we must,*/
589 {
590 double x1, y1;
591
592 if (FnctData(x,y)->m_sTopLen != -1) /* is it already in the array */
593 return(FnctData(x,y)->m_dFnVal);
594
595 /* not in the array, create new array element */
596 x1 = m_pLimits[0]+m_dDx*x;
597 y1 = m_pLimits[2]+m_dDy*y;
598 FnctData(x,y)->m_sTopLen = 0;
599 FnctData(x,y)->m_sBotLen = 0;
600 FnctData(x,y)->m_sRightLen = 0;
601 FnctData(x,y)->m_sLeftLen = 0;
602 return (FnctData(x,y)->m_dFnVal = (*m_pFieldFcn)(x1, y1,m_pFieldFcnData));
603 }
604
605 inline void ccontour::set_planes(const std::vector<double>& vPlanes)
606 {
607 // cleaning memory
608 CleanMemory();
609
610 m_vPlanes = vPlanes;
611 }
612
613 inline void ccontour::set_field_fcn(double (*_pFieldFcn)(double, double,void*),void* aData)
614 {
615 m_pFieldFcnData = aData;
616 m_pFieldFcn=_pFieldFcn;
617 }
618
619 inline void ccontour::set_first_grid(int iCol, int iRow)
620 {
621 m_iColFir=mx<int>(iCol,2);
622 m_iRowFir=mx<int>(iRow,2);
623 }
624
625 inline void ccontour::set_secondary_grid(int iCol, int iRow)
626 {
627 // cleaning work matrices if allocated
628 CleanMemory();
629
630 m_iColSec=mx<int>(iCol,2);
631 m_iRowSec=mx<int>(iRow,2);
632 }
633
634 inline void ccontour::set_limits(double pLimits[4])
635 {
636 _ASSERT_(pLimits[0]<pLimits[1],"ccontour::set_limits");
637 _ASSERT_(pLimits[2]<pLimits[3],"ccontour::set_limits");
638 for (int i=0;i<4;i++)
639 {
640 m_pLimits[i]=pLimits[i];
641 }
642 }
643
644 inline void ccontour::get_limits(double pLimits[4])
645 {
646 for (int i=0;i<4;i++)
647 {
648 pLimits[i]=m_pLimits[i];
649 }
650 }
651
652 //G.Barrand : from .h to .cxx to avoid _ASSERT_ in .h
653 inline double ccontour::get_plane(unsigned int i) const {
654 /*_ASSERT_(i>=0);*/
655 _ASSERT_(i<m_vPlanes.size(),"ccontour::get_plane");
656 return m_vPlanes[i];
657 }
658
659 inline double ccontour::get_yi(int i) const {
660 if(i<0) ::printf("ccontour::get_yi : %d\n",i);
661 _ASSERT_(i>=0,"ccontour::get_yi");
662 return m_pLimits[2] + i/(m_iColSec+1)*(m_pLimits[3]-m_pLimits[2])/(double)( m_iRowSec );
663 }
664
665 }
666
667 #endif