diff --git a/PointFeature.cpp b/PointFeature.cpp index a4957b5..3a182b0 100644 --- a/PointFeature.cpp +++ b/PointFeature.cpp @@ -301,7 +301,6 @@ static bool ExtractStatFromSF( const CCVector3& queryPoint, assert(false); return false; } - std::numeric_limits::quiet_NaN(); //spherical neighborhood extraction structure CCLib::DgmOctree::NearestNeighboursSphericalSearchStruct nNSS; @@ -437,6 +436,66 @@ static bool ExtractStatFromSF( const CCVector3& queryPoint, return true; } +static bool PrepareOctree(ccPointCloud* sourceCloud, CCLib::GenericProgressCallback* progressCb = nullptr) +{ + if (!sourceCloud) + { + //invalid input parameters + assert(false); + return false; + } + + ccOctree::Shared octree = sourceCloud->getOctree(); + if (!octree) + { + ccLog::Print(QString("Computing octree of cloud %1 (%2 points)").arg(sourceCloud->getName()).arg(sourceCloud->size())); + octree = sourceCloud->computeOctree(progressCb); + if (!octree) + { + ccLog::Warning("Failed to compute octree"); + return nullptr; + } + } + + return true; +} + +static CCLib::ScalarField* PrepareSF(const CorePoints& corePoints, const char* resultSFName) +{ + if (!corePoints.cloud || !resultSFName) + { + //invalid input parameters + assert(false); + return nullptr; + } + + CCLib::ScalarField* resultSF = nullptr; + int sfIdx = corePoints.cloud->getScalarFieldIndexByName(resultSFName); + if (sfIdx >= 0) + { + resultSF = corePoints.cloud->getScalarField(sfIdx); + } + else + { + ccScalarField* newSF = new ccScalarField(resultSFName); + if (!newSF->resizeSafe(corePoints.cloud->size())) + { + ccLog::Warning("Not enough memory"); + newSF->release(); + return nullptr; + } + corePoints.cloud->addScalarField(newSF); + + resultSF = newSF; + + } + + assert(resultSF); + resultSF->fill(NAN_VALUE); + + return resultSF; +} + static CCLib::ScalarField* ExtractStat( const CorePoints& corePoints, ccPointCloud* sourceCloud, const IScalarFieldWrapper* sourceField, @@ -580,11 +639,13 @@ bool PointFeature::prepare( const CorePoints& corePoints, { //invalid input assert(false); + error = "internal error (no input core points)"; return false; } //look for the source field - QSharedPointer field1 = retrieveField(cloud1, error); + assert(!field1); + field1 = retrieveField(cloud1, error); if (!field1) { //error should be up to date @@ -601,12 +662,12 @@ bool PointFeature::prepare( const CorePoints& corePoints, return false; } - QSharedPointer field2; if (cloud2) { //no need to compute the second scalar field if no MATH operation has to be performed?! if (op != Feature::NO_OPERATION) { + assert(!field2); field2 = retrieveField(cloud2, error); if (!field2) { @@ -630,10 +691,21 @@ bool PointFeature::prepare( const CorePoints& corePoints, } resultSFName += "@" + QString::number(scale); - CCLib::ScalarField* statSF1 = ExtractStat(corePoints, cloud1, field1.data(), scale, stat, qPrintable(resultSFName), progressCb); + //prepare the octree + //if (!PrepareOctree(cloud1, progressCb)) + //{ + // error = "Failed to compute octree (not enough memory?)"; + // return false; + //} + + //and the scalar fielda + assert(!statSF1); + statSF1 = PrepareSF(corePoints, qPrintable(resultSFName)); + //CCLib::ScalarField* statSF1 = ExtractStat(corePoints, cloud1, field1.data(), scale, stat, qPrintable(resultSFName), progressCb); if (!statSF1) { - error = QString("Failed to extract stat. from field '%1' @ scale %2").arg(field1->getName()).arg(scale); + //error = QString("Failed to extract stat. from field '%1' @ scale %2").arg(field1->getName()).arg(scale); + error = QString("Failed to prepare scalar field for field '%1' @ scale %2").arg(field1->getName()).arg(scale); return false; } sourceName = statSF1->getName(); @@ -641,8 +713,18 @@ bool PointFeature::prepare( const CorePoints& corePoints, if (cloud2 && field2 && op != Feature::NO_OPERATION) { QString resultSFName2 = cloud2Label + "." + field2->getName() + QString("_") + Feature::StatToString(stat) + "@" + QString::number(scale); - int sfIndex2 = corePoints.cloud->getScalarFieldIndexByName(qPrintable(resultSFName2)); - CCLib::ScalarField* statSF2 = ExtractStat(corePoints, cloud2, field2.data(), scale, stat, qPrintable(resultSFName2), progressCb); + keepStatSF2 = (corePoints.cloud->getScalarFieldIndexByName(qPrintable(resultSFName2)) >= 0); //we remember that the scalar field was already existing! + + //prepare the octree + //if (!PrepareOctree(cloud2, progressCb)) + //{ + // error = "Failed to compute octree (not enough memory?)"; + // return false; + //} + + assert(!statSF2); + statSF2 = PrepareSF(corePoints, qPrintable(resultSFName2)); + //statSF2 = ExtractStat(corePoints, cloud2, field2.data(), scale, stat, qPrintable(resultSFName2), progressCb); if (!statSF2) { error = QString("Failed to extract stat. from field '%1' @ scale %2").arg(field2->getName()).arg(scale); @@ -650,18 +732,18 @@ bool PointFeature::prepare( const CorePoints& corePoints, } //now perform the math operation - if (!PerformMathOp(statSF1, statSF2, op)) - { - error = "Failed to perform the MATH operation"; - return false; - } + //if (!PerformMathOp(statSF1, statSF2, 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); - } + //if (sfIndex2 < 0) + //{ + // //release some memory + // sfIndex2 = corePoints.cloud->getScalarFieldIndexByName(qPrintable(resultSFName2)); + // corePoints.cloud->deleteScalarField(sfIndex2); + //} } return true; @@ -692,10 +774,6 @@ bool PointFeature::prepare( const CorePoints& corePoints, //build the final SF name QString resultSFName = /*cloud1Label + "." + */field1->getName(); - //if (cloud2 && field2 && op != Feature::NO_OPERATION) - //{ - // resultSFName += QString("_") + Feature::OpToString(op) + "_" + field2->getName(); - //} //retrieve/create a SF to host the result CCLib::ScalarField* resultSF = nullptr; @@ -733,24 +811,203 @@ bool PointFeature::prepare( const CorePoints& corePoints, sourceName = resultSF->getName(); - //if (cloud2 && field2 && op != Feature::NO_OPERATION) - //{ - // //now perform the math operation - // if (!PerformMathOp(*field1, *field2, 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 true; } } +bool PointFeature::computeStat(const CCLib::DgmOctree::NeighboursSet& pointsInNeighbourhood, const QSharedPointer& sourceField, double& outputValue) const +{ + outputValue = std::numeric_limits::quiet_NaN(); + + if (!sourceField || stat == Feature::NO_STAT) + { + //invalid input parameters + assert(false); + return false; + } + + size_t kNN = pointsInNeighbourhood.size(); + if (kNN == 0) + { + assert(false); + return false; + } + + //specific case + if (stat == Feature::RANGE) + { + double minValue = 0; + double maxValue = 0; + + for (size_t k = 0; k < kNN; ++k) + { + unsigned index = pointsInNeighbourhood[k].pointIndex; + double v = sourceField->pointValue(index); + + //track min and max values + if (k != 0) + { + if (v < minValue) + minValue = v; + else if (v > maxValue) + maxValue = v; + } + else + { + minValue = maxValue = v; + } + } + + outputValue = maxValue - minValue; + return true; + } + else + { + bool withSums = (stat == Feature::MEAN || stat == Feature::STD); + bool storeValues = (stat == Feature::MODE || stat == Feature::SKEW); + double sum = 0.0; + double sum2 = 0.0; + + CCLib::WeibullDistribution::ScalarContainer values; + if (storeValues) + { + try + { + values.resize(kNN); + } + catch (const std::bad_alloc&) + { + ccLog::Warning("Not enough memory"); + return false; + } + } + + for (unsigned k = 0; k < kNN; ++k) + { + unsigned index = pointsInNeighbourhood[k].pointIndex; + double v = sourceField->pointValue(index); + + if (withSums) + { + //compute average and std. dev. + sum += v; + sum2 += v * v; + } + + if (storeValues) + { + values[k] = static_cast(v); + } + } + + switch (stat) + { + case Feature::MEAN: + { + outputValue = sum / kNN; + } + break; + + case Feature::MODE: + { + CCLib::WeibullDistribution w; + w.computeParameters(values); + outputValue = w.computeMode(); + } + break; + + case Feature::STD: + { + outputValue = sqrt(std::abs(sum2 * kNN - sum * sum)) / kNN; + } + break; + + case Feature::RANGE: + { + //we can't be here + assert(false); + } + return false; + + case Feature::SKEW: + { + CCLib::WeibullDistribution w; + w.computeParameters(values); + outputValue = w.computeSkewness(); + } + break; + + default: + { + ccLog::Warning("Unhandled STAT measure"); + assert(false); + } + return false; + + } + } + + return true; +} + +bool PointFeature::finish(const CorePoints& corePoints, QString& error) +{ + if (!scaled()) + { + //nothing to do + return true; + } + + if (!corePoints.cloud) + { + //invalid input + assert(false); + error = "internal error (no input core points)"; + return false; + } + + bool success = true; + + if (statSF1) + { + statSF1->computeMinAndMax(); + } + + if (statSF2) + { + //now perform the math operation + if (op != Feature::NO_OPERATION) + { + if (!PerformMathOp(statSF1, statSF2, op)) + { + error = "Failed to perform the MATH operation"; + success = false; + } + } + + if (keepStatSF2) + { + statSF2->computeMinAndMax(); + } + else + { + int sfIndex2 = corePoints.cloud->getScalarFieldIndexByName(statSF2->getName()); + if (sfIndex2 >= 0) + { + corePoints.cloud->deleteScalarField(sfIndex2); + } + else + { + assert(false); + statSF2->release(); + } + statSF2 = nullptr; + } + } + + return success; +} + QString PointFeature::toString() const { //default keyword otherwise diff --git a/PointFeature.h b/PointFeature.h index 898d13b..9734b9b 100644 --- a/PointFeature.h +++ b/PointFeature.h @@ -31,6 +31,8 @@ namespace masc { public: //PointFeatureType + typedef QSharedPointer Shared; + enum PointFeatureType { Invalid = 0 @@ -142,6 +144,11 @@ namespace masc PointFeature(PointFeatureType p_type) : type(p_type) , sourceSFIndex(-1) + , field1(nullptr) + , field2(nullptr) + , statSF1(nullptr) + , statSF2(nullptr) + , keepStatSF2(false) { //auomatically set the right source for specific features switch (type) @@ -188,6 +195,12 @@ namespace masc //! Returns the descriptor for this particular feature virtual QString toString() const override; + //! Finishes the feature preparation (update the scalar field, etc.) + bool finish(const CorePoints& corePoints, QString& error); + + //! Compute the associated 'stat' on a set of points (and with a given field) + bool computeStat(const CCLib::DgmOctree::NeighboursSet& pointsInNeighbourhood, const QSharedPointer& sourceField, double& outputValue) const; + protected: //methods //! Returns the 'source' field from a given cloud @@ -202,5 +215,10 @@ namespace masc //! Source scalar field index (if the feature source is 'ScalarField') int sourceSFIndex; + + //! For scaled features + QSharedPointer field1, field2; + CCLib::ScalarField *statSF1, *statSF2; + bool keepStatSF2; }; } diff --git a/q3DMASCTools.cpp b/q3DMASCTools.cpp index c85c9a5..dff3975 100644 --- a/q3DMASCTools.cpp +++ b/q3DMASCTools.cpp @@ -792,6 +792,12 @@ CCLib::ScalarField* Tools::RetrieveSF(const ccPointCloud* cloud, const QString& } } +struct FeaturesAndScales +{ + std::vector scales; + std::vector features; +}; + bool Tools::PrepareFeatures(const CorePoints& corePoints, Feature::Set& features, QString& error, CCLib::GenericProgressCallback* progressCb/*=nullptr*/) { if (features.empty() || !corePoints.origin) @@ -800,6 +806,9 @@ bool Tools::PrepareFeatures(const CorePoints& corePoints, Feature::Set& features assert(false); return false; } + + //gather all the scales that need to be extracted + QMap cloudsWithScaledFeatures; for (const Feature::Shared& feature : features) { @@ -816,9 +825,179 @@ bool Tools::PrepareFeatures(const CorePoints& corePoints, Feature::Set& features //something failed (error should be up to date) return false; } + + if (feature->getType() == Feature::Type::PointFeature && feature->scaled()) + { + try + { + //build the scaled feature list attached to the first cloud + if (feature->cloud1) + { + FeaturesAndScales& fas = cloudsWithScaledFeatures[feature->cloud1]; + fas.features.push_back(qSharedPointerCast(feature)); + if (std::find(fas.scales.begin(), fas.scales.end(), feature->scale) == fas.scales.end()) + { + fas.scales.push_back(feature->scale); + } + } + + //build the scaled feature list attached to the second cloud (if any) + if (feature->cloud2 && feature->cloud2 != feature->cloud1 && feature->op != Feature::NO_OPERATION) + { + FeaturesAndScales& fas = cloudsWithScaledFeatures[feature->cloud2]; + fas.features.push_back(qSharedPointerCast(feature)); + if (std::find(fas.scales.begin(), fas.scales.end(), feature->scale) == fas.scales.end()) + { + fas.scales.push_back(feature->scale); + } + } + } + catch (const std::bad_alloc&) + { + error = "Not enough memory"; + return false; + } + } } - return true; + bool success = true; + + //if we have scaled features + if (!cloudsWithScaledFeatures.empty()) + { + for (QMap::iterator it = cloudsWithScaledFeatures.begin(); it != cloudsWithScaledFeatures.end(); ++it) + { + FeaturesAndScales& fas = it.value(); + ccPointCloud* sourceCloud = it.key(); + + //sort the scales + std::sort(fas.scales.begin(), fas.scales.end()); + + //get the octree + ccOctree::Shared octree = sourceCloud->getOctree(); + if (!octree) + { + ccLog::Print(QString("Computing octree of cloud %1 (%2 points)").arg(sourceCloud->getName()).arg(sourceCloud->size())); + octree = sourceCloud->computeOctree(progressCb); + if (!octree) + { + error = "Failed to compute octree (not enough memory?)"; + return false; + } + } + + //now extract the neighborhoods from the biggest to the smallest scale + double largetScale = fas.scales.back(); + PointCoordinateType largestRadius = static_cast(largetScale / 2); //scale is the diameter! + unsigned char octreeLevel = octree->findBestLevelForAGivenNeighbourhoodSizeExtraction(largestRadius); + + unsigned pointCount = corePoints.size(); + if (progressCb) + { + progressCb->setInfo(qPrintable(QString("Computing fields for cloud %1\n(core points: %2)").arg(sourceCloud->getName()).arg(pointCount))); + } + ccLog::Print(QString("Computing fields for cloud %1 (core points: %2)").arg(sourceCloud->getName()).arg(pointCount)); + CCLib::NormalizedProgress nProgress(progressCb, pointCount); + + for (unsigned i = 0; i < pointCount; ++i) + { + //spherical neighborhood extraction structure + CCLib::DgmOctree::NearestNeighboursSphericalSearchStruct nNSS; + { + nNSS.level = octreeLevel; + nNSS.queryPoint = *corePoints.cloud->getPoint(i); + nNSS.prepare(largestRadius, octree->getCellSize(nNSS.level)); + octree->getTheCellPosWhichIncludesThePoint(&nNSS.queryPoint, nNSS.cellPos, nNSS.level); + octree->computeCellCenter(nNSS.cellPos, nNSS.level, nNSS.cellCenter); + } + + //we extract the point's neighbors + unsigned kNN = octree->findNeighborsInASphereStartingFromCell(nNSS, largestRadius, true); + if (kNN == 0) + { + //nothing todo + continue; + } + nNSS.pointsInNeighbourhood.resize(kNN); + + //for each scale (from the largest to the smallest) + for (size_t scaleIndex = 0; scaleIndex < fas.scales.size(); ++scaleIndex) + { + if (scaleIndex != 0) + { + double radius = fas.scales[fas.scales.size() - 1 - scaleIndex] / 2; //scale is the diameter! + double sqRadius = radius * radius; + //remove the farthest points + for (; kNN > 0; --kNN) + { + if (nNSS.pointsInNeighbourhood[kNN - 1].squareDistd <= sqRadius) + { + break; + } + } + + if (kNN == 0) + { + //no need to go further + break; + } + nNSS.pointsInNeighbourhood.resize(kNN); + } + + double outputValue = 0; + for (PointFeature::Shared& feature : fas.features) + { + if (feature->cloud1 == sourceCloud && feature->statSF1 && feature->field1) + { + if (!feature->computeStat(nNSS.pointsInNeighbourhood, feature->field1, outputValue)) + { + //an error occurred + success = false; + break; + } + + ScalarType v1 = static_cast(outputValue); + feature->statSF1->setValue(i, v1); + } + + if (feature->cloud2 == sourceCloud &&feature->statSF2 && feature->field2) + { + assert(feature->op != Feature::NO_OPERATION); + if (!feature->computeStat(nNSS.pointsInNeighbourhood, feature->field2, outputValue)) + { + //an error occurred + success = false; + break; + } + + ScalarType v2 = static_cast(outputValue); + feature->statSF2->setValue(i, v2); + } + } + + if (!success) + { + break; + } + + if (progressCb && !nProgress.oneStep()) + { + //process cancelled by the user + ccLog::Warning("Process cancelled"); + error = true; + break; + } + + } //for each scale + + } //for each point + + } //for each cloud + + //now we can end + } + + return success; } bool Tools::RandomSubset(ccPointCloud* cloud, float ratio, CCLib::ReferenceCloud* inRatioSubset, CCLib::ReferenceCloud* outRatioSubset)