diff --git a/GammaRay.pro b/GammaRay.pro index a779a1df..ec80f659 100644 --- a/GammaRay.pro +++ b/GammaRay.pro @@ -36,6 +36,7 @@ win32 { SOURCES += main.cpp\ dialogs/choosevariabledialog.cpp \ dialogs/contactanalysisdialog.cpp \ + dialogs/driftanalysisdialog.cpp \ dialogs/faciestransitionmatrixoptionsdialog.cpp \ dialogs/gridrepositiondialog.cpp \ dialogs/listbuilderdialog.cpp \ @@ -47,16 +48,21 @@ SOURCES += main.cpp\ domain/auxiliary/verticalproportioncurvemaker.cpp \ domain/section.cpp \ domain/verticalproportioncurve.cpp \ + geometry/boundingbox.cpp \ geometry/intersectionfinder.cpp \ geometry/quadrilateral.cpp \ geometry/triangle.cpp \ geostats/contactanalysis.cpp \ + geostats/driftanalysis.cpp \ geostats/mcmcdataimputation.cpp \ + geostats/quadratic3dtrendmodelfitting.cpp \ geostats/searchannulus.cpp \ geostats/searchannulusstratigraphic.cpp \ + geostats/searchbox.cpp \ geostats/searchsphericalshell.cpp \ geostats/searchverticaldumbbell.cpp \ geostats/searchwasher.cpp \ + gslib/gslibparameterfiles/kt3dtrendmodelparameters.cpp \ gslib/gslibparams/gslibpardir.cpp \ gslib/gslibparams/widgets/widgetgslibpardir.cpp \ mainwindow.cpp \ @@ -316,6 +322,7 @@ SOURCES += main.cpp\ HEADERS += mainwindow.h \ dialogs/choosevariabledialog.h \ dialogs/contactanalysisdialog.h \ + dialogs/driftanalysisdialog.h \ dialogs/faciestransitionmatrixoptionsdialog.h \ dialogs/gridrepositiondialog.h \ dialogs/listbuilderdialog.h \ @@ -332,16 +339,21 @@ HEADERS += mainwindow.h \ domain/projectroot.h \ domain/section.h \ domain/verticalproportioncurve.h \ + geometry/boundingbox.h \ geometry/intersectionfinder.h \ geometry/quadrilateral.h \ geometry/triangle.h \ geostats/contactanalysis.h \ + geostats/driftanalysis.h \ geostats/mcmcdataimputation.h \ + geostats/quadratic3dtrendmodelfitting.h \ geostats/searchannulus.h \ geostats/searchannulusstratigraphic.h \ + geostats/searchbox.h \ geostats/searchsphericalshell.h \ geostats/searchverticaldumbbell.h \ geostats/searchwasher.h \ + gslib/gslibparameterfiles/kt3dtrendmodelparameters.h \ gslib/gslibparams/gslibpardir.h \ gslib/gslibparams/widgets/widgetgslibpardir.h \ util.h \ @@ -604,6 +616,7 @@ HEADERS += mainwindow.h \ FORMS += mainwindow.ui \ dialogs/choosevariabledialog.ui \ dialogs/contactanalysisdialog.ui \ + dialogs/driftanalysisdialog.ui \ dialogs/faciestransitionmatrixoptionsdialog.ui \ dialogs/gridrepositiondialog.ui \ dialogs/listbuilderdialog.ui \ @@ -868,7 +881,7 @@ win32 { # The application version # Don't forget to update the Util::importSettingsFromPreviousVersion() method to # enable the import of registry/user settings of previous versions. -VERSION = 6.20 +VERSION = 6.22 # Define a preprocessor macro so we can get the application version in application code. DEFINES += APP_VERSION=\\\"$$VERSION\\\" diff --git a/README.md b/README.md index b8a27bc8..9ef681ff 100644 --- a/README.md +++ b/README.md @@ -15,6 +15,7 @@ If you enjoyed this project, you might also enjoy GeostatsPy: https://github.com Python script to convert Eclipse grids to Paraview-compatible VTU format: https://github.com/BinWang0213/PyGRDECL VERSION HISTORY:
+   Version 6.22 - Drift analysis and drift model fitting to data.
   Version 6.20 - Contact Analysis.
   Version 6.18 - MCRFSim execution in batch/unattended mode; dependencies upgrades (VTK, ITK, Boost, C++14, ...).
   Version 6.17 - Transiography and MCRFSim for Bayesian approach; some fixes and improvements.
diff --git a/art/iconsHD/calc32.png b/art/iconsHD/calc32.png new file mode 100644 index 00000000..a3f557b7 Binary files /dev/null and b/art/iconsHD/calc32.png differ diff --git a/dialogs/cokrigingdialog.cpp b/dialogs/cokrigingdialog.cpp index cc7cb566..bf18de45 100644 --- a/dialogs/cokrigingdialog.cpp +++ b/dialogs/cokrigingdialog.cpp @@ -326,8 +326,8 @@ void CokrigingDialog::onParametersCokb3d() //load the data psInputData->loadData(); //init max with min and min with max ;-) - double minAll = std::numeric_limits::max(); - double maxAll = std::numeric_limits::min(); + double minAll = std::numeric_limits::max(); + double maxAll = -std::numeric_limits::max(); //gather the indexes of the selected variables QVector selectedColumns; selectedColumns.append( m_inputPrimVarSelector->getSelectedVariableGEOEASIndex() - 1 ); @@ -531,8 +531,8 @@ void CokrigingDialog::onParametersNewcokb3d() //load the data psInputData->loadData(); //init max with min and min with max ;-) - double minAll = std::numeric_limits::max(); - double maxAll = std::numeric_limits::min(); + double minAll = std::numeric_limits::max(); + double maxAll = -std::numeric_limits::max(); //gather the indexes of the selected variables QVector selectedColumns; selectedColumns.append( m_inputPrimVarSelector->getSelectedVariableGEOEASIndex() - 1 ); diff --git a/dialogs/contactanalysisdialog.cpp b/dialogs/contactanalysisdialog.cpp index 24e894ad..779c321d 100644 --- a/dialogs/contactanalysisdialog.cpp +++ b/dialogs/contactanalysisdialog.cpp @@ -153,7 +153,7 @@ void ContactAnalysisDialog::onProceed() LineChartWidget* lcw = new LineChartWidget(this); lcw->setSharedYaxis( true ); lcw->setChartTitle( "Contact analysis of " + m_dataFile->getName() ); - lcw->setData( chartData, 0, + lcw->setData( chartData, 0, true, {{1, nameDomain1 + " (Domain 1)"}, {2, nameDomain2 + " (Domain 2)"}}, {{2, "mean " + m_attributeGrade->getName() }}, // shared y-axis is enabled, set only the last one {{1, colorDomain1} , {2, colorDomain2}} ); diff --git a/dialogs/driftanalysisdialog.cpp b/dialogs/driftanalysisdialog.cpp new file mode 100644 index 00000000..e9a1b499 --- /dev/null +++ b/dialogs/driftanalysisdialog.cpp @@ -0,0 +1,395 @@ +#include "driftanalysisdialog.h" +#include "ui_driftanalysisdialog.h" + +#include "domain/attribute.h" +#include "domain/application.h" +#include "domain/datafile.h" +#include "geostats/driftanalysis.h" +#include "geostats/quadratic3dtrendmodelfitting.h" +#include "widgets/linechartwidget.h" +#include "viewer3d/view3dcolortables.h" +#include "gslib/gslibparametersdialog.h" +#include "gslib/gslibparameterfiles/kt3dtrendmodelparameters.h" +#include "geometry/boundingbox.h" + +#include +#include +#include +#include +#include + +DriftAnalysisDialog::DriftAnalysisDialog(DataFile *dataFile, Attribute *attribute, QWidget *parent) : + QDialog(parent), + ui(new Ui::DriftAnalysisDialog), + m_dataFile(dataFile), + m_attribute(attribute), + m_lastNameForDriftVariable("drift_model") +{ + ui->setupUi(this); + + this->setWindowTitle("Drift Analysis"); + + ui->lblDataFile->setText("" + m_dataFile->getName() + ""); + ui->lblAttribute->setText("" + m_attribute->getName() + ""); + + m_chartWestEast = new LineChartWidget(); + m_chartWestEast ->setLegendVisible( false ); + m_chartSouthNorth = new LineChartWidget(); + m_chartSouthNorth->setLegendVisible( false ); + m_chartVertical = new LineChartWidget( nullptr, LineChartWidget::ZoomDirection::HORIZONTAL ); + m_chartVertical ->setLegendVisible( false ); + + ui->grpWEdrift->layout()->addWidget( m_chartWestEast ); + ui->grpSNdrift->layout()->addWidget( m_chartSouthNorth ); + ui->grpVerticalDrift->layout()->addWidget( m_chartVertical ); + + //if data set is not 3D, hide the vertical drift chart + if( ! m_dataFile->isTridimensional() ) + ui->grpVerticalDrift->hide(); + + ui->grpFitTrendModel->hide(); +} + +DriftAnalysisDialog::~DriftAnalysisDialog() +{ + delete ui; +} + +void DriftAnalysisDialog::performDriftAnalysis( DriftAnalysis& driftAnalysis, bool clear_curves ) +{ + if( ! driftAnalysis.run() ) { + + Application::instance()->logError("DriftAnalysisDialog::performDriftAnalysis(): failed execution:"); + Application::instance()->logError(" " + driftAnalysis.getLastError()); + QMessageBox::critical( this, "Error", "Drift analysis failed. Further details in the message panel." ); + + } else { //if the drift analysis completed successfully + + //shortcut for the STL's not-a-number value. + const double NaN = std::numeric_limits::quiet_NaN(); + + //get the drift analysis results + std::vector< std::pair< DriftAnalysis::coordX, DriftAnalysis::Mean > > + resultsWestEast = driftAnalysis.getResultsWestEast(); + std::vector< std::pair< DriftAnalysis::coordY, DriftAnalysis::Mean > > + resultsSouthNorth = driftAnalysis.getResultsSouthNorth(); + std::vector< std::pair< DriftAnalysis::coordZ, DriftAnalysis::Mean > > + resultsVertical = driftAnalysis.getResultsVertical(); + + //prepare data points for the three chart plotting + //outer vector: the series of multivariate data + //inner vectors: each multivariate datum (each element is a X, Y1, Y2,... value). + std::vector< std::vector > chartDataWestEast; + std::vector< std::vector > chartDataSouthNorth; + std::vector< std::vector > chartDataVertical; + + //traverse the results to fill the three chart data vectors + for( const std::pair< DriftAnalysis::coordX, DriftAnalysis::Mean >& result : resultsWestEast ) + chartDataWestEast .push_back( { result.first, result.second } ); + for( const std::pair< DriftAnalysis::coordY, DriftAnalysis::Mean >& result : resultsSouthNorth ) + chartDataSouthNorth.push_back( { result.first, result.second } ); + for( const std::pair< DriftAnalysis::coordZ, DriftAnalysis::Mean >& result : resultsVertical ) + chartDataVertical .push_back( { result.second, result.first } ); + + //set some properties of the domain categories relevant to make + //the chart informative + QColor colorWestEast = QColorConstants::Red; + QColor colorSouthNorth = QColorConstants::DarkGreen; + QColor colorVertical = QColorConstants::DarkBlue; + + QPen style; + if( ! clear_curves ) //user wants to display drift model along observed drift + style.setStyle( Qt::DashLine ); + + //display the results index 0: X values; 1: Y values + m_chartWestEast->setData( chartDataWestEast, 0, clear_curves, + {{}}, + {{1, "mean " + m_attribute->getName() }}, + {{1, colorWestEast}}, + {{1, style}} ); + m_chartWestEast->setXaxisCaption( "Easting" ); + m_chartSouthNorth->setData( chartDataSouthNorth, 0, clear_curves, + {{}}, + {{1, "mean " + m_attribute->getName() }}, + {{1, colorSouthNorth}}, + {{1, style}} ); + m_chartSouthNorth->setXaxisCaption( "Northing" ); + m_chartVertical->setData( chartDataVertical, 0, clear_curves, + {{}}, + {{1, "Z" }}, + {{1, colorVertical}}, + {{1, style}} ); + m_chartVertical->setXaxisCaption( "mean " + m_attribute->getName() ); + + } + +} + +void DriftAnalysisDialog::displayParamaters(const Quad3DTrendModelFittingAuxDefs::Parameters &model_parameters) +{ + //a lambda to compute the decimal logarithm of the number that returns 0.0 if value is 0.0 and inverts the sign + //of negative numbers, that is, does not return a NaN nor throws an exception. + auto lambdaQuietLog10 = [](double value){ + if( Util::almostEqual2sComplement( value, 0.0, 1 ) ) + return 0.0; + else + return std::log10( std::abs( value ) ); + }; + + //find the magnitudes of x, y and z of the dataset + BoundingBox bbox = m_dataFile->getBoundingBox(); + double magX = lambdaQuietLog10( std::max( std::abs( bbox.m_minX), std::abs( bbox.m_maxX ) ) ); + double magY = lambdaQuietLog10( std::max( std::abs( bbox.m_minY), std::abs( bbox.m_maxY ) ) ); + double magZ = lambdaQuietLog10( std::max( std::abs( bbox.m_minZ), std::abs( bbox.m_maxZ ) ) ); + + //find the magnitudes of each of the nine terms of the trend model + std::vector mag_term(9); + mag_term[0] = magX * 2 + lambdaQuietLog10( model_parameters[0] ); + mag_term[1] = magX + magY + lambdaQuietLog10( model_parameters[1] ); + mag_term[2] = magX + magZ + lambdaQuietLog10( model_parameters[2] ); + mag_term[3] = magY * 2 + lambdaQuietLog10( model_parameters[3] ); + mag_term[4] = magY + magZ + lambdaQuietLog10( model_parameters[4] ); + mag_term[5] = magZ * 2 + lambdaQuietLog10( model_parameters[5] ); + mag_term[6] = magX + lambdaQuietLog10( model_parameters[6] ); + mag_term[7] = magY + lambdaQuietLog10( model_parameters[7] ); + mag_term[8] = magZ + lambdaQuietLog10( model_parameters[8] ); + + //zero-out magnitudes whose parameters are zero so they are rendered + //with neutral color to convey its low impact on the trend model + for( uint i = 0; i < Quad3DTrendModelFittingAuxDefs::N_PARAMS; i++ ) + if( Util::almostEqual2sComplement( model_parameters[i], 0.0, 1) ) + mag_term[i] = 0.0; + + //find the highest and lowest magnitude values for color table + double absmax = -std::numeric_limits::max(); + for( int i = 0; i < Quad3DTrendModelFittingAuxDefs::N_PARAMS; i++) + if( std::abs(mag_term[i]) > absmax ) absmax = std::abs(mag_term[i]); + + //a lambda to automatize the generation of HTML to display a background color proportional to the magnitude of each model term + //if sign_value is negative, then value is converted to -value before retrieving the color form the color table. + auto lambdaMakeBGColor = [absmax](double value, double sign_value){ + if( sign_value >= 0 ) + return " bgcolor='" + Util::getHTMLColorFromValue( std::abs(value), ColorTable::SEISMIC, -absmax, absmax ) + "'"; + else + return " bgcolor='" + Util::getHTMLColorFromValue( -std::abs(value), ColorTable::SEISMIC, -absmax, absmax ) + "'"; + }; + + //a lambda to automatize the generation of HTML to render text in a contrasting color with respect to the background color + // scale_value: value to be used do determine color with respect to the color scale + // face_vale: value to be printed to screen + auto lambdaMakeFontColor = [absmax](double scale_value, double face_value){ + if( face_value >= 0 ) + return Util::fontColorTag( QString::number(face_value), Util::getColorFromValue( std::abs(scale_value), ColorTable::SEISMIC, -absmax, absmax ) ); + else + return Util::fontColorTag( QString::number(face_value), Util::getColorFromValue( -std::abs(scale_value), ColorTable::SEISMIC, -absmax, absmax ) ); + }; + + const QString btdSansColor = ""; + const QString btdAvecColor = "" + lambdaMakeFontColor(mag_term[0], model_parameters.a) + etd + + btdSansColor + "x2 + " + etd; + output += btdAvecColor + lambdaMakeBGColor( mag_term[1], model_parameters.b) + ">" + lambdaMakeFontColor(mag_term[1], model_parameters.b) + etd + + btdSansColor + "xy + " + etd; + output += btdAvecColor + lambdaMakeBGColor( mag_term[2], model_parameters.c) + ">" + lambdaMakeFontColor(mag_term[2], model_parameters.c) + etd + + btdSansColor + "xz + " + etd; + output += btdAvecColor + lambdaMakeBGColor( mag_term[3], model_parameters.d) + ">" + lambdaMakeFontColor(mag_term[3], model_parameters.d) + etd + + btdSansColor + "y2 + " + etd; + output += btdAvecColor + lambdaMakeBGColor( mag_term[4], model_parameters.e) + ">" + lambdaMakeFontColor(mag_term[4], model_parameters.e) + etd + + btdSansColor + "yz + " + etd; + output += btdAvecColor + lambdaMakeBGColor( mag_term[5], model_parameters.f) + ">" + lambdaMakeFontColor(mag_term[5], model_parameters.f) + etd + + btdSansColor + "z2 + " + etd; + output += btdAvecColor + lambdaMakeBGColor( mag_term[6], model_parameters.g) + ">" + lambdaMakeFontColor(mag_term[6], model_parameters.g) + etd + + btdSansColor + "x + " + etd; + output += btdAvecColor + lambdaMakeBGColor( mag_term[7], model_parameters.h) + ">" + lambdaMakeFontColor(mag_term[7], model_parameters.h) + etd + + btdSansColor + "y + " + etd; + output += btdAvecColor + lambdaMakeBGColor( mag_term[8], model_parameters.i) + ">" + lambdaMakeFontColor(mag_term[8], model_parameters.i) + etd + + btdSansColor + "z" + etd; + output += ""; + + ui->lblTrendModel->setText( output ); + +} + +void DriftAnalysisDialog::onRun() +{ + DriftAnalysis driftAnalysis; + driftAnalysis.setAttribute( m_attribute ); + driftAnalysis.setInputDataFile( m_dataFile ); + driftAnalysis.setNumberOfSteps( ui->spinNumberOfSteps->value() ); + + performDriftAnalysis( driftAnalysis ); +} + +void DriftAnalysisDialog::onFitTrendModel() +{ + //toggle the frame with the trend model fitting controls + ui->grpFitTrendModel->setVisible( ! ui->grpFitTrendModel->isVisible() ); +} + +void DriftAnalysisDialog::onRunFitTrendModel() +{ + //fit a trend model to the data + Quadratic3DTrendModelFitting q3dtmf( m_dataFile, m_attribute ); + Quad3DTrendModelFittingAuxDefs::Parameters modelParameters = + q3dtmf.processWithNonLinearLeastSquares(); + + //update the trend model disaply in this dialog + displayParamaters( modelParameters ); + + //keeps the last trend model fit for further use while another one is not fit + m_lastFitDriftModelParameters.reset( new Quad3DTrendModelFittingAuxDefs::Parameters( modelParameters ) ); +} + +void DriftAnalysisDialog::onSaveNewVariableWithDriftModel() +{ + if( ! m_lastFitDriftModelParameters ) { + Application::instance()->logError("DriftAnalysisDialog::onSaveNewVariableWithDriftModel(): null drift model. Run model fitting at least once.", true); + return; + } + + // make sure the data is loaded from filesystem + m_dataFile->loadData(); + + //present a variable naming dialog to the user + { + bool ok; + QString var_name = QInputDialog::getText(this, "Name the new variable", + "New variable with the trend model evaluated in data locations: ", + QLineEdit::Normal, m_lastNameForDriftVariable, &ok); + if( ! ok ) + return; + else + m_lastNameForDriftVariable = var_name; + } + + //allocate container for the drift values + std::vector drift_values( m_dataFile->getDataLineCount() ); + + //for each data sample + for( uint64_t iRow = 0; iRow < m_dataFile->getDataLineCount(); iRow++){ + + //get the spatial location of the current sample + double x, y, z; + m_dataFile->getDataSpatialLocation( iRow, x, y, z ); + + //evaluate the trend model at the current sample location + double drift_value = + m_lastFitDriftModelParameters->a * x * x + + m_lastFitDriftModelParameters->b * x * y + + m_lastFitDriftModelParameters->c * x * z + + m_lastFitDriftModelParameters->d * y * y + + m_lastFitDriftModelParameters->e * y * z + + m_lastFitDriftModelParameters->f * z * z + + m_lastFitDriftModelParameters->g * x + + m_lastFitDriftModelParameters->h * y + + m_lastFitDriftModelParameters->i * z ; + + //store the trend model value at the current sample location + drift_values[iRow] = drift_value; + } + + //add or update the trend model values in the data files + uint driftVatGEOEASindex = m_dataFile->getFieldGEOEASIndex( m_lastNameForDriftVariable ); + uint drift_col_index = 9999; + if( driftVatGEOEASindex == 0 ) { //the variable does not exist + //adds the drift values as a new attribute to the data set file + drift_col_index = m_dataFile->addNewDataColumn( m_lastNameForDriftVariable, drift_values ); + } else { //drift variable already exists in the file + drift_col_index = driftVatGEOEASindex - 1; + //update the data file with the recomputed drift model values + for( uint iRow = 0; iRow < m_dataFile->getDataLineCount(); iRow++ ) + m_dataFile->setData( iRow, drift_col_index, drift_values[iRow] ); + m_dataFile->writeToFS(); + } + + //gets the new Attribute object associated with the newly added variable with the trend model values + Attribute* driftVariable = m_dataFile->getAttributeFromGEOEASIndex( drift_col_index + 1 ); + + //sanity check + if( ! driftVariable ){ + Application::instance()->logError( "DriftAnalysisDialog::onSaveNewVariableWithDriftModel(): failed to retrieve the Attribute object" + " associated with the newly computed trend model values." ); + return; + } + + //run drift analysis on the drift itself so the use can assess the quality of the fit + Attribute* tmp = m_attribute; + { + m_attribute = driftVariable; + DriftAnalysis driftAnalysis; + driftAnalysis.setAttribute( driftVariable ); + driftAnalysis.setInputDataFile( m_dataFile ); + driftAnalysis.setNumberOfSteps( ui->spinNumberOfSteps->value() ); + performDriftAnalysis( driftAnalysis, false ); + } + m_attribute = tmp; +} + +void DriftAnalysisDialog::onCopyTrendModelAsCalcScript(){ + + if( ! m_lastFitDriftModelParameters ) { + Application::instance()->logError("DriftAnalysisDialog::onCopyTrendModelAsCalcScript(): null drift model. Run model fitting at least once.", true); + return; + } + + // assemble calculator script + QString script; + script += " [output variable] := (" + QString::number(m_lastFitDriftModelParameters->a) + ") * X_ * X_ + \n"; + script += " (" + QString::number(m_lastFitDriftModelParameters->b) + ") * X_ * Y_ + \n"; + script += " (" + QString::number(m_lastFitDriftModelParameters->c) + ") * X_ * Z_ + \n"; + script += " (" + QString::number(m_lastFitDriftModelParameters->d) + ") * Y_ * Y_ + \n"; + script += " (" + QString::number(m_lastFitDriftModelParameters->e) + ") * Y_ * Z_ + \n"; + script += " (" + QString::number(m_lastFitDriftModelParameters->f) + ") * Z_ * Z_ + \n"; + script += " (" + QString::number(m_lastFitDriftModelParameters->g) + ") * X_ + \n"; + script += " (" + QString::number(m_lastFitDriftModelParameters->h) + ") * Y_ + \n"; + script += " (" + QString::number(m_lastFitDriftModelParameters->i) + ") * Z_ ; \n"; + + // copy script to clipboard (CTRL+C) + QClipboard *clipboard = QGuiApplication::clipboard(); + QString originalText = clipboard->text(); + clipboard->setText(script); + QMessageBox::information(this, "Information", "Calculator script copied to the clipboard."); +} + +void DriftAnalysisDialog::onEditTrendModelParameters() +{ + if( ! m_lastFitDriftModelParameters ) { + Application::instance()->logError("DriftAnalysisDialog::onEditTrendModelParameters(): null drift model. Run model fitting at least once.", true); + return; + } + + // open a parameter editor dialog + Kt3dTrendModelParameters kt3dTrendModelParameters; + kt3dTrendModelParameters.setA( m_lastFitDriftModelParameters->a ); + kt3dTrendModelParameters.setB( m_lastFitDriftModelParameters->b ); + kt3dTrendModelParameters.setC( m_lastFitDriftModelParameters->c ); + kt3dTrendModelParameters.setD( m_lastFitDriftModelParameters->d ); + kt3dTrendModelParameters.setE( m_lastFitDriftModelParameters->e ); + kt3dTrendModelParameters.setF( m_lastFitDriftModelParameters->f ); + kt3dTrendModelParameters.setG( m_lastFitDriftModelParameters->g ); + kt3dTrendModelParameters.setH( m_lastFitDriftModelParameters->h ); + kt3dTrendModelParameters.setI( m_lastFitDriftModelParameters->i ); + GSLibParametersDialog gpd( &kt3dTrendModelParameters, this ); + gpd.setWindowTitle( "Edit kt3d's trend model parameters" ); + + //show the dialog modally an treat the user response + int ret = gpd.exec(); + if( ret != QDialog::Accepted ) + return; + + //if user didn't dismiss the dialog, assign the entered values to the model + m_lastFitDriftModelParameters->a = kt3dTrendModelParameters.getA(); + m_lastFitDriftModelParameters->b = kt3dTrendModelParameters.getB(); + m_lastFitDriftModelParameters->c = kt3dTrendModelParameters.getC(); + m_lastFitDriftModelParameters->d = kt3dTrendModelParameters.getD(); + m_lastFitDriftModelParameters->e = kt3dTrendModelParameters.getE(); + m_lastFitDriftModelParameters->f = kt3dTrendModelParameters.getF(); + m_lastFitDriftModelParameters->g = kt3dTrendModelParameters.getG(); + m_lastFitDriftModelParameters->h = kt3dTrendModelParameters.getH(); + m_lastFitDriftModelParameters->i = kt3dTrendModelParameters.getI(); + + //update the model display in this dialog + displayParamaters( *m_lastFitDriftModelParameters ); +} diff --git a/dialogs/driftanalysisdialog.h b/dialogs/driftanalysisdialog.h new file mode 100644 index 00000000..67b9ca67 --- /dev/null +++ b/dialogs/driftanalysisdialog.h @@ -0,0 +1,58 @@ +#ifndef DRIFTANALYSISDIALOG_H +#define DRIFTANALYSISDIALOG_H + +#include + +class DataFile; +class Attribute; +class LineChartWidget; + +class DriftAnalysis; + +namespace Quad3DTrendModelFittingAuxDefs{ + class Parameters; +} + +namespace Ui { + class DriftAnalysisDialog; +} + +class DriftAnalysisDialog : public QDialog +{ + Q_OBJECT + +public: + explicit DriftAnalysisDialog(DataFile* dataFile, Attribute* attribute, QWidget *parent = nullptr); + ~DriftAnalysisDialog(); + +private: + Ui::DriftAnalysisDialog *ui; + DataFile* m_dataFile; + Attribute* m_attribute; + LineChartWidget* m_chartWestEast; + LineChartWidget* m_chartSouthNorth; + LineChartWidget* m_chartVertical; + std::unique_ptr m_lastFitDriftModelParameters; + QString m_lastNameForDriftVariable; + + /** + * Runs the drift analysis algorithm passed as parameter and displays its results in the dialog's charts. + * @param clear_curves If true, clears all current curves, otherwise the resulting curves pile up. + */ + void performDriftAnalysis( DriftAnalysis& driftAnalysis, bool clear_curves = true ); + + /** + * Updates the label showing the trend model formula with the passed model parameters. + */ + void displayParamaters(const Quad3DTrendModelFittingAuxDefs::Parameters& model_parameters); + +private Q_SLOTS: + void onRun(); + void onFitTrendModel(); + void onRunFitTrendModel(); + void onSaveNewVariableWithDriftModel(); + void onCopyTrendModelAsCalcScript(); + void onEditTrendModelParameters(); +}; + +#endif // DRIFTANALYSISDIALOG_H diff --git a/dialogs/driftanalysisdialog.ui b/dialogs/driftanalysisdialog.ui new file mode 100644 index 00000000..2acad471 --- /dev/null +++ b/dialogs/driftanalysisdialog.ui @@ -0,0 +1,463 @@ + + + DriftAnalysisDialog + + + + 0 + 0 + 663 + 453 + + + + Dialog + + + + + + QLayout::SetDefaultConstraint + + + + + + 0 + 0 + + + + data file + + + + + + + + 0 + 0 + + + + variable + + + + + + + + 0 + 0 + + + + Drift analysis on: + + + + + + + + 0 + 0 + + + + variable: + + + + + + + number of steps: + + + + + + + 5 + + + 999 + + + 20 + + + + + + + + + + 0 + 0 + + + + Qt::Horizontal + + + + + 0 + 0 + + + + W-E drift + + + + 0 + + + 0 + + + 0 + + + 0 + + + 0 + + + + + + + 0 + 0 + + + + S-N drift + + + + 0 + + + 0 + + + 0 + + + 0 + + + 0 + + + + + + + 0 + 0 + + + + vertical drift + + + + 0 + + + 0 + + + 0 + + + 0 + + + 0 + + + + + + + + + + + Run + + + true + + + + + + + Fit trend model + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + Dismiss + + + + + + + + + + 0 + 0 + + + + Fit a trend model for kriging with a trend (kt3d) + + + + 6 + + + 1 + + + + + Save new variable with evaluation of drift model. + + + + + + + :/icons32/save32:/icons32/save32 + + + + + + + Copy trend model as calculator script. + + + + + + + :/icons32/calc32:/icons32/calc32 + + + + + + + fit a trend model + + + + + + + :/icons32/play32:/icons32/play32 + + + + + + + + 10 + + + + <html><head/><body><p>Trend model:x<sup>2</sup></p></body></html> + + + + + + + Qt::Horizontal + + + + 40 + 20 + + + + + + + + Edit trend model parameters + + + + + + + :/icons32/setting32:/icons32/setting32 + + + + + + + + + + + + + + btnDismiss + clicked() + DriftAnalysisDialog + reject() + + + 652 + 362 + + + 391 + 337 + + + + + btnRun + clicked() + DriftAnalysisDialog + onRun() + + + 84 + 362 + + + 115 + 399 + + + + + btnFitTrendModel + clicked() + DriftAnalysisDialog + onFitTrendModel() + + + 170 + 362 + + + 190 + 399 + + + + + btnRunFitTrendModel + clicked() + DriftAnalysisDialog + onRunFitTrendModel() + + + 46 + 416 + + + 467 + 397 + + + + + btnSaveDrift + clicked() + DriftAnalysisDialog + onSaveNewVariableWithDriftModel() + + + 80 + 416 + + + 549 + 509 + + + + + btnCopyModelAsCalcScript + clicked() + DriftAnalysisDialog + onCopyTrendModelAsCalcScript() + + + 148 + 416 + + + 131 + 449 + + + + + btnEditTrendModel + clicked() + DriftAnalysisDialog + onEditTrendModelParameters() + + + 102 + 411 + + + 171 + 450 + + + + + + onRun() + onFitTrendModel() + onRunFitTrendModel() + onSaveNewVariableWithDriftModel() + onCopyTrendModelAsCalcScript() + onEditTrendModelParameters() + + diff --git a/dialogs/gridrepositiondialog.cpp b/dialogs/gridrepositiondialog.cpp index a51a750b..d772ce71 100644 --- a/dialogs/gridrepositiondialog.cpp +++ b/dialogs/gridrepositiondialog.cpp @@ -22,9 +22,9 @@ GridRepositionDialog::GridRepositionDialog(CartesianGrid *cg, QWidget *parent) : ui->spinLLB_J->setValue( 0 ); ui->spinLLB_K->setValue( 0 ); - ui->dblSpinLLB_X->setRange( std::numeric_limits::min(), std::numeric_limits::max() ); - ui->dblSpinLLB_Y->setRange( std::numeric_limits::min(), std::numeric_limits::max() ); - ui->dblSpinLLB_Z->setRange( std::numeric_limits::min(), std::numeric_limits::max() ); + ui->dblSpinLLB_X->setRange( -std::numeric_limits::max(), std::numeric_limits::max() ); + ui->dblSpinLLB_Y->setRange( -std::numeric_limits::max(), std::numeric_limits::max() ); + ui->dblSpinLLB_Z->setRange( -std::numeric_limits::max(), std::numeric_limits::max() ); cg->getCellLocation( 0, 0, 0, x, y, z ); ui->dblSpinLLB_X->setValue( x ); ui->dblSpinLLB_Y->setValue( y ); @@ -37,9 +37,9 @@ GridRepositionDialog::GridRepositionDialog(CartesianGrid *cg, QWidget *parent) : ui->spinURT_J->setValue( cg->getNJ()-1 ); ui->spinURT_K->setValue( cg->getNK()-1 ); - ui->dblSpinURT_X->setRange( std::numeric_limits::min(), std::numeric_limits::max() ); - ui->dblSpinURT_Y->setRange( std::numeric_limits::min(), std::numeric_limits::max() ); - ui->dblSpinURT_Z->setRange( std::numeric_limits::min(), std::numeric_limits::max() ); + ui->dblSpinURT_X->setRange( -std::numeric_limits::max(), std::numeric_limits::max() ); + ui->dblSpinURT_Y->setRange( -std::numeric_limits::max(), std::numeric_limits::max() ); + ui->dblSpinURT_Z->setRange( -std::numeric_limits::max(), std::numeric_limits::max() ); cg->getCellLocation( cg->getNI()-1 , cg->getNJ()-1 , cg->getNK()-1 , x, y, z ); ui->dblSpinURT_X->setValue( x ); ui->dblSpinURT_Y->setValue( y ); diff --git a/dialogs/sisimdialog.cpp b/dialogs/sisimdialog.cpp index 66bdd5ac..1c50219c 100644 --- a/dialogs/sisimdialog.cpp +++ b/dialogs/sisimdialog.cpp @@ -442,8 +442,8 @@ void SisimDialog::onConfigureAndRun() data_max += fabs( data_max/100.0 ); //get min and max of secondary variable (for sisim_gs) - double secData_min = std::numeric_limits::max(); - double secData_max = std::numeric_limits::min(); + double secData_min = std::numeric_limits::max(); + double secData_max = -std::numeric_limits::max(); if( sisimProgram == "sisim_gs" ){ DataFile* secondaryDataFile = static_cast( m_SoftDataSetSelector->getSelectedFile() ); if( secondaryDataFile ){ diff --git a/docs/GammaRayManual.docx b/docs/GammaRayManual.docx index e488f616..55eecb07 100755 Binary files a/docs/GammaRayManual.docx and b/docs/GammaRayManual.docx differ diff --git a/domain/cartesiangrid.cpp b/domain/cartesiangrid.cpp index bb1fb0b5..05fd97e3 100644 --- a/domain/cartesiangrid.cpp +++ b/domain/cartesiangrid.cpp @@ -13,6 +13,7 @@ #include "geogrid.h" #include "domain/section.h" #include "domain/pointset.h" +#include "geometry/boundingbox.h" #include "spectral/spectral.h" //eigen third party library @@ -902,7 +903,7 @@ void CartesianGrid::reposition(uint llb_I, uint llb_J, uint llb_K, updateMetaDataFile(); } -double CartesianGrid::getDataSpatialLocation(uint line, CartesianCoord whichCoord) +double CartesianGrid::getDataSpatialLocation(uint line, CartesianCoord whichCoord) const { uint i, j, k; double x, y, z; @@ -916,14 +917,14 @@ double CartesianGrid::getDataSpatialLocation(uint line, CartesianCoord whichCoor } } -void CartesianGrid::getDataSpatialLocation(uint line, double &x, double &y, double &z) +void CartesianGrid::getDataSpatialLocation(uint line, double &x, double &y, double &z) const { uint i, j, k; indexToIJK( line, i, j, k ); IJKtoXYZ( i, j, k, x, y, z); } -bool CartesianGrid::isTridimensional() +bool CartesianGrid::isTridimensional() const { return m_nK > 1; } @@ -940,3 +941,14 @@ void CartesianGrid::probe(double pickedX, double pickedY, double pickedZ, Attrib else Application::instance()->logWarn("CartesianGrid::probe(): picked attribute not passed (check probe() caller code) or not displayed." ); } + +BoundingBox CartesianGrid::getBoundingBox() const +{ + double minX = _x0 - _dx / 2; + double minY = _y0 - _dy / 2; + double minZ = _z0 - _dz / 2; + double maxX = minX + _dx * m_nI; + double maxY = minY + _dy * m_nJ; + double maxZ = minZ + _dz * m_nK; + return BoundingBox(minX,minY,minZ,maxX,maxY,maxZ); +} diff --git a/domain/cartesiangrid.h b/domain/cartesiangrid.h index 0f29b46a..9030614f 100644 --- a/domain/cartesiangrid.h +++ b/domain/cartesiangrid.h @@ -133,10 +133,11 @@ class CartesianGrid : public GridFile, public IJAbstractCartesianGrid /** Cartesian grids never have declustering weights. At least they are not supposed to be. */ virtual Attribute* getVariableOfWeight( Attribute* /*at*/ ) { return nullptr; } virtual bool isRegular() { return true; } - virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ); - virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ); - virtual bool isTridimensional(); + virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ) const; + virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ) const; + virtual bool isTridimensional() const; virtual void probe( double pickedX, double pickedY, double pickedZ, Attribute* targetAttribute ); + virtual BoundingBox getBoundingBox( ) const; // File interface public: diff --git a/domain/datafile.cpp b/domain/datafile.cpp index 2d0450b9..69de090c 100644 --- a/domain/datafile.cpp +++ b/domain/datafile.cpp @@ -29,6 +29,7 @@ #include "algorithms/ialgorithmdatasource.h" #include "calculator/icalcproperty.h" #include "geogrid.h" +#include "geometry/boundingbox.h" /****************************** THE DATASOURCE INTERFACE TO THE ALGORITHM CLASSES * ****************************/ @@ -734,7 +735,7 @@ void DataFile::setCategorical(uint variableIndex, const CategoryDefinition *cd) void DataFile::probe(double pickedX, double pickedY, double pickedZ, Attribute *targetAttribute) { Application::instance()->logWarn( "DataFile::probe(): Default implementation called. Objects of type " - + getTypeName() + " must implement one according to their particular characteristics. Probe information is limited." ); + + getTypeName() + " must implement one according to their particular characteristics. Probe information is limited." ); } void DataFile::replacePhysicalFile(const QString from_file_path) diff --git a/domain/datafile.h b/domain/datafile.h index 0fc3359a..e1a1d810 100644 --- a/domain/datafile.h +++ b/domain/datafile.h @@ -14,6 +14,7 @@ class Attribute; class UnivariateCategoryClassification; class CategoryDefinition; class IAlgorithmDataSource; +class BoundingBox; enum class CartesianCoord : int { X, @@ -334,15 +335,15 @@ class DataFile : public File, public ICalcPropertyCollection /** * Returns one of the spatial coordinates (x, y or z) of the data value given its line number. */ - virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ) = 0; + virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ) const = 0; /** * Returns all the spatial coordinates (x, y and z as output parameters) of the data value given its line number. */ - virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ) = 0; + virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ) const = 0; /** Returns whether this data set is tridimensional. */ - virtual bool isTridimensional() = 0; + virtual bool isTridimensional() const = 0; /** Removes one data line from the internal data array. * It is necessary to call writeToFS() to commit the change to filesystem. @@ -446,6 +447,9 @@ class DataFile : public File, public ICalcPropertyCollection */ virtual bool isGridded() const = 0; + /** Returns the spatial extent of the object represented by this data file. */ + virtual BoundingBox getBoundingBox( ) const = 0; + //File interface virtual void deleteFromFS(); virtual void writeToFS(); diff --git a/domain/geogrid.cpp b/domain/geogrid.cpp index abceac7d..de2aa5d6 100644 --- a/domain/geogrid.cpp +++ b/domain/geogrid.cpp @@ -12,6 +12,7 @@ #include "domain/project.h" #include "geometry/vector3d.h" #include "geometry/face3d.h" +#include "geometry/boundingbox.h" #include "exceptions/invalidmethodexception.h" #include @@ -228,13 +229,13 @@ GeoGrid::GeoGrid(QString path, std::vector zones) : } -void GeoGrid::getBoundingBox(uint cellIndex, - double & minX, double & minY, double & minZ, - double & maxX, double & maxY, double & maxZ ) +void GeoGrid::getCellBoundingBox(uint cellIndex, + double & minX, double & minY, double & minZ, + double & maxX, double & maxY, double & maxZ ) const { //initialize the results to ensure the returned extrema are those of the cell. - minX = minY = minZ = std::numeric_limits::max(); - maxX = maxY = maxZ = std::numeric_limits::min(); + minX = minY = minZ = std::numeric_limits::max(); + maxX = maxY = maxZ = -std::numeric_limits::max(); //Get the cell CellDefRecordPtr cellDef = m_cellDefsPart.at( cellIndex ); //for each of the eight vertexes of the cell @@ -442,7 +443,7 @@ uint GeoGrid::getMeshNumberOfVertexes() return m_vertexesPart.size(); } -void GeoGrid::getMeshVertexLocation(uint index, double & x, double & y, double & z) +void GeoGrid::getMeshVertexLocation(uint index, double & x, double & y, double & z) const { x = m_vertexesPart[index]->X; y = m_vertexesPart[index]->Y; @@ -455,7 +456,7 @@ uint GeoGrid::getMeshNumberOfCells() return m_cellDefsPart.size(); } -void GeoGrid::getMeshCellDefinition(uint index, uint (&vIds)[8]) +void GeoGrid::getMeshCellDefinition(uint index, uint (&vIds)[8]) const { vIds[0] = m_cellDefsPart[index]->vId[0]; vIds[1] = m_cellDefsPart[index]->vId[1]; @@ -862,7 +863,7 @@ CartesianGrid *GeoGrid::getUnderlyingCartesianGrid() return nullptr; } -Hexahedron GeoGrid::makeHexahedron( uint cellIndex ) +Hexahedron GeoGrid::makeHexahedron( uint cellIndex ) const { Hexahedron hexa; uint vertexIndexes[8]; @@ -1186,7 +1187,7 @@ bool GeoGrid::XYZtoIJK( double x, double y, double z, uint& i, uint& j, uint& k return false; } -double GeoGrid::getDataSpatialLocation(uint line, CartesianCoord whichCoord) +double GeoGrid::getDataSpatialLocation(uint line, CartesianCoord whichCoord) const { uint i, j, k; double x, y, z; @@ -1200,7 +1201,7 @@ double GeoGrid::getDataSpatialLocation(uint line, CartesianCoord whichCoord) } } -void GeoGrid::getDataSpatialLocation(uint line, double &x, double &y, double &z) +void GeoGrid::getDataSpatialLocation(uint line, double &x, double &y, double &z) const { uint i, j, k; indexToIJK( line, i, j, k ); @@ -1233,6 +1234,38 @@ void GeoGrid::freeLoadedData() DataFile::freeLoadedData(); } +BoundingBox GeoGrid::getBoundingBox() const +{ + if( m_vertexesPart.empty() || m_cellDefsPart.empty() ){ + Application::instance()->logError("GeoGrid::getBoundingBox(): mesh data not loaded." + " Maybe a prior call to loadMesh() is missing. "); + return BoundingBox(); + } else { + double minX = std::numeric_limits::max(); + double minY = std::numeric_limits::max(); + double minZ = std::numeric_limits::max(); + double maxX = -std::numeric_limits::max(); + double maxY = -std::numeric_limits::max(); + double maxZ = -std::numeric_limits::max(); + //loop to update the limits of the bounding box + for( int iDataLine = 0; iDataLine < getDataLineCount(); ++iDataLine ){ + + double cellMinX, cellMinY, cellMinZ; + double cellMaxX, cellMaxY, cellMaxZ; + + getCellBoundingBox( iDataLine, cellMinX, cellMinY, cellMinZ, cellMaxX, cellMaxY, cellMaxZ ); + + if( cellMaxX > maxX ) maxX = cellMaxX; + if( cellMaxY > maxY ) maxY = cellMaxY; + if( cellMaxZ > maxZ ) maxZ = cellMaxZ; + if( cellMinX < minX ) minX = cellMinX; + if( cellMinY < minY ) minY = cellMinY; + if( cellMinZ < minZ ) minZ = cellMinZ; + } + return BoundingBox( minX, minY, minZ, maxX, maxY, maxZ ); + } +} + bool GeoGrid::canHaveMetaData() { return true; diff --git a/domain/geogrid.h b/domain/geogrid.h index 5326d985..96f2eaff 100644 --- a/domain/geogrid.h +++ b/domain/geogrid.h @@ -98,9 +98,9 @@ class GeoGrid : public GridFile /** * Returns (via output parameters) the bounding box of a cell given its cell index. */ - void getBoundingBox(uint cellIndex, - double& minX, double& minY, double& minZ, - double& maxX, double& maxY, double& maxZ ); + void getCellBoundingBox(uint cellIndex, + double& minX, double& minY, double& minZ, + double& maxX, double& maxY, double& maxZ ) const; /** Returns the path to the file that stores the geometry data, that is, * the contents of the m_vertexesPart and m_cellDefsPart members. @@ -132,13 +132,13 @@ class GeoGrid : public GridFile uint getMeshNumberOfVertexes(); /** Returns, via output parameters, the coordinates of a mesh vertex given its id (index). */ - void getMeshVertexLocation( uint index, double& x, double& y, double& z ); + void getMeshVertexLocation( uint index, double& x, double& y, double& z ) const; /** Returns the number of cells in the GeoGrid's mesh. */ uint getMeshNumberOfCells(); /** Returns, via output parameter, the indexes of mesh vertexes of a cell given its id (index). */ - void getMeshCellDefinition( uint index, uint (&vIds)[8] ); + void getMeshCellDefinition( uint index, uint (&vIds)[8] ) const; /** Returns a new point set by transforming the input point set from XYZ space to UVW space. * Returns null pointer if unfold fails for any reason. @@ -192,7 +192,7 @@ class GeoGrid : public GridFile /** * Creates a Hexahedron object from a cell's geometry data given the cell's index. */ - Hexahedron makeHexahedron( uint cellIndex ); + Hexahedron makeHexahedron( uint cellIndex ) const; /** * Exports this GeoGrid as ASCII Eclipse Grid (.grdecl format). @@ -220,13 +220,14 @@ class GeoGrid : public GridFile /** GeoGrids never have declustering weights. At least they are not supposed to have. */ virtual Attribute* getVariableOfWeight( Attribute* /*at*/ ) { return nullptr; } virtual bool isRegular() { return false; } - virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ); - virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ); + virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ) const; + virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ) const; /** GeoGrids are assumed to be always 3D. */ - virtual bool isTridimensional(){ return true; } + virtual bool isTridimensional() const { return true; } /** NOTE: override the default counting-only behavior of DataFile::getProportion(). */ virtual double getProportion(int variableIndex, double value0, double value1 ); virtual void freeLoadedData(); + virtual BoundingBox getBoundingBox( ) const; // File interface public: diff --git a/domain/pointset.cpp b/domain/pointset.cpp index 5143755e..fd9874a2 100644 --- a/domain/pointset.cpp +++ b/domain/pointset.cpp @@ -12,6 +12,7 @@ #include #include "viewer3d/view3dviewdata.h" #include "viewer3d/view3dbuilders.h" +#include "geometry/boundingbox.h" PointSet::PointSet( QString path ) : DataFile( path ) @@ -156,30 +157,30 @@ void PointSet::deleteVariable(uint columnToDelete) DataFile::deleteVariable( columnToDelete ); } -double PointSet::getDataSpatialLocation(uint line, CartesianCoord whichCoord) +double PointSet::getDataSpatialLocation(uint line, CartesianCoord whichCoord) const { switch ( whichCoord ) { - case CartesianCoord::X: return data( line, _x_field_index - 1 ); //x,y,z is in data file directly - case CartesianCoord::Y: return data( line, _y_field_index - 1 ); //x,y,z is in data file directly + case CartesianCoord::X: return dataConst( line, _x_field_index - 1 ); //x,y,z is in data file directly + case CartesianCoord::Y: return dataConst( line, _y_field_index - 1 ); //x,y,z is in data file directly case CartesianCoord::Z: if( isTridimensional() ) - return data( line, _z_field_index - 1 ); //x,y,z is in data file directly + return dataConst( line, _z_field_index - 1 ); //x,y,z is in data file directly else return 0.0; //returns z=0.0 for datasets in 2D. } } -void PointSet::getDataSpatialLocation(uint line, double &x, double &y, double &z) +void PointSet::getDataSpatialLocation(uint line, double &x, double &y, double &z) const { - x = data( line, _x_field_index - 1 ); //x,y,z is in data file directly - y = data( line, _y_field_index - 1 ); //x,y,z is in data file directly + x = dataConst( line, _x_field_index - 1 ); //x,y,z is in data file directly + y = dataConst( line, _y_field_index - 1 ); //x,y,z is in data file directly if( isTridimensional() ) - z = data( line, _z_field_index - 1 ); //x,y,z is in data file directly + z = dataConst( line, _z_field_index - 1 ); //x,y,z is in data file directly else z = 0.0; //returns z=0.0 for datasets in 2D. } -bool PointSet::isTridimensional() +bool PointSet::isTridimensional() const { return is3D(); } @@ -207,6 +208,37 @@ bool PointSet::getCenter(double &x, double &y, double &z) const } } +BoundingBox PointSet::getBoundingBox() const +{ + if( getDataLineCount() == 0){ + Application::instance()->logError("PointSet::getBoundingBox(): data not loaded." + " Maybe a prior call to readFromFS() is missing. "); + return BoundingBox(); + } else { + double minX = std::numeric_limits::max(); + double minY = std::numeric_limits::max(); + double minZ = std::numeric_limits::max(); + double maxX = -std::numeric_limits::max(); + double maxY = -std::numeric_limits::max(); + double maxZ = -std::numeric_limits::max(); + //loop to update the limits of the bounding box + for( int iDataLine = 0; iDataLine < getDataLineCount(); ++iDataLine ){ + double xCoord = dataConst( iDataLine, _x_field_index - 1 ) ; + double yCoord = dataConst( iDataLine, _y_field_index - 1 ) ; + double zCoord = 0.0; + if (is3D()) + zCoord = dataConst( iDataLine, _z_field_index - 1 ) ; + if( xCoord > maxX ) maxX = xCoord; + if( yCoord > maxY ) maxY = yCoord; + if( zCoord > maxZ ) maxZ = zCoord; + if( xCoord < minX ) minX = xCoord; + if( yCoord < minY ) minY = yCoord; + if( zCoord < minZ ) minZ = zCoord; + } + return BoundingBox( minX, minY, minZ, maxX, maxY, maxZ ); + } +} + void PointSet::setInfoFromMetadataFile() { QString md_file_path( this->_path ); diff --git a/domain/pointset.h b/domain/pointset.h index 842a2c93..d56182d5 100644 --- a/domain/pointset.h +++ b/domain/pointset.h @@ -93,11 +93,12 @@ class PointSet : public DataFile Attribute* getVariableOfWeight( Attribute* weight ); virtual void deleteVariable( uint columnToDelete ); virtual bool isRegular() { return false; } - virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ); - virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ); - virtual bool isTridimensional(); + virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ) const; + virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ) const; + virtual bool isTridimensional() const; virtual bool getCenter( double& x, double& y, double& z ) const; virtual bool isGridded() const override { return false; } + virtual BoundingBox getBoundingBox( ) const override; // File interface public: diff --git a/domain/segmentset.cpp b/domain/segmentset.cpp index ea17f049..7455f4f9 100644 --- a/domain/segmentset.cpp +++ b/domain/segmentset.cpp @@ -4,6 +4,7 @@ #include "domain/attribute.h" #include "domain/application.h" #include "domain/project.h" +#include "geometry/boundingbox.h" #include "util.h" #include #include @@ -512,25 +513,25 @@ File *SegmentSet::duplicatePhysicalFiles(const QString new_file_name) return newSS; } -double SegmentSet::getDataSpatialLocation(uint line, CartesianCoord whichCoord) +double SegmentSet::getDataSpatialLocation(uint line, CartesianCoord whichCoord) const { switch ( whichCoord ) { - case CartesianCoord::X: return ( data( line, _x_field_index - 1 ) + data( line, _x_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly - case CartesianCoord::Y: return ( data( line, _y_field_index - 1 ) + data( line, _y_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly + case CartesianCoord::X: return ( dataConst( line, _x_field_index - 1 ) + dataConst( line, _x_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly + case CartesianCoord::Y: return ( dataConst( line, _y_field_index - 1 ) + dataConst( line, _y_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly case CartesianCoord::Z: if( isTridimensional() ) - return ( data( line, _z_field_index - 1 ) + data( line, _z_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly + return ( dataConst( line, _z_field_index - 1 ) + dataConst( line, _z_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly else return 0.0; //returns z=0.0 for datasets in 2D. } } -void SegmentSet::getDataSpatialLocation(uint line, double &x, double &y, double &z) +void SegmentSet::getDataSpatialLocation(uint line, double &x, double &y, double &z) const { - x = ( data( line, _x_field_index - 1 ) + data( line, _x_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly - y = ( data( line, _y_field_index - 1 ) + data( line, _y_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly + x = ( dataConst( line, _x_field_index - 1 ) + dataConst( line, _x_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly + y = ( dataConst( line, _y_field_index - 1 ) + dataConst( line, _y_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly if( isTridimensional() ) - z = ( data( line, _z_field_index - 1 ) + data( line, _z_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly + z = ( dataConst( line, _z_field_index - 1 ) + dataConst( line, _z_final_field_index - 1 ) ) / 2.0; //x,y,z is in data file directly else z = 0.0; //returns z=0.0 for datasets in 2D. } @@ -581,6 +582,44 @@ bool SegmentSet::getCenter(double &x, double &y, double &z) const } } +BoundingBox SegmentSet::getBoundingBox() const +{ + if( getDataLineCount() == 0){ + Application::instance()->logError("SegmentSet::getBoundingBox(): data not loaded." + " Maybe a prior call to readFromFS() is missing. "); + return BoundingBox(); + } else { + double minX = std::numeric_limits::max(); + double minY = std::numeric_limits::max(); + double minZ = std::numeric_limits::max(); + double maxX = -std::numeric_limits::max(); + double maxY = -std::numeric_limits::max(); + double maxZ = -std::numeric_limits::max(); + //loop to update the limits of the bounding box + for( int iDataLine = 0; iDataLine < getDataLineCount(); ++iDataLine ){ + double xCoord = dataConst( iDataLine, _x_field_index - 1 ) ; + double yCoord = dataConst( iDataLine, _y_field_index - 1 ) ; + double zCoord = dataConst( iDataLine, _z_field_index - 1 ) ; + if( xCoord > maxX ) maxX = xCoord; + if( yCoord > maxY ) maxY = yCoord; + if( zCoord > maxZ ) maxZ = zCoord; + if( xCoord < minX ) minX = xCoord; + if( yCoord < minY ) minY = yCoord; + if( zCoord < minZ ) minZ = zCoord; + double xCoord2 = dataConst( iDataLine, _x_final_field_index - 1 ) ; + double yCoord2 = dataConst( iDataLine, _y_final_field_index - 1 ) ; + double zCoord2 = dataConst( iDataLine, _z_final_field_index - 1 ) ; + if( xCoord2 > maxX ) maxX = xCoord2; + if( yCoord2 > maxY ) maxY = yCoord2; + if( zCoord2 > maxZ ) maxZ = zCoord2; + if( xCoord2 < minX ) minX = xCoord2; + if( yCoord2 < minY ) minY = yCoord2; + if( zCoord2 < minZ ) minZ = zCoord2; + } + return BoundingBox( minX, minY, minZ, maxX, maxY, maxZ ); + } +} + void SegmentSet::setInfoFromOtherPointSet(PointSet *otherPS) { Q_UNUSED( otherPS ); diff --git a/domain/segmentset.h b/domain/segmentset.h index 5e3ff027..367ae7ab 100644 --- a/domain/segmentset.h +++ b/domain/segmentset.h @@ -144,11 +144,12 @@ class SegmentSet : public PointSet // DataFile interface public: /** NOTE: this returns the middle point of the segment. */ - virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ); - virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ); + virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ) const; + virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ) const; /** NOTE: override the default counting-only behavior of DataFile::getProportion(). */ virtual double getProportion(int variableIndex, double value0, double value1 ); virtual bool getCenter( double& x, double& y, double& z ) const; + virtual BoundingBox getBoundingBox( ) const; // PointSet interface public: diff --git a/domain/verticalproportioncurve.cpp b/domain/verticalproportioncurve.cpp index 15605840..8e8f0af9 100644 --- a/domain/verticalproportioncurve.cpp +++ b/domain/verticalproportioncurve.cpp @@ -8,6 +8,7 @@ #include "domain/project.h" #include "domain/objectgroup.h" #include "domain/attribute.h" +#include "geometry/boundingbox.h" VerticalProportionCurve::VerticalProportionCurve(QString path, QString associatedCategoryDefinitionName) : DataFile( path ), @@ -385,14 +386,14 @@ bool VerticalProportionCurve::isRegular() { " a spatial object." ); } -double VerticalProportionCurve::getDataSpatialLocation(uint line, CartesianCoord whichCoord) +double VerticalProportionCurve::getDataSpatialLocation(uint line, CartesianCoord whichCoord) const { Q_UNUSED( line ) Q_UNUSED( whichCoord ) assert( false && "VerticalProportionCurve::getDataSpatialLocation(): a VerticalProportionCurve is not a spatial object." ); } -void VerticalProportionCurve::getDataSpatialLocation(uint line, double &x, double &y, double &z) +void VerticalProportionCurve::getDataSpatialLocation(uint line, double &x, double &y, double &z) const { Q_UNUSED( line ) Q_UNUSED( x ) @@ -401,7 +402,7 @@ void VerticalProportionCurve::getDataSpatialLocation(uint line, double &x, doubl assert( false && "VerticalProportionCurve::getDataSpatialLocation(): a VerticalProportionCurve is not a spatial object." ); } -bool VerticalProportionCurve::isTridimensional() +bool VerticalProportionCurve::isTridimensional() const { assert( false && "VerticalProportionCurve::isTridimensional(): a VerticalProportionCurve is not a spatial object." ); } @@ -457,6 +458,11 @@ void VerticalProportionCurve::readFromFS() _data.clear(); } +BoundingBox VerticalProportionCurve::getBoundingBox() const +{ + assert( false && "VerticalProportionCurve::getBoundingBox(): a VerticalProportionCurve is not a spatial object." ); +} + void VerticalProportionCurve::getSpatialAndTopologicalCoordinates(int iRecord, double &x, double &y, double &z, int &i, int &j, int &k) { assert( false && "VerticalProportionCurve::getSpatialAndTopologicalCoordinates(): a VerticalProportionCurve is not a spatial object." ); diff --git a/domain/verticalproportioncurve.h b/domain/verticalproportioncurve.h index f85d8962..8c9b25ab 100644 --- a/domain/verticalproportioncurve.h +++ b/domain/verticalproportioncurve.h @@ -145,12 +145,13 @@ class VerticalProportionCurve : public DataFile Attribute* getVariableOfWeight( Attribute* weight ); virtual void deleteVariable( uint columnToDelete ); virtual bool isRegular(); - virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ); - virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ); - virtual bool isTridimensional(); + virtual double getDataSpatialLocation( uint line, CartesianCoord whichCoord ) const; + virtual void getDataSpatialLocation( uint line, double& x, double& y, double& z ) const; + virtual bool isTridimensional() const; virtual void writeToFS(); virtual void readFromFS(); virtual bool isGridded() const override { return false; } + virtual BoundingBox getBoundingBox( ) const override; // ICalcPropertyCollection interface public: diff --git a/geometry/boundingbox.cpp b/geometry/boundingbox.cpp new file mode 100644 index 00000000..57e632b8 --- /dev/null +++ b/geometry/boundingbox.cpp @@ -0,0 +1,14 @@ +#include "boundingbox.h" + +BoundingBox::BoundingBox() +{} + +BoundingBox::BoundingBox(double minX, double minY, double minZ, + double maxX, double maxY, double maxZ) : + m_minX(minX), + m_minY(minY), + m_minZ(minZ), + m_maxX(maxX), + m_maxY(maxY), + m_maxZ(maxZ) +{} diff --git a/geometry/boundingbox.h b/geometry/boundingbox.h new file mode 100644 index 00000000..72dab177 --- /dev/null +++ b/geometry/boundingbox.h @@ -0,0 +1,27 @@ +#ifndef BOUNDINGBOX_H +#define BOUNDINGBOX_H + + +/** + * The BoundingBox class represents the spatial extent of some geometric object. + */ +class BoundingBox +{ +public: + BoundingBox(); + + BoundingBox(double minX, double minY, double minZ, + double maxX, double maxY, double maxZ); + + double getXsize() const { return m_maxX - m_minX; } + double getYsize() const { return m_maxY - m_minY; } + double getZsize() const { return m_maxZ - m_minZ; } + + double getCenterX() const { return (m_maxX+m_minX)/2; } + double getCenterY() const { return (m_maxY+m_minY)/2; } + double getCenterZ() const { return (m_maxZ+m_minZ)/2; } + + double m_minX, m_minY, m_minZ, m_maxX, m_maxY, m_maxZ; +}; + +#endif // BOUNDINGBOX_H diff --git a/geostats/contactanalysis.cpp b/geostats/contactanalysis.cpp index d3406bd5..897d8393 100644 --- a/geostats/contactanalysis.cpp +++ b/geostats/contactanalysis.cpp @@ -374,7 +374,7 @@ bool ContactAnalysis::run() double maxBBdY = -std::numeric_limits::max(); uint maxIndex = gg->getDataLineCount(); for( uint i = 0; i < maxIndex; i++ ){ - gg->getBoundingBox( i, bbMinX, bbMinY, bbMinZ, bbMaxX, bbMaxY, bbMaxZ ); + gg->getCellBoundingBox( i, bbMinX, bbMinY, bbMinZ, bbMaxX, bbMaxY, bbMaxZ ); double bbDx = bbMaxX - bbMinX; double bbDy = bbMaxY - bbMinY; if( bbDx > maxBBdX ) maxBBdX = bbDx; diff --git a/geostats/driftanalysis.cpp b/geostats/driftanalysis.cpp new file mode 100644 index 00000000..093d134b --- /dev/null +++ b/geostats/driftanalysis.cpp @@ -0,0 +1,283 @@ +#include "driftanalysis.h" + +#include "domain/attribute.h" +#include "domain/cartesiangrid.h" +#include "domain/datafile.h" +#include "domain/geogrid.h" +#include "domain/pointset.h" +#include "domain/segmentset.h" +#include "geometry/boundingbox.h" +#include "geostats/searchbox.h" +#include "geostats/searchstrategy.h" +#include "spatialindex/spatialindex.h" + +#include +#include + +DriftAnalysis::DriftAnalysis() : + m_inputDataFile(nullptr), + m_attribute(nullptr), + m_NumberOfSteps(2), + m_lastError(""), + m_inputDataType( DataSetType::UNDEFINED ) +{} + +// getters and setters +DataFile *DriftAnalysis::getInputDataFile() const{return m_inputDataFile;} +void DriftAnalysis::setInputDataFile(DataFile *inputDataFile){ m_inputDataFile = inputDataFile;} +Attribute *DriftAnalysis::getAttribute() const{ return m_attribute;} +void DriftAnalysis::setAttribute(Attribute *attribute){ m_attribute = attribute;} +uint16_t DriftAnalysis::getNumberOfSteps() const{ return m_NumberOfSteps;} +void DriftAnalysis::setNumberOfSteps(const uint16_t &numberOfSteps){ m_NumberOfSteps = numberOfSteps;} +QString DriftAnalysis::getLastError() const{ return m_lastError;} +std::vector > DriftAnalysis::getResultsWestEast() const{ return m_resultsWestEast;} +std::vector > DriftAnalysis::getResultsSouthNorth() const{ return m_resultsSouthNorth;} +std::vector > DriftAnalysis::getResultsVertical() const{ return m_resultsVertical;} + +bool DriftAnalysis::run() +{ + //assuming the algorithm will finish normally + m_lastError = ""; + + //check whether everything is ok to run + if( ! isOKtoRun() ) + return false; + + //make sure the results vectors are empty + m_resultsWestEast.clear(); + m_resultsSouthNorth.clear(); + m_resultsVertical.clear(); + + //load the input data file + m_inputDataFile->loadData(); + + //get the data file column index corresponding to the attribute whose drift is to be analysed + uint indexOfVariable = m_inputDataFile->getFieldGEOEASIndex( m_attribute->getName() ) - 1; + if( indexOfVariable > m_inputDataFile->getDataColumnCount() ){ + m_lastError = "Error getting the data file column index of the categorical attribute with the domains."; + return false; + } + + //determine the type of the input data set once (avoid repetitive calls to the slow File::getFileType()) + //define a cell object that represents the current simulation cell. Also construct a data cell concrete object + //depending on the input data set type. + if( m_inputDataFile->getFileType() == "POINTSET" ) + m_inputDataType = DataSetType::POINTSET; + else if( m_inputDataFile->getFileType() == "CARTESIANGRID" ) + m_inputDataType = DataSetType::CARTESIANGRID; + else if( m_inputDataFile->getFileType() == "GEOGRID" ) + m_inputDataType = DataSetType::GEOGRID; + else if( m_inputDataFile->getFileType() == "SEGMENTSET" ) + m_inputDataType = DataSetType::SEGMENTSET; + else { + m_inputDataType = DataSetType::UNDEFINED; + m_lastError = "Input data type not currently supported:" + m_inputDataFile->getFileType(); + return false; + } + + //build the spatial index + SpatialIndex spatialIndex; + { + ////////////////////////////////// + QProgressDialog progressDialog; + progressDialog.show(); + progressDialog.setLabelText("Building spatial index..."); + progressDialog.setMinimum(0); + progressDialog.setValue(0); + progressDialog.setMaximum( 0 ); + QApplication::processEvents(); + ////////////////////////////////// + //the spatial index is filled differently depending on the type of the input data set + if( m_inputDataFile->getFileType() == "POINTSET" ){ + PointSet* psAspect = dynamic_cast( m_inputDataFile ); + spatialIndex.fill( psAspect, 0.1 ); + } else if ( m_inputDataFile->getFileType() == "SEGMENTSET") { + SegmentSet* ssAspect = dynamic_cast( m_inputDataFile ); + spatialIndex.fill( ssAspect, 0.1 ); //use cell size as tolerance + } else if ( m_inputDataFile->getFileType() == "CARTESIANGRID") { + CartesianGrid* cgAspect = dynamic_cast( m_inputDataFile ); + spatialIndex.fill( cgAspect ); + } else if ( m_inputDataFile->getFileType() == "GEOGRID") { + GeoGrid* ggAspect = dynamic_cast( m_inputDataFile ); + spatialIndex.fillWithCenters( ggAspect, 0.0001 ); + } else { + m_lastError = "Internal error building spatial index: input data of type " + m_inputDataFile->getFileType() + " are not currently supported."; + return false; + } + } + + //determine the total number of processing steps (#of lags X #of samples) + int total_steps = m_NumberOfSteps * ( m_inputDataFile->isTridimensional() ? 3 : 2 ); + + //configure and display a progress bar for the simulation task + ////////////////////////////////// + QProgressDialog progressDialog; + progressDialog.show(); + progressDialog.setLabelText("Running contact analysis..."); + progressDialog.setMinimum(0); + progressDialog.setValue(0); + progressDialog.setMaximum( total_steps ); + ///////////////////////////////// + + //this is to keep track of processing progress to update the progress bar + int total_done_so_far = 0; + + //get the spatial extent of the data file + BoundingBox bbox = m_inputDataFile->getBoundingBox(); + + //get the sizes of the slices the three axes + double sliceSizeWE = bbox.getXsize() / m_NumberOfSteps; + double sliceSizeSN = bbox.getYsize() / m_NumberOfSteps; + double sliceSizeVertical = bbox.getZsize() / m_NumberOfSteps; + + //get the input data file's NDV as floating point number + double dummyValue = m_inputDataFile->getNoDataValueAsDouble(); + bool hasNDV = m_inputDataFile->hasNoDataValue(); + + //compute drift along West->East diretcion (x axis). + for( int iSlice = 0; iSlice < m_NumberOfSteps; iSlice++, total_done_so_far++ ){ + + //define the X interval of the current X slice + double sliceMin = bbox.m_minX + iSlice * sliceSizeWE; + double sliceMax = sliceMin + sliceSizeWE; + + //build a bounding box representing the current X slice + BoundingBox bbox( sliceMin, -std::numeric_limits::max(), -std::numeric_limits::max(), + sliceMax, std::numeric_limits::max(), std::numeric_limits::max() ); + + //do a spatial search, fetching the indexes of the samples lying within the current X slice + QList samplesIndexes = spatialIndex.getWithinBoundingBox( bbox ); + + //for each sample found inside the current X slice + double sum = 0.0; + uint64_t count = 0; + for( uint sampleIndex : samplesIndexes ){ + double sampleValue = m_inputDataFile->dataConst( sampleIndex, indexOfVariable ); + if( ! hasNDV || ! Util::almostEqual2sComplement( sampleValue, dummyValue, 1 ) ){ + sum += sampleValue; + count++; + } + } + + //compute the mean of the attribute within the current X slice + double mean = std::numeric_limits::quiet_NaN(); + if( count > 0 ) + mean = sum / count; + + //store the result for the current X slice + m_resultsWestEast.push_back( { bbox.getCenterX(), mean } ); + + progressDialog.setValue( total_done_so_far ); //update the progress bar + QApplication::processEvents(); //let Qt repaint its widgets + } // for each X slice + + //compute drift along South->North diretcion (y axis). + for( int iSlice = 0; iSlice < m_NumberOfSteps; iSlice++, total_done_so_far++ ){ + + //define the Y interval of the current Y slice + double sliceMin = bbox.m_minY + iSlice * sliceSizeSN; + double sliceMax = sliceMin + sliceSizeSN; + + //build a bounding box representing the current Y slice + BoundingBox bbox( -std::numeric_limits::max(), sliceMin, -std::numeric_limits::max(), + std::numeric_limits::max(), sliceMax, std::numeric_limits::max() ); + + //do a spatial search, fetching the indexes of the samples lying within the current Y slice + QList samplesIndexes = spatialIndex.getWithinBoundingBox( bbox ); + + //for each sample found inside the current Y slice + double sum = 0.0; + uint64_t count = 0; + for( uint sampleIndex : samplesIndexes ){ + double sampleValue = m_inputDataFile->dataConst( sampleIndex, indexOfVariable ); + if( ! Util::almostEqual2sComplement( sampleValue, dummyValue, 1 ) ){ + sum += sampleValue; + count++; + } + } + + //compute the mean of the attribute within the current Y slice + double mean = std::numeric_limits::quiet_NaN(); + if( count > 0 ) + mean = sum / count; + + //store the result for the current Y slice + m_resultsSouthNorth.push_back( { bbox.getCenterY(), mean } ); + + progressDialog.setValue( total_done_so_far ); //update the progress bar + QApplication::processEvents(); //let Qt repaint its widgets + } // for each Y slice + + //compute drift along vertical diretcion (z axis). + if( m_inputDataFile->isTridimensional() ) { + for( int iSlice = 0; iSlice < m_NumberOfSteps; iSlice++, total_done_so_far++ ){ + + //define the Z interval of the current Z slice + double sliceMin = bbox.m_minZ + iSlice * sliceSizeVertical; + double sliceMax = sliceMin + sliceSizeVertical; + + //build a bounding box representing the current Z slice + BoundingBox bbox( -std::numeric_limits::max(), -std::numeric_limits::max(), sliceMin, + std::numeric_limits::max(), std::numeric_limits::max(), sliceMax ); + + //do a spatial search, fetching the indexes of the samples lying within the current Z slice + QList samplesIndexes = spatialIndex.getWithinBoundingBox( bbox ); + + //for each sample found inside the current Z slice + double sum = 0.0; + uint64_t count = 0; + for( uint sampleIndex : samplesIndexes ){ + double sampleValue = m_inputDataFile->dataConst( sampleIndex, indexOfVariable ); + if( ! Util::almostEqual2sComplement( sampleValue, dummyValue, 1 ) ){ + sum += sampleValue; + count++; + } + } + + //compute the mean of the attribute within the current Z slice + double mean = std::numeric_limits::quiet_NaN(); + if( count > 0 ) + mean = sum / count; + + //store the result for the current Z slice + m_resultsVertical.push_back( { bbox.getCenterZ(), mean } ); + + progressDialog.setValue( total_done_so_far ); //update the progress bar + QApplication::processEvents(); //let Qt repaint its widgets + } // for each Z slice + } // if data file is 3D + + + return true; +} + +bool DriftAnalysis::isOKtoRun() +{ + //assuming everything is ok + m_lastError = ""; + + if( ! m_attribute ){ + m_lastError = "Attribute not provided."; + return false; + } + + if( ! m_inputDataFile ){ + m_lastError = "Input data file not provided."; + return false; + } else { + DataFile* dataFile = dynamic_cast( m_attribute->getContainingFile() ); + if( ! dataFile ) { + m_lastError = "Make sure the attributes have a valid parent data file."; + return false; + } + } + + if( m_NumberOfSteps < 2 ){ + m_lastError = "Number of steps must be 2 or more."; + return false; + } + + return true; +} + + diff --git a/geostats/driftanalysis.h b/geostats/driftanalysis.h new file mode 100644 index 00000000..fa2fba9c --- /dev/null +++ b/geostats/driftanalysis.h @@ -0,0 +1,85 @@ +#ifndef DRIFTANALYSIS_H +#define DRIFTANALYSIS_H + +#include "util.h" + +#include + +#include + +class DataFile; +class Attribute; + +class DriftAnalysis +{ +public: + + typedef double coordX; + typedef double coordY; + typedef double coordZ; + typedef double Mean; + + DriftAnalysis(); + + //@{ + /** The input dataset. */ + DataFile *getInputDataFile() const; + void setInputDataFile(DataFile *inputDataFile); + //@} + + //@{ + /** The input variable. */ + Attribute *getAttribute() const; + void setAttribute(Attribute *attribute); + //@} + + //@{ + /** The number of steps in which divide the south-north, west-east and vertical extents. + * The greater this number, the greater the resolution of the output data. + */ + uint16_t getNumberOfSteps() const; + void setNumberOfSteps(const uint16_t &numberOfSteps); + //@} + + /** The input dataset. */ + QString getLastError() const; + + //@{ + /** The results of the drift analysis. */ + std::vector< std::pair< DriftAnalysis::coordX, DriftAnalysis::Mean > > getResultsWestEast() const; + std::vector< std::pair< DriftAnalysis::coordY, DriftAnalysis::Mean > > getResultsSouthNorth() const; + std::vector< std::pair< DriftAnalysis::coordZ, DriftAnalysis::Mean > > getResultsVertical() const; + //@} + + /** Executes the algorithm. isOKtoRun() is called to check whether all parameters are set properly. + * @return True if the algorithm finishes without issues (check getLastError() otherwise). + */ + bool run(); + +private: + + /** + * Verifies whether everything is correct to execute the algorithm. + * If any, the reason for last check failure must be checked with getLastError(). + * This is called by run(). + * @return True if the algorithm is all set to start. + */ + bool isOKtoRun(); + + //-------------agorithm parameters------------- + DataFile* m_inputDataFile; + Attribute* m_attribute; + uint16_t m_NumberOfSteps; + //--------------------------------------------- + + QString m_lastError; + DataSetType m_inputDataType; + + //stores the drift analysis results, that is, pairs of X/Y/Z values and mean of + //the input variable in the three directions (vertical will be empty for 2D datasets). + std::vector< std::pair< DriftAnalysis::coordX, DriftAnalysis::Mean > > m_resultsWestEast; + std::vector< std::pair< DriftAnalysis::coordY, DriftAnalysis::Mean > > m_resultsSouthNorth; + std::vector< std::pair< DriftAnalysis::coordZ, DriftAnalysis::Mean > > m_resultsVertical; +}; + +#endif // DRIFTANALYSIS_H diff --git a/geostats/quadratic3dtrendmodelfitting.cpp b/geostats/quadratic3dtrendmodelfitting.cpp new file mode 100644 index 00000000..594a9a39 --- /dev/null +++ b/geostats/quadratic3dtrendmodelfitting.cpp @@ -0,0 +1,777 @@ +#include "quadratic3dtrendmodelfitting.h" + +#include "util.h" +#include "mainwindow.h" +#include "dialogs/emptydialog.h" +#include "domain/application.h" +#include "domain/attribute.h" +#include "domain/datafile.h" +#include "geometry/boundingbox.h" + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +std::vector< double > Quadratic3DTrendModelFitting::s_objectiveFunctionValues; + +////////////////////////////////////////CLASS FOR THE GENETIC ALGORITHM////////////////////////////////////////// + +class Individual{ +public: + //constructors + Individual() : + m_parameters(), + m_fValue( std::numeric_limits::max() ) + {} + Individual( const Quad3DTrendModelFittingAuxDefs::Parameters& pparameters ) : + m_parameters( pparameters ), + m_fValue( std::numeric_limits::max() ) + {} + Individual( const Individual& otherIndividual ) : + m_parameters( otherIndividual.m_parameters ), + m_fValue( otherIndividual.m_fValue ) + {} + + //genetic operators + std::pair crossOver( const Individual& otherIndividual, + int pointOfCrossOver ) const { + Individual child1, child2; + for( int iParameter = 0; iParameter < Quad3DTrendModelFittingAuxDefs::N_PARAMS; ++iParameter ){ + if( iParameter < pointOfCrossOver ){ + child1.m_parameters[iParameter] = m_parameters[iParameter]; + child2.m_parameters[iParameter] = otherIndividual.m_parameters[iParameter]; + } else { + child1.m_parameters[iParameter] = otherIndividual.m_parameters[iParameter]; + child2.m_parameters[iParameter] = m_parameters[iParameter]; + } + } + return { child1, child2 }; + } + void mutate( double mutationRate, + const Quad3DTrendModelFittingAuxDefs::Parameters& lowBoundaries, + const Quad3DTrendModelFittingAuxDefs::Parameters& highBoundaries ){ + //compute the mutation probability for a single gene (parameter) + double probOfMutation = 1.0 / Quad3DTrendModelFittingAuxDefs::N_PARAMS * mutationRate; + //traverse all genes (parameters) + for( int iPar = 0; iPar < Quad3DTrendModelFittingAuxDefs::N_PARAMS; ++iPar ){ + //draw a value between 0.0 and 1.0 from an uniform distribution + double p = std::rand() / (double)RAND_MAX; + //if a mutation is due... + if( p < probOfMutation ) { + //perform mutation by randomly sorting a value within the domain. + double LO = lowBoundaries[iPar]; + double HI = highBoundaries[iPar]; + m_parameters[iPar] = LO + std::rand() / (RAND_MAX/(HI-LO)); + } + } + } + + //attribution operator + Individual& operator=( const Individual& otherIndividual ){ + m_parameters = otherIndividual.m_parameters; + m_fValue = otherIndividual.m_fValue; + return *this; + } + + //comparison operator + bool operator<( const Individual& otherIndividual ) const { + return m_fValue < otherIndividual.m_fValue; + } + + //member variables + Quad3DTrendModelFittingAuxDefs::Parameters m_parameters; //the genes (model parameters) of this individual + double m_fValue; //objective function value corresponding to this individual +}; +typedef Individual Solution; //make a synonym just for code readbility +///////////////////////////////////////////////////////////////////////////////////////////////////////// + +/** + * The code for multithreaded evaluation of the objective function for a range of individuals (a set of variogam parameters) + * in the Genetic Algorithm. + * @param quad3DTrendModelFitRef Reference to the Quadratic3DTrendModelFitting object so its objective function can be called. + * @param iIndividual_initial First individual index to process. + * @param iIndividual_final Last individual index to process. If final_i==initial_i, this thread processes only one individual. + * OUTPUT PARAMETER: + * @param population The set of individuals to receive the evaluation of the objective function. + */ +void taskEvaluateObjetiveInRangeOfIndividualsForGenetic( + const Quadratic3DTrendModelFitting& quad3DTrendModelFitRef, + int iIndividual_initial, + int iIndividual_final, + std::vector< Individual >& population //-->output parameter + ) { + for( int iInd = iIndividual_initial; iInd <= iIndividual_final; ++iInd ){ + Individual& ind = population[iInd]; + ind.m_fValue = quad3DTrendModelFitRef.objectiveFunction( ind.m_parameters ); + } +} + + +Quadratic3DTrendModelFitting::Quadratic3DTrendModelFitting( DataFile* dataFile, Attribute* attribute ) : + m_attribute(attribute), m_dataFile(dataFile) +{} + +double Quadratic3DTrendModelFitting::objectiveFunction( const Quad3DTrendModelFittingAuxDefs::Parameters& parameters ) const +{ + //the sum to compute the error mean to be returned + double sum = 0.0; + + //get the no-data value in numeric form to improve performance + double NDV = m_dataFile->getNoDataValueAsDouble(); + bool hasNDV = m_dataFile->hasNoDataValue(); + + //if data is 2D, then the z-bearing terms are invariant. + constexpr int IS_3D = 1; + constexpr int IS_2D = 0; + int is3D = m_dataFile->isTridimensional() ? IS_3D : IS_2D; + + uint16_t indexOfDependentVariable = m_attribute->getAttributeGEOEASgivenIndex()-1; + + //for each observation + uint64_t count = 0; + for( int iRow = 0; iRow < m_dataFile->getDataLineCount(); iRow++ ){ + //get the observed value + double observedValue = m_dataFile->dataConst( iRow, indexOfDependentVariable ); + //if observed value is valid (not no-data-value) + if( ! hasNDV || ! Util::almostEqual2sComplement( observedValue, NDV, 1 ) ){ + // get the xyz location of the current sample. + double x, y, z; + m_dataFile->getDataSpatialLocation( iRow, x, y, z ); + // use the current trend model (defined by its parameters) + // to predict the value at the data location + double predictedValue; + switch (is3D) { //if is slower than a switch + case IS_3D: + predictedValue = parameters.a * x * x + + parameters.b * x * y + + parameters.c * x * z + + parameters.d * y * y + + parameters.e * y * z + + parameters.f * z * z + + parameters.g * x + + parameters.h * y + + parameters.i * z ; + break; + case IS_2D: //z-bearing terms are invariant in 2D datasets + predictedValue = parameters.a * x * x + + parameters.b * x * y + + parameters.d * y * y + + parameters.g * x + + parameters.h * y; + } + //accumulate the observed - predicted differences (prediction errors) + double error = observedValue - predictedValue; + sum += std::abs(error); + count++; + } + } + return sum / count; +} + +void Quadratic3DTrendModelFitting::initDomain( Quad3DTrendModelFittingAuxDefs::ParametersDomain& domain, + double searchWindowSize ) const +{ + for( int i = 0; i < Quad3DTrendModelFittingAuxDefs::N_PARAMS; ++i){ + //define the domain boundaries + domain.min[i] = -searchWindowSize; + domain.max[i] = searchWindowSize; + } +} + +Quad3DTrendModelFittingAuxDefs::Parameters Quadratic3DTrendModelFitting::processWithGenetic( + uint16_t nThreads, uint32_t seed, uint16_t maxNumberOfGenerations, + uint16_t nPopulationSize, uint16_t nSelectionSize, double probabilityOfCrossOver, + uint8_t pointOfCrossover, double mutationRate, double searchWindowSize, double windowWindowShiftThreshold) const +{ + + //clear the collected objective function values. + s_objectiveFunctionValues.clear(); + + //Intialize the random number generator with the same seed + std::srand (seed); + + //sanity checks + if( nSelectionSize >= nPopulationSize ){ + QMessageBox::critical( Application::instance()->getMainWindow(), + "Error", + "Quadratic3DTrendModelFitting::processWithGenetic(): Selection pool size must be less than population size."); + return Quad3DTrendModelFittingAuxDefs::Parameters(); + } + if( nPopulationSize % 2 + nSelectionSize % 2 ){ + QMessageBox::critical( Application::instance()->getMainWindow(), + "Error", + "Quadratic3DTrendModelFitting::processWithGenetic(): Sizes must be even numbers."); + return Quad3DTrendModelFittingAuxDefs::Parameters(); + } + if( pointOfCrossover >= Quad3DTrendModelFittingAuxDefs::N_PARAMS ){ + QMessageBox::critical( Application::instance()->getMainWindow(), + "Error", + "Quadratic3DTrendModelFitting::processWithGenetic(): Point of crossover must be less than the number of parameters."); + return Quad3DTrendModelFittingAuxDefs::Parameters(); + } + + // Fetch data from the data source. + m_dataFile->loadData(); + + //Initialize the optimization domain (boundary conditions) + Quad3DTrendModelFittingAuxDefs::ParametersDomain domain; + Quad3DTrendModelFittingAuxDefs::Parameters domainCenter = {0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0}; //init domain center at origin of the 9-D space + initDomain( domain, searchWindowSize ); + + //initialize the best fit paramaters in the domain center + Quad3DTrendModelFittingAuxDefs::Parameters gbest_pw = domainCenter; + + //if data is 2D, then the z-bearing terms are invariant. + constexpr int IS_3D = 1; + constexpr int IS_2D = 0; + int dimensionality = m_dataFile->isTridimensional() ? IS_3D : IS_2D; + + //loop for the domain window shift if a solution is near the current domain boundaries + //executes at least once + while(true){ + + //=========================================THE GENETIC ALGORITHM================================================== + + //distribute as evenly as possible (load balance) the starting + //points (by their indexes) amongst the threads. + std::vector< std::pair< int, int > > individualsIndexesRanges = + Util::generateSubRanges( 0, nPopulationSize-1, nThreads ); + + //sanity check + assert( individualsIndexesRanges.size() == nThreads && "Quadratic3DTrendModelFitting::processWithGenetic(): " + "number of threads different from individual index ranges. " + " This is likely a bug in Util::generateSubRanges() function." ); + + QProgressDialog progressDialog; + progressDialog.setRange(0, maxNumberOfGenerations); + progressDialog.setValue( 0 ); + progressDialog.show(); + progressDialog.setLabelText("Genetic Algorithm in progress..."); + + //the main genetic algorithm algorithm loop + std::vector< Individual > population; + for( int iGen = 0; iGen < maxNumberOfGenerations; ++iGen ){ + + //the current best solution found so far is kept as part of the population + if( ! gbest_pw.isTrivial() ) { + Individual currentBest( gbest_pw ); + population.push_back( currentBest ); + } + + //complete the population with randomly generated individuals. + while( population.size() < nPopulationSize ){ + //create a set of genes (parameters) + Quad3DTrendModelFittingAuxDefs::Parameters pw; + //randomize the individual's position in the domain. + for( int i = 0; i < Quad3DTrendModelFittingAuxDefs::N_PARAMS; ++i ){ + double LO = domain.min[i]; + double HI = domain.max[i]; + pw[i] = LO + std::rand() / (RAND_MAX/(HI-LO)); + } + //zero-out z-bearing parameters if data set is 2D + if( dimensionality == IS_2D ) + pw.zeroOutZBearingTerms(); + //create an individual with the random genes + Individual ind( pw ); + //add the new individual to the population + population.push_back( ind ); + } + + //create and start the threads. Each thread evaluates the objective function for a series of individuals. + std::thread threads[nThreads]; + unsigned int iThread = 0; + for( const std::pair< int, int >& individualsIndexesRange : individualsIndexesRanges ) { + threads[iThread] = std::thread( taskEvaluateObjetiveInRangeOfIndividualsForGenetic, + std::cref(*this), + individualsIndexesRange.first, + individualsIndexesRange.second, + std::ref( population ) //--> OUTPUT PARAMETER + ); + ++iThread; + } //for each thread (ranges of starting points) + + //wait for the threads to finish. + for( unsigned int iThread = 0; iThread < nThreads; ++iThread) + threads[iThread].join(); + + //sort the population in ascending order (lower value == better fitness) + std::sort( population.begin(), population.end() ); + + //collect the iteration's best objective function value + s_objectiveFunctionValues.push_back( population[0].m_fValue ); + + //clip the population (the excessive worst fit individuals die) + while( population.size() > nPopulationSize ) + population.pop_back(); + + //perform selection by binary tournament + std::vector< Individual > selection; + for( uint iSel = 0; iSel < nSelectionSize; ++iSel ){ + //perform binary tournament + std::vector< Individual > tournament; + { + //draw two different individuals at random from the population for the tournament. + int tournCandidate1 = std::rand() / (double)RAND_MAX * ( population.size() - 1 ); + int tournCandidate2 = tournCandidate1; + while( tournCandidate2 == tournCandidate1 ) + tournCandidate2 = std::rand() / (double)RAND_MAX * ( population.size() - 1 ); + //add the participants in the tournament + tournament.push_back( population[tournCandidate1] ); + tournament.push_back( population[tournCandidate2] ); + //sort the binary tournament + std::sort( tournament.begin(), tournament.end()); + } + //add the best of tournament to the selection pool + selection.push_back( tournament.front() ); + } + + //perform crossover and mutation on the selected individuals + std::vector< Individual > nextGen; + while( ! selection.empty() ){ + //draw two different selected individuals at random for crossover. + int parentIndex1 = std::rand() / (double)RAND_MAX * ( selection.size() - 1 ); + int parentIndex2 = parentIndex1; + while( parentIndex2 == parentIndex1 ) + parentIndex2 = std::rand() / (double)RAND_MAX * ( selection.size() - 1 ); + Individual parent1 = selection[ parentIndex1 ]; + Individual parent2 = selection[ parentIndex2 ]; + selection.erase( selection.begin() + parentIndex1 ); + selection.erase( selection.begin() + parentIndex2 ); + //draw a value between 0.0 and 1.0 from an uniform distribution + double p = std::rand() / (double)RAND_MAX; + //if crossover is due... + if( p < probabilityOfCrossOver ){ + //crossover + std::pair< Individual, Individual> offspring = parent1.crossOver( parent2, pointOfCrossover ); + Individual child1 = offspring.first; + Individual child2 = offspring.second; + //mutate all + child1.mutate( mutationRate, domain.min, domain.max ); + child2.mutate( mutationRate, domain.min, domain.max ); + parent1.mutate( mutationRate, domain.min, domain.max ); + parent2.mutate( mutationRate, domain.min, domain.max ); + if( dimensionality == IS_2D ){ //zero-out z-bearing parameters (genes) + child1.m_parameters.zeroOutZBearingTerms(); + child2.m_parameters.zeroOutZBearingTerms(); + parent1.m_parameters.zeroOutZBearingTerms(); + parent2.m_parameters.zeroOutZBearingTerms(); + } + //add them to the next generation pool + nextGen.push_back( child1 ); + nextGen.push_back( child2 ); + nextGen.push_back( parent1 ); + nextGen.push_back( parent2 ); + } else { //no crossover took place + //simply mutate and insert the parents into the next generation pool + parent1.mutate( mutationRate, domain.min, domain.max ); + parent2.mutate( mutationRate, domain.min, domain.max ); + if( dimensionality == IS_2D ){ //zero-out z-bearing parameters (genes) + parent1.m_parameters.zeroOutZBearingTerms(); + parent2.m_parameters.zeroOutZBearingTerms(); + } + nextGen.push_back( parent1 ); + nextGen.push_back( parent2 ); + } + } + + //make the next generation + population = nextGen; + + //update progress bar + progressDialog.setValue( iGen ); + QApplication::processEvents(); // let Qt update the UI + + } //main algorithm loop + + //=====================================GET RESULTS======================================== + progressDialog.hide(); + + //evaluate the individuals of final population + for( uint iInd = 0; iInd < population.size(); ++iInd ){ + Individual& ind = population[iInd]; + ind.m_fValue = objectiveFunction( ind.m_parameters ); + } + + //sort the population in ascending order (lower value == better fitness) + std::sort( population.begin(), population.end() ); + + //get the parameters of the best individual (set of parameters) + gbest_pw = population[0].m_parameters; + + //if the solution is too close to the current domain boundaries, then + //it is possible that the actual solution lies outside, then the program + //shifts the domain so the current solution lies in its center and rerun + //the genetic algorithm within the new search domain + if( domain.isNearBoundary( gbest_pw, windowWindowShiftThreshold ) ) + domain.centerAt( gbest_pw, searchWindowSize ); + else + break; //terminates the search if the solution is not near the domain boundaries + + } //loop for domain window shift if a solution is near the current domain boundaries + + // Display the results in a window. + //displayResults( variogramStructures, inputFFTimagPhase, inputVarmap, false ); + showObjectiveFunctionEvolution(); + + //return the fitted model + return gbest_pw; +} + +//locally defined data structure to store observation data in a form compatible with GSL's +//non-linear least squares solver. +struct data { + size_t count; + double* x; + double* y; + double* z; + double* value; +}; + +//locally defined C-style function for use with GSL's non-linear least squares solver. +//this function returns the resuiduals (predicted - observed) in residuals_at_data_locations. +int cost_f (const gsl_vector* model_parameters_x, + void* data, + gsl_vector* residuals_at_data_locations) { + + // get data information + size_t data_count = static_cast(data)->count; + double *x = static_cast(data)->x; + double *y = static_cast(data)->y; + double *z = static_cast(data)->z; + double *data_value = static_cast(data)->value; + + // get trend model parameters + double p0 = gsl_vector_get (model_parameters_x, 0); + double p1 = gsl_vector_get (model_parameters_x, 1); + double p2 = gsl_vector_get (model_parameters_x, 2); + double p3 = gsl_vector_get (model_parameters_x, 3); + double p4 = gsl_vector_get (model_parameters_x, 4); + double p5 = gsl_vector_get (model_parameters_x, 5); + double p6 = gsl_vector_get (model_parameters_x, 6); + double p7 = gsl_vector_get (model_parameters_x, 7); + double p8 = gsl_vector_get (model_parameters_x, 8); + + // for each observation + for (size_t i = 0; i < data_count; i++) + { + + // evaluate the trend model at the data location + double predicted_value_i = + p0 * x[i]*x[i] + + p1 * x[i]*y[i] + + p2 * x[i]*z[i] + + p3 * y[i]*y[i] + + p4 * y[i]*z[i] + + p5 * z[i]*z[i] + + p6 * x[i] + + p7 * y[i] + + p8 * z[i]; + + // compute and return the error between measurement and predicted. + gsl_vector_set (residuals_at_data_locations, i, predicted_value_i - data_value[i]); + } + + return GSL_SUCCESS; +} + + +//locally defined C-style function to compute the 1st derivative of +//the cost function in the form of a Jacobian matrix compatible with GSL's +//non-linear least-squares solver +int cost_df ( const gsl_vector * model_parameters_x, + void *data, + gsl_matrix * J ) { + + // get data information + size_t data_count = static_cast(data)->count; + double *x = static_cast(data)->x; + double *y = static_cast(data)->y; + double *z = static_cast(data)->z; + + for (size_t i = 0; i < data_count; i++) { + // An element of the Jacobian matrix: J(i,j) = ∂fi / ∂xj, + // where fi = (Yi - yi)/sigma[i], + // Yi = predicted value at x,y,z , + // yi = data value at x,y,z, + // xj = the nine trend model parameters p0...p8 + // f(x,y,z) = p0*x² + p1*x*y + p2*x*z + p3*y² + p4*y*z + p5*z² + p6*x + p7*y + p8*z + gsl_matrix_set (J, i, 0, x[i]*x[i]); // ∂f/∂p0 + gsl_matrix_set (J, i, 1, x[i]*y[i]); // ∂f/∂p1 + gsl_matrix_set (J, i, 2, x[i]*z[i]); // ∂f/∂p2 + gsl_matrix_set (J, i, 3, y[i]*y[i]); // ∂f/∂p3 + gsl_matrix_set (J, i, 4, y[i]*z[i]); // ∂f/∂p4 + gsl_matrix_set (J, i, 5, z[i]*z[i]); // ∂f/∂p5 + gsl_matrix_set (J, i, 6, x[i]); // ∂f/∂p6 + gsl_matrix_set (J, i, 7, y[i]); // ∂f/∂p7 + gsl_matrix_set (J, i, 8, z[i]); // ∂f/∂p8 + } + + return GSL_SUCCESS; +} + + +// locally defined C-style function that is called for each outer loop iteration of the non-linear +// least squares algorithm of GSL. +// this is usally useful to keep track of progress or to record an execution log +void iteration_callback( const size_t iter, + void *params, + const gsl_multifit_nlinear_workspace *w ) +{ +// gsl_vector *f = gsl_multifit_nlinear_residual(w); +// gsl_vector *x = gsl_multifit_nlinear_position(w); +// double rcond; + +// /* compute reciprocal condition number of J(x) */ +// gsl_multifit_nlinear_rcond(&rcond, w); + +// fprintf(stderr, "iter %2zu: A = %.4f, lambda = %.4f, b = %.4f, cond(J) = %8.4f, |f(x)| = %.4f\n", +// iter, +// gsl_vector_get(x, 0), +// gsl_vector_get(x, 1), +// gsl_vector_get(x, 2), +// 1.0 / rcond, +// gsl_blas_dnrm2(f)); + + //update the progress bar (if set). + if( s_progressDiag_for_iteration_callback ) { + s_progressDiag_for_iteration_callback->setValue( iter ); + QApplication::processEvents(); //let Qt repaint widgets + } +} + +Quad3DTrendModelFittingAuxDefs::Parameters Quadratic3DTrendModelFitting::processWithNonLinearLeastSquares() const +{ + // load dataset + m_dataFile->loadData(); + + // set whether the data is 2D or 3D in Cartesian space + bool is3D = m_dataFile->isTridimensional(); + + // determine the number of valid samples (no NDV) + size_t number_of_valid_samples = 0; + { + if( ! m_dataFile->hasNoDataValue() ){ + number_of_valid_samples = m_dataFile->getDataLineCount(); + } else { + double NDV = m_dataFile->getNoDataValueAsDouble(); + int iColumn = m_attribute->getAttributeGEOEASgivenIndex() - 1; + for( int iRow = 0; iRow < m_dataFile->getDataLineCount(); iRow++){ + double sample_value = m_dataFile->data( iRow, iColumn ); + if( ! Util::almostEqual2sComplement( NDV, sample_value, 1 )) + number_of_valid_samples++; + } + } + } + + // select non-linear least squares method (so far GSL only supports TRS). + // TRS = Trust Region Subproblem. This means that in a small 9-D cube around a given + // x,y,z location, the trend model is assumed near-constant value. + const gsl_multifit_nlinear_type* T = gsl_multifit_nlinear_trust; + + // set the sizes of the problem + const size_t number_of_samples = number_of_valid_samples; + constexpr size_t number_of_parameters = Quad3DTrendModelFittingAuxDefs::N_PARAMS; + + // allocate data structures for algorithm operation + double* x = new double[number_of_samples]; + double* y = new double[number_of_samples]; + double* z = new double[number_of_samples]; + double* data_values = new double[number_of_samples]; + double* weights = new double[number_of_samples]; + struct data dataSet = { number_of_samples, x, y, z, data_values }; + + // GSL documentation is not very clear as to whatever this actually does but it is necessary. + gsl_rng_env_setup(); + gsl_rng* trusted_region = gsl_rng_alloc( gsl_rng_default ); //the returned pointer is not being used by this client code + + // define the function to be minimized + gsl_multifit_nlinear_fdf fdf; + fdf.f = cost_f; + fdf.df = cost_df; // set to NULL for finite-difference Jacobian + fdf.fvv = NULL; // not using geodesic acceleration + fdf.n = number_of_samples; + fdf.p = number_of_parameters; + fdf.params = &dataSet; + + // read the data to be fitted into the arrays + { + bool hasNDV = m_dataFile->hasNoDataValue(); + double NDV = m_dataFile->getNoDataValueAsDouble(); + int iColumn = m_attribute->getAttributeGEOEASgivenIndex() - 1; + size_t iDataArray = 0; + for( int iDataFileRow = 0; iDataFileRow < m_dataFile->getDataLineCount(); iDataFileRow++){ + double sample_value = m_dataFile->data( iDataFileRow, iColumn ); + if( ! hasNDV || ! Util::almostEqual2sComplement( NDV, sample_value, 1 )) { + data_values[iDataArray] = sample_value; + m_dataFile->getDataSpatialLocation( iDataFileRow, x[iDataArray], y[iDataArray], z[iDataArray] ); + weights[iDataArray] = 1.0; + iDataArray++; + } + } + } + + // from this point on, the data file's contents are no longer needed, so deallocate them + // to free up memory space + m_dataFile->freeLoadedData(); + + // allocate workspace with default parameters + gsl_multifit_nlinear_parameters fdf_params = gsl_multifit_nlinear_default_parameters(); + gsl_multifit_nlinear_workspace* workspace = gsl_multifit_nlinear_alloc (T, &fdf_params, number_of_samples, number_of_parameters); + + // initialize solver with starting point and weights + // 0 1 2 3 4 5 6 7 8 + // x² + xy + xz + y²+ yz + z² + x + y + z + double params_init[9] = { 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 }; /* starting parameter values */ + if( ! is3D ) //zero-out z-bearing terms + params_init[2] = params_init[4] = params_init[5] = params_init[8] = 0.0; + gsl_vector_view parameters = gsl_vector_view_array (params_init, number_of_parameters); + gsl_vector_view wts = gsl_vector_view_array (weights, number_of_samples); + gsl_multifit_nlinear_winit (¶meters.vector, &wts.vector, &fdf, workspace); + + // compute initial cost function + gsl_vector* f = gsl_multifit_nlinear_residual(workspace); + double chisq0; + gsl_blas_ddot(f, f, &chisq0); + + // perform the non-linear least squares of 100 iterations + int info; + constexpr int max_iterations = 100; + constexpr double xtol = 1e-8; //parameter epsilon (e.g. 1.0 + epsilon == 1.0) + constexpr double gtol = 1e-8; //1st derivative epsilon + constexpr double ftol = 0.0; //2nd derivative epsilon (not using geodesic acceleration) + int status; + { + QProgressDialog progressDialog; + progressDialog.setRange(0, max_iterations); + progressDialog.setValue( 0 ); + progressDialog.show(); + progressDialog.setLabelText("Non-linear least squares in progress (max. " + QString::number(max_iterations) + " iterations) ..."); + s_progressDiag_for_iteration_callback = &progressDialog; + status = gsl_multifit_nlinear_driver(max_iterations, xtol, gtol, ftol, iteration_callback, NULL, &info, workspace); + s_progressDiag_for_iteration_callback = nullptr; //avoid dangling pointer after progressDialog goes out of scope + } + + // get the final Jacobian matrix + gsl_matrix* J = gsl_multifit_nlinear_jac(workspace); + + // compute the final covariance matrix of best fit parameters to get fit errors + gsl_matrix* covar = gsl_matrix_alloc (number_of_parameters, number_of_parameters); + gsl_multifit_nlinear_covar (J, 0.0, covar); + + // compute final cost value + double chisq; + gsl_blas_ddot(f, f, &chisq); + + // output execution summary + Application::instance()->logInfo( "summary from method '" + + QString(gsl_multifit_nlinear_name(workspace)) + "/" + + QString(gsl_multifit_nlinear_trs_name(workspace)) + "'" ); + Application::instance()->logInfo( "number of iterations: " + + QString::number(gsl_multifit_nlinear_niter(workspace)) ); + Application::instance()->logInfo( "function evaluations: " + + QString::number(fdf.nevalf) ); + Application::instance()->logInfo( "Jacobian evaluations: " + + QString::number(fdf.nevaldf) ); + Application::instance()->logInfo( "reason for stopping:: " + + (info == 1 ? QString("small step size") : QString("small gradient")) ); + Application::instance()->logInfo( "initial |f(x)| = " + + QString::number(sqrt(chisq0)) ); + Application::instance()->logInfo( "final |f(x)| = " + + QString::number(sqrt(chisq)) ); + + // output the best fit trend model parameters with their fit errors + { + double dof = number_of_samples - number_of_parameters; + double c = GSL_MAX_DBL(1, sqrt(chisq / dof)); + + Application::instance()->logInfo( "chisq/dof = " + + QString::number(chisq / dof) ); + + // some lambda functions to serve as "macros" + auto FIT = [workspace](size_t i){ return gsl_vector_get(workspace->x, i); }; + auto ERR = [covar](size_t i){ return std::sqrt(gsl_matrix_get(covar,i,i)); }; + + //output the parameter values of the best fit trend model. + for(size_t i = 0; i < Quad3DTrendModelFittingAuxDefs::N_PARAMS; i++ ) + Application::instance()->logInfo( "p" + QString::number(i) + " = " + + QString::number(FIT(i)) + " +/- " + + QString::number(c*ERR(i)) ); + } + + // output final execution status (e.g. whether it completed by convergence or + // maximum number of iterations has been reached) + Application::instance()->logInfo( "status = " + + QString( gsl_strerror(status) ) ); + + // free up memory + gsl_multifit_nlinear_free (workspace); + gsl_matrix_free (covar); + gsl_rng_free (trusted_region); + delete[] x; + delete[] y; + delete[] z; + delete[] data_values; + delete[] weights; + + // return best-fit trend model + auto get_par = [workspace](size_t i){ return gsl_vector_get(workspace->x, i); }; + return Quad3DTrendModelFittingAuxDefs::Parameters( get_par(0), + get_par(1), + get_par(2), + get_par(3), + get_par(4), + get_par(5), + get_par(6), + get_par(7), + get_par(8) ); +} + +void Quadratic3DTrendModelFitting::showObjectiveFunctionEvolution() const +{ + //load the x,y data for the chart + QtCharts::QLineSeries *chartSeries = new QtCharts::QLineSeries(); + double max = std::numeric_limits::lowest(); + for(uint i = 0; i < s_objectiveFunctionValues.size(); ++i){ + chartSeries->append( i+1, s_objectiveFunctionValues[i] ); + if( s_objectiveFunctionValues[i] > max ) + max = s_objectiveFunctionValues[i]; + } + + //create a new chart object + QtCharts::QChart *objFuncValuesChart = new QtCharts::QChart(); + { + objFuncValuesChart->addSeries( chartSeries ); + objFuncValuesChart->axisX( chartSeries ); + + QtCharts::QValueAxis* axisX = new QtCharts::QValueAxis(); + axisX->setLabelFormat("%i"); + objFuncValuesChart->setAxisX(axisX, chartSeries); + + QtCharts::QValueAxis *axisY = new QtCharts::QValueAxis(); + axisY->setLabelFormat("%3.2f"); + axisY->setRange( 0.0, max ); + objFuncValuesChart->setAxisY(axisY, chartSeries); + + objFuncValuesChart->legend()->hide(); + } + + //create the chart dialog + EmptyDialog* ed = new EmptyDialog( Application::instance()->getMainWindow() ); + QtCharts::QChartView* chartView = new QtCharts::QChartView( objFuncValuesChart ); + ed->addWidget( chartView ); + ed->setWindowTitle( "Objective function with iterations." ); + ed->show(); +} diff --git a/geostats/quadratic3dtrendmodelfitting.h b/geostats/quadratic3dtrendmodelfitting.h new file mode 100644 index 00000000..643793b1 --- /dev/null +++ b/geostats/quadratic3dtrendmodelfitting.h @@ -0,0 +1,171 @@ +#ifndef QUADRATIC3DTRENDMODELFITTING_H +#define QUADRATIC3DTRENDMODELFITTING_H + +#include +#include +#include +#include + +class DataFile; +class Attribute; +class QProgressDialog; + +/** + * A pointer to a progress dialog so a GSL's C-style callback function can access it + * to display algorithm progress. + */ +static QProgressDialog* s_progressDiag_for_iteration_callback = nullptr; + +namespace Quad3DTrendModelFittingAuxDefs { + + /** The number of model parameters. */ + constexpr static int N_PARAMS = 9; + + /** A shortcut to a standard NaN value. */ + constexpr static double NaN = std::numeric_limits::quiet_NaN(); + + /** Data structure: the nine model parameters to find. */ + struct Parameters { + Parameters():a(NaN),b(NaN),c(NaN),d(NaN),e(NaN),f(NaN),g(NaN),h(NaN),i(NaN){} //all-NaN intialization + Parameters(double a, double b, double c, double d, double e, double f, double g, double h, double i): + a(a),b(b),c(c),d(d),e(e),f(f),g(g),h(h),i(i){} //valid values intialization + double a, b, c, d, e, f, g, h, i; //the individual parameters + //xx,xy,xz,yy,yz,zz, x, y, z //the model terms corresponding to the parameters + double& operator [](int index){ //this allows the set to be used as an array (l-value). + switch(index){ + case 0: return a; case 1: return b; case 2: return c; + case 3: return d; case 4: return e; case 5: return f; + case 6: return g; case 7: return h; case 8: return i; + default: assert( false && "Quadratic3DTrendModelFitting::Parameters: index out of bounds (must be between 0 and 8)."); + } + } + double operator [](int index) const { //this allows the set to be used as an array (r-value). + switch(index){ + case 0: return a; case 1: return b; case 2: return c; + case 3: return d; case 4: return e; case 5: return f; + case 6: return g; case 7: return h; case 8: return i; + default: assert( false && "Quadratic3DTrendModelFitting::Parameters: index out of bounds (must be between 0 and 8)."); + } + } + bool isTrivial() const { + return a==0.0 && b==0.0 && c==0.0 && d==0.0 && e==0.0 && f==0.0 && g==0.0 && h==0.0 && i==0.0; + } + void zeroOutZBearingTerms(){ c = e = f = i = 0.0; } //this is mostly used with 2D datasets (z-less data) + }; + + /** Data structure: the parameters domain for the optimization methods (boundary conditions). */ + struct ParametersDomain { + Parameters min; + Parameters max; + bool isNearBoundary( const Parameters& parametersToTest, double threshold ) const { + for( int i = 0; i < N_PARAMS; i++ ){ + if( std::abs(parametersToTest[i] - min[i]) <= threshold ) + return true; + else if( std::abs(parametersToTest[i] - max[i]) <= threshold ) + return true; + } + return false; + } + void centerAt( const Parameters& targetParameters, double domainSize ){ + double halfSize = domainSize/2; + for( int i = 0; i < N_PARAMS; i++ ){ + min[i] = targetParameters[i] - halfSize; + max[i] = targetParameters[i] + halfSize; + } + } + }; +} + +/** + * The Quadratic3DTrendModelFitting class comprises the functionalities to fit a 3D quadratic trend model + * to a data set. The trend model has nine terms in the form: + * f(x,y,z) = a*x*x + b*x*y + c*x*z + d*y*y + e*y*z + f*z*z + g*x + h*y + i*z + * where: a...i are the model parameters. + * This form follows the one expected in the configuration of the kt3d GSLib program for kriging with + * a trend model. The model predicts the value of a continuous variable given a location in space. + * The fitting process searchs for a set of a...i parameters such that the differences between f(x,y,z) + * and the observed values at data locations are minimized. + */ +class Quadratic3DTrendModelFitting +{ +public: + + /** Constructor. */ + Quadratic3DTrendModelFitting(DataFile* dataFile, Attribute* attribute); + + + + /** The objective function for an optimization process to find a minimum (ideally the global one). + * @param parameters The parameters of the model to test. + * @return The sum of the differences between all the observed and predicted values. + */ + double objectiveFunction (const Quad3DTrendModelFittingAuxDefs::Parameters& parameters) const; + + + /** Method called by a method of optimization to initialize the + * parameter domain (boundary conditions) and the set of parameters. + * All non-const method parameters are output parameters. + * @param domain The domain boundaries to initialize as two sets of + * parameters: with the minima and with the maxima. + * @param searchWindowSize Defines the size of the parameter domain. E.g. if set to 10.00, then each of the nine coefficients + * will be restricted to an initial interval between -10.00 and 10.00. + */ + void initDomain( Quad3DTrendModelFittingAuxDefs::ParametersDomain& domain, + double searchWindowSize) const; + + /** Performs the trend model fitting using a Genetic Algorithm as optimization method. + *...................................Global Parameters.................................... + * @param nThreads Number of parallel execution threads. + * @param seed Seed for the the random number generator. + *...................................GA Parameters................................. + * @param maxNumberOfGenerations number of generations (iterations). + * @param nPopulationSize Number of individuals (solutions). Mas be an even number. + * @param nSelectionSize The size of the selection pool (must be < nPopulationSize and be an even number). + * @param probabilityOfCrossOver The probability of crossover (between 0.0 and 1.0); + * @param pointOfCrossover The parameter index from which crossover takes place (must be less than the total + * number of parameters per individual). + * @param mutationRate Mutation rate means how many paramaters are expected to change per mutation + * the probability of any parameter parameter (gene) to be changed is 1/nParameters * mutationRate + * thus, 1.0 means that one gene will surely be mutated per mutation on average. Fractionary + * values are possible. 0.0 means no mutation will take place. + * @param searchWindowSize Defines the size of the parameter domain. E.g. if set to 10.00, then each of the nine coefficients + * will be restricted to an initial interval between -10.00 and 10.00. + * @param windowWindowShiftThreshold If any of the coefficients of the optimal model lies close to the boundaries of the domain + * then the algorithm defines a new parameter domain centered on it and resumes search. This + * value controls how close a parameter must be to a domain boundary to enable a new search. + * @returns The set of parameters of the fit trend model. An all-NaN set of parameters is returned if something goes + * wrong during processing (error messages will be displayed). + */ + Quad3DTrendModelFittingAuxDefs::Parameters processWithGenetic( + uint16_t nThreads, + uint32_t seed, + uint16_t maxNumberOfGenerations, + uint16_t nPopulationSize, + uint16_t nSelectionSize, + double probabilityOfCrossOver, + uint8_t pointOfCrossover, + double mutationRate, + double searchWindowSize, + double windowWindowShiftThreshold + ) const; + + Quad3DTrendModelFittingAuxDefs::Parameters processWithNonLinearLeastSquares() const; + + +private: + Attribute* m_attribute; + DataFile* m_dataFile; + + /** The objective function values collected during the last execution + * of an optimization method. + */ + static std::vector< double > s_objectiveFunctionValues; + + /** + * Displays a dialog with a chart showing the evolution of the objective function + * value versus iterations of the optimization method. + */ + void showObjectiveFunctionEvolution( ) const; +}; + +#endif // QUADRATIC3DTRENDMODELFITTING_H diff --git a/geostats/searchbox.cpp b/geostats/searchbox.cpp new file mode 100644 index 00000000..e9ff1e3b --- /dev/null +++ b/geostats/searchbox.cpp @@ -0,0 +1,31 @@ +#include "searchbox.h" + +SearchBox::SearchBox( const BoundingBox& boundingBox ) : + m_boundingBox( boundingBox ) +{ + +} + +void SearchBox::getBBox(double centerX, double centerY, double centerZ, + double &minX, double &minY, double &minZ, + double &maxX, double &maxY, double &maxZ) const +{ + double semiWidth = m_boundingBox.getXsize() / 2; + double semiLength = m_boundingBox.getYsize() / 2; + double semiHeight = m_boundingBox.getZsize() / 2; + minX = centerX - semiWidth; + maxX = centerX + semiWidth; + minY = centerY - semiLength; + maxY = centerY + semiLength; + minZ = centerZ - semiHeight; + maxZ = centerZ + semiHeight; +} + +bool SearchBox::isInside(double centerX, double centerY, double centerZ, + double x, double y, double z) const +{ + double minX, minY, minZ, maxX, maxY, maxZ; + //the geometry for inside test of a search box is exactly its bounding box + getBBox( centerX, centerY, centerZ, minX, minY, minZ, maxX, maxY, maxZ ); + return ( x >= minX && x <= maxX && y >= minY && y <= maxY && z >= minZ && z <= maxZ ); +} diff --git a/geostats/searchbox.h b/geostats/searchbox.h new file mode 100644 index 00000000..1db1e49d --- /dev/null +++ b/geostats/searchbox.h @@ -0,0 +1,39 @@ +#ifndef SEARCHBOX_H +#define SEARCHBOX_H + +#include "geometry/boundingbox.h" +#include "geostats/searchneighborhood.h" + +/** + * The SearchBox class represents a search neighborhood with the shape of a 3D box. + */ +class SearchBox : public SearchNeighborhood +{ +public: + + SearchBox( const BoundingBox& boundingBox ); + + /** Move constructor. */ + SearchBox( SearchBox&& right_hand_side ); + + + ////-------------SearchNeighborhood interface------------------------------------------ + virtual void getBBox( double centerX, double centerY, double centerZ, + double& minX, double& minY, double& minZ, + double& maxX, double& maxY, double& maxZ ) const; + + virtual bool isInside(double centerX, double centerY, double centerZ, + double x, double y, double z ) const; + + virtual bool hasSpatialFiltering() const { return false; } + + virtual void performSpatialFilter( double centerX, double centerY, double centerZ, + std::vector< IndexedSpatialLocationPtr >& samplesLocations, + const SearchStrategy& parentSearchStrategy ) const {} + ////------------------------------------------------------------------------------------- + + + BoundingBox m_boundingBox; +}; + +#endif // SEARCHBOX_H diff --git a/gslib/gslibparameterfiles/kt3dtrendmodelparameters.cpp b/gslib/gslibparameterfiles/kt3dtrendmodelparameters.cpp new file mode 100644 index 00000000..ef750eca --- /dev/null +++ b/gslib/gslibparameterfiles/kt3dtrendmodelparameters.cpp @@ -0,0 +1,116 @@ +#include "kt3dtrendmodelparameters.h" +#include "gslib/gslibparameterfiles/gslibparamtypes.h" + +Kt3dTrendModelParameters::Kt3dTrendModelParameters() : + GSLibParameterFile({ + " -x²", // 0 + " -xy", // 1 + " -xz", // 2 + " -y²", // 3 + " -yz", // 4 + " -z²", // 5 + " -x", // 6 + " -y", // 7 + " -z", // 8 + }) +{ + getParameter(0)->_value = 0.0; + getParameter(1)->_value = 0.0; + getParameter(2)->_value = 0.0; + getParameter(3)->_value = 0.0; + getParameter(4)->_value = 0.0; + getParameter(5)->_value = 0.0; + getParameter(6)->_value = 0.0; + getParameter(7)->_value = 0.0; + getParameter(8)->_value = 0.0; +} + +double Kt3dTrendModelParameters::getA() +{ + return getParameter(0)->_value; +} + +double Kt3dTrendModelParameters::getB() +{ + return getParameter(1)->_value; +} + +double Kt3dTrendModelParameters::getC() +{ + return getParameter(2)->_value; +} + +double Kt3dTrendModelParameters::getD() +{ + return getParameter(3)->_value; +} + +double Kt3dTrendModelParameters::getE() +{ + return getParameter(4)->_value; +} + +double Kt3dTrendModelParameters::getF() +{ + return getParameter(5)->_value; +} + +double Kt3dTrendModelParameters::getG() +{ + return getParameter(6)->_value; +} + +double Kt3dTrendModelParameters::getH() +{ + return getParameter(7)->_value; +} + +double Kt3dTrendModelParameters::getI() +{ + return getParameter(8)->_value; +} + +void Kt3dTrendModelParameters::setA(double value) +{ + getParameter(0)->_value = value; +} + +void Kt3dTrendModelParameters::setB(double value) +{ + getParameter(1)->_value = value; +} + +void Kt3dTrendModelParameters::setC(double value) +{ + getParameter(2)->_value = value; +} + +void Kt3dTrendModelParameters::setD(double value) +{ + getParameter(3)->_value = value; +} + +void Kt3dTrendModelParameters::setE(double value) +{ + getParameter(4)->_value = value; +} + +void Kt3dTrendModelParameters::setF(double value) +{ + getParameter(5)->_value = value; +} + +void Kt3dTrendModelParameters::setG(double value) +{ + getParameter(6)->_value = value; +} + +void Kt3dTrendModelParameters::setH(double value) +{ + getParameter(7)->_value = value; +} + +void Kt3dTrendModelParameters::setI(double value) +{ + getParameter(8)->_value = value; +} diff --git a/gslib/gslibparameterfiles/kt3dtrendmodelparameters.h b/gslib/gslibparameterfiles/kt3dtrendmodelparameters.h new file mode 100644 index 00000000..5d829494 --- /dev/null +++ b/gslib/gslibparameterfiles/kt3dtrendmodelparameters.h @@ -0,0 +1,36 @@ +#ifndef KT3DTRENDMODELPARAMETERS_H +#define KT3DTRENDMODELPARAMETERS_H + +#include "gslib/gslibparameterfiles/gslibparameterfile.h" + +/** Specialization of GSLibParameterFile used to store kt3d's trend model's nine parameters: + * Ax² + Bxy + Cxz + Dy² + Eyz + Fz² + Gx + Hy + Iz + */ +class Kt3dTrendModelParameters : public GSLibParameterFile +{ +public: + + Kt3dTrendModelParameters(); + + double getA(); + double getB(); + double getC(); + double getD(); + double getE(); + double getF(); + double getG(); + double getH(); + double getI(); + + void setA( double value ); + void setB( double value ); + void setC( double value ); + void setD( double value ); + void setE( double value ); + void setF( double value ); + void setG( double value ); + void setH( double value ); + void setI( double value ); +}; + +#endif // KT3DTRENDMODELPARAMETERS_H diff --git a/imagejockey/imagejockeyutils.cpp b/imagejockey/imagejockeyutils.cpp index fbed9d62..d079e956 100644 --- a/imagejockey/imagejockeyutils.cpp +++ b/imagejockey/imagejockeyutils.cpp @@ -769,8 +769,8 @@ spectral::array ImageJockeyUtils::interpolateNullValuesShepard(const spectral::a //populate the VTK collections above double x, y, z; - double minValue = std::numeric_limits::max(); - double maxValue = std::numeric_limits::min(); + double minValue = std::numeric_limits::max(); + double maxValue = -std::numeric_limits::max(); for( int k = 0; k < nK; ++k ) for( int j = 0; j < nJ; ++j ) for( int i = 0; i < nI; ++i ){ diff --git a/mainwindow.cpp b/mainwindow.cpp index d4a35dae..d8b7fb2e 100644 --- a/mainwindow.cpp +++ b/mainwindow.cpp @@ -93,6 +93,7 @@ #include "dialogs/subgriddialog.h" #include "dialogs/gridrepositiondialog.h" #include "dialogs/contactanalysisdialog.h" +#include "dialogs/driftanalysisdialog.h" #include "vertpropcurves/verticalproportioncurvedialog.h" #include "viewer3d/view3dwidget.h" #include "imagejockey/imagejockeydialog.h" @@ -633,10 +634,11 @@ void MainWindow::onProjectContextMenu(const QPoint &mouse_location) makeMenuFlipData(); _projectContextMenu->addMenu( m_subMenuFlipData ); } - if( Util::isIn( parent_file->getFileType(), {"POINTSET","CARTESIANGRID","SEGMENTSET"} ) ) - _projectContextMenu->addAction("Delete variable", this, SLOT(onDeleteVariable())); if( _right_clicked_attribute->isCategorical() ) _projectContextMenu->addAction("Make facies transition matrix", this, SLOT(onMakeFaciesTransitionMatrix())); + _projectContextMenu->addAction("Drift analysis", this, SLOT(onPerformDriftAnalysis())); + if( Util::isIn( parent_file->getFileType(), {"POINTSET","CARTESIANGRID","SEGMENTSET"} ) ) + _projectContextMenu->addAction("Delete variable", this, SLOT(onDeleteVariable())); } //two items were selected. The context menu depends on the combination of items. } else if ( selected_indexes.size() == 2 ) { @@ -1509,6 +1511,17 @@ void MainWindow::onPerformContactAnaysis() } } +void MainWindow::onPerformDriftAnalysis() +{ + DataFile* df = dynamic_cast( _right_clicked_attribute->getContainingFile() ); + if( ! df ){ + QMessageBox::critical(this, "Error", "File is not a data file."); + return; + } + DriftAnalysisDialog* dag = new DriftAnalysisDialog(df, _right_clicked_attribute, this); + dag->show(); +} + void MainWindow::openTransiographyBayesian() { TransiogramBandDialog* tbd = new TransiogramBandDialog( nullptr, nullptr, this ); diff --git a/mainwindow.h b/mainwindow.h index 16f25b9e..d8510163 100644 --- a/mainwindow.h +++ b/mainwindow.h @@ -253,6 +253,7 @@ private slots: void onReviewTransiogramBand(); void onRepositionGrid(); void onPerformContactAnaysis(); + void onPerformDriftAnalysis(); private: /** diff --git a/resources.qrc b/resources.qrc index 339ccb94..c4747b42 100644 --- a/resources.qrc +++ b/resources.qrc @@ -139,5 +139,6 @@ art/iconsHD/v3Druler.png art/iconsHD/v3DprojPar32.png art/iconsHD/v3DprojPer32.png + art/iconsHD/calc32.png diff --git a/spatialindex/spatialindex.cpp b/spatialindex/spatialindex.cpp index 388b9ac3..b55fd8f8 100644 --- a/spatialindex/spatialindex.cpp +++ b/spatialindex/spatialindex.cpp @@ -2,6 +2,7 @@ #include "domain/pointset.h" #include "domain/application.h" +#include "geometry/boundingbox.h" #include "geostats/searchellipsoid.h" #include "geostats/datacell.h" #include "geostats/searchstrategy.h" @@ -118,7 +119,7 @@ void SpatialIndex::fillWithBBoxes(GeoGrid *gg) for( uint iLine = 0; iLine < totlines; ++iLine){ //get the cell's bounding box (each line corresponds to a cell) double minX, minY, minZ, maxX, maxY, maxZ; - gg->getBoundingBox( iLine, minX, minY, minZ, maxX, maxY, maxZ ); + gg->getCellBoundingBox( iLine, minX, minY, minZ, maxX, maxY, maxZ ); //make the bounding box object Box box( Point3D(minX, minY, minZ), @@ -404,6 +405,33 @@ QList SpatialIndex::getNearestWithinGenericRTreeBased(const DataCell& data return result; } +QList SpatialIndex::getWithinBoundingBox(const BoundingBox &bbox) const +{ + assert( m_dataFile && "SpatialIndexPoints::getWithinBoundingBox(): No data file. Make sure there's a call to DataSet::fill() prior to making queries."); + + QList result; + + double maxX = bbox.m_maxX; + double maxY = bbox.m_maxY; + double maxZ = bbox.m_maxZ; + double minX = bbox.m_minX; + double minY = bbox.m_minY; + double minZ = bbox.m_minZ; + Box searchBB( Point3D( minX, minY, minZ ), + Point3D( maxX, maxY, maxZ )); + + //Get all the points within the bounding box of the search neighborhood. + std::vector pointsInSearchBB; + m_rtree.query( bgi::intersects( searchBB ), std::back_inserter( pointsInSearchBB ) ); + + //Get all the row indexes of the samples that fall inside the passed bounding box. + std::vector::iterator it = pointsInSearchBB.begin(); + for(; it != pointsInSearchBB.end(); ++it) + result.push_back( (*it).second ); + + return result; +} + inline bool operator< (const BoxAndDataIndexAndDistance& boxAndDataIndexAndDistance1, const BoxAndDataIndexAndDistance& boxAndDataIndexAndDistance2) { return boxAndDataIndexAndDistance1.second < boxAndDataIndexAndDistance2.second ; diff --git a/spatialindex/spatialindex.h b/spatialindex/spatialindex.h index 66843ec7..a0e14291 100644 --- a/spatialindex/spatialindex.h +++ b/spatialindex/spatialindex.h @@ -14,6 +14,7 @@ class DataFile; class GeoGrid; class SegmentSet; class GridCell; +class BoundingBox; namespace bg = boost::geometry; namespace bgi = boost::geometry::index; @@ -92,6 +93,13 @@ class SpatialIndex QList getNearestWithinGenericRTreeBased(const DataCell& dataCell, const SearchStrategy & searchStrategy ) const; + /** + * Returns the indexes of all the data lines within the given bounding box. + * The indexes are the data line indexes (file data lines) of the DataFile + * used fill the index. May return an empty list. + */ + QList getWithinBoundingBox( const BoundingBox& bbox ) const; + /** * Does the same as getNearestWithinGenericRTreeBased() but is tuned for large, high-density data sets. * It may run slower for smaller data sets than the former, though. diff --git a/util.cpp b/util.cpp index 13228e13..f2a8b450 100644 --- a/util.cpp +++ b/util.cpp @@ -1089,8 +1089,8 @@ void Util::importSettingsFromPreviousVersion() QSettings currentSettings; //The list of previous versions (order from latest to oldest version is advised) QStringList previousVersions; - previousVersions << "6.18" << "6.17" << "6.16" << "6.14" << "6.12" << "6.9" << "6.7" << "6.6" << "6.5" << "6.3" << "6.2" - << "6.1" << "6.0" << "5.7.1" + previousVersions << "6.20" << "6.18" << "6.17" << "6.16" << "6.14" << "6.12" << "6.9" << "6.7" << "6.6" + << "6.5" << "6.3" << "6.2" << "6.1" << "6.0" << "5.7.1" << "5.7" << "5.5" << "5.3" << "5.1" << "5.0" << "4.9" << "4.7" << "4.5.1" << "4.5" << "4.3.3" << "4.3" << "4.0" << "3.8" << "3.6.1" << "3.6" << "3.5" << "3.2" << "3.0" << "2.7.2" << "2.7.1" << "2.7" << "2.5.1" << "2.5" << "2.4" << "2.3" << "2.2" << "2.1" diff --git a/widgets/linechartwidget.cpp b/widgets/linechartwidget.cpp index 5ddceb44..3369b7f5 100644 --- a/widgets/linechartwidget.cpp +++ b/widgets/linechartwidget.cpp @@ -5,12 +5,49 @@ #include #include +#include #include -LineChartWidget::LineChartWidget(QWidget *parent) : + +/** Deriving QChartView to support a custom mouse wheel zoom mechanics. + * changed from the original from https://stackoverflow.com/users/6622587/eyllanesc (https://stackoverflow.com/a/48626725/2153955) + */ +class MyChartView : public QtCharts::QChartView { + +public: + MyChartView( QtCharts::QChart* chart, + LineChartWidget::ZoomDirection zoomDirection = LineChartWidget::ZoomDirection::VERTICAL ) + : QtCharts::QChartView( chart ), + mZoomDirection( zoomDirection ) { } + +private: + qreal mFactor=1.0; + LineChartWidget::ZoomDirection mZoomDirection; + +protected: + void wheelEvent(QWheelEvent *event) Q_DECL_OVERRIDE { + chart()->zoomReset(); + + mFactor *= event->angleDelta().y() > 0 ? 0.5 : 2; + + QRectF rect = chart()->plotArea(); + QPointF c = chart()->plotArea().center(); + if( mZoomDirection == LineChartWidget::ZoomDirection::VERTICAL ) + rect.setHeight(mFactor*rect.height()); + else + rect.setWidth(mFactor*rect.width()); + rect.moveCenter(c); + chart()->zoomIn(rect); + + QChartView::wheelEvent(event); + } +}; + +LineChartWidget::LineChartWidget(QWidget *parent, ZoomDirection zoomDirection) : QWidget(parent), ui(new Ui::LineChartWidget), - m_sharedYaxis(false) + m_sharedYaxis(false), + m_ZoomDirection( zoomDirection ) { ui->setupUi(this); @@ -30,7 +67,7 @@ LineChartWidget::LineChartWidget(QWidget *parent) : m_axisX->setTickCount(11); //create and adds the chart view widget - m_chartView = new QtCharts::QChartView( m_chart ); + m_chartView = new MyChartView( m_chart, zoomDirection ); m_chartView->setRenderHint(QPainter::Antialiasing); this->layout()->addWidget( m_chartView ); } @@ -42,9 +79,11 @@ LineChartWidget::~LineChartWidget() void LineChartWidget::setData(const std::vector > &data, int indexForXAxis, + bool clearCurrentCurves, const std::map& yVariablesCaptions, const std::map& yVariablesYaxisTitles, - const std::map& yVariablesColors) + const std::map& yVariablesColors, + const std::map& yVariablesStyles ) { // Does nothing if the input data table is empty. if( data.empty() ) @@ -55,7 +94,8 @@ void LineChartWidget::setData(const std::vector > &data, double maxX = -std::numeric_limits::max(); // Clears all current data series. - m_chart->removeAllSeries(); + if( clearCurrentCurves ) + m_chart->removeAllSeries(); // Initialize the limits for the Y axes (one per dependent variable or a global Y axis). double minY = std::numeric_limits::max(); @@ -83,6 +123,11 @@ void LineChartWidget::setData(const std::vector > &data, series->setName( yVariablesCaptions.at( iSeries ) ); } catch (...){} + // tries to set a line style for the current series (user may not have set one) + try { + series->setPen( yVariablesStyles.at( iSeries ) ); + } catch (...){} + // tries to set a color for the current series (user may not have set one) try { series->setColor( yVariablesColors.at( iSeries ) ); @@ -109,6 +154,9 @@ void LineChartWidget::setData(const std::vector > &data, //create the single Y axis or a new Y axis per dependent variable. if( ! m_sharedYaxis || ! axisY ){ + //remove any previously added Y axes + for( QtCharts::QAbstractAxis* axis : m_chart->axes( Qt::Vertical ) ) + m_chart->removeAxis( axis ); axisY = new QtCharts::QValueAxis(); axisY->setTickCount(11); m_chart->addAxis(axisY, Qt::AlignLeft); @@ -116,6 +164,10 @@ void LineChartWidget::setData(const std::vector > &data, axisY->setRange( minY, maxY ); series->attachAxis( axisY ); + //make sure the curves get rescaled accordingly + for( QtCharts::QAbstractSeries* a_series : m_chart->series() ) + a_series->attachAxis( axisY ); + // tries to set a Y-axis title for the current series (user may not have set one) try { axisY->setTitleText( yVariablesYaxisTitles.at( iSeries ) ); @@ -138,3 +190,11 @@ void LineChartWidget::setXaxisCaption(const QString caption) { m_axisX->setTitleText( caption ); } + +void LineChartWidget::setLegendVisible(const bool value) +{ + if( value ) + m_chart->legend()->show(); + else + m_chart->legend()->hide(); +} diff --git a/widgets/linechartwidget.h b/widgets/linechartwidget.h index bc36b56c..48afc099 100644 --- a/widgets/linechartwidget.h +++ b/widgets/linechartwidget.h @@ -20,12 +20,17 @@ class LineChartWidget : public QWidget public: + enum class ZoomDirection : unsigned int{ + VERTICAL, + HORIZONTAL + }; + /** * Displays a multivariate line chart with a shared X-axis. * The Y-axis can be shared or one per Y variable * @see setSharedYaxis() */ - explicit LineChartWidget(QWidget *parent = nullptr); + explicit LineChartWidget( QWidget *parent = nullptr, ZoomDirection zoomDirection = ZoomDirection::VERTICAL ); ~LineChartWidget(); /** @@ -35,6 +40,7 @@ class LineChartWidget : public QWidget * @param indexForXAxis The index of element in the inner vectors corresponding to the * independent variable, that is, the values in the X-axis. The other * values will be considered Y variables. + * @param clearCurrentCurves If true, the chart's current curves are cleared. * @param yVariablesCaptions Sets the captions for each curve legend. If a caption is not set for a given curve, * the legend will be displayed with a blank caption. The uint8_t value is the index of * the Y value in the inner vectors in data vector. @@ -46,12 +52,18 @@ class LineChartWidget : public QWidget * @param yVariablesColors Sets the colors for each curve. If a color is not set for a given curve, * it'll will be displayed with a default color. The uint8_t value is the index of * the Y value in the inner vectors in data vector. + * @param yVariablesStyles Sets the line styles for each curve. If a line style is not set for a given curve, + * it'll will be displayed with a default style. The uint8_t value is the index of + * the Y value in the inner vectors in data vector. Note, the color in yVariablesColors, + * if set, overrides the style's color. */ void setData( const std::vector< std::vector< double > >& data, int indexForXAxis, + bool clearCurrentCurves = true, const std::map& yVariablesCaptions = {}, const std::map& yVariablesYaxisTitles = {}, - const std::map& yVariablesColors = {} ); + const std::map& yVariablesColors = {}, + const std::map& yVariablesStyles = {} ); /** * Enables/disables whether all Y series share the same vertical axis in the chart. @@ -60,7 +72,7 @@ class LineChartWidget : public QWidget void setSharedYaxis( bool value ){ m_sharedYaxis = value; } /** - * Sets the chart's title. The new title becomes effective after the next call to setData(). + * Sets the chart's title. Takes effect upon calling. */ void setChartTitle( const QString chartTitle ); @@ -69,6 +81,11 @@ class LineChartWidget : public QWidget */ void setXaxisCaption( const QString caption ); + /** + * Shows or hides the chart legend (default is true). Takes effect upon calling. + */ + void setLegendVisible( const bool value ); + private: Ui::LineChartWidget *ui; @@ -79,6 +96,7 @@ class LineChartWidget : public QWidget bool m_sharedYaxis; + ZoomDirection m_ZoomDirection; }; #endif // LINECHARTWIDGET_H