Compute features at various scales in a smarter way (work in progress)

PointFeature.cpp

@@ -301,7 +301,6 @@ static bool ExtractStatFromSF(	const CCVector3& queryPoint,
 		assert(false);
 		return false;
 	}
-	std::numeric_limits<double>::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<IScalarFieldWrapper> 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<IScalarFieldWrapper> 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<IScalarFieldWrapper>& sourceField, double& outputValue) const
+{
+	outputValue = std::numeric_limits<double>::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<ScalarType>(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

PointFeature.h

@@ -31,6 +31,8 @@ namespace masc
 	{
 	public:	//PointFeatureType
 
+		typedef QSharedPointer<PointFeature> 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<IScalarFieldWrapper>& 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<IScalarFieldWrapper> field1, field2;
+		CCLib::ScalarField *statSF1, *statSF2;
+		bool keepStatSF2;
 	};
 }

q3DMASCTools.cpp

@@ -792,6 +792,12 @@ CCLib::ScalarField* Tools::RetrieveSF(const ccPointCloud* cloud, const QString&
 	}
 }
 
+struct FeaturesAndScales
+{
+	std::vector<double> scales;
+	std::vector<PointFeature::Shared> 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<ccPointCloud*, FeaturesAndScales> 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<PointFeature>(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<PointFeature>(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<ccPointCloud*, FeaturesAndScales>::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<PointCoordinateType>(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<ScalarType>(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<ScalarType>(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)