diff --git a/Features.h b/Features.h index a1ffc08..bcd2823 100644 --- a/Features.h +++ b/Features.h @@ -36,6 +36,8 @@ public: //PointFeatureType { Invalid = 0 , Intensity + , X + , Y , Z , NbRet , RetNb @@ -59,6 +61,10 @@ public: //PointFeatureType return "Invalid"; case Intensity: return "INT"; + case X: + return "X"; + case Y: + return "Y"; case Z: return "Z"; case NbRet: @@ -92,11 +98,15 @@ public: //PointFeatureType return "Invalid"; } - static PointFeatureType FromString(QString token) + static inline PointFeatureType FromString(const QString& token) { return FromUpperString(token.toUpper()); } + static PointFeatureType FromUpperString(const QString& token) { - token = token.toUpper(); if (token == "INT") return Intensity; + else if (token == "X") + return X; + else if (token == "Y") + return Y; else if (token == "Z") return Z; else if (token == "NBRET") @@ -113,15 +123,15 @@ public: //PointFeatureType return B; else if (token == "NIR") return NIR; - else if (token == "DIPANG") + else if (token == "NORMDIPANG") return DipAng; - else if (token == "DIPDIR") + else if (token == "NORMDIPDIR") return DipDir; else if (token == "M3C2") return M3C2; else if (token == "PCV") return PCV; - else if (token == "SF#") + else if (token.startsWith("SF")) return SF; assert(false); @@ -134,10 +144,43 @@ public: //methods virtual Type getType() override { return Type::PointFeature; } //! Default constructor - PointFeature(ccPointCloud* p_cloud, PointFeatureType p_type, Source p_source, QString p_name) - : Feature(p_cloud, p_source, p_name) + PointFeature(PointFeatureType p_type, ccPointCloud* p_cloud = nullptr) + : Feature(p_cloud) , type(p_type) - {} + { + //auomatically set the right source for specific features + switch (type) + { + case Z: + source = DimZ; + sourceName = "X"; + break; + case Z: + source = DimZ; + sourceName = "Y"; + break; + case Z: + source = DimZ; + sourceName = "Z"; + break; + case R: + source = Red; + sourceName = "Red"; + break; + case G: + source = Green; + sourceName = "Green"; + break; + case B: + source = Blue; + sourceName = "Blue"; + break; + default: + source = ScalarField; + //sourceName --> TBD later + break; + } + } //! Returns the descriptor for this particular feature virtual QString toString() const override @@ -253,10 +296,9 @@ public: //NeighborhoodFeatureType return "Invalid"; } - static NeighborhoodFeatureType FromString(QString token) + static inline NeighborhoodFeatureType FromString(const QString& token) { return FromUpperString(token.toUpper()); } + static NeighborhoodFeatureType FromUpperString(const QString& token) { - token = token.toUpper(); - if (token == "PCA1") return PCA1; else if (token == "PCA2") @@ -302,11 +344,11 @@ public: //methods virtual Type getType() override { return Type::NeighborhoodFeature; } //! Default constructor - NeighborhoodFeature(ccPointCloud* p_cloud, NeighborhoodFeatureType p_type, Source p_source, QString p_name, double p_scale = std::numeric_limits::quiet_NaN()) - : Feature(p_cloud, p_source, p_name) - , scale(p_scale) + NeighborhoodFeature(NeighborhoodFeatureType p_type, ccPointCloud* p_cloud = nullptr) + : Feature(p_cloud) , type(p_type) - {} + { + } //! Returns the descriptor for this particular feature virtual QString toString() const override @@ -317,9 +359,6 @@ public: //methods public: //members - //! Neighborhood scale (diameter) - double scale; - //! Neighborhood feature type /** \warning different from the feature type **/ @@ -330,8 +369,8 @@ public: //members struct TwoCloudsFeature : public Feature { //! Default constructor - TwoCloudsFeature(ccPointCloud* p_cloud, ccPointCloud* p_otherCloud, Source p_source, QString p_name, double p_scale = std::numeric_limits::quiet_NaN()) - : Feature(p_cloud, p_source, p_name) + TwoCloudsFeature(ccPointCloud* p_cloud = nullptr, ccPointCloud* p_otherCloud = nullptr, double p_scale = std::numeric_limits::quiet_NaN()) + : Feature(p_cloud, p_scale) , otherCloud(p_otherCloud) { assert(otherCloud); @@ -370,10 +409,9 @@ public: //ContextBasedFeatureType return "Invalid"; } - static ContextBasedFeatureType FromString(QString token) + static inline ContextBasedFeatureType FromString(const QString& token) { return FromUpperString(token.toUpper()); } + static ContextBasedFeatureType FromUpperString(const QString& token) { - token = token.toUpper(); - if (token.startsWith("DZ")) return DZ; else if (token.startsWith("DH")) @@ -389,18 +427,15 @@ public: //methods virtual Type getType() override { return Type::ContextBasedFeature; } //! Default constructor - ContextBasedFeature(ccPointCloud* p_cloud, - ccPointCloud* p_otherCloud, - ContextBasedFeatureType p_type, - int p_kNN, - double p_scale, - int p_ctxClassLabel, - Source p_source, - QString p_name) - : TwoCloudsFeature(p_cloud, p_otherCloud, p_source, p_name) + ContextBasedFeature(ContextBasedFeatureType p_type, + ccPointCloud* p_cloud = nullptr, + ccPointCloud* p_otherCloud = nullptr, + int p_kNN = 0, + double p_scale = std::numeric_limits::quiet_NaN(), + int p_ctxClassLabel = 0) + : TwoCloudsFeature(p_cloud, p_otherCloud, p_scale) , type(p_type) , kNN(p_kNN) - , scale(p_scale) , ctxClassLabel(p_ctxClassLabel) {} @@ -418,8 +453,6 @@ public: //methods //Number of neighbors int kNN; - //! Scale (optional) - double scale; //! Context class (label) int ctxClassLabel; }; @@ -451,10 +484,9 @@ public: //DualCloudFeatureType return "Invalid"; } - static DualCloudFeatureType FromString(QString token) + static inline DualCloudFeatureType FromString(const QString& token) { return FromUpperString(token.toUpper()); } + static DualCloudFeatureType FromUpperString(const QString& token) { - token = token.toUpper(); - if (token == "IDIFF") return IDIFF; @@ -468,16 +500,13 @@ public: //methods virtual Type getType() override { return Type::DualCloudFeature; } //! Default constructor - DualCloudFeature( ccPointCloud* p_cloud, - ccPointCloud* p_otherCloud, - DualCloudFeatureType p_type, - double p_scale, - Source p_source, - QString p_name - ) - : TwoCloudsFeature(p_cloud, p_otherCloud, p_source, p_name) + DualCloudFeature( DualCloudFeatureType p_type, + ccPointCloud* p_cloud = nullptr, + ccPointCloud* p_otherCloud = nullptr, + double p_scale = std::numeric_limits::quiet_NaN() + ) + : TwoCloudsFeature(p_cloud, p_otherCloud, p_scale) , type(p_type) - , scale(p_scale) {} //! Returns the descriptor for this particular feature @@ -491,9 +520,6 @@ public: //methods /** \warning different from the feature type **/ DualCloudFeatureType type; - - //! Scale (optional) - double scale; }; struct Scales @@ -504,23 +530,122 @@ struct Scales struct FeatureRule { - typedef std::vector Set; + typedef QSharedPointer Shared; + typedef std::vector Set; enum Stat { - NO_STAT, MEAN, STD, MIN, MAX + NO_STAT, + MEAN, + MODE, //number with the highest frequency + STD, + RANGE, + SKEW //(SKEW = (MEAN - MODE)/STD) }; Stat stat = NO_STAT; //only considered if a scale is defined - + + static QString StatToString(Stat stat) + { + switch (stat) + { + case MEAN: + return "MEAN"; + case MODE: + return "MODE"; + case STD: + return "STD"; + case RANGE: + return "RANGE"; + case SKEW: + return "SKEW"; + default: + break; + }; + return QString(); + } + enum Operation { NO_OPERATION, MINUS, PLUS, DIVIDE, MULTIPLY }; Operation op = NO_OPERATION; //only considered if 2 clouds are defined + static QString OpToString(Operation op) + { + switch (op) + { + case MINUS: + return "MINUS"; + case PLUS: + return "PLUS"; + case DIVIDE: + return "DIVIDE"; + case MULTIPLY: + return "MULTIPLY"; + default: + break; + }; + return QString(); + } + Scales::Shared scales; - Feature::Type featureType = Feature::Type::Invalid; Feature::Shared feature; ccPointCloud* cloud1 = nullptr; ccPointCloud* cloud2 = nullptr; + + //! Source scalar field index (if the feature source is 'ScalarField') + int sourceSFIndex = -1; + + //! Checks the rule validity + bool checkValidity(QString &error) const + { + int cloudCount = (cloud1 ? (cloud2 ? 2 : 1) : 0); + + if (feature == nullptr) //feature object should have already been instantiated + { + error = "feature rule has no associated feature"; + return false; + } + if (scales != nullptr && scales->values.empty()) + { + error = "invalid scales definition"; + return false; + } + if (stat != FeatureRule::NO_STAT) + { + if (feature->getType() != Feature::Type::PointFeature) + { + 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"; + return false; + } + } + if (stat != FeatureRule::NO_OPERATION) + { + if (feature->getType() == Feature::Type::DualCloudFeature) + { + error = "math operation can't be defined on dual-cloud features"; + return false; + } + if (cloudCount < 2) + { + error = "at least two clouds are required to apply math operations"; + return false; + } + } + if (feature->getType() == Feature::Type::DualCloudFeature || feature->getType() == Feature::Type::ContextBasedFeature) + { + if (cloudCount < 2) + { + error = "at least two clouds are required to compute dual-cloud or context-based features"; + return false; + } + } + + return true; + } }; diff --git a/FeaturesInterface.h b/FeaturesInterface.h index 519f426..570bf7d 100644 --- a/FeaturesInterface.h +++ b/FeaturesInterface.h @@ -58,9 +58,10 @@ struct Feature }; //! Default constructor - Feature(ccPointCloud* p_cloud, Source p_source, QString p_sourceName) + Feature(ccPointCloud* p_cloud = nullptr, double p_scale = std::numeric_limits::quiet_NaN(), Source p_source = ScalarField, QString p_sourceName = QString()) : cloud(p_cloud) , source(p_source) + , scale(p_scale) , sourceName(p_sourceName) { assert(cloud); @@ -68,8 +69,10 @@ struct Feature //! Associated cloud ccPointCloud* cloud; + //! Scale (diameter) + double scale; //! Values source Source source; - //! Feature source name (mandatory for scalar fields) + //! Feature source name (mandatory for scalar fields if the SF index is not set) QString sourceName; }; diff --git a/q3DMASCClassifier.cpp b/q3DMASCClassifier.cpp index 8d7a969..f0d6aa9 100644 --- a/q3DMASCClassifier.cpp +++ b/q3DMASCClassifier.cpp @@ -55,6 +55,7 @@ static QSharedPointer GetSource(const Feature::Shared& f, c { case Feature::ScalarField: { + assert(!f->sourceName.isEmpty()); int sfIdx = cloud->getScalarFieldIndexByName(qPrintable(f->sourceName)); if (sfIdx >= 0) { @@ -63,6 +64,7 @@ static QSharedPointer GetSource(const Feature::Shared& f, c else { ccLog::Warning(QObject::tr("Internal error: unknwon scalar field '%1'").arg(f->sourceName)); + return QSharedPointer(nullptr); } } break; diff --git a/q3DMASCTools.cpp b/q3DMASCTools.cpp index d96b8c8..9dcd637 100644 --- a/q3DMASCTools.cpp +++ b/q3DMASCTools.cpp @@ -17,8 +17,17 @@ #include "q3DMASCTools.h" +//Local +#include "ScalarFieldWrappers.h" + //qCC_io #include +#include +//qCC_db +#include + +//CCLib +#include //Qt #include @@ -46,10 +55,11 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea Scales::Shared scales(new Scales); assert(features.empty()); - QMap > clouds; + QMap > clouds; QTextStream stream(&file); - while (true) + + for (int lineNumber = 0; ; ++lineNumber) { QString line = stream.readLine(); if (line.isNull()) @@ -57,6 +67,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea //eof break; } + ++lineNumber; if (line.startsWith("#")) { @@ -64,7 +75,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea continue; } - //strip out the potential comment at the end of the mine as well + //strip out the potential comment at the end of the line as well int commentIndex = line.indexOf('#'); if (commentIndex >= 0) line = line.left(commentIndex); @@ -76,7 +87,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea QStringList tokens = command.split(':'); if (tokens.size() != 2) { - ccLog::Warning("Malformed file: expecting 2 tokens after 'cloud:'"); + ccLog::Warning("Malformed file: expecting 2 tokens after 'cloud:' on line #" + QString::number(lineNumber)); return false; } QString pcName = tokens[0]; @@ -100,7 +111,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea delete object; return false; } - clouds.insert(pcName, QScopedPointer(static_cast(object))); + clouds.insert(pcName, QSharedPointer(static_cast(object))); } } else if (upperLine.startsWith("SCALES:")) //scales @@ -109,7 +120,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea QStringList tokens = command.split(';'); if (tokens.empty()) { - ccLog::Warning("Malformed file: expecting at least one token after 'scales:'"); + ccLog::Warning("Malformed file: expecting at least one token after 'scales:' on line #" + QString::number(lineNumber)); return false; } @@ -132,12 +143,12 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea double stop = subTokens[2].toDouble(ok + 2); if (!ok[0] || !ok[1] || !ok[2]) { - ccLog::Warning(QString("Malformed file: invalid values in scales range (%1)").arg(token)); + ccLog::Warning(QString("Malformed file: invalid values in scales range (%1) on line #%2").arg(token).arg(lineNumber)); return false; } if (stop < start || step <= 1.0-6) { - ccLog::Warning(QString("Malformed file: invalid range (%1)").arg(token)); + 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) @@ -151,7 +162,7 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea double v = token.toDouble(&ok); if (!ok) { - ccLog::Warning(QString("Malformed file: invalid scale value (%1)").arg(token)); + ccLog::Warning(QString("Malformed file: invalid scale value (%1) on line #%2").arg(token).arg(lineNumber)); return false; } scales->values.push_back(v); @@ -170,17 +181,238 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea QStringList tokens = command.split('_'); if (tokens.empty()) { - ccLog::Warning("Malformed file: expecting at least one token after 'feature:'"); + ccLog::Warning("Malformed file: expecting at least one token after 'feature:' on line #" + QString::number(lineNumber)); return false; } + FeatureRule::Shared rule(new FeatureRule); + + //read the type QString typeStr = tokens[0].toUpper(); - //TODO - + { + for (int iteration = 0; iteration < 1; ++iteration) //fake loop for easy break + { + PointFeature::PointFeatureType pointFeatureType = PointFeature::FromUpperString(typeStr); + if (pointFeatureType != PointFeature::Invalid) + { + //we have a point feature + PointFeature::Shared pointFeature(new PointFeature(pointFeatureType)); + + //specific case: 'SF#' + if (pointFeatureType == PointFeature::SF) + { + QString sfIndexStr = typeStr.mid(2); + bool ok = true; + int sfIndex = sfIndexStr.toInt(&ok); + if (!ok) + { + ccLog::Warning(QString("Malformed file: expecting a valid integer value after 'SF' on line #%1").arg(lineNumber)); + return false; + } + pointFeature->sourceSFIndex = sfIndex; + } + + rule->feature = pointFeature; + break; + } + NeighborhoodFeature::NeighborhoodFeatureType neighborhoodFeatureType = NeighborhoodFeature::FromUpperString(typeStr); + if (neighborhoodFeatureType != NeighborhoodFeature::Invalid) + { + //we have a neighborhood feature + rule->feature = NeighborhoodFeature::Shared(new NeighborhoodFeature(neighborhoodFeatureType)); + break; + } + ContextBasedFeature::ContextBasedFeatureType contextBasedFeatureType = ContextBasedFeature::FromUpperString(typeStr); + if (contextBasedFeatureType != ContextBasedFeature::Invalid) + { + //we have a context-based feature + rule->feature = ContextBasedFeature::Shared(new ContextBasedFeature(contextBasedFeatureType)); + break; + } + DualCloudFeature::DualCloudFeatureType dualCloudFeatureType = DualCloudFeature::FromUpperString(typeStr); + if (dualCloudFeatureType != DualCloudFeature::Invalid) + { + //we have a dual cloud feature + rule->feature = DualCloudFeature::Shared(new DualCloudFeature(dualCloudFeatureType)); + break; + } + + ccLog::Warning(QString("Malformed file: unrecognized token '%1' after 'feature:' on line #%2").arg(typeStr).arg(lineNumber)); + return false; + } + } + assert(rule->feature); + + //read the scales + { + QString scaleStr = tokens[1].toUpper(); + if (!scaleStr.startsWith("SC")) + { + ccLog::Warning(QString("Malformed file: unrecognized token '%1' (expecting the scale descriptor 'SC...' on line #%2").arg(typeStr).arg(lineNumber)); + return false; + } + + if (scaleStr == "SC0") + { + //no scale + } + else if (scaleStr == "SCX") + { + //all scales + rule->scales = scales; + } + else + { + //read the specific scale index + QString scaleStr = scaleStr.mid(2); + bool ok = true; + double scale = scaleStr.toDouble(&ok); + if (!ok) + { + ccLog::Warning(QString("Malformed file: expecting a valid number after 'SC:' on line #%1").arg(lineNumber)); + return false; + } + rule->scales = Scales::Shared(new Scales); + rule->scales->values.resize(1); + rule->scales->values.front() = scale; + } + } + + //process the next tokens (may not be ordered) + int cloudCount = 0; + bool statDefined = false; + bool mathDefined = false; + for (int i = 2; i < tokens.size(); ++i) + { + QString token = tokens[i].toUpper(); + + //is the token a 'stat' one? + if (!statDefined) + { + if (token == "MEAN") + { + rule->stat = FeatureRule::MEAN; + statDefined = true; + } + else if (token == "MODE") + { + rule->stat = FeatureRule::MODE; + statDefined = true; + } + else if (token == "STD") + { + rule->stat = FeatureRule::STD; + statDefined = true; + } + else if (token == "RANGE") + { + rule->stat = FeatureRule::RANGE; + statDefined = true; + } + else if (token == "SKEW") + { + rule->stat = FeatureRule::SKEW; + statDefined = true; + } + + if (statDefined) + { + continue; + } + } + + //is the token a cloud name? + if (cloudCount < 2) + { + bool cloudNameMatches = false; + for (QMap >::const_iterator it = clouds.begin(); it != clouds.end(); ++it) + { + if (it.key().toUpper() == token) + { + if (cloudCount == 0) + rule->cloud1 = it.value().data(); + else + rule->cloud2 = it.value().data(); + ++cloudCount; + cloudNameMatches = true; + break; + } + } + + if (cloudNameMatches) + { + continue; + } + } + + //is the token a 'math' one? + if (cloudCount == 2 && rule->feature->getType() != Feature::Type::DualCloudFeature && !mathDefined) + { + if (token == "MINUS") + { + rule->op = FeatureRule::MINUS; + mathDefined = true; + } + else if (token == "PLUS") + { + rule->op = FeatureRule::PLUS; + mathDefined = true; + } + else if (token == "DIVIDE") + { + rule->op = FeatureRule::DIVIDE; + mathDefined = true; + } + else if (token == "MULTIPLY") + { + rule->op = FeatureRule::MULTIPLY; + mathDefined = true; + } + + if (mathDefined) + { + continue; + } + } + + //is the token a 'context' descriptor? + if (rule->feature->getType() == Feature::Type::ContextBasedFeature && token.startsWith("CTX")) + { + //read the context label + QString ctxLabelStr = token.mid(2); + bool ok = true; + int ctxLabel = ctxLabelStr.toInt(&ok); + if (!ok) + { + ccLog::Warning(QString("Malformed file: expecting a valid integer value after 'CTX' on line #%1").arg(lineNumber)); + return false; + } + static_cast(rule->feature.data())->ctxClassLabel = ctxLabel; + continue; + } + + //if we are here, it means we couldn't find a correspondance for the current token + ccLog::Warning(QString("Malformed file: unrecognized or unexpected token '%1' on line #%2").arg(token).arg(lineNumber)); + return false; + } + + //now check the consistency of the rule + assert(rule && rule->feature); + + QString errorMessage; + bool ruleIsValid = rule->checkValidity(errorMessage); + if (!ruleIsValid) + { + ccLog::Warning("Malformed feature: " + errorMessage + QString("(line %1)").arg(lineNumber)); + return false; + } + + //otherwise save it + features.push_back(rule); } else { - ccLog::Warning("Unrecognized token/line: " + (line.length() < 10 ? line : line.left(10) + "...")); + ccLog::Warning(QString("Line #%1: unrecognized token/command: ").arg(lineNumber) + (line.length() < 10 ? line : line.left(10) + "...")); return false; } } @@ -188,6 +420,555 @@ bool Tools::LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, Fea return true; } +static CCLib::ScalarField* RetrieveSF(const ccPointCloud* cloud, const QString& sfName, bool caseSensitive = true) +{ + if (!cloud) + { + assert(false); + return nullptr; + } + int sfIdx = -1; + if (caseSensitive) + { + sfIdx = cloud->getScalarFieldIndexByName(qPrintable(sfName)); + } + else + { + QString sfNameUpper = sfName.toUpper(); + for (unsigned i = 0; i < cloud->getNumberOfScalarFields(); ++i) + { + if (QString(cloud->getScalarField(i)->getName()).toUpper() == sfNameUpper) + { + sfIdx = static_cast(i); + break; + } + } + } + + if (sfIdx >= 0) + { + return cloud->getScalarField(sfIdx); + } + else + { + return nullptr; + } +} + +static const char* s_echoRatioSFName = "EchoRat"; +static const char* s_NIRSFName = "NIR"; +static const char* s_M3C2SFName = "M3C2 distance"; +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) +{ + QString sfName; + switch (featureType) + { + case PointFeature::Intensity: + { + CCLib::ScalarField* sf = RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_INTENSITY], false); + if (!sf) + { + error = "Cloud has no 'intensity' scalar field"; + return nullptr; + } + return sf; + } + case PointFeature::X: + case PointFeature::Y: + case PointFeature::Z: + //not a ScalarField source + error = "Internal error (source is not a scalar field)"; + return nullptr; + case PointFeature::NbRet: + { + CCLib::ScalarField* sf = RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_NUMBER_OF_RETURNS], false); + if (!sf) + { + error = "Cloud has no 'number of returns' scalar field"; + return nullptr; + } + return sf; + } + case PointFeature::RetNb: + { + CCLib::ScalarField* sf = RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_RETURN_NUMBER], false); + if (!sf) + { + error = "Cloud has no 'return number' scalar field"; + return nullptr; + } + return 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 + CCLib::ScalarField* numberOfRetSF = RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_NUMBER_OF_RETURNS], false); + if (!numberOfRetSF) + { + error = "Can't compute the 'echo ratio' field: no 'Number of Return' SF available"; + return nullptr; + } + CCLib::ScalarField* retNumberSF = RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_RETURN_NUMBER], false); + if (!retNumberSF) + { + error = "Can't compute the 'echo ratio' field: no 'Return number' SF available"; + return nullptr; + } + if (retNumberSF->size() != numberOfRetSF->size() || retNumberSF->size() != cloud->size()) + { + 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; + } + case PointFeature::R: + case PointFeature::G: + case PointFeature::B: + //not a ScalarField source + error = "Internal error (source is not a scalar field)"; + return nullptr; + case PointFeature::NIR: + { + CCLib::ScalarField* sf = RetrieveSF(cloud, s_NIRSFName, false); + if (!sf) + { + error = "Cloud has no 'NIR' scalar field"; + return nullptr; + } + return 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! + 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; + } + case PointFeature::M3C2: + { + CCLib::ScalarField* sf = RetrieveSF(cloud, s_M3C2SFName, true); + if (!sf) + { + error = "Cloud has no 'm3c2 distance' scalar field"; + return nullptr; + } + return sf; + } + case PointFeature::PCV: + { + CCLib::ScalarField* sf = RetrieveSF(cloud, s_PCVSFName, true); + if (!sf) + { + error = "Cloud has no 'PCV/Illuminance' scalar field"; + return nullptr; + } + return sf; + } + case PointFeature::SF: + if (sourceSFIndex < 0 || sourceSFIndex >= static_cast(cloud->getNumberOfScalarFields())) + { + error = QString("Can't retrieve the specified SF: invalid index (%1)").arg(sourceSFIndex); + return nullptr; + } + return cloud->getScalarField(sourceSFIndex); + default: + break; + } + + error = "Unhandled feature type"; + return nullptr; +} + +static bool ExtractStatFromSF( const CCLib::DgmOctree::octreeCell& cell, + void** additionalParameters, + CCLib::NormalizedProgress* nProgress = nullptr) +{ + //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(); + + //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; + } + + for (unsigned i = 0; i < n; ++i) + { + //retrieve the points around the current cell point + cell.points->getPoint(i, nNSS.queryPoint); + + //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; + + for (unsigned k = 0; k < kNN; ++k) + { + unsigned index = nNSS.pointsInNeighbourhood[k].pointIndex; + ScalarType v = inputSF->getValue(index); + + //track min and max values + if (k != 0) + { + if (v < minValue) + minValue = v; + else if (v > maxValue) + maxValue = v; + } + else + { + minValue = maxValue = v; + } + + //compute average and std. dev. + sum += v; + sum2 += static_cast(v) * v; + + if (withMode) + { + if (modeCounter.contains(v)) + { + ++modeCounter[v]; + } + else + { + modeCounter[v] = 1; + } + } + } + + 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) + { + if (it.value() > maxCounter) + { + maxCounter = it.value(); + mode = it.key(); + } + } + } + + ScalarType outValue = NAN_VALUE; + switch (stat) + { + 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: + { + 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; + } + default: + assert(false); + break; + } + resultSF->setValue(cell.points->getPointGlobalIndex(i), outValue); + + if (nProgress && !nProgress->oneStep()) + { + return false; + } + } + + return true; +} + +static CCLib::ScalarField* ExtractStat( ccPointCloud* cloud, + CCLib::ScalarField* sf, + double scale, + FeatureRule::Stat stat, + CCLib::GenericProgressCallback* progressCb = nullptr) +{ + if (!cloud || !sf || scale <= 0.0 || stat == FeatureRule::NO_STAT) + { + //invalid input parameters + assert(false); + return nullptr; + } + + ccOctree::Shared octree = cloud->getOctree(); + if (!octree) + { + octree = cloud->computeOctree(progressCb); + if (!octree) + { + ccLog::Warning("Failed to compute octree"); + return nullptr; + } + } + + CCLib::ScalarField* resultSF = nullptr; + QString resultSFName = sf->getName() + QString("_") + FeatureRule::StatToString(stat) + "_" + QString::number(scale); + int sfIdx = cloud->getScalarFieldIndexByName(qPrintable(resultSFName)); + if (sfIdx >= 0) + { + resultSF = cloud->getScalarField(sfIdx); + } + else + { + resultSF = new ccScalarField(qPrintable(resultSFName)); + if (!resultSF->reserveSafe(cloud->size())) + { + ccLog::Warning("Not enough memory"); + resultSF->release(); + return nullptr; + } + } + resultSF->fill(NAN_VALUE); + + 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) + }; + + if (octree->executeFunctionForAllCellsAtLevel( octreeLevel, + ExtractStatFromSF, + additionalParameters, + true, + progressCb, + qPrintable(QString("Extract stat @ scale %1").arg(scale))) == 0) + { + //something went wrong + ccLog::Warning("Process failed"); + resultSF->release(); + return nullptr; + } + + resultSF->computeMinAndMax(); + cloud->addScalarField(static_cast(resultSF)); + + return resultSF; +} + + +static bool PreparePointBasedFeature(FeatureRule& rule, QString& error) +{ + assert(rule.feature && rule.feature->getType() == Feature::Type::PointFeature); + + PointFeature* feature = static_cast(rule.feature.data()); + + std::vector preparedFeatures; + + //look for the source field (and compute it if necessary) + CCLib::ScalarField* sf1 = RetrieveOrComputeSF(feature->type, rule.sourceSFIndex, rule.cloud1, error); + if (!sf1) + { + //error should be up to date + return false; + } + + CCLib::ScalarField* sf2 = nullptr; + if (rule.cloud2 && rule.op != FeatureRule::NO_OPERATION) + { + sf2 = RetrieveOrComputeSF(feature->type, rule.sourceSFIndex, rule.cloud2, error); + if (!sf2) + { + //error should be up to date + 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) + { + CCLib::ScalarField* statSF1 = ExtractStat(rule.cloud1, sf1, s, rule.stat); + if (!statSF1) + { + 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) + { + ccLog::Warning(QString("Failed to extract stat. from sf '%1' @ scale %2").arg(sf2->getName()).arg(s)); + return false; + } + PointFeature::Shared f2(new PointFeature(*feature)); + f2->cloud = rule.cloud2; + f2->sourceName = statSF2->getName(); + preparedFeatures.push_back(f2); + } + } + } + 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: + } + +} + +bool Tools::PrepareFeatures(const FeatureRule::Set& rules, Feature::Set& features, QString& error) +{ + for (const FeatureRule::Shared& rule : rules) + { + QString errorMessage("invalid pointer"); + if (!rule || !rule->checkValidity(errorMessage)) + { + error = "Invalid rule/feature: " + error; + return false; + } + + if () + } + + return true; +} + bool Tools::RandomSubset(ccPointCloud* cloud, float ratio, CCLib::ReferenceCloud* inRatioSubset, CCLib::ReferenceCloud* outRatioSubset) { if (!cloud) @@ -272,4 +1053,3 @@ bool Tools::RandomSubset(ccPointCloud* cloud, float ratio, CCLib::ReferenceCloud return true; } - diff --git a/q3DMASCTools.h b/q3DMASCTools.h index 48123ad..f6a2e61 100644 --- a/q3DMASCTools.h +++ b/q3DMASCTools.h @@ -37,6 +37,8 @@ namespace masc static bool LoadFile(QString filename, ccPointCloud* pc1, ccPointCloud* pc2, FeatureRule::Set& features); + static bool PrepareFeatures(const FeatureRule::Set& rules, Feature::Set& features, QString& error); + static bool RandomSubset(ccPointCloud* cloud, float ratio, CCLib::ReferenceCloud* inRatioSubset, CCLib::ReferenceCloud* outRatioSubset); };