00001 /*************************************************************************** 00002 ** 00003 ** This file is part of QGpCoreTools. 00004 ** 00005 ** This library is free software; you can redistribute it and/or 00006 ** modify it under the terms of the GNU Lesser General Public 00007 ** License as published by the Free Software Foundation; either 00008 ** version 2.1 of the License, or (at your option) any later version. 00009 ** 00010 ** This file is distributed in the hope that it will be useful, but WITHOUT 00011 ** ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 00012 ** FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 00013 ** License for more details. 00014 ** 00015 ** You should have received a copy of the GNU Lesser General Public 00016 ** License along with this library; if not, write to the Free Software 00017 ** Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00018 ** 00019 ** See http://www.geopsy.org for more information. 00020 ** 00021 ** Created : 2009-02-15 00022 ** Authors : 00023 ** Marc Wathelet 00024 ** Marc Wathelet (LGIT, Grenoble, France) 00025 ** 00026 ***************************************************************************/ 00027 00028 #ifndef GRID3D_H 00029 #define GRID3D_H 00030 00031 #include <math.h> 00032 00033 #include "XMLClass.h" 00034 #include "Grid2D.h" 00035 #include "Point.h" 00036 #include "Translations.h" 00037 #include "CoreApplication.h" 00038 #include "Plane.h" 00039 #include "QGpCoreToolsDLLExport.h" 00040 00041 namespace QGpCoreTools { 00042 00043 template <class ValueType> 00044 class QGPCORETOOLS_EXPORT Grid3D : public XMLClass 00045 { 00046 TRANSLATIONS( "Grid3D" ) 00047 public: 00048 Grid3D(); 00049 Grid3D(int nx, int ny, int nz); 00050 Grid3D(const Grid3D& m); 00051 virtual ~Grid3D(); 00052 00053 virtual const QString& xml_tagName() const {return xmlGrid3DTag;} 00054 static const QString xmlGrid3DTag; 00055 00056 int nx() const {return _nx;} 00057 int ny() const {return _ny;} 00058 int nz() const {return _nz;} 00059 void clear(); 00060 void init(int nx, int ny, int nz); 00061 void init(ValueType value); 00062 inline void init(int nx, int ny, int nz, ValueType value); 00063 00064 ValueType value(int ix, int iy, int iz) const {return _values[ix+_nx*(iy+_ny*iz)];} 00065 void setValue(int ix, int iy, int iz, double val) {_values[ix+_nx*(iy+_ny*iz)]=val;} 00066 ValueType * valuePointer(int ix, int iy, int iz) {return _values+ix+_nx*(iy+_ny*iz);} 00067 const ValueType * valuePointer(int ix, int iy, int iz) const {return _values+ix+_nx*(iy+_ny*iz);} 00068 00069 void setOrigin(const Point& p) {_origin=p;} 00070 const Point& origin() const {return _origin;} 00071 void setDeltaX(double dx) {_delta.setX(dx); _invDelta.setX(1.0/dx);} 00072 void setDeltaY(double dy) {_delta.setY(dy); _invDelta.setY(1.0/dy);} 00073 void setDeltaZ(double dz) {_delta.setZ(dz); _invDelta.setZ(1.0/dz);} 00074 double deltaX() const {return _delta.x();} 00075 double deltaY() const {return _delta.y();} 00076 double deltaZ() const {return _delta.z();} 00077 double x(int ix) const {return _origin.x()+ix*_delta.x();} 00078 double y(int iy) const {return _origin.y()+iy*_delta.y();} 00079 double z(int iz) const {return _origin.z()+iz*_delta.z();} 00080 double east(int ix) const {return x(ix)-_delta.x()*0.5;} 00081 double west(int ix) const {return x(ix)+_delta.x()*0.5;} 00082 double south(int iy) const {return y(iy)-_delta.y()*0.5;} 00083 double north(int iy) const {return y(iy)+_delta.y()*0.5;} 00084 double top(int iz) const {return z(iz)-_delta.z()*0.5;} 00085 double bottom(int iz) const {return z(iz)+_delta.z()*0.5;} 00086 Point coordinates(int ix, int iy, int iz) const; 00087 int indexOfX(double x) const {return (int) round((x - _origin.x()) * _invDelta.x());} 00088 int indexOfY(double y) const {return (int) round((y - _origin.y()) * _invDelta.y());} 00089 int indexOfZ(double z) const {return (int) round((z - _origin.z()) * _invDelta.z());} 00090 00091 ValueType minimumValue() const; 00092 ValueType maximumValue() const; 00093 00094 Grid2D<ValueType> * crossSectionX(double x) const; 00095 Grid2D<ValueType> * crossSectionY(double y) const; 00096 Grid2D<ValueType> * crossSectionZ(double z) const; 00097 protected: 00098 virtual void xml_writeProperties(XML_WRITEPROPERTIES_ARGS) const; 00099 virtual bool xml_setProperty(XML_SETPROPERTY_ARGS); 00100 virtual XMLMember xml_member(XML_MEMBER_ARGS); 00101 virtual void xml_writeBinaryData(XML_WRITEBINARYDATA_ARGS) const; 00102 virtual bool xml_setBinaryData(XML_SETBINARYDATA_ARGS); 00103 private: 00104 int _nx, _ny, _nz; 00105 Point _origin, _delta, _invDelta; 00106 ValueType *_values; 00107 }; 00108 00109 00110 // Restore from a QTextStream 00111 template <class ValueType> 00112 QGPCORETOOLS_EXPORT QTextStream& operator<<(QTextStream& s, Grid3D<ValueType>& grd); 00113 template <class ValueType> 00114 QGPCORETOOLS_EXPORT QTextStream& operator>>(QTextStream& s, Grid3D<ValueType>& grd); 00115 00116 template <class ValueType> 00117 const QString Grid3D<ValueType>::xmlGrid3DTag="Grid3D"; 00118 00119 template <class ValueType> 00120 Grid3D<ValueType>::Grid3D() 00121 { 00122 TRACE; 00123 _nx=_ny=_nz=0; 00124 _delta=Point(1., 1. ); 00125 _invDelta=Point(1., 1. ); 00126 _origin=Point(0., 0. ); 00127 _values=0; 00128 } 00129 00130 template <class ValueType> 00131 Grid3D<ValueType>::Grid3D(int nx, int ny, int nz) 00132 { 00133 TRACE; 00134 _values=0; 00135 _delta=Point(1., 1., 1. ); 00136 _invDelta=Point(1., 1., 1. ); 00137 _origin=Point(0., 0., 0. ); 00138 init(nx, ny, nz); 00139 } 00140 00141 template <class ValueType> 00142 Grid3D<ValueType>::Grid3D(const Grid3D& m) : XMLClass() 00143 { 00144 TRACE; 00145 _values=0; 00146 _delta=m._delta; 00147 _invDelta=m._invDelta; 00148 _origin=m._origin; 00149 init(m._nx, m._ny, m._nz); 00150 int n=_ny * _nx * _nz; 00151 for(int i=0; i < n; i++ ) 00152 _values[ i ]=m._values[ i ]; 00153 } 00154 00155 template <class ValueType> 00156 Grid3D<ValueType>::~Grid3D() 00157 { 00158 TRACE; 00159 if(_values) delete [] _values; 00160 } 00161 00162 template <class ValueType> 00163 void Grid3D<ValueType>::init(int nx, int ny, int nz) 00164 { 00165 TRACE; 00166 if(_values) delete [] _values; 00167 _nx=nx; 00168 _ny=ny; 00169 _nz=nz; 00170 _values=new ValueType[ _nz * _ny * _nx ]; 00171 } 00172 00173 template <class ValueType> 00174 void Grid3D<ValueType>::init(ValueType value) 00175 { 00176 TRACE; 00177 int n=_nx * _ny * _nz; 00178 for(int i=0; i < n; i++ ) _values[ i ]=value; 00179 } 00180 00181 template <class ValueType> 00182 inline void Grid3D<ValueType>::init(int nx, int ny, int nz, ValueType value) 00183 { 00184 init(nx, ny, nz); 00185 init(value); 00186 } 00187 00188 template <class ValueType> 00189 ValueType Grid3D<ValueType>::minimumValue() const 00190 { 00191 TRACE; 00192 int n=_nx*_ny*_nz; 00193 if(n<1) return 0; 00194 ValueType val=_values[0]; 00195 for(int i=1; i<n; i++) { 00196 if(_values[i]<val) val=_values[i]; 00197 } 00198 return val; 00199 } 00200 00201 template <class ValueType> 00202 ValueType Grid3D<ValueType>::maximumValue() const 00203 { 00204 TRACE; 00205 int n=_nx*_ny*_nz; 00206 if(n<1) return 0; 00207 ValueType val=_values[0]; 00208 for(int i=1; i<n; i++) { 00209 if(_values[i]>val) val=_values[i]; 00210 } 00211 return val; 00212 } 00213 00214 template <class ValueType> 00215 Point Grid3D<ValueType>::coordinates(int ix, int iy, int iz) const 00216 { 00217 TRACE; 00218 double x=_origin.x() + (double) ix * _delta.x(); 00219 double y=_origin.y() + (double) iy * _delta.y(); 00220 double z=_origin.z() + (double) iz * _delta.z(); 00221 return (Point(x, y, z)); 00222 } 00223 00224 template <class ValueType> 00225 Grid2D<ValueType> * Grid3D<ValueType>::crossSectionX(double x) const 00226 { 00227 int ix=indexOfX(x); 00228 if(ix<0 || ix>=_nx) { 00229 return 0; 00230 } 00231 Grid2D<ValueType> * section=new Grid2D<ValueType>(_ny, _nz); 00232 section->setDeltaX(_delta.y()); 00233 section->setDeltaY(_delta.z()); 00234 section->setOrigin(Point2D(_origin.y(), _origin.z())); 00235 ValueType * gridPtr=_values; 00236 ValueType * sectionPtr=section->valuePointer(0, 0); 00237 gridPtr+=ix; 00238 for(int iz=0; iz<_nz; iz++) { 00239 for(int iy=0; iy<_ny; iy++) { 00240 *(sectionPtr++)=*gridPtr; 00241 gridPtr+=_nx; 00242 } 00243 } 00244 return section; 00245 } 00246 00247 template <class ValueType> 00248 Grid2D<ValueType> * Grid3D<ValueType>::crossSectionY(double y) const 00249 { 00250 int iy=indexOfY(y); 00251 if(iy<0 || iy>=_ny) { 00252 return 0; 00253 } 00254 Grid2D<ValueType> * section=new Grid2D<ValueType>(_nx, _nz); 00255 section->setDeltaX(_delta.x()); 00256 section->setDeltaY(_delta.z()); 00257 section->setOrigin(Point2D(_origin.x(), _origin.z())); 00258 ValueType * gridPtr=_values; 00259 ValueType * sectionPtr=section->valuePointer(0, 0); 00260 gridPtr+=iy*_nx; 00261 int d=(_ny-1)*_nx; 00262 for(int iz=0; iz<_nz; iz++) { 00263 for(int ix=0; ix<_nx; ix++) { 00264 *(sectionPtr++)=*(gridPtr++); 00265 } 00266 gridPtr+=d; 00267 } 00268 return section; 00269 } 00270 00271 template <class ValueType> 00272 Grid2D<ValueType> * Grid3D<ValueType>::crossSectionZ(double z) const 00273 { 00274 int iz=indexOfZ(z); 00275 if(iz<0 || iz>=_nz) { 00276 return 0; 00277 } 00278 Grid2D<ValueType> * section=new Grid2D<ValueType>(_nx, _ny); 00279 section->setDeltaX(_delta.x()); 00280 section->setDeltaY(_delta.y()); 00281 section->setOrigin(Point2D(_origin.x(), _origin.y())); 00282 ValueType * gridPtr=_values; 00283 ValueType * sectionPtr=section->valuePointer(0, 0); 00284 gridPtr+=iz*_nx*_ny; 00285 for(int iy=0; iy<_ny; iy++) { 00286 for(int ix=0; ix<_nx; ix++) { 00287 *(sectionPtr++)=*(gridPtr++); 00288 } 00289 } 00290 return section; 00291 } 00292 00293 template <class ValueType> 00294 void Grid3D<ValueType>::xml_writeProperties(XML_WRITEPROPERTIES_ARGS) const 00295 { 00296 TRACE; 00297 writeProperty(s, "nx", _nx); 00298 writeProperty(s, "ny", _ny); 00299 writeProperty(s, "nz", _nz); 00300 writeProperty(s, "origin", _origin.toString()); 00301 writeProperty(s, "delta", _delta.toString()); 00302 writeBinaryData(s, context); 00303 } 00304 00305 template <class ValueType> 00306 void Grid3D<ValueType>::xml_writeBinaryData(XML_WRITEBINARYDATA_ARGS) const 00307 { 00308 TRACE; 00309 Q_UNUSED(context); 00310 s << _nx; 00311 s << _ny; 00312 s << _nz; 00313 s.writeRawData((const char *)_values, sizeof(double)*_nx * _ny * _nz); 00314 } 00315 00316 template <class ValueType> 00317 bool Grid3D<ValueType>::xml_setBinaryData(XML_SETBINARYDATA_ARGS) 00318 { 00319 TRACE; 00320 Q_UNUSED(context); 00321 int nx, ny, nz; 00322 s >> nx; 00323 s >> ny; 00324 s >> nz; 00325 if(nx!=_nx || ny!=ny || nz!=nz) { 00326 App::stream() << tr( "Grid3D size in binary file: %1x%2x%3\n" 00327 " in XML data : %4x%5x%6" ) 00328 .arg(nx).arg(ny).arg(nz).arg(_nx).arg(_ny).arg(_nz) << endl; 00329 _nx=nx; 00330 _ny=ny; 00331 _nz=nz; 00332 } 00333 init(_nx, _ny, _nz); 00334 s.readRawData((char *)_values, sizeof(double)*_nx * _ny * _nz); 00335 return true; 00336 } 00337 00338 template <class ValueType> 00339 XMLMember Grid3D<ValueType>::xml_member(XML_MEMBER_ARGS) 00340 { 00341 TRACE; 00342 Q_UNUSED(attributes); 00343 Q_UNUSED(context); 00344 if(tag=="origin" ) return XMLMember(0); 00345 else if(tag=="delta" ) return XMLMember(1); 00346 else if(tag=="nx" ) return XMLMember(2); 00347 else if(tag=="ny" ) return XMLMember(3); 00348 else if(tag=="nz" ) return XMLMember(4); 00349 return XMLMember(XMLMember::Unknown); 00350 } 00351 00352 template <class ValueType> 00353 bool Grid3D<ValueType>::xml_setProperty(XML_SETPROPERTY_ARGS) 00354 { 00355 TRACE; 00356 Q_UNUSED(attributes); 00357 Q_UNUSED(context); 00358 switch (memberID) { 00359 case 0: _origin.fromString(content); return true; 00360 case 1: 00361 _delta.fromString(content); 00362 _invDelta.setX(1.0/_delta.x()); 00363 _invDelta.setY(1.0/_delta.y()); 00364 _invDelta.setZ(1.0/_delta.z()); 00365 return true; 00366 case 2: _nx=content.toInt(); return true; 00367 case 3: _ny=content.toInt(); return true; 00368 case 4: _nz=content.toInt(); return true; 00369 } 00370 return false; 00371 } 00372 00373 template <class ValueType> 00374 QTextStream& operator<<(QTextStream& s, const Grid3D<double>& grd) 00375 { 00376 TRACE; 00377 s << "x y z val" << ::endl; 00378 const double * ptr=grd.valuePointer(0, 0, 0); 00379 for(int iz=0;iz < grd.nz();iz++ ) { 00380 for(int iy=0;iy < grd.ny();iy++ ) { 00381 for(int ix=0;ix < grd.nx();ix++ ) 00382 s << grd.x(ix) << " " << grd.y(iy) << " " << grd.z(iz) << " " << *(ptr++ ) << ::endl; 00383 } 00384 } 00385 return s; 00386 } 00387 00388 template <class ValueType> 00389 QTextStream& operator>>(QTextStream& s, Grid3D<ValueType>& grd) 00390 { 00391 TRACE; 00392 QString str; 00393 bool ok; 00394 int i; 00395 int line=0; 00396 str=s.readLine().trimmed(); line++; 00397 while(str.isEmpty() || str[0]=='#') { // Accept blank lines and comments only before grid definition 00398 str=s.readLine().trimmed(); line++; 00399 } 00400 if(str=="x y z val" ) { 00401 QVector<Point> points; 00402 QVector<ValueType> values; 00403 Point p; 00404 ValueType val; 00405 const QChar * ptr; 00406 StringSection f; 00407 bool ok=true; 00408 while( !s.atEnd()) { 00409 str=s.readLine().trimmed(); line++; 00410 if(str.isEmpty() || str[0]=='#') break; 00411 StringSection pStr(str); 00412 ptr=0; 00413 f=pStr.nextField(ptr); 00414 if(f.isValid()) p.setX(f.toDouble()); else ok=false; 00415 f=pStr.nextField(ptr); 00416 if(f.isValid()) p.setY(f.toDouble()); else ok=false; 00417 f=pStr.nextField(ptr); 00418 if(f.isValid()) p.setZ(f.toDouble()); else ok=false; 00419 f=pStr.nextField(ptr); 00420 if(f.isValid()) val=f.toDouble(); else ok=false; 00421 if( !ok) { 00422 App::stream() << QCoreApplication::translate( "Grid3D", "Wrong format at line %1, expected: float value" ).arg(line) << endl; 00423 return s; 00424 } 00425 points.append(p); 00426 values.append(val); 00427 } 00428 // Deduce nx, ny, nz, delta and origin from x, y and z vectors 00429 QMap<double,int> xMap, yMap, zMap; 00430 QMap<double,int>::iterator itMap; 00431 for(QVector<Point>::iterator it=points.begin(); it!=points.end(); it++ ) { 00432 itMap=xMap.find(it->x()); 00433 if(itMap==xMap.end()) { 00434 xMap.insert(it->x(),0); 00435 } 00436 itMap=yMap.find(it->y()); 00437 if(itMap==yMap.end()) { 00438 yMap.insert(it->y(),0); 00439 } 00440 itMap=zMap.find(it->z()); 00441 if(itMap==zMap.end()) { 00442 zMap.insert(it->z(),0); 00443 } 00444 } 00445 if(xMap.isEmpty() || yMap.isEmpty() || zMap.isEmpty()) { 00446 return s; 00447 } 00448 grd.setOrigin(Point(xMap.begin().value(), yMap.begin().value(), zMap.begin().value())); 00449 grd.setDeltaX(((--xMap.end()).value() - xMap.begin().value())/xMap.count()); 00450 grd.setDeltaY(((--yMap.end()).value() - yMap.begin().value())/yMap.count()); 00451 grd.setDeltaZ(((--zMap.end()).value() - zMap.begin().value())/zMap.count()); 00452 grd.init(xMap.count(), yMap.count(), zMap.count(), 0.0); 00453 // Set index in x, y and z maps 00454 int i; 00455 i=0; 00456 for(QMap<double,int>::iterator it=xMap.begin();it!=xMap.end(); it++, i++) { 00457 grd.setX(i, it.key()); 00458 it.value()=i; 00459 //printf("x[%i]=%lf\n",i,it.key()); 00460 } 00461 i=0; 00462 for(QMap<double,int>::iterator it=yMap.begin();it!=yMap.end(); it++, i++) { 00463 grd.setY(i, it.key()); 00464 it.value()=i; 00465 //printf("y[%i]=%lf\n",i,it.key()); 00466 } 00467 i=0; 00468 for(QMap<double,int>::iterator it=zMap.begin();it!=zMap.end(); it++, i++) { 00469 grd.setZ(i, it.key()); 00470 it.value()=i; 00471 //printf("z[%i]=%lf\n",i,it.key()); 00472 } 00473 // Set values 00474 for(int i=points.count()-1; i>=0; i--) { 00475 Point& p=points.at(i); 00476 grd.setValue(xMap[p.x()], yMap[p.y()], zMap[p.z()], values.at(i)); 00477 } 00478 } else { 00479 App::stream() << QCoreApplication::translate( "Grid3D", "Wrong format at line %1, expected: values" ).arg(line) << endl; 00480 } 00481 return s; 00482 } 00483 00484 } // namespace QGpCoreTools 00485 00486 #endif // GRID3D.H
