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[4]) 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[4]) 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