diff --git a/CorePoints.cpp b/CorePoints.cpp new file mode 100644 index 0000000..0dd7f75 --- /dev/null +++ b/CorePoints.cpp @@ -0,0 +1,112 @@ +//########################################################################## +//# # +//# CLOUDCOMPARE PLUGIN: q3DMASC # +//# # +//# This program is free software; you can redistribute it and/or modify # +//# it under the terms of the GNU General Public License as published by # +//# the Free Software Foundation; version 2 or later of the License. # +//# # +//# This program is distributed in the hope that it will be useful, # +//# but WITHOUT ANY WARRANTY; without even the implied warranty of # +//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +//# GNU General Public License for more details. # +//# # +//# COPYRIGHT: Dimitri Lague / CNRS / UEB # +//# # +//########################################################################## + +#include "CorePoints.h" + +//qCC_db +#include + +//CCLib +#include + +//system +#include + +using namespace masc; + +bool CorePoints::prepare(CCLib::GenericProgressCallback* progressCb/*=nullptr*/) +{ + if (!origin) + { + assert(false); + return false; + } + + if (selection) + { + //nothing to do + return true; + } + + //now we can compute the subsampled version + CCLib::ReferenceCloud* ref = nullptr; + switch (selectionMethod) + { + case SPATIAL: + { + //we'll need an octree + if (!origin->getOctree()) + { + if (!origin->computeOctree(progressCb)) + { + ccLog::Warning("[CorePoints::prepare] Failed to compute the octree"); + return false; + } + } + + CCLib::CloudSamplingTools::SFModulationParams modParams; + modParams.enabled = false; + ref = CCLib::CloudSamplingTools::resampleCloudSpatially( + origin, + static_cast(selectionParam), + modParams, + origin->getOctree().data(), + progressCb); + + break; + } + + case RANDOM: + { + if (selectionParam <= 0.0 || selectionParam >= 1.0) + { + ccLog::Warning("[CorePoints::prepare] Random subsampling ration must be between 0 and 1 (excluded)"); + return false; + } + int targetCount = static_cast(origin->size() * selectionParam); + ref = CCLib::CloudSamplingTools::subsampleCloudRandomly(origin, targetCount, progressCb); + break; + } + + case NONE: + //nothing to do + cloud = origin; + return true; + + default: + assert(false); + break; + } + + //store the references + if (!ref) + { + ccLog::Warning("[CorePoints::prepare] Failed to subsampled the origin cloud"); + return false; + } + selection.reset(ref); + + //and create the subsampled version of the cloud + cloud = origin->partialClone(ref); + if (!cloud) + { + ccLog::Warning("[CorePoints::prepare] Failed to subsampled the origin cloud (not enough memory)"); + return false; + } + + return true; +} diff --git a/CorePoints.h b/CorePoints.h new file mode 100644 index 0000000..658ae8c --- /dev/null +++ b/CorePoints.h @@ -0,0 +1,57 @@ +#pragma once + +//########################################################################## +//# # +//# CLOUDCOMPARE PLUGIN: q3DMASC # +//# # +//# This program is free software; you can redistribute it and/or modify # +//# it under the terms of the GNU General Public License as published by # +//# the Free Software Foundation; version 2 or later of the License. # +//# # +//# This program is distributed in the hope that it will be useful, # +//# but WITHOUT ANY WARRANTY; without even the implied warranty of # +//# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +//# GNU General Public License for more details. # +//# # +//# COPYRIGHT: Dimitri Lague / CNRS / UEB # +//# # +//########################################################################## + +//qCC_db +#include + +//CCLib +#include +#include + +//Qt +#include + +//! 3DMASC classifier +namespace masc +{ + //! Core points descriptor + struct CorePoints + { + //origin cloud + ccPointCloud* origin = nullptr; + + //core points cloud + ccPointCloud* cloud = nullptr; + + //! Return the size + inline unsigned size() const { return (cloud ? cloud->size() : 0); } + //! Return the point index + inline unsigned originIndex(unsigned i) const { return selection ? selection->getPointGlobalIndex(i) : i; } + + //selection (if any) + QSharedPointer selection; + enum SubSamplingMethod { NONE, RANDOM, SPATIAL }; + SubSamplingMethod selectionMethod = NONE; + double selectionParam = std::numeric_limits::quiet_NaN(); + + //! Prepares the selection (must be called once) + bool prepare(CCLib::GenericProgressCallback* progressCb = nullptr); + }; + +}; //namespace masc diff --git a/Features.h b/Features.h index bcd2823..c2f0cf7 100644 --- a/Features.h +++ b/Features.h @@ -19,6 +19,7 @@ //Local #include "FeaturesInterface.h" +#include "CorePoints.h" //qCC_db #include @@ -151,12 +152,12 @@ public: //methods //auomatically set the right source for specific features switch (type) { - case Z: - source = DimZ; + case X: + source = DimX; sourceName = "X"; break; - case Z: - source = DimZ; + case Y: + source = DimY; sourceName = "Y"; break; case Z: @@ -597,7 +598,7 @@ struct FeatureRule int sourceSFIndex = -1; //! Checks the rule validity - bool checkValidity(QString &error) const + bool checkValidity(/*const masc::CorePoints& corePoints, */QString &error) const { int cloudCount = (cloud1 ? (cloud2 ? 2 : 1) : 0); @@ -606,29 +607,50 @@ struct FeatureRule error = "feature rule has no associated feature"; return false; } - if (scales != nullptr && scales->values.empty()) + if (scales != nullptr) { - error = "invalid scales definition"; - return false; + if (scales->values.empty()) + { + error = "invalid scales definition"; + return false; + } + if (stat == FeatureRule::NO_STAT) + { + error = "scaled features need a STAT measure to be defined"; + return false; + } + } + else //no scales + { + //if (corePoints.origin != cloud1) + //{ + // error = "feature with no scale must be computed/extracted from the core points origin cloud"; + // return false; + //} } if (stat != FeatureRule::NO_STAT) { if (feature->getType() != Feature::Type::PointFeature) { - error = "stat. measures can only be defined on Point features"; + error = "STAT measures can only be defined on Point features"; return false; } if (!scales) { - error = "stat. measures need at least one scale to be defined"; + error = "STAT measures need at least one scale to be defined"; return false; } } - if (stat != FeatureRule::NO_OPERATION) + if (op != FeatureRule::NO_OPERATION) { + if (!scales) + { + error = "math operations can't be defined on scale-less features (SC0)"; + return false; + } if (feature->getType() == Feature::Type::DualCloudFeature) { - error = "math operation can't be defined on dual-cloud features"; + error = "math operations can't be defined on dual-cloud features"; return false; } if (cloudCount < 2) diff --git a/FeaturesInterface.h b/FeaturesInterface.h index 570bf7d..406f78b 100644 --- a/FeaturesInterface.h +++ b/FeaturesInterface.h @@ -63,9 +63,7 @@ struct Feature , source(p_source) , scale(p_scale) , sourceName(p_sourceName) - { - assert(cloud); - } + {} //! Associated cloud ccPointCloud* cloud; diff --git a/ScalarFieldWrappers.h b/ScalarFieldWrappers.h index 2f68bd6..fcd1bbd 100644 --- a/ScalarFieldWrappers.h +++ b/ScalarFieldWrappers.h @@ -27,6 +27,8 @@ class IScalarFieldWrapper public: virtual double pointValue(unsigned index) const = 0; virtual bool isValid() const = 0; + virtual QString getName() const = 0; + virtual unsigned size() const = 0; }; class ScalarFieldWrapper : public IScalarFieldWrapper @@ -38,11 +40,64 @@ public: virtual inline double pointValue(unsigned index) const override { return m_sf->at(index); } virtual inline bool isValid() const { return m_sf != nullptr; } + virtual inline QString getName() const { return m_sf->getName(); } + virtual unsigned size() const { return m_sf->size(); } protected: CCLib::ScalarField* m_sf; }; +class ScalarFieldRatioWrapper : public IScalarFieldWrapper +{ +public: + ScalarFieldRatioWrapper(CCLib::ScalarField* sfp, CCLib::ScalarField* sfq, QString name) + : m_sfp(sfp) + , m_sfq(sfq) + , m_name(name) + {} + + virtual inline double pointValue(unsigned index) const override + { + ScalarType p = m_sfp->getValue(index); + ScalarType q = m_sfq->getValue(index); + ScalarType ratio = (std::abs(q) > std::numeric_limits::epsilon() ? p / q : NAN_VALUE); + return ratio; + } + virtual inline bool isValid() const { return (m_sfp != nullptr && m_sfq != nullptr); } + virtual inline QString getName() const { return m_name; } + virtual inline unsigned size() const { return std::min(m_sfp->size(), m_sfq->size()); } + +protected: + CCLib::ScalarField *m_sfp, *m_sfq; + QString m_name; +}; + +class NormDipAndDipDirFieldWrapper : public IScalarFieldWrapper +{ +public: + enum Mode { Dip = 0, DipDir = 1 }; + + NormDipAndDipDirFieldWrapper(ccPointCloud* cloud, Mode mode) + : m_cloud(cloud) + , m_mode(mode) + {} + + virtual double pointValue(unsigned index) const override + { + const CCVector3& N = m_cloud->getPointNormal(index); + PointCoordinateType dip_deg, dipDir_deg; + ccNormalVectors::ConvertNormalToDipAndDipDir(N, dip_deg, dipDir_deg); + return (m_mode == Dip ? dip_deg : dipDir_deg); + } + virtual inline bool isValid() const { return m_cloud != nullptr && m_cloud->hasNormals(); } + virtual inline QString getName() const { static const char s_names[][14] = { "Norm dip", "Norm dip dir." }; return s_names[m_mode]; } + virtual inline unsigned size() const { return m_cloud->size(); } + +protected: + ccPointCloud* m_cloud; + Mode m_mode; +}; + class DimScalarFieldWrapper : public IScalarFieldWrapper { public: @@ -55,6 +110,8 @@ public: virtual inline double pointValue(unsigned index) const override { return m_cloud->getPoint(index)->u[m_dim]; } virtual inline bool isValid() const { return m_cloud != nullptr; } + virtual inline QString getName() const { static const char s_names[][5] = { "DimX", "DimY", "DimZ" }; return s_names[m_dim]; } + virtual inline unsigned size() const { return m_cloud->size(); } protected: ccPointCloud* m_cloud; @@ -73,6 +130,8 @@ public: virtual inline double pointValue(unsigned index) const override { return m_cloud->getPointColor(index).rgb[m_band]; } virtual inline bool isValid() const { return m_cloud != nullptr && m_cloud->hasColors(); } + virtual inline QString getName() const { static const char s_names[][6] = { "Red", "Green", "Blue" }; return s_names[m_band]; } + virtual inline unsigned size() const { return m_cloud->size(); } protected: ccPointCloud* m_cloud; diff --git a/q3DMASC.cpp b/q3DMASC.cpp index 237ce01..b3de451 100644 --- a/q3DMASC.cpp +++ b/q3DMASC.cpp @@ -25,11 +25,13 @@ //qCC_db #include +#include //Qt #include #include #include +#include q3DMASCPlugin::q3DMASCPlugin(QObject* parent/*=0*/) : QObject(parent) @@ -49,7 +51,8 @@ void q3DMASCPlugin::onNewSelection(const ccHObject::Container& selectedEntities) if (m_trainAction) { - m_trainAction->setEnabled(m_app && m_app->dbRootObject() && m_app->dbRootObject()->getChildrenNumber() != 0); //need some loaded entities to train the classifier! + //m_trainAction->setEnabled(m_app && m_app->dbRootObject() && m_app->dbRootObject()->getChildrenNumber() != 0); //need some loaded entities to train the classifier! + m_trainAction->setEnabled(true); } m_selectedEntities = selectedEntities; @@ -114,14 +117,6 @@ void q3DMASCPlugin::doTrainAction() //ccPointCloud* cloud1 = static_cast(m_selectedEntities[0]); //ccPointCloud* cloud2 = static_cast(m_selectedEntities[1]); - if (m_selectedEntities.empty() || !m_selectedEntities.front()->isA(CC_TYPES::POINT_CLOUD)) - { - m_app->dispToConsole("Select one and only one point cloud!", ccMainAppInterface::ERR_CONSOLE_MESSAGE); - return; - } - - ccPointCloud* cloud = static_cast(m_selectedEntities.front()); - masc::TrainParameters params; if (params.testDataRatio < 0 || params.testDataRatio > 0.99f) { @@ -130,33 +125,129 @@ void q3DMASCPlugin::doTrainAction() } Feature::Set features; +#if 0 + if (m_selectedEntities.empty() || !m_selectedEntities.front()->isA(CC_TYPES::POINT_CLOUD)) { - features.push_back(Feature::Shared(new PointFeature(cloud, PointFeature::Z, Feature::DimZ, "Z"))); - features.push_back(Feature::Shared(new PointFeature(cloud, PointFeature::Intensity, Feature::ScalarField, "Intensity"))); + m_app->dispToConsole("Select one and only one point cloud!", ccMainAppInterface::ERR_CONSOLE_MESSAGE); + return; } - QString outputFilename = QCoreApplication::applicationDirPath() + "/classifier.yaml"; + ccPointCloud* cloud = static_cast(m_selectedEntities.front()); + + //features + { + Feature::Shared featureZ(new PointFeature(PointFeature::Z, cloud)); + features.push_back(featureZ); + + Feature::Shared featureIntensity(new PointFeature(PointFeature::Intensity, cloud)); + features.push_back(featureIntensity); + } +#else + QString inputFilename; + { + QSettings settings; + settings.beginGroup("3DMASC"); + QString inputPath = settings.value("FilePath", QCoreApplication::applicationDirPath()).toString(); + inputFilename = QFileDialog::getOpenFileName(m_app->getMainWindow(), "Load 3DMASC script file", inputPath, "*.txt"); + if (inputFilename.isNull()) + { + //process cancelled by the user + return; + } + settings.setValue("FilePath", QFileInfo(inputFilename).absolutePath()); + settings.endGroup(); + } + + FeatureRule::Set rules; + std::vector loadedClouds; + masc::CorePoints corePoints; + if (!masc::Tools::LoadFile(inputFilename, rules, loadedClouds, corePoints)) + { + while (!loadedClouds.empty()) + { + delete loadedClouds.back(); + loadedClouds.pop_back(); + } + return; + } + + //add the loaded clouds to the main DB (so that we don't need to handle them anymore) + ccHObject* group = new ccHObject("3DMASC"); + for (ccPointCloud* pc : loadedClouds) + { + group->addChild(pc); + } + + ccProgressDialog pDlg(true, m_app->getMainWindow()); + if (!corePoints.prepare(&pDlg)) + { + m_app->dispToConsole("Failed to compute/prepare the core points!", ccMainAppInterface::ERR_CONSOLE_MESSAGE); + delete group; + return; + } + if (corePoints.cloud != corePoints.origin) + { + //auto-hide the other clouds + for (ccPointCloud* pc : loadedClouds) + { + pc->setEnabled(false); + } + //set an explicit name for the core points + QString corePointsName = corePoints.origin->getName(); + switch (corePoints.selectionMethod) + { + case masc::CorePoints::NONE: + break; + case masc::CorePoints::RANDOM: + corePointsName += "_SS_Random@" + QString::number(corePoints.selectionParam); + break; + case masc::CorePoints::SPATIAL: + corePointsName += "_SS_Spatial@" + QString::number(corePoints.selectionParam); + break; + default: + assert(false); + } + corePoints.cloud->setName(QString("Core points (%1)").arg(corePointsName)); + group->addChild(corePoints.cloud); + } + m_app->addToDB(group); + QCoreApplication::processEvents(); + + pDlg.setAutoClose(false); //we don't want the progress dialog to 'pop' for each feature + QString error; + if (!masc::Tools::PrepareFeatures(rules, corePoints, features, error, &pDlg)) + { + m_app->dispToConsole(error, ccMainAppInterface::ERR_CONSOLE_MESSAGE); + delete group; + return; + } + pDlg.close(); + QCoreApplication::processEvents(); + pDlg.setAutoClose(true); //restore the default behavior of the progress dialog + +#endif //randomly select the training points - QScopedPointer trainSubset(new CCLib::ReferenceCloud(cloud)); - QScopedPointer testSubset(new CCLib::ReferenceCloud(cloud)); - if (!masc::Tools::RandomSubset(cloud, params.testDataRatio, trainSubset.data(), testSubset.data())) + QScopedPointer trainSubset(new CCLib::ReferenceCloud(corePoints.cloud)); + QScopedPointer testSubset(new CCLib::ReferenceCloud(corePoints.cloud)); + if (!masc::Tools::RandomSubset(corePoints.cloud, params.testDataRatio, trainSubset.data(), testSubset.data())) { m_app->dispToConsole("Not enough memory", ccMainAppInterface::ERR_CONSOLE_MESSAGE); return; } masc::Classifier classifier; - if (QFile(outputFilename).exists()) - { - if (!classifier.fromFile(outputFilename, m_app->getMainWindow())) - { - m_app->dispToConsole("Failed to load previous classifier file", ccMainAppInterface::ERR_CONSOLE_MESSAGE); - return; - } - m_app->dispToConsole("Previous classifier loaded", ccMainAppInterface::WRN_CONSOLE_MESSAGE); - } - else + //QString outputFilename = QCoreApplication::applicationDirPath() + "/classifier.yaml"; + //if (QFile(outputFilename).exists()) + //{ + // if (!classifier.fromFile(outputFilename, m_app->getMainWindow())) + // { + // m_app->dispToConsole("Failed to load previous classifier file", ccMainAppInterface::ERR_CONSOLE_MESSAGE); + // return; + // } + // m_app->dispToConsole("Previous classifier loaded", ccMainAppInterface::WRN_CONSOLE_MESSAGE); + //} + //else { QString errorMessage; if (!classifier.train(params.rt, features, errorMessage, trainSubset.data(), m_app->getMainWindow())) @@ -165,20 +256,38 @@ void q3DMASCPlugin::doTrainAction() return; } + QString outputFilename; + { + QSettings settings; + settings.beginGroup("3DMASC"); + QString outputPath = settings.value("FilePath", QCoreApplication::applicationDirPath()).toString(); + outputFilename = QFileDialog::getSaveFileName(m_app->getMainWindow(), "Save 3DMASC classifier", outputPath, "*.yaml"); + if (outputFilename.isNull()) + { + //process cancelled by the user + return; + } + settings.setValue("FilePath", QFileInfo(outputFilename).absolutePath()); + settings.endGroup(); + } + //save the classifier classifier.toFile(outputFilename, m_app->getMainWindow()); m_app->dispToConsole("Classifier succesfully created", ccMainAppInterface::WRN_CONSOLE_MESSAGE); } - masc::Classifier::AccuracyMetrics metrics; - QString errorMessage; - if (!classifier.evaluate(features, testSubset.data(), metrics, errorMessage, m_app->getMainWindow())) + //test classifier { - m_app->dispToConsole(errorMessage, ccMainAppInterface::ERR_CONSOLE_MESSAGE); - return; - } + masc::Classifier::AccuracyMetrics metrics; + QString errorMessage; + if (!classifier.evaluate(features, testSubset.data(), metrics, errorMessage, m_app->getMainWindow())) + { + m_app->dispToConsole(errorMessage, ccMainAppInterface::ERR_CONSOLE_MESSAGE); + return; + } - m_app->dispToConsole(QString("Correct = %1 / %2 --> accuracy = %3").arg(metrics.goodGuess).arg(metrics.sampleCount).arg(metrics.ratio), ccMainAppInterface::STD_CONSOLE_MESSAGE); + m_app->dispToConsole(QString("Correct = %1 / %2 --> accuracy = %3").arg(metrics.goodGuess).arg(metrics.sampleCount).arg(metrics.ratio), ccMainAppInterface::STD_CONSOLE_MESSAGE); + } } void q3DMASCPlugin::registerCommands(ccCommandLineInterface* cmd) diff --git a/q3DMASCTools.cpp b/q3DMASCTools.cpp index 9dcd637..7287415 100644 --- a/q3DMASCTools.cpp +++ b/q3DMASCTools.cpp @@ -25,40 +25,43 @@ #include //qCC_db #include - -//CCLib -#include +#include //Qt #include #include +#include +#include //system #include using namespace masc; -bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, FeatureRule::Set& features) +bool Tools::LoadFile( QString filename, + FeatureRule::Set& features, + std::vector& loadedClouds, + CorePoints& corePoints) { - QFile file(filename); - if (!file.exists()) + QFileInfo fi(filename); + if (!fi.exists()) { ccLog::Warning(QString("Can't find file '%1'").arg(filename)); return false; } + + QFile file(filename); if (!file.open(QFile::Text | QFile::ReadOnly)) { ccLog::Warning(QString("Can't open file '%1'").arg(filename)); return false; } - Scales::Shared scales(new Scales); assert(features.empty()); - - QMap > clouds; + Scales::Shared scales(new Scales); + QMap clouds; QTextStream stream(&file); - for (int lineNumber = 0; ; ++lineNumber) { QString line = stream.readLine(); @@ -84,14 +87,14 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea if (upperLine.startsWith("CLOUD:")) //clouds { QString command = line.mid(6); - QStringList tokens = command.split(':'); + QStringList tokens = command.split('='); if (tokens.size() != 2) { ccLog::Warning("Malformed file: expecting 2 tokens after 'cloud:' on line #" + QString::number(lineNumber)); return false; } - QString pcName = tokens[0]; - QString pcFilename = tokens[1]; + QString pcName = tokens[0].trimmed(); + QString pcFilename = fi.absoluteDir().absoluteFilePath(tokens[1].trimmed()); //try to open the cloud { FileIOFilter::LoadParameters parameters; @@ -100,18 +103,92 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea ccHObject* object = FileIOFilter::LoadFromFile(pcFilename, parameters, error); if (error != CC_FERR_NO_ERROR || !object) { - ccLog::Warning("Failed to open the file (see console)"); + //error message already issued if (object) delete object; return false; } - if (!object->isA(CC_TYPES::POINT_CLOUD)) + ccHObject::Container cloudsInFile; + object->filterChildren(cloudsInFile, false, CC_TYPES::POINT_CLOUD, true); + if (cloudsInFile.empty()) { ccLog::Warning("File doesn't contain a single cloud"); delete object; return false; } - clouds.insert(pcName, QSharedPointer(static_cast(object))); + else if (cloudsInFile.size() > 1) + { + ccLog::Warning("File contains more than one cloud, only the first one will be kept"); + } + ccPointCloud* pc = static_cast(cloudsInFile.front()); + for (size_t i = 1; i < cloudsInFile.size(); ++i) + { + delete cloudsInFile[i]; + } + if (pc->getParent()) + pc->getParent()->detachChild(pc); + pc->setName(pcName); + clouds.insert(pcName, pc); + loadedClouds.push_back(pc); + } + } + else if (upperLine.startsWith("CORE_POINTS:")) //core points + { + if (corePoints.origin) + { + ccLog::Warning("Malformed file: can't declare core points twice! (line #" + QString::number(lineNumber) + ")"); + return false; + } + QString command = line.mid(12); + QStringList tokens = command.split('_'); + if (tokens.empty()) + { + ccLog::Warning("Malformed file: expecting tokens after 'core_points:' on line #" + QString::number(lineNumber)); + return false; + } + QString pcName = tokens[0].trimmed(); + if (!clouds.contains(pcName)) + { + ccLog::Warning(QString("Malformed file: unknown cloud '%1' on line #%2 (make sure it is declared before the core points)").arg(pcName).arg(lineNumber)); + return false; + } + corePoints.origin = clouds[pcName]; + + //should we sub-sample the origin cloud? + if (tokens.size() > 1) + { + if (tokens[1].toUpper() == "SS") + { + if (tokens.size() < 3) + { + ccLog::Warning("Malformed file: missing token after 'SS' on line #" + QString::number(lineNumber)); + return false; + } + QString options = tokens[2]; + if (options.startsWith('R')) + { + corePoints.selectionMethod = CorePoints::RANDOM; + } + else if (options.startsWith('S')) + { + corePoints.selectionMethod = CorePoints::SPATIAL; + } + else + { + ccLog::Warning("Malformed file: unknown option after 'SS' on line #" + QString::number(lineNumber)); + return false; + } + + //read the subsampling parameter (ignore the first character) + bool ok = false; + corePoints.selectionParam = options.mid(1).toDouble(&ok); + if (!ok) + { + ccLog::Warning("Malformed file: expecting a number after 'SS_X' on line #" + QString::number(lineNumber)); + return false; + } + + } //end of subsampling options } } else if (upperLine.startsWith("SCALES:")) //scales @@ -131,14 +208,14 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea if (token.contains(':')) { //it's probably a range - QStringList subTokens = token.split(':'); + QStringList subTokens = token.trimmed().split(':'); if (subTokens.size() != 3) { ccLog::Warning(QString("Malformed file: expecting 3 tokens for a range of scales (%1)").arg(token)); return false; } bool ok[3] = { true, true, true }; - double start = subTokens[0].toDouble(ok); + double start = subTokens[0].trimmed().toDouble(ok); double step = subTokens[1].toDouble(ok + 1); double stop = subTokens[2].toDouble(ok + 2); if (!ok[0] || !ok[1] || !ok[2]) @@ -151,7 +228,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea ccLog::Warning(QString("Malformed file: invalid range (%1) on line #%2").arg(token).arg(lineNumber)); return false; } - for (double v = start; v <= stop + 1.0 - 6; v += step) + for (double v = start; v <= stop + 1.0e-6; v += step) { scales->values.push_back(v); } @@ -159,7 +236,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea else { bool ok = true; - double v = token.toDouble(&ok); + double v = token.trimmed().toDouble(&ok); if (!ok) { ccLog::Warning(QString("Malformed file: invalid scale value (%1) on line #%2").arg(token).arg(lineNumber)); @@ -188,7 +265,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea FeatureRule::Shared rule(new FeatureRule); //read the type - QString typeStr = tokens[0].toUpper(); + QString typeStr = tokens[0].trimmed().toUpper(); { for (int iteration = 0; iteration < 1; ++iteration) //fake loop for easy break { @@ -209,7 +286,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea ccLog::Warning(QString("Malformed file: expecting a valid integer value after 'SF' on line #%1").arg(lineNumber)); return false; } - pointFeature->sourceSFIndex = sfIndex; + rule->sourceSFIndex = sfIndex; } rule->feature = pointFeature; @@ -264,7 +341,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea else { //read the specific scale index - QString scaleStr = scaleStr.mid(2); + scaleStr = scaleStr.mid(2); bool ok = true; double scale = scaleStr.toDouble(&ok); if (!ok) @@ -284,7 +361,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea bool mathDefined = false; for (int i = 2; i < tokens.size(); ++i) { - QString token = tokens[i].toUpper(); + QString token = tokens[i].trimmed().toUpper(); //is the token a 'stat' one? if (!statDefined) @@ -325,14 +402,15 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea if (cloudCount < 2) { bool cloudNameMatches = false; - for (QMap >::const_iterator it = clouds.begin(); it != clouds.end(); ++it) + for (QMap::const_iterator it = clouds.begin(); it != clouds.end(); ++it) { - if (it.key().toUpper() == token) + QString key = it.key().toUpper(); + if (key == token) { if (cloudCount == 0) - rule->cloud1 = it.value().data(); + rule->cloud1 = it.value(); else - rule->cloud2 = it.value().data(); + rule->cloud2 = it.value(); ++cloudCount; cloudNameMatches = true; break; @@ -400,7 +478,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea assert(rule && rule->feature); QString errorMessage; - bool ruleIsValid = rule->checkValidity(errorMessage); + bool ruleIsValid = rule->checkValidity(/*corePoints, */errorMessage); if (!ruleIsValid) { ccLog::Warning("Malformed feature: " + errorMessage + QString("(line %1)").arg(lineNumber)); @@ -462,7 +540,7 @@ static const char* s_PCVSFName = "Illuminance (PCV)"; static const char* s_normDipSFName = "Norm dip"; static const char* s_normDipDirSFName = "Norm dip dir."; -static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType featureType, int sourceSFIndex, ccPointCloud* cloud, QString& error) +static QSharedPointer RetrieveField(PointFeature::PointFeatureType featureType, int sourceSFIndex, ccPointCloud* cloud, QString& error) { QString sfName; switch (featureType) @@ -475,14 +553,14 @@ static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType fe error = "Cloud has no 'intensity' scalar field"; return nullptr; } - return sf; + return QSharedPointer(new ScalarFieldWrapper(sf)); } case PointFeature::X: + return QSharedPointer(new DimScalarFieldWrapper(cloud, DimScalarFieldWrapper::DimX)); case PointFeature::Y: + return QSharedPointer(new DimScalarFieldWrapper(cloud, DimScalarFieldWrapper::DimY)); case PointFeature::Z: - //not a ScalarField source - error = "Internal error (source is not a scalar field)"; - return nullptr; + return QSharedPointer(new DimScalarFieldWrapper(cloud, DimScalarFieldWrapper::DimZ)); case PointFeature::NbRet: { CCLib::ScalarField* sf = RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_NUMBER_OF_RETURNS], false); @@ -491,7 +569,7 @@ static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType fe error = "Cloud has no 'number of returns' scalar field"; return nullptr; } - return sf; + return QSharedPointer(new ScalarFieldWrapper(sf)); } case PointFeature::RetNb: { @@ -501,17 +579,11 @@ static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType fe error = "Cloud has no 'return number' scalar field"; return nullptr; } - return sf; + return QSharedPointer(new ScalarFieldWrapper(sf)); } case PointFeature::EchoRat: { - CCLib::ScalarField* _echoRatioSF = RetrieveSF(cloud, s_echoRatioSFName, true); - if (_echoRatioSF) - { - //SF was already computed? - return _echoRatioSF; - } - //otherwise we need to compute it + //retrieve the two scalar fields 'p/q' CCLib::ScalarField* numberOfRetSF = RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_NUMBER_OF_RETURNS], false); if (!numberOfRetSF) { @@ -529,31 +601,14 @@ static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType fe error = "Internal error (inconsistent scalar fields)"; return nullptr; } - ccScalarField* echoRatioSF = new ccScalarField(s_echoRatioSFName); - if (!echoRatioSF->reserveSafe(retNumberSF->size())) - { - error = "Not enough memory"; - echoRatioSF->release(); - return nullptr; - } - - for (unsigned i = 0; i < cloud->size(); ++i) - { - ScalarType p = retNumberSF->getValue(i); - ScalarType q = numberOfRetSF->getValue(i); - ScalarType ratio = (std::abs(q) > std::numeric_limits::epsilon() ? p / q : NAN_VALUE); - echoRatioSF->addElement(ratio); - } - echoRatioSF->computeMinAndMax(); - cloud->addScalarField(echoRatioSF); - return echoRatioSF; + return QSharedPointer(new ScalarFieldRatioWrapper(retNumberSF, numberOfRetSF, "EchoRat")); } case PointFeature::R: + return QSharedPointer(new ColorScalarFieldWrapper(cloud, ColorScalarFieldWrapper::Red)); case PointFeature::G: + return QSharedPointer(new ColorScalarFieldWrapper(cloud, ColorScalarFieldWrapper::Green)); case PointFeature::B: - //not a ScalarField source - error = "Internal error (source is not a scalar field)"; - return nullptr; + return QSharedPointer(new ColorScalarFieldWrapper(cloud, ColorScalarFieldWrapper::Blue)); case PointFeature::NIR: { CCLib::ScalarField* sf = RetrieveSF(cloud, s_NIRSFName, false); @@ -562,46 +617,18 @@ static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType fe error = "Cloud has no 'NIR' scalar field"; return nullptr; } - return sf; + return QSharedPointer(new ScalarFieldWrapper(sf)); } case PointFeature::DipAng: case PointFeature::DipDir: { - CCLib::ScalarField* _dipSF = RetrieveSF(cloud, (featureType == PointFeature::DipAng ? s_normDipSFName : s_normDipDirSFName), true); - if (_dipSF) - { - //SF was already computed? - return _dipSF; - } - //otherwise we need to compute it - - static const char* s_normDipSFName = "Norm dip"; - static const char* s_normDipDirSFName = "Norm dip dir."; - //we need normals to cumpute Dip and Dip Dir. angles! + //we need normals to compute the dip and dip direction! if (!cloud->hasNormals()) { error = "Cloud has no normals: can't compute dip or dip dir. angles"; return nullptr; } - - ccScalarField* dipSF = new ccScalarField(featureType == PointFeature::DipAng ? s_normDipSFName : s_normDipDirSFName); - if (!dipSF->reserveSafe(cloud->size())) - { - error = "Not enough memory"; - dipSF->release(); - return nullptr; - } - - for (unsigned i = 0; i < cloud->size(); ++i) - { - const CCVector3& N = cloud->getPointNormal(i); - PointCoordinateType dip_deg, dipDir_deg; - ccNormalVectors::ConvertNormalToDipAndDipDir(N, dip_deg, dipDir_deg); - dipSF->addElement(static_cast(featureType == PointFeature::DipAng ? dip_deg : dipDir_deg)); - } - dipSF->computeMinAndMax(); - cloud->addScalarField(dipSF); - return dipSF; + return QSharedPointer(new NormDipAndDipDirFieldWrapper(cloud, featureType == PointFeature::DipAng ? NormDipAndDipDirFieldWrapper::Dip : NormDipAndDipDirFieldWrapper::DipDir)); } case PointFeature::M3C2: { @@ -611,7 +638,7 @@ static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType fe error = "Cloud has no 'm3c2 distance' scalar field"; return nullptr; } - return sf; + return QSharedPointer(new ScalarFieldWrapper(sf)); } case PointFeature::PCV: { @@ -621,7 +648,7 @@ static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType fe error = "Cloud has no 'PCV/Illuminance' scalar field"; return nullptr; } - return sf; + return QSharedPointer(new ScalarFieldWrapper(sf)); } case PointFeature::SF: if (sourceSFIndex < 0 || sourceSFIndex >= static_cast(cloud->getNumberOfScalarFields())) @@ -629,7 +656,7 @@ static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType fe error = QString("Can't retrieve the specified SF: invalid index (%1)").arg(sourceSFIndex); return nullptr; } - return cloud->getScalarField(sourceSFIndex); + return QSharedPointer(new ScalarFieldWrapper(cloud->getScalarField(sourceSFIndex))); default: break; } @@ -638,71 +665,48 @@ static CCLib::ScalarField* RetrieveOrComputeSF(PointFeature::PointFeatureType fe return nullptr; } -static bool ExtractStatFromSF( const CCLib::DgmOctree::octreeCell& cell, - void** additionalParameters, - CCLib::NormalizedProgress* nProgress = nullptr) +static bool ExtractStatFromSF( const CCVector3& queryPoint, + const CCLib::DgmOctree* octree, + unsigned char octreeLevel, + FeatureRule::Stat stat, + const IScalarFieldWrapper& inputField, + PointCoordinateType radius, + double& outputValue) { - //additional parameters - FeatureRule::Stat stat = *reinterpret_cast (additionalParameters[0]); - CCLib::ScalarField* inputSF = reinterpret_cast (additionalParameters[1]); - CCLib::ScalarField* resultSF = reinterpret_cast (additionalParameters[2]); - PointCoordinateType radius = *reinterpret_cast(additionalParameters[3]); - assert(inputSF && resultSF); - - //number of points inside the current cell - unsigned n = cell.points->size(); + if (!octree) + { + assert(false); + return false; + } + outputValue = std::numeric_limits::quiet_NaN(); //spherical neighborhood extraction structure CCLib::DgmOctree::NearestNeighboursSphericalSearchStruct nNSS; - nNSS.level = cell.level; - nNSS.prepare(radius, cell.parentOctree->getCellSize(nNSS.level)); - cell.parentOctree->getCellPos(cell.truncatedCode, cell.level, nNSS.cellPos, true); - cell.parentOctree->computeCellCenter(nNSS.cellPos, cell.level, nNSS.cellCenter); - - //we already know the points inside the current cell { - try - { - nNSS.pointsInNeighbourhood.resize(n); - } - catch (.../*const std::bad_alloc&*/) //out of memory - { - return false; - } - CCLib::DgmOctree::NeighboursSet::iterator it = nNSS.pointsInNeighbourhood.begin(); - for (unsigned j = 0; j < n; ++j, ++it) - { - it->point = cell.points->getPointPersistentPtr(j); - it->pointIndex = cell.points->getPointGlobalIndex(j); - } - nNSS.alreadyVisitedNeighbourhoodSize = 1; + nNSS.level = octreeLevel; + nNSS.queryPoint = queryPoint; + nNSS.prepare(radius, octree->getCellSize(nNSS.level)); + octree->getTheCellPosWhichIncludesThePoint(&nNSS.queryPoint, nNSS.cellPos, nNSS.level); + octree->computeCellCenter(nNSS.cellPos, nNSS.level, nNSS.cellCenter); } - for (unsigned i = 0; i < n; ++i) + //we extract the point's neighbors + unsigned kNN = octree->findNeighborsInASphereStartingFromCell(nNSS, radius, true); + if (kNN == 0) { - //retrieve the points around the current cell point - cell.points->getPoint(i, nNSS.queryPoint); + return true; + } - //we extract the point's neighbors - //warning: there may be more points at the end of nNSS.pointsInNeighbourhood than the actual nearest neighbors (k)! - unsigned kNN = cell.parentOctree->findNeighborsInASphereStartingFromCell(nNSS, radius, true); - if (kNN == 0) - { - assert(false); - continue; - } - - double sum = 0.0; - double sum2 = 0.0; - ScalarType minValue = 0; - ScalarType maxValue = 0; - bool withMode = (stat == FeatureRule::MODE || stat == FeatureRule::SKEW); - QMap modeCounter; + //specific case + if (stat == FeatureRule::RANGE) + { + double minValue = 0; + double maxValue = 0; for (unsigned k = 0; k < kNN; ++k) { unsigned index = nNSS.pointsInNeighbourhood[k].pointIndex; - ScalarType v = inputSF->getValue(index); + double v = inputField.pointValue(index); //track min and max values if (k != 0) @@ -716,96 +720,114 @@ static bool ExtractStatFromSF( const CCLib::DgmOctree::octreeCell& cell, { minValue = maxValue = v; } + } + outputValue = maxValue - minValue; + return true; + } + + bool withSums = (stat == FeatureRule::MEAN || stat == FeatureRule::STD || stat == FeatureRule::SKEW); + bool withMode = (stat == FeatureRule::MODE || stat == FeatureRule::SKEW); + double sum = 0.0; + double sum2 = 0.0; + QMap modeCounter; + + for (unsigned k = 0; k < kNN; ++k) + { + unsigned index = nNSS.pointsInNeighbourhood[k].pointIndex; + double v = inputField.pointValue(index); + + if (withSums) + { //compute average and std. dev. sum += v; - sum2 += static_cast(v) * v; - - if (withMode) - { - if (modeCounter.contains(v)) - { - ++modeCounter[v]; - } - else - { - modeCounter[v] = 1; - } - } + sum2 += v * v; } - double mode = NAN_VALUE; if (withMode) { - int maxCounter = 0; - //look for the value with the highest frequency - for (QMap::const_iterator it = modeCounter.begin(); it != modeCounter.end(); ++it) + //store the number of occurences of each value + //DGM TODO: it would be better with a custom 'resolution' if the field is not an integer one + float vf = static_cast(v); + if (modeCounter.contains(vf)) { - if (it.value() > maxCounter) - { - maxCounter = it.value(); - mode = it.key(); - } + ++modeCounter[vf]; + } + else + { + modeCounter[vf] = 1; } } + } - ScalarType outValue = NAN_VALUE; - switch (stat) + double mode = std::numeric_limits::quiet_NaN(); + if (withMode) + { + //look for the value with the highest frequency + unsigned maxCounter = 0; + for (QMap::const_iterator it = modeCounter.begin(); it != modeCounter.end(); ++it) { - case FeatureRule::MEAN: - outValue = static_cast(sum / kNN); - break; - case FeatureRule::MODE: - outValue = static_cast(mode); - break; - case FeatureRule::STD: - outValue = static_cast(sqrt(std::abs(sum2 * kNN - sum * sum)) / kNN); - break; - case FeatureRule::RANGE: - outValue = maxValue - minValue; - break; - case FeatureRule::SKEW: + if (it.value() > maxCounter) { - double mean = sum / kNN; - double std = sqrt(std::abs(sum2 / kNN - mean * mean)); - if (std > std::numeric_limits::epsilon()) //arbitrary epsilon - { - outValue = static_cast((mean - mode) / std); - } - break; + maxCounter = it.value(); + mode = it.key(); } - default: - assert(false); - break; } - resultSF->setValue(cell.points->getPointGlobalIndex(i), outValue); + } - if (nProgress && !nProgress->oneStep()) + switch (stat) + { + case FeatureRule::MEAN: + outputValue = sum / kNN; + break; + case FeatureRule::MODE: + outputValue = mode; + break; + case FeatureRule::STD: + outputValue = sqrt(std::abs(sum2 * kNN - sum * sum)) / kNN; + break; + case FeatureRule::RANGE: + //we can't be here + assert(false); + return false; + case FeatureRule::SKEW: + { + double mean = sum / kNN; + double std = sqrt(std::abs(sum2 / kNN - mean * mean)); + if (std > std::numeric_limits::epsilon()) //arbitrary epsilon { - return false; + outputValue = (mean - mode) / std; } + break; + } + default: + ccLog::Warning("Unhandled STAT measure"); + assert(false); + return false; } return true; } -static CCLib::ScalarField* ExtractStat( ccPointCloud* cloud, - CCLib::ScalarField* sf, +static CCLib::ScalarField* ExtractStat( const CorePoints& corePoints, + ccPointCloud* sourceCloud, + const IScalarFieldWrapper* sourceField, double scale, FeatureRule::Stat stat, + const char* resultSFName, CCLib::GenericProgressCallback* progressCb = nullptr) { - if (!cloud || !sf || scale <= 0.0 || stat == FeatureRule::NO_STAT) + if (!corePoints.cloud || !sourceCloud || !sourceField || scale <= 0.0 || stat == FeatureRule::NO_STAT || !resultSFName) { //invalid input parameters assert(false); return nullptr; } - ccOctree::Shared octree = cloud->getOctree(); + ccOctree::Shared octree = sourceCloud->getOctree(); if (!octree) { - octree = cloud->computeOctree(progressCb); + octree = sourceCloud->computeOctree(progressCb); if (!octree) { ccLog::Warning("Failed to compute octree"); @@ -814,16 +836,15 @@ static CCLib::ScalarField* ExtractStat( ccPointCloud* cloud, } CCLib::ScalarField* resultSF = nullptr; - QString resultSFName = sf->getName() + QString("_") + FeatureRule::StatToString(stat) + "_" + QString::number(scale); - int sfIdx = cloud->getScalarFieldIndexByName(qPrintable(resultSFName)); + int sfIdx = corePoints.cloud->getScalarFieldIndexByName(resultSFName); if (sfIdx >= 0) { - resultSF = cloud->getScalarField(sfIdx); + resultSF = corePoints.cloud->getScalarField(sfIdx); } else { - resultSF = new ccScalarField(qPrintable(resultSFName)); - if (!resultSF->reserveSafe(cloud->size())) + resultSF = new ccScalarField(resultSFName); + if (!resultSF->resizeSafe(corePoints.cloud->size())) { ccLog::Warning("Not enough memory"); resultSF->release(); @@ -835,135 +856,331 @@ static CCLib::ScalarField* ExtractStat( ccPointCloud* cloud, PointCoordinateType radius = static_cast(scale / 2); unsigned char octreeLevel = octree->findBestLevelForAGivenNeighbourhoodSizeExtraction(radius); //scale is the diameter! - //additionnal parameters - void* additionalParameters[] = { static_cast(&stat), - static_cast(&sf), - static_cast(&resultSF), - static_cast(&radius) - }; + unsigned pointCount = corePoints.size(); + progressCb->setInfo(qPrintable(QString("Computing field: %1\n(core points: %2)").arg(resultSFName).arg(pointCount))); + CCLib::NormalizedProgress nProgress(progressCb, pointCount); - if (octree->executeFunctionForAllCellsAtLevel( octreeLevel, - ExtractStatFromSF, - additionalParameters, - true, - progressCb, - qPrintable(QString("Extract stat @ scale %1").arg(scale))) == 0) + for (unsigned i = 0; i < pointCount; ++i) { - //something went wrong - ccLog::Warning("Process failed"); - resultSF->release(); - return nullptr; + double outputValue = 0; + if (!ExtractStatFromSF( *corePoints.cloud->getPoint(i), + octree.data(), + octreeLevel, + stat, + *sourceField, + radius, + outputValue)) + { + //unexpected error + resultSF->release(); + return nullptr; + } + + ScalarType v = static_cast(outputValue); + resultSF->setValue(i, v); + + if (progressCb && !nProgress.oneStep()) + { + //process cancelled by the user + ccLog::Warning("Process cancelled"); + resultSF->release(); + return nullptr; + } } resultSF->computeMinAndMax(); - cloud->addScalarField(static_cast(resultSF)); + int newSFIdx = corePoints.cloud->addScalarField(static_cast(resultSF)); + //update display + if (corePoints.cloud->getDisplay()) + { + corePoints.cloud->setCurrentDisplayedScalarField(newSFIdx); + corePoints.cloud->getDisplay()->redraw(); + } return resultSF; } - -static bool PreparePointBasedFeature(FeatureRule& rule, QString& error) +static bool PerformMathOp(CCLib::ScalarField* sf1, const CCLib::ScalarField* sf2, FeatureRule::Operation op) { - assert(rule.feature && rule.feature->getType() == Feature::Type::PointFeature); + if (!sf1 || !sf2 || sf1->size() != sf2->size() || op == FeatureRule::NO_OPERATION) + { + //invalid input parameters + return false; + } - PointFeature* feature = static_cast(rule.feature.data()); + for (unsigned i = 0; i < sf1->size(); ++i) + { + ScalarType s1 = sf1->getValue(i); + ScalarType s2 = sf2->getValue(i); + ScalarType s = NAN_VALUE; + switch (op) + { + case FeatureRule::MINUS: + s = s1 - s2; + break; + case FeatureRule::PLUS: + s = s1 + s2; + break; + case FeatureRule::DIVIDE: + if (std::abs(s2) > std::numeric_limits::epsilon()) + s = s1 / s2; + break; + case FeatureRule::MULTIPLY: + s = s1 * s2; + break; + default: + assert(false); + break; + } + sf1->setValue(i, s); + } + sf1->computeMinAndMax(); - std::vector preparedFeatures; + return true; +} - //look for the source field (and compute it if necessary) - CCLib::ScalarField* sf1 = RetrieveOrComputeSF(feature->type, rule.sourceSFIndex, rule.cloud1, error); - if (!sf1) +static Feature::Shared PreparePointBasedFeature(const FeatureRule& rule, + double scale, + const CorePoints& corePoints, + QString& error, + CCLib::GenericProgressCallback* progressCb = nullptr) +{ + if (!rule.cloud1 || !rule.feature || rule.feature->getType() != Feature::Type::PointFeature || !corePoints.cloud) + { + //invalid input + assert(false); + return false; + } + PointFeature::PointFeatureType featureType = static_cast(rule.feature.data())->type; + + //look for the source field + QSharedPointer field1 = RetrieveField(featureType, rule.sourceSFIndex, rule.cloud1, error); + if (!field1) { //error should be up to date return false; } - CCLib::ScalarField* sf2 = nullptr; - if (rule.cloud2 && rule.op != FeatureRule::NO_OPERATION) + //shall we extract a statistical measure? (= scaled feature) + if (std::isfinite(scale)) { - sf2 = RetrieveOrComputeSF(feature->type, rule.sourceSFIndex, rule.cloud2, error); - if (!sf2) + if (rule.stat == FeatureRule::NO_STAT) { - //error should be up to date + assert(false); + ccLog::Warning("Scaled features (SCx) must have an associated STAT measure"); return false; } - } - //shall we extract a statistical measure? - if (rule.scales && rule.stat != FeatureRule::NO_STAT) - { - //duplicate the feature for each scale - for (double s : rule.scales->values) + QSharedPointer field2; + if (rule.cloud2) { - CCLib::ScalarField* statSF1 = ExtractStat(rule.cloud1, sf1, s, rule.stat); - if (!statSF1) + //no need to compute the second scalar field if no MATH operation has to be performed?! + if (rule.op != FeatureRule::NO_OPERATION) { - ccLog::Warning(QString("Failed to extract stat. from sf '%1' @ scale %2").arg(sf1->getName()).arg(s)); - return false; - } - PointFeature::Shared f1(new PointFeature(*feature)); - f1->cloud = rule.cloud1; - f1->sourceName = statSF1->getName(); - preparedFeatures.push_back(f1); - - if (rule.cloud2 && sf2) - { - assert(rule.op != FeatureRule::NO_OPERATION); - CCLib::ScalarField* statSF2 = ExtractStat(rule.cloud2, sf2, s, rule.stat); - if (!statSF2) + field2 = RetrieveField(featureType, rule.sourceSFIndex, rule.cloud2, error); + if (!field2) { - ccLog::Warning(QString("Failed to extract stat. from sf '%1' @ scale %2").arg(sf2->getName()).arg(s)); + //error should be up to date return false; } - PointFeature::Shared f2(new PointFeature(*feature)); - f2->cloud = rule.cloud2; - f2->sourceName = statSF2->getName(); - preparedFeatures.push_back(f2); + } + else + { + assert(false); + ccLog::Warning("Feature has a second cloud associated but no MATH operation is defined"); } } - } - else - { - //only one version of the main feature - feature->cloud = rule.cloud1; - feature->sourceName = sf1->getName(); - preparedFeatures.push_back(rule.feature); - } - switch (feature->type) - { - case PointFeature::Intensity: - case PointFeature::X: - case PointFeature::Y: - case PointFeature::Z: - case PointFeature::NbRet: - case PointFeature::RetNb: - case PointFeature::EchoRat: - case PointFeature::R: - case PointFeature::G: - case PointFeature::B: - case PointFeature::NIR: - case PointFeature::DipAng: - case PointFeature::DipDir: - case PointFeature::M3C2: - case PointFeature::PCV: - case PointFeature::SF: - } + //build the final SF name + QString resultSFName = rule.cloud1->getName() + "." + field1->getName() + QString("_") + FeatureRule::StatToString(rule.stat); + if (field2 && rule.op != FeatureRule::NO_OPERATION) + { + //include the math operation as well if necessary! + resultSFName += "_" + FeatureRule::OpToString(rule.op) + "_" + rule.cloud2->getName() + "." + field2->getName() + QString("_") + FeatureRule::StatToString(rule.stat); + } + resultSFName += "@" + QString::number(scale); + CCLib::ScalarField* statSF1 = ExtractStat(corePoints, rule.cloud1, field1.data(), scale, rule.stat, qPrintable(resultSFName), progressCb); + if (!statSF1) + { + error = QString("Failed to extract stat. from field '%1' @ scale %2").arg(field1->getName()).arg(scale); + return false; + } + + PointFeature::Shared feature(new PointFeature(*static_cast(rule.feature.data()))); + feature->cloud = corePoints.cloud; + feature->sourceName = statSF1->getName(); + feature->scale = scale; + + if (rule.cloud2 && field2 && rule.op != FeatureRule::NO_OPERATION) + { + QString resultSFName2 = rule.cloud2->getName() + "." + field2->getName() + QString("_") + FeatureRule::StatToString(rule.stat) + "@" + QString::number(scale); + int sfIndex2 = corePoints.cloud->getScalarFieldIndexByName(qPrintable(resultSFName2)); + CCLib::ScalarField* statSF2 = ExtractStat(corePoints, rule.cloud2, field2.data(), scale, rule.stat, qPrintable(resultSFName2), progressCb); + if (!statSF2) + { + error = QString("Failed to extract stat. from field '%1' @ scale %2").arg(field2->getName()).arg(scale); + return false; + } + + //now perform the math operation + if (!PerformMathOp(statSF1, statSF2, rule.op)) + { + error = "Failed to perform the MATH operation"; + return false; + } + + if (sfIndex2 < 0) + { + //release some memory + sfIndex2 = corePoints.cloud->getScalarFieldIndexByName(qPrintable(resultSFName2)); + corePoints.cloud->deleteScalarField(sfIndex2); + } + } + + return feature; + } + else //non scaled feature + { + if (rule.cloud1 != corePoints.cloud && rule.cloud1 != corePoints.origin) + { + assert(false); + error = "Scale-less features (SC0) can only be defined on the core points (origin) cloud"; + return false; + } + + if (rule.cloud2) + { + if (rule.op != FeatureRule::NO_OPERATION) + { + assert(false); + ccLog::Warning("MATH operations cannot be performed on scale-less features (SC0)"); + return false; + } + else + { + assert(false); + ccLog::Warning("Feature has a second cloud associated but no MATH operation is defined"); + } + } + + //build the final SF name + QString resultSFName = /*rule.cloud1->getName() + "." + */field1->getName(); + //if (rule.cloud2 && field2 && rule.op != FeatureRule::NO_OPERATION) + //{ + // resultSFName += QString("_") + FeatureRule::OpToString(rule.op) + "_" + field2->getName(); + //} + + //retrieve/create a SF to host the result + CCLib::ScalarField* resultSF = nullptr; + int sfIdx = corePoints.cloud->getScalarFieldIndexByName(qPrintable(resultSFName)); + if (sfIdx >= 0) + { + //reuse the existing field + resultSF = corePoints.cloud->getScalarField(sfIdx); + } + else + { + //copy the SF1 field + resultSF = new ccScalarField(qPrintable(resultSFName)); + if (!resultSF->resizeSafe(corePoints.cloud->size())) + { + error = "Not enough memory"; + resultSF->release(); + return false; + } + + //copy the values + for (unsigned i = 0; i < corePoints.size(); ++i) + { + resultSF->setValue(i, field1->pointValue(corePoints.originIndex(i))); + } + resultSF->computeMinAndMax(); + int newSFIdx = corePoints.cloud->addScalarField(static_cast(resultSF)); + //update display + if (corePoints.cloud->getDisplay()) + { + corePoints.cloud->setCurrentDisplayedScalarField(newSFIdx); + corePoints.cloud->getDisplay()->redraw(); + } + } + + rule.feature->cloud = corePoints.cloud; + rule.feature->sourceName = resultSF->getName(); + rule.feature->scale = scale; + + //if (rule.cloud2 && field2 && rule.op != FeatureRule::NO_OPERATION) + //{ + // //now perform the math operation + // if (!PerformMathOp(*field1, *field2, rule.op, resultSF)) + // { + // error = "Failed to perform the MATH operation"; + // return false; + // } + + // //sf2 is held by the second cloud for now + // //sf2->release(); + // //sf2 = nullptr; + //} + + return rule.feature; + } } -bool Tools::PrepareFeatures(const FeatureRule::Set& rules, Feature::Set& features, QString& error) +bool Tools::PrepareFeatures(const FeatureRule::Set& rules, const CorePoints& corePoints, Feature::Set& features, QString& error, CCLib::GenericProgressCallback* progressCb/*=nullptr*/) { + if (rules.empty() || !corePoints.origin) + { + //invalid input parameters + assert(false); + return false; + } + for (const FeatureRule::Shared& rule : rules) { QString errorMessage("invalid pointer"); - if (!rule || !rule->checkValidity(errorMessage)) + if (!rule || !rule->checkValidity(/*corePoints, */errorMessage)) { error = "Invalid rule/feature: " + error; return false; } - if () + size_t scaleCount = (rule->scales ? rule->scales->values.size(): 1); + for (size_t i = 0; i < scaleCount; ++i) + { + //retrieve the right scale + double scale = std::numeric_limits::quiet_NaN(); + if (rule->scales) + { + scale = rule->scales->values[i]; + } + + Feature::Shared preparedFeature; + + //we will prepare the different versions of the feature (one per scale, etc.) + //depending on the feature type + switch (rule->feature->getType()) + { + case Feature::Type::PointFeature: + { + //Point feature + preparedFeature = PreparePointBasedFeature(*rule, scale, corePoints, error, progressCb); + break; + } + default: + assert(false); + break; + } + + if (!preparedFeature) + { + //something failed (error should be up to date) + return false; + } + + //otherwise add the new feature + features.push_back(preparedFeature); + } } return true; diff --git a/q3DMASCTools.h b/q3DMASCTools.h index f6a2e61..9bafbf6 100644 --- a/q3DMASCTools.h +++ b/q3DMASCTools.h @@ -19,12 +19,13 @@ //Local #include "Features.h" - -//qCC_db -#include +#include "CorePoints.h" //CCLib -#include +#include + +//qCC_db +class ccPointCloud; class QWidget; @@ -35,9 +36,9 @@ namespace masc { public: - static bool LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, FeatureRule::Set& features); + static bool LoadFile(QString filename, FeatureRule::Set& features, std::vector& loadedClouds, CorePoints& corePoints); - static bool PrepareFeatures(const FeatureRule::Set& rules, Feature::Set& features, QString& error); + static bool PrepareFeatures(const FeatureRule::Set& rules, const CorePoints& corePoints, Feature::Set& features, QString& error, CCLib::GenericProgressCallback* progressCb = nullptr); static bool RandomSubset(ccPointCloud* cloud, float ratio, CCLib::ReferenceCloud* inRatioSubset, CCLib::ReferenceCloud* outRatioSubset); };