| 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 unsigned int m_iColSec; // secondary grid, number of columns
172 unsigned int m_iRowSec; // sec 172 unsigned 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 (unsigned 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 for (unsigned int i=0;i<m_iColSec+1;i++)
267 { 267 {
268 if (m_ppFnData[i]) 268 if (m_ppFnData[i])
269 delete[] (m_ppFnData[i]); 269 delete[] (m_ppFnData[i]);
270 } 270 }
271 delete[] m_ppFnData; 271 delete[] m_ppFnData;
272 m_ppFnData=NULL; 272 m_ppFnData=NULL;
273 } 273 }
274 } 274 }
275 275
276 inline void ccontour::generate() 276 inline void ccontour::generate()
277 { 277 {
278 278
279 int i, j; 279 int i, j;
280 int x3, x4, y3, y4, x, y, oldx3, xlow; 280 int x3, x4, y3, y4, x, y, oldx3, xlow;
281 const unsigned int cols=m_iColSec+1; 281 const unsigned int cols=m_iColSec+1;
282 const unsigned int rows=m_iRowSec+1; 282 const unsigned int rows=m_iRowSec+1;
283 //double xoff,yoff; 283 //double xoff,yoff;
284 284
285 // Initialize memroy if needed 285 // Initialize memroy if needed
286 InitMemory(); 286 InitMemory();
287 287
288 m_dDx = (m_pLimits[1]-m_pLimits[0])/(double) 288 m_dDx = (m_pLimits[1]-m_pLimits[0])/(double)(m_iColSec);
289 //xoff = m_pLimits[0]; 289 //xoff = m_pLimits[0];
290 m_dDy = (m_pLimits[3]-m_pLimits[2])/(double) 290 m_dDy = (m_pLimits[3]-m_pLimits[2])/(double)(m_iRowSec);
291 //yoff = m_pLimits[2]; 291 //yoff = m_pLimits[2];
292 292
293 xlow = 0; 293 xlow = 0;
294 oldx3 = 0; 294 oldx3 = 0;
295 x3 = (cols-1)/m_iRowFir; 295 x3 = (cols-1)/m_iRowFir;
296 x4 = ( 2*(cols-1) )/m_iRowFir; 296 x4 = ( 2*(cols-1) )/m_iRowFir;
297 for (x = oldx3; x <= x4; x++) 297 for (x = oldx3; x <= x4; x++)
298 { /* allocate new columns needed 298 { /* allocate new columns needed
299 */ 299 */
300 if (x >= (int)cols) 300 if (x >= (int)cols)
301 break; 301 break;
302 if (m_ppFnData[x]==NULL) 302 if (m_ppFnData[x]==NULL)
303 m_ppFnData[x] = new CFnStr[rows]; 303 m_ppFnData[x] = new CFnStr[rows];
304 304
305 for (y = 0; y < (int)rows; y++) 305 for (y = 0; y < (int)rows; y++)
306 FnctData(x,y)->m_sTopLen = -1; 306 FnctData(x,y)->m_sTopLen = -1;
307 } 307 }
308 308
309 y4 = 0; 309 y4 = 0;
310 for (j = 0; j < m_iColFir; j++) 310 for (j = 0; j < m_iColFir; j++)
311 { 311 {
312 y3 = y4; 312 y3 = y4;
313 y4 = ((j+1)*(rows-1))/m_iColFir; 313 y4 = ((j+1)*(rows-1))/m_iColFir;
314 Cntr1(oldx3, x3, y3, y4); 314 Cntr1(oldx3, x3, y3, y4);
315 } 315 }
316 316
317 for (i = 1; i < m_iRowFir; i++) 317 for (i = 1; i < m_iRowFir; i++)
318 { 318 {
319 y4 = 0; 319 y4 = 0;
320 for (j = 0; j < m_iColFir; j++) 320 for (j = 0; j < m_iColFir; j++)
321 { 321 {
322 y3 = y4; 322 y3 = y4;
323 y4 = ((j+1)*(rows-1))/m_iColFir; 323 y4 = ((j+1)*(rows-1))/m_iColFir;
324 Cntr1(x3, x4, y3, y4); 324 Cntr1(x3, x4, y3, y4);
325 } 325 }
326 326
327 y4 = 0; 327 y4 = 0;
328 for (j = 0; j < m_iColFir; j++) 328 for (j = 0; j < m_iColFir; j++)
329 { 329 {
330 y3 = y4; 330 y3 = y4;
331 y4 = ((j+1)*(rows-1))/m_iColFir; 331 y4 = ((j+1)*(rows-1))/m_iColFir;
332 Pass2(oldx3,x3,y3,y4); 332 Pass2(oldx3,x3,y3,y4);
333 } 333 }
334 334
335 if (i < (m_iRowFir-1)) 335 if (i < (m_iRowFir-1))
336 { /* re-use columns no longer needed */ 336 { /* re-use columns no longer needed */
337 oldx3 = x3; 337 oldx3 = x3;
338 x3 = x4; 338 x3 = x4;
339 x4 = ((i+2)*(cols-1))/m_iRowFir; 339 x4 = ((i+2)*(cols-1))/m_iRowFir;
340 for (x = x3+1; x <= x4; x++) 340 for (x = x3+1; x <= x4; x++)
341 { 341 {
342 if (xlow < oldx3) 342 if (xlow < oldx3)
343 { 343 {
344 if (m_ppFnData[x]) 344 if (m_ppFnData[x])
345 delete[] m_ppFnData[x]; 345 delete[] m_ppFnData[x];
346 m_ppFnData[x] = m_ppFnData[xlow]; 346 m_ppFnData[x] = m_ppFnData[xlow];
347 m_ppFnData[ xlow++ ] = NULL; 347 m_ppFnData[ xlow++ ] = NULL;
348 } 348 }
349 else 349 else
350 if (m_ppFnData[x]==NULL) 350 if (m_ppFnData[x]==NULL)
351 m_ppFnData[x] = new CFnStr[rows]; 351 m_ppFnData[x] = new CFnStr[rows];
352 352
353 for (y = 0; y < (int)rows; y++) 353 for (y = 0; y < (int)rows; y++)
354 FnctData(x,y)->m_sTopLen = -1; 354 FnctData(x,y)->m_sTopLen = -1;
355 } 355 }
356 } 356 }
357 } 357 }
358 358
359 y4 = 0; 359 y4 = 0;
360 for (j = 0; j < m_iColFir; j++) 360 for (j = 0; j < m_iColFir; j++)
361 { 361 {
362 y3 = y4; 362 y3 = y4;
363 y4 = ((j+1)*(rows-1))/m_iColFir; 363 y4 = ((j+1)*(rows-1))/m_iColFir;
364 Pass2(x3,x4,y3,y4); 364 Pass2(x3,x4,y3,y4);
365 } 365 }
366 } 366 }
367 367
368 inline void ccontour::Cntr1(int x1, int x2, in 368 inline void ccontour::Cntr1(int x1, int x2, int y1, int y2)
369 { 369 {
370 double f11, f12, f21, f22, f33; 370 double f11, f12, f21, f22, f33;
371 int x3, y3, i, j; 371 int x3, y3, i, j;
372 372
373 if ((x1 == x2) || (y1 == y2)) /* if not a re 373 if ((x1 == x2) || (y1 == y2)) /* if not a real cell, punt */
374 return; 374 return;
375 f11 = Field(x1, y1); 375 f11 = Field(x1, y1);
376 f12 = Field(x1, y2); 376 f12 = Field(x1, y2);
377 f21 = Field(x2, y1); 377 f21 = Field(x2, y1);
378 f22 = Field(x2, y2); 378 f22 = Field(x2, y2);
379 if ((x2 > x1+1) || (y2 > y1+1)) { /* is cell 379 if ((x2 > x1+1) || (y2 > y1+1)) { /* is cell divisible? */
380 x3 = (x1+x2)/2; 380 x3 = (x1+x2)/2;
381 y3 = (y1+y2)/2; 381 y3 = (y1+y2)/2;
382 f33 = Field(x3, y3); 382 f33 = Field(x3, y3);
383 i = j = 0; 383 i = j = 0;
384 if (f33 < f11) i++; else if (f33 > f11) j+ 384 if (f33 < f11) i++; else if (f33 > f11) j++;
385 if (f33 < f12) i++; else if (f33 > f12) j+ 385 if (f33 < f12) i++; else if (f33 > f12) j++;
386 if (f33 < f21) i++; else if (f33 > f21) j+ 386 if (f33 < f21) i++; else if (f33 > f21) j++;
387 if (f33 < f22) i++; else if (f33 > f22) j+ 387 if (f33 < f22) i++; else if (f33 > f22) j++;
388 if ((i > 2) || (j > 2)) /* should we divid 388 if ((i > 2) || (j > 2)) /* should we divide cell? */
389 { 389 {
390 /* subdivide cell */ 390 /* subdivide cell */
391 Cntr1(x1, x3, y1, y3); 391 Cntr1(x1, x3, y1, y3);
392 Cntr1(x3, x2, y1, y3); 392 Cntr1(x3, x2, y1, y3);
393 Cntr1(x1, x3, y3, y2); 393 Cntr1(x1, x3, y3, y2);
394 Cntr1(x3, x2, y3, y2); 394 Cntr1(x3, x2, y3, y2);
395 return; 395 return;
396 } 396 }
397 } 397 }
398 /* install cell in array */ 398 /* install cell in array */
399 FnctData(x1,y2)->m_sBotLen = FnctData(x1,y1) 399 FnctData(x1,y2)->m_sBotLen = FnctData(x1,y1)->m_sTopLen = x2-x1;
400 FnctData(x2,y1)->m_sLeftLen = FnctData(x1,y1 400 FnctData(x2,y1)->m_sLeftLen = FnctData(x1,y1)->m_sRightLen = y2-y1;
401 } 401 }
402 402
403 inline void ccontour::Pass2(int x1, int x2, in 403 inline void ccontour::Pass2(int x1, int x2, int y1, int y2)
404 { 404 {
405 int left = 0, right = 0, top = 0, bot = 0,ol 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, x 406 double yy0 = 0, yy1 = 0, xx0 = 0, xx1 = 0, xx3, yy3;
407 double v, f11, f12, f21, f22, f33, fold, fne 407 double v, f11, f12, f21, f22, f33, fold, fnew, f;
408 double xoff=m_pLimits[0]; 408 double xoff=m_pLimits[0];
409 double yoff=m_pLimits[2]; 409 double yoff=m_pLimits[2];
410 410
411 if ((x1 == x2) || (y1 == y2)) /* if not a re 411 if ((x1 == x2) || (y1 == y2)) /* if not a real cell, punt */
412 return; 412 return;
413 f11 = FnctData(x1,y1)->m_dFnVal; 413 f11 = FnctData(x1,y1)->m_dFnVal;
414 f12 = FnctData(x1,y2)->m_dFnVal; 414 f12 = FnctData(x1,y2)->m_dFnVal;
415 f21 = FnctData(x2,y1)->m_dFnVal; 415 f21 = FnctData(x2,y1)->m_dFnVal;
416 f22 = FnctData(x2,y2)->m_dFnVal; 416 f22 = FnctData(x2,y2)->m_dFnVal;
417 if ((x2 > x1+1) || (y2 > y1+1)) /* is cell d 417 if ((x2 > x1+1) || (y2 > y1+1)) /* is cell divisible? */
418 { 418 {
419 x3 = (x1+x2)/2; 419 x3 = (x1+x2)/2;
420 y3 = (y1+y2)/2; 420 y3 = (y1+y2)/2;
421 f33 = FnctData(x3, y3)->m_dFnVal; 421 f33 = FnctData(x3, y3)->m_dFnVal;
422 i = j = 0; 422 i = j = 0;
423 if (f33 < f11) i++; else if (f33 > f11) j+ 423 if (f33 < f11) i++; else if (f33 > f11) j++;
424 if (f33 < f12) i++; else if (f33 > f12) j+ 424 if (f33 < f12) i++; else if (f33 > f12) j++;
425 if (f33 < f21) i++; else if (f33 > f21) j+ 425 if (f33 < f21) i++; else if (f33 > f21) j++;
426 if (f33 < f22) i++; else if (f33 > f22) j+ 426 if (f33 < f22) i++; else if (f33 > f22) j++;
427 if ((i > 2) || (j > 2)) /* should we divid 427 if ((i > 2) || (j > 2)) /* should we divide cell? */
428 { 428 {
429 /* subdivide cell */ 429 /* subdivide cell */
430 Pass2(x1, x3, y1, y3); 430 Pass2(x1, x3, y1, y3);
431 Pass2(x3, x2, y1, y3); 431 Pass2(x3, x2, y1, y3);
432 Pass2(x1, x3, y3, y2); 432 Pass2(x1, x3, y3, y2);
433 Pass2(x3, x2, y3, y2); 433 Pass2(x3, x2, y3, y2);
434 return; 434 return;
435 } 435 }
436 } 436 }
437 437
438 for (i = 0; i < (int)m_vPlanes.size(); i++) 438 for (i = 0; i < (int)m_vPlanes.size(); i++)
439 { 439 {
440 v = m_vPlanes[i]; 440 v = m_vPlanes[i];
441 j = 0; 441 j = 0;
442 if (f21 > v) j++; 442 if (f21 > v) j++;
443 if (f11 > v) j |= 2; 443 if (f11 > v) j |= 2;
444 if (f22 > v) j |= 4; 444 if (f22 > v) j |= 4;
445 if (f12 > v) j |= 010; 445 if (f12 > v) j |= 010;
446 if ((f11 > v) ^ (f12 > v)) 446 if ((f11 > v) ^ (f12 > v))
447 { 447 {
448 if ((FnctData(x1,y1)->m_sLeftLen != 0) & 448 if ((FnctData(x1,y1)->m_sLeftLen != 0) &&
449 (FnctData(x1,y1)->m_sLeftLen < FnctDat 449 (FnctData(x1,y1)->m_sLeftLen < FnctData(x1,y1)->m_sRightLen))
450 { 450 {
451 old = y1; 451 old = y1;
452 fold = f11; 452 fold = f11;
453 while (1) 453 while (1)
454 { 454 {
455 iNew = old+FnctData(x1,old)->m_sLeft 455 iNew = old+FnctData(x1,old)->m_sLeftLen;
456 fnew = FnctData(x1,iNew)->m_dFnVal; 456 fnew = FnctData(x1,iNew)->m_dFnVal;
457 if ((fnew > v) ^ (fold > v)) 457 if ((fnew > v) ^ (fold > v))
458 break; 458 break;
459 old = iNew; 459 old = iNew;
460 fold = fnew; 460 fold = fnew;
461 } 461 }
462 yy0 = ((old-y1)+(iNew-old)*(v-fold)/(f 462 yy0 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
463 } 463 }
464 else 464 else
465 yy0 = (v-f11)/(f12-f11); 465 yy0 = (v-f11)/(f12-f11);
466 466
467 left = (int)(y1+(y2-y1)*yy0+0.5); 467 left = (int)(y1+(y2-y1)*yy0+0.5);
468 } 468 }
469 if ((f21 > v) ^ (f22 > v)) 469 if ((f21 > v) ^ (f22 > v))
470 { 470 {
471 if ((FnctData(x2,y1)->m_sRightLen != 0) 471 if ((FnctData(x2,y1)->m_sRightLen != 0) &&
472 (FnctData(x2,y1)->m_sRightLen < FnctDa 472 (FnctData(x2,y1)->m_sRightLen < FnctData(x2,y1)->m_sLeftLen))
473 { 473 {
474 old = y1; 474 old = y1;
475 fold = f21; 475 fold = f21;
476 while (1) 476 while (1)
477 { 477 {
478 iNew = old+FnctData(x2,old)->m_sRigh 478 iNew = old+FnctData(x2,old)->m_sRightLen;
479 fnew = FnctData(x2,iNew)->m_dFnVal; 479 fnew = FnctData(x2,iNew)->m_dFnVal;
480 if ((fnew > v) ^ (fold > v)) 480 if ((fnew > v) ^ (fold > v))
481 break; 481 break;
482 old = iNew; 482 old = iNew;
483 fold = fnew; 483 fold = fnew;
484 } 484 }
485 yy1 = ((old-y1)+(iNew-old)*(v-fold)/(f 485 yy1 = ((old-y1)+(iNew-old)*(v-fold)/(fnew-fold))/(y2-y1);
486 } 486 }
487 else 487 else
488 yy1 = (v-f21)/(f22-f21); 488 yy1 = (v-f21)/(f22-f21);
489 489
490 right = (int)(y1+(y2-y1)*yy1+0.5); 490 right = (int)(y1+(y2-y1)*yy1+0.5);
491 } 491 }
492 if ((f21 > v) ^ (f11 > v)) 492 if ((f21 > v) ^ (f11 > v))
493 { 493 {
494 if ((FnctData(x1,y1)->m_sBotLen != 0) && 494 if ((FnctData(x1,y1)->m_sBotLen != 0) &&
495 (FnctData(x1,y1)->m_sBotLen < FnctData 495 (FnctData(x1,y1)->m_sBotLen < FnctData(x1,y1)->m_sTopLen)) {
496 old = x1; 496 old = x1;
497 fold = f11; 497 fold = f11;
498 while (1) { 498 while (1) {
499 iNew = old+FnctData(old,y1)->m_sBotL 499 iNew = old+FnctData(old,y1)->m_sBotLen;
500 fnew = FnctData(iNew,y1)->m_dFnVal; 500 fnew = FnctData(iNew,y1)->m_dFnVal;
501 if ((fnew > v) ^ (fold > v)) 501 if ((fnew > v) ^ (fold > v))
502 break; 502 break;
503 old = iNew; 503 old = iNew;
504 fold = fnew; 504 fold = fnew;
505 } 505 }
506 xx0 = ((old-x1)+(iNew-old)*(v-fold)/(f 506 xx0 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
507 } 507 }
508 else 508 else
509 xx0 = (v-f11)/(f21-f11); 509 xx0 = (v-f11)/(f21-f11);
510 510
511 bot = (int)(x1+(x2-x1)*xx0+0.5); 511 bot = (int)(x1+(x2-x1)*xx0+0.5);
512 } 512 }
513 if ((f22 > v) ^ (f12 > v)) 513 if ((f22 > v) ^ (f12 > v))
514 { 514 {
515 if ((FnctData(x1,y2)->m_sTopLen != 0) && 515 if ((FnctData(x1,y2)->m_sTopLen != 0) &&
516 (FnctData(x1,y2)->m_sTopLen < FnctData 516 (FnctData(x1,y2)->m_sTopLen < FnctData(x1,y2)->m_sBotLen)) {
517 old = x1; 517 old = x1;
518 fold = f12; 518 fold = f12;
519 while (1) { 519 while (1) {
520 iNew = old+FnctData(old,y2)->m_sTopL 520 iNew = old+FnctData(old,y2)->m_sTopLen;
521 fnew = FnctData(iNew,y2)->m_dFnVal; 521 fnew = FnctData(iNew,y2)->m_dFnVal;
522 if ((fnew > v) ^ (fold > v)) 522 if ((fnew > v) ^ (fold > v))
523 break; 523 break;
524 old = iNew; 524 old = iNew;
525 fold = fnew; 525 fold = fnew;
526 } 526 }
527 xx1 = ((old-x1)+(iNew-old)*(v-fold)/(f 527 xx1 = ((old-x1)+(iNew-old)*(v-fold)/(fnew-fold))/(x2-x1);
528 } 528 }
529 else 529 else
530 xx1 = (v-f12)/(f22-f12); 530 xx1 = (v-f12)/(f22-f12);
531 531
532 top = (int)(x1+(x2-x1)*xx1+0.5); 532 top = (int)(x1+(x2-x1)*xx1+0.5);
533 } 533 }
534 534
535 switch (j) 535 switch (j)
536 { 536 {
537 case 7: 537 case 7:
538 case 010: 538 case 010:
539 ExportLine(i,x1,left,top,y2); 539 ExportLine(i,x1,left,top,y2);
540 break; 540 break;
541 case 5: 541 case 5:
542 case 012: 542 case 012:
543 ExportLine(i,bot,y1,top,y2); 543 ExportLine(i,bot,y1,top,y2);
544 break; 544 break;
545 case 2: 545 case 2:
546 case 015: 546 case 015:
547 ExportLine(i,x1,left,bot,y1); 547 ExportLine(i,x1,left,bot,y1);
548 break; 548 break;
549 case 4: 549 case 4:
550 case 013: 550 case 013:
551 ExportLine(i,top,y2,x2,right); 551 ExportLine(i,top,y2,x2,right);
552 break; 552 break;
553 case 3: 553 case 3:
554 case 014: 554 case 014:
555 ExportLine(i,x1,left,x2,right); 555 ExportLine(i,x1,left,x2,right);
556 break; 556 break;
557 case 1: 557 case 1:
558 case 016: 558 case 016:
559 ExportLine(i,bot,y1,x2,right); 559 ExportLine(i,bot,y1,x2,right);
560 break; 560 break;
561 case 0: 561 case 0:
562 case 017: 562 case 017:
563 break; 563 break;
564 case 6: 564 case 6:
565 case 011: 565 case 011:
566 yy3 = (xx0*(yy1-yy0)+yy0)/(1.0-(xx1-xx0) 566 yy3 = (xx0*(yy1-yy0)+yy0)/(1.0-(xx1-xx0)*(yy1-yy0));
567 xx3 = yy3*(xx1-xx0)+xx0; 567 xx3 = yy3*(xx1-xx0)+xx0;
568 xx3 = x1+xx3*(x2-x1); 568 xx3 = x1+xx3*(x2-x1);
569 yy3 = y1+yy3*(y2-y1); 569 yy3 = y1+yy3*(y2-y1);
570 xx3 = xoff+xx3*m_dDx; 570 xx3 = xoff+xx3*m_dDx;
571 yy3 = yoff+yy3*m_dDy; 571 yy3 = yoff+yy3*m_dDy;
572 f = (*m_pFieldFcn)(xx3, yy3,m_pFieldFcnD 572 f = (*m_pFieldFcn)(xx3, yy3,m_pFieldFcnData);
573 if (f == v) { 573 if (f == v) {
574 ExportLine(i,bot,y1,top,y2); 574 ExportLine(i,bot,y1,top,y2);
575 ExportLine(i,x1,left,x2,right); 575 ExportLine(i,x1,left,x2,right);
576 } else 576 } else
577 if (((f > v) && (f22 > v)) || ((f < v) 577 if (((f > v) && (f22 > v)) || ((f < v) && (f22 < v))) {
578 ExportLine(i,x1,left,top,y2); 578 ExportLine(i,x1,left,top,y2);
579 ExportLine(i,bot,y1,x2,right); 579 ExportLine(i,bot,y1,x2,right);
580 } else { 580 } else {
581 ExportLine(i,x1,left,bot,y1); 581 ExportLine(i,x1,left,bot,y1);
582 ExportLine(i,top,y2,x2,right); 582 ExportLine(i,top,y2,x2,right);
583 } 583 }
584 } 584 }
585 } 585 }
586 } 586 }
587 587
588 inline double ccontour::Field(int x, int y) / 588 inline double ccontour::Field(int x, int y) /* evaluate funct if we must,*/
589 { 589 {
590 double x1, y1; 590 double x1, y1;
591 591
592 if (FnctData(x,y)->m_sTopLen != -1) /* is i 592 if (FnctData(x,y)->m_sTopLen != -1) /* is it already in the array */
593 return(FnctData(x,y)->m_dFnVal); 593 return(FnctData(x,y)->m_dFnVal);
594 594
595 /* not in the array, create new array elemen 595 /* not in the array, create new array element */
596 x1 = m_pLimits[0]+m_dDx*x; 596 x1 = m_pLimits[0]+m_dDx*x;
597 y1 = m_pLimits[2]+m_dDy*y; 597 y1 = m_pLimits[2]+m_dDy*y;
598 FnctData(x,y)->m_sTopLen = 0; 598 FnctData(x,y)->m_sTopLen = 0;
599 FnctData(x,y)->m_sBotLen = 0; 599 FnctData(x,y)->m_sBotLen = 0;
600 FnctData(x,y)->m_sRightLen = 0; 600 FnctData(x,y)->m_sRightLen = 0;
601 FnctData(x,y)->m_sLeftLen = 0; 601 FnctData(x,y)->m_sLeftLen = 0;
602 return (FnctData(x,y)->m_dFnVal = (*m_pField 602 return (FnctData(x,y)->m_dFnVal = (*m_pFieldFcn)(x1, y1,m_pFieldFcnData));
603 } 603 }
604 604
605 inline void ccontour::set_planes(const std::ve 605 inline void ccontour::set_planes(const std::vector<double>& vPlanes)
606 { 606 {
607 // cleaning memory 607 // cleaning memory
608 CleanMemory(); 608 CleanMemory();
609 609
610 m_vPlanes = vPlanes; 610 m_vPlanes = vPlanes;
611 } 611 }
612 612
613 inline void ccontour::set_field_fcn(double (*_ 613 inline void ccontour::set_field_fcn(double (*_pFieldFcn)(double, double,void*),void* aData)
614 { 614 {
615 m_pFieldFcnData = aData; 615 m_pFieldFcnData = aData;
616 m_pFieldFcn=_pFieldFcn; 616 m_pFieldFcn=_pFieldFcn;
617 } 617 }
618 618
619 inline void ccontour::set_first_grid(int iCol, 619 inline void ccontour::set_first_grid(int iCol, int iRow)
620 { 620 {
621 m_iColFir=mx<int>(iCol,2); 621 m_iColFir=mx<int>(iCol,2);
622 m_iRowFir=mx<int>(iRow,2); 622 m_iRowFir=mx<int>(iRow,2);
623 } 623 }
624 624
625 inline void ccontour::set_secondary_grid(int i 625 inline void ccontour::set_secondary_grid(int iCol, int iRow)
626 { 626 {
627 // cleaning work matrices if allocated 627 // cleaning work matrices if allocated
628 CleanMemory(); 628 CleanMemory();
629 629
630 m_iColSec=mx<int>(iCol,2); 630 m_iColSec=mx<int>(iCol,2);
631 m_iRowSec=mx<int>(iRow,2); 631 m_iRowSec=mx<int>(iRow,2);
632 } 632 }
633 633
634 inline void ccontour::set_limits(double pLimit 634 inline void ccontour::set_limits(double pLimits[4])
635 { 635 {
636 _ASSERT_(pLimits[0]<pLimits[1],"ccontour::se 636 _ASSERT_(pLimits[0]<pLimits[1],"ccontour::set_limits");
637 _ASSERT_(pLimits[2]<pLimits[3],"ccontour::se 637 _ASSERT_(pLimits[2]<pLimits[3],"ccontour::set_limits");
638 for (int i=0;i<4;i++) 638 for (int i=0;i<4;i++)
639 { 639 {
640 m_pLimits[i]=pLimits[i]; 640 m_pLimits[i]=pLimits[i];
641 } 641 }
642 } 642 }
643 643
644 inline void ccontour::get_limits(double pLimit 644 inline void ccontour::get_limits(double pLimits[4])
645 { 645 {
646 for (int i=0;i<4;i++) 646 for (int i=0;i<4;i++)
647 { 647 {
648 pLimits[i]=m_pLimits[i]; 648 pLimits[i]=m_pLimits[i];
649 } 649 }
650 } 650 }
651 651
652 //G.Barrand : from .h to .cxx to avoid _ASSERT 652 //G.Barrand : from .h to .cxx to avoid _ASSERT_ in .h
653 inline double ccontour::get_plane(unsigned int 653 inline double ccontour::get_plane(unsigned int i) const {
654 /*_ASSERT_(i>=0);*/ 654 /*_ASSERT_(i>=0);*/
655 _ASSERT_(i<m_vPlanes.size(),"ccontour::get_p 655 _ASSERT_(i<m_vPlanes.size(),"ccontour::get_plane");
656 return m_vPlanes[i]; 656 return m_vPlanes[i];
657 } 657 }
658 658
659 inline double ccontour::get_yi(int i) const { 659 inline double ccontour::get_yi(int i) const {
660 if(i<0) ::printf("ccontour::get_yi : %d\n",i 660 if(i<0) ::printf("ccontour::get_yi : %d\n",i);
661 _ASSERT_(i>=0,"ccontour::get_yi"); 661 _ASSERT_(i>=0,"ccontour::get_yi");
662 return m_pLimits[2] + i/(m_iColSec+1)*(m_pL 662 return m_pLimits[2] + i/(m_iColSec+1)*(m_pLimits[3]-m_pLimits[2])/(double)( m_iRowSec );
663 } 663 }
664 664
665 } 665 }
666 666
667 #endif 667 #endif