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