PointFeature.cpp

//##########################################################################
//#                                                                        #
//#                     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 "PointFeature.h"

//Local
#include "q3DMASCTools.h"

#if defined(_OPENMP)
#include <omp.h>
#endif

//qPDALIO
#ifdef PLUGIN_IO_QPDAL
#include "../../../core/IO/qPDALIO/include/LASFields.h"
#else
#include "../../../core/IO/qLASIO/include/LasDetails.h"
enum LAS_FIELDS {
	LAS_X = 0,
	LAS_Y = 1,
	LAS_Z = 2,
	LAS_INTENSITY = 3,
	LAS_RETURN_NUMBER = 4,
	LAS_NUMBER_OF_RETURNS = 5,
	LAS_SCAN_DIRECTION = 6,
	LAS_FLIGHT_LINE_EDGE = 7,
	LAS_CLASSIFICATION = 8,
	LAS_SCAN_ANGLE_RANK = 9,
	LAS_USER_DATA = 10,
	LAS_POINT_SOURCE_ID = 11,
	LAS_RED = 12,
	LAS_GREEN = 13,
	LAS_BLUE = 14,
	LAS_TIME = 15,
	LAS_EXTRA = 16,
	//Sub fields
	LAS_CLASSIF_VALUE = 17,
	LAS_CLASSIF_SYNTHETIC = 18,
	LAS_CLASSIF_KEYPOINT = 19,
	LAS_CLASSIF_WITHHELD = 20,
	LAS_CLASSIF_OVERLAP = 21,
	//Invald flag
	LAS_INVALID = 255
};

constexpr const char* LAS_FIELD_NAMES[22] = {"X",
									"Y",
									"Z",
									LasNames::Intensity,
									LasNames::ReturnNumber,
									LasNames::NumberOfReturns,
									LasNames::ScanDirectionFlag,
									LasNames::EdgeOfFlightLine,
									LasNames::Classification,
									LasNames::ScanAngleRank,
									LasNames::UserData,
									LasNames::PointSourceId,
									"Red",
									"Green",
									"Blue",
									LasNames::GpsTime,
									"extra",
									"[Classif] Value",
									"[Classif] Synthetic flag",
									"[Classif] Key-point flag",
									"[Classif] Withheld flag",
									"[Classif] Overlap flag",
};
#endif

//qCC_db
#include <ccPointCloud.h>
#include <ccScalarField.h>

//CCLib
#include <WeibullDistribution.h>

//system
#include <assert.h>

//Qt
#include <QCoreApplication>
#include <QMutex>

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.";

using namespace masc;

bool PointFeature::checkValidity(QString corePointRole, QString &error) const
{
	if (!Feature::checkValidity(corePointRole, error))
	{
		return false;
	}

	if (type == Invalid)
	{
		assert(false);
		error = "invalid feature type";
		return false;
	}

	assert(cloud1);

	if (!corePointRole.isEmpty() && !scaled() && cloud1Label != corePointRole) //in some cases, we don't know the role of the core points yet!
	{
		error = "Scale-less features can only be computed on the core points / classified cloud";
		return false;
	}

	if (scaled() && stat == NO_STAT)
	{
		error = "scaled point features need a STAT measure to be defined";
		return false;
	}

	if (op != NO_OPERATION && !cloud2)
	{
		error = "math operations require two clouds";
		return false;
	}
	
	switch (type)
	{
	case PointFeature::Intensity:
	{
		if (cloud1->getScalarFieldIndexByName(LAS_FIELD_NAMES[LAS_INTENSITY]) < 0)
		{
			error = QString("Cloud %0 has no '%1' scalar field").arg(cloud1->getName()).arg(LAS_FIELD_NAMES[LAS_INTENSITY]);
			return false;
		}
		return true;
	}
	case PointFeature::X:
	case PointFeature::Y:
	case PointFeature::Z:
		return true;
	case PointFeature::NbRet:
	{
		if (cloud1->getScalarFieldIndexByName(LAS_FIELD_NAMES[LAS_NUMBER_OF_RETURNS]) < 0)
		{
			error = QString("Cloud %0 has no '%1' scalar field").arg(cloud1->getName()).arg(LAS_FIELD_NAMES[LAS_NUMBER_OF_RETURNS]);
			return false;
		}
		return true;
	}
	case PointFeature::RetNb:
	{
		if (cloud1->getScalarFieldIndexByName(LAS_FIELD_NAMES[LAS_RETURN_NUMBER]) < 0)
		{
			error = QString("Cloud %0 has no '%1' scalar field").arg(cloud1->getName()).arg(LAS_FIELD_NAMES[LAS_RETURN_NUMBER]);
			return false;
		}
		return true;
	}
	case PointFeature::EchoRat:
	{
		if (cloud1->getScalarFieldIndexByName(LAS_FIELD_NAMES[LAS_NUMBER_OF_RETURNS]) < 0)
		{
			error = QString("Cloud %0 has no '%1' scalar field").arg(cloud1->getName()).arg(LAS_FIELD_NAMES[LAS_NUMBER_OF_RETURNS]);
			return false;
		}
		if (cloud1->getScalarFieldIndexByName(LAS_FIELD_NAMES[LAS_RETURN_NUMBER]) < 0)
		{
			error = QString("Cloud %0 has no '%1' scalar field").arg(cloud1->getName()).arg(LAS_FIELD_NAMES[LAS_RETURN_NUMBER]);
			return false;
		}
		return true;
	}
	case PointFeature::R:
	case PointFeature::G:
	case PointFeature::B:
		if (!cloud1->hasColors())
		{
			error = QString("Cloud %0 has no RGB color").arg(cloud1->getName());
			return false;
		}
		return true;
	case PointFeature::NIR:
	{
		if (cloud1->getScalarFieldIndexByName(s_NIRSFName) < 0)
		{
			error = QString("Cloud %0 has no '%1' scalar field").arg(cloud1->getName()).arg(s_NIRSFName);
			return false;
		}
		return true;
	}
	case PointFeature::Dip:
	case PointFeature::DipDir:
	{
		if (!cloud1->hasNormals())
		{
			error = QString("Cloud %0 has no normals").arg(cloud1->getName());
			return false;
		}
		return true;
	}
	case PointFeature::M3C2:
	{
		if (cloud1->getScalarFieldIndexByName(s_M3C2SFName) < 0)
		{
			error = QString("Cloud %0 has no '%1' scalar field").arg(cloud1->getName()).arg(s_M3C2SFName);
			return false;
		}
		return true;
	}
	case PointFeature::PCV:
	{
		if (cloud1->getScalarFieldIndexByName(s_PCVSFName) < 0)
		{
			error = QString("Cloud %0 has no '%1' scalar field").arg(cloud1->getName()).arg(s_PCVSFName);
			return false;
		}
		return true;
	}
	case PointFeature::SF:
		if (sourceSFIndex >= static_cast<int>(cloud1->getNumberOfScalarFields()))
		{
			error = QString("Cloud %0 has no scalar field #%1").arg(cloud1->getName()).arg(sourceSFIndex);
			return false;
		}
		return true;
	default:
		break;
	}

	return true;
}

IScalarFieldWrapper::Shared PointFeature::retrieveField(ccPointCloud* cloud, QString& error)
{
	if (!cloud)
	{
		assert(false);
		return IScalarFieldWrapper::Shared(nullptr);
	}
	
	switch (type)
	{
	case PointFeature::Intensity:
	{
		CCCoreLib::ScalarField* sf = Tools::RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_INTENSITY], false);
		if (!sf)
		{
			error = "Cloud has no 'intensity' scalar field";
			return nullptr;
		}
		return IScalarFieldWrapper::Shared(new ScalarFieldWrapper(sf));
	}
	case PointFeature::X:
		return IScalarFieldWrapper::Shared(new DimScalarFieldWrapper(cloud, DimScalarFieldWrapper::DimX));
	case PointFeature::Y:
		return IScalarFieldWrapper::Shared(new DimScalarFieldWrapper(cloud, DimScalarFieldWrapper::DimY));
	case PointFeature::Z:
		return IScalarFieldWrapper::Shared(new DimScalarFieldWrapper(cloud, DimScalarFieldWrapper::DimZ));
	case PointFeature::NbRet:
	{
		CCCoreLib::ScalarField* sf = Tools::RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_NUMBER_OF_RETURNS], false);
		if (!sf)
		{
			error = "Cloud has no 'number of returns' scalar field";
			return nullptr;
		}
		return IScalarFieldWrapper::Shared(new ScalarFieldWrapper(sf));
	}
	case PointFeature::RetNb:
	{
		CCCoreLib::ScalarField* sf = Tools::RetrieveSF(cloud, LAS_FIELD_NAMES[LAS_RETURN_NUMBER], false);
		if (!sf)
		{
			error = "Cloud has no 'return number' scalar field";
			return nullptr;
		}
		return IScalarFieldWrapper::Shared(new ScalarFieldWrapper(sf));
	}
	case PointFeature::EchoRat:
	{
		//retrieve the two scalar fields 'p/q'
		CCCoreLib::ScalarField* numberOfRetSF = Tools::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;
		}
		CCCoreLib::ScalarField* retNumberSF = Tools::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;
		}
		return IScalarFieldWrapper::Shared(new ScalarFieldRatioWrapper(retNumberSF, numberOfRetSF, "EchoRat"));
	}
	case PointFeature::R:
		return IScalarFieldWrapper::Shared(new ColorScalarFieldWrapper(cloud, ColorScalarFieldWrapper::Red));
	case PointFeature::G:
		return IScalarFieldWrapper::Shared(new ColorScalarFieldWrapper(cloud, ColorScalarFieldWrapper::Green));
	case PointFeature::B:
		return IScalarFieldWrapper::Shared(new ColorScalarFieldWrapper(cloud, ColorScalarFieldWrapper::Blue));
	case PointFeature::NIR:
	{
		CCCoreLib::ScalarField* sf = Tools::RetrieveSF(cloud, s_NIRSFName, false);
		if (!sf)
		{
			error = "Cloud has no 'NIR' scalar field";
			return nullptr;
		}
		return IScalarFieldWrapper::Shared(new ScalarFieldWrapper(sf));
	}
	case PointFeature::Dip:
	case PointFeature::DipDir:
	{
		//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;
		}
		return IScalarFieldWrapper::Shared(new NormDipAndDipDirFieldWrapper(cloud, type == PointFeature::Dip ? NormDipAndDipDirFieldWrapper::Dip : NormDipAndDipDirFieldWrapper::DipDir));
	}
	case PointFeature::M3C2:
	{
		CCCoreLib::ScalarField* sf = Tools::RetrieveSF(cloud, s_M3C2SFName, true);
		if (!sf)
		{
			error = "Cloud has no 'm3c2 distance' scalar field";
			return nullptr;
		}
		return IScalarFieldWrapper::Shared(new ScalarFieldWrapper(sf));
	}
	case PointFeature::PCV:
	{
		CCCoreLib::ScalarField* sf = Tools::RetrieveSF(cloud, s_PCVSFName, true);
		if (!sf)
		{
			error = "Cloud has no 'PCV/Illuminance' scalar field";
			return nullptr;
		}
		return IScalarFieldWrapper::Shared(new ScalarFieldWrapper(sf));
	}
	case PointFeature::SF:
		if (sourceSFIndex < 0 || sourceSFIndex >= static_cast<int>(cloud->getNumberOfScalarFields()))
		{
			error = QString("Can't retrieve the specified SF: invalid index (%1)").arg(sourceSFIndex);
			return nullptr;
		}
		return IScalarFieldWrapper::Shared(new ScalarFieldWrapper(cloud->getScalarField(sourceSFIndex)));
	default:
		break;
	}

	error = "Unhandled feature type";
	return nullptr;
}

static bool ComputeMathOpWithNearestNeighbor(	const CorePoints& corePoints,
												const IScalarFieldWrapper& field1,
												CCCoreLib::ScalarField* outSF,
												ccPointCloud& cloud2,
												const IScalarFieldWrapper& field2,
												masc::Feature::Operation op,
												QString& error,
												CCCoreLib::GenericProgressCallback* progressCb = nullptr)
{
	if (op == masc::Feature::NO_OPERATION || !outSF || outSF->size() != corePoints.size())
	{
		//invalid input parameters
		assert(false);
		error = "invalid input parameters";
		return false;
	}
	
	ccOctree::Shared octree = cloud2.getOctree();
	if (!octree)
	{
		octree = cloud2.computeOctree(progressCb);
		if (!octree)
		{
			error = "failed to compute octree on cloud " + cloud2.getName();
			return false;
		}
	}

	//now extract the neighborhoods
	unsigned char octreeLevel = octree->findBestLevelForAGivenPopulationPerCell(3);
	ccLog::Print(QString("[Initial octree level] level = %1").arg(octreeLevel));

	unsigned pointCount = corePoints.size();
	QString logMessage = QString("Extracting %1 core points nearest neighbors in cloud %2").arg(pointCount).arg(cloud2.getName());
	if (progressCb)
	{
		progressCb->setMethodTitle("Compute math operation");
		progressCb->setInfo(qPrintable(logMessage));
	}
	ccLog::Print(logMessage);
	CCCoreLib::NormalizedProgress nProgress(progressCb, pointCount);

	double meanNeighborhoodSize = 0;
	int tenth = pointCount / 10;
	error.clear();
	bool cancelled = false;
#ifndef _DEBUG
#if defined(_OPENMP)
#pragma omp parallel for num_threads(std::max(1, omp_get_max_threads() - 2))
#endif
#endif
	for (int i = 0; i < static_cast<int>(pointCount); ++i)
	{
	if (!cancelled)
	{
		const CCVector3* P = corePoints.cloud->getPoint(i);
		CCCoreLib::ReferenceCloud Yk(&cloud2);
		double maxSquareDist = 0;

		ScalarType s = CCCoreLib::NAN_VALUE;

		int neighborhoodSize = 0;
		if (octree->findPointNeighbourhood(P, &Yk, 1, octreeLevel, maxSquareDist, 0.0, &neighborhoodSize) >= 1)
		{
			double s1 = field1.pointValue(corePoints.originIndex(i));
			double s2 = field2.pointValue(Yk.getPointGlobalIndex(0));
			s = masc::Feature::PerformMathOp(s1, s2, op);
		}

		outSF->setValue(i, s);

		if (i && (i % tenth) == 0)
		{
			double density = meanNeighborhoodSize / tenth;
			if (density < 1.1)
			{
				if (octreeLevel + 1 < CCCoreLib::DgmOctree::MAX_OCTREE_LEVEL)
					++octreeLevel;
			}
			else while (density > 2.9)
			{
				if (octreeLevel <= 5)
					break;
				--octreeLevel;
				density /= 2.0;
			}
			ccLog::Print(QString("[Adaptative octree level] Mean neighborhood size: %1 --> new level = %2").arg(meanNeighborhoodSize / tenth).arg(octreeLevel));
			meanNeighborhoodSize = 0;
		}
		else
		{
			meanNeighborhoodSize += neighborhoodSize;
		}

		if (progressCb)
		{
			cancelled = !nProgress.oneStep();
			if (cancelled)
			{
				//process cancelled by the user
				error = "[Point feature] Process cancelled";
			}
		}
	}
	}

	outSF->computeMinAndMax();

	if (progressCb)
	{
		progressCb->stop();
	}

	return error.isEmpty();
}

bool PointFeature::prepare(	const CorePoints& corePoints,
							QString& error,
							CCCoreLib::GenericProgressCallback* progressCb/*=nullptr*/,
							SFCollector* generatedScalarFields/*=nullptr*/)
{
	if (!cloud1 || !corePoints.cloud)
	{
		//invalid input
		assert(false);
		error = "internal error (no input core points)";
		return false;
	}

	//look for the source field
	assert(!field1);
	field1 = retrieveField(cloud1, error);
	if (!field1)
	{
		//error should be up to date
		return false;
	}

	assert(!field2);
	if (cloud2)
	{
		//no need to compute the second scalar field if no MATH operation has to be performed?!
		if (op != Feature::NO_OPERATION)
		{
			field2 = retrieveField(cloud2, error);
			if (!field2)
			{
				//error should be up to date
				return false;
			}
		}
		else
		{
			assert(false);
			error = "Feature has a second cloud associated but no MATH operation is defined";
			return false;
		}
	}

	bool isScaled = scaled();

	//build the final SF name
	QString resultSF1Name = field1->getName();
	if (cloud2 || corePoints.role != cloud1Label)
	{
		resultSF1Name += "_" + cloud1Label;
	}

	if (isScaled)
	{
		//shall we extract a statistical measure? (mandatory for scaled feature)
		if (stat == Feature::NO_STAT)
		{
			assert(false);
			error = "Scaled features (SCx) must have an associated STAT measure";
			return false;
		}
		resultSF1Name += QString("_") + Feature::StatToString(stat);
	}
	else //not scaled
	{
		if (cloud1 != corePoints.cloud && cloud1 != corePoints.origin)
		{
			assert(false);
			error = "Scale-less features (SC0) can only be defined on the core points (origin) cloud";
			return false;
		}
	}

	if (field2 && op != Feature::NO_OPERATION)
	{
		//include the math operation as well if necessary!
		resultSF1Name += "_" + Feature::OpToString(op) + "_" + field2->getName() + "_" + cloud2Label;
		if (isScaled)
		{
			assert(stat != Feature::NO_STAT);
			resultSF1Name += QString("_") + Feature::StatToString(stat);
		}
	}

	if (isScaled)
	{
		resultSF1Name += "@" + QString::number(scale);

		//prepare the corresponding scalar field
		sf1WasAlreadyExisting = CheckSFExistence(corePoints.cloud, qPrintable(resultSF1Name));
		if (sf1WasAlreadyExisting)
		{
			// if the SF exists, it is not added to generatedScalarFields
			statSF1 = PrepareSF(corePoints.cloud, qPrintable(resultSF1Name), generatedScalarFields, SFCollector::ALWAYS_KEEP);
			if (generatedScalarFields->scalarFields.contains(statSF1)) // i.e. the SF is existing but was not present at the startup of the plugin
				generatedScalarFields->setBehavior(statSF1, SFCollector::CAN_REMOVE);
		}
		else
			statSF1 = PrepareSF(corePoints.cloud, qPrintable(resultSF1Name), generatedScalarFields, SFCollector::CAN_REMOVE);
		if (!statSF1)
		{
			error = QString("Failed to prepare scalar field for field '%1' @ scale %2").arg(field1->getName()).arg(scale);
			return false;
		}
		source.name = statSF1->getName();

		if (field2 && op != Feature::NO_OPERATION && !sf1WasAlreadyExisting) // nothing to do if statSF1 was already there
		{
			QString resultSF2Name = field2->getName() + QString("_") + cloud2Label + "_" + Feature::StatToString(stat) + "@" + QString::number(scale);
			//keepStatSF2 = (corePoints.cloud->getScalarFieldIndexByName(qPrintable(resultSFName2)) >= 0); //we remember that the scalar field was already existing!

			assert(!statSF2);
			sf2WasAlreadyExisting = CheckSFExistence(corePoints.cloud, qPrintable(resultSF2Name));
			if (sf2WasAlreadyExisting)
				statSF2 = PrepareSF(corePoints.cloud, qPrintable(resultSF2Name), generatedScalarFields, SFCollector::ALWAYS_KEEP);
			else
				statSF2 = PrepareSF(corePoints.cloud, qPrintable(resultSF2Name), generatedScalarFields, SFCollector::ALWAYS_REMOVE);
			if (!statSF2)
			{
				error = QString("Failed to prepare scalar field for field '%1' @ scale %2").arg(field2->getName()).arg(scale);
				return false;
			}
		}

		return true;
	}
	else //non scaled feature
	{
		assert(cloud1 == corePoints.cloud || cloud1 == corePoints.origin);

		//retrieve/create a SF to host the result
		int sfIdx = corePoints.cloud->getScalarFieldIndexByName(qPrintable(resultSF1Name));

		CCCoreLib::ScalarField* resultSF = nullptr;
		if (sfIdx >= 0)
		{
			//reuse the existing field
			resultSF = corePoints.cloud->getScalarField(sfIdx);
		}
		else
		{
			//copy the SF1 field
			resultSF = new ccScalarField(qPrintable(resultSF1Name));
			if (!resultSF->resizeSafe(corePoints.cloud->size()))
			{
				error = "Not enough memory";
				resultSF->release();
				return false;
			}

			if (op == NO_OPERATION)
			{
				//simply copy the values
				for (unsigned i = 0; i < corePoints.size(); ++i)
				{
					resultSF->setValue(i, field1->pointValue(corePoints.originIndex(i)));
				}
				resultSF->computeMinAndMax();
			}
			else if (field2)
			{
				if (!ComputeMathOpWithNearestNeighbor(	corePoints,
														*field1,
														resultSF,
														*cloud2,
														*field2,
														op,
														error,
														progressCb)
					)
				{
					error = "Failed to perform the MATH operation (" + error + ")";
					resultSF->release();
					return false;
				}
			}

			int newSFIdx = corePoints.cloud->addScalarField(static_cast<ccScalarField*>(resultSF));
			if (generatedScalarFields)
			{
				//track the generated scalar-field
				generatedScalarFields->push(corePoints.cloud, resultSF, SFCollector::CAN_REMOVE);
			}

			corePoints.cloud->setCurrentDisplayedScalarField(newSFIdx);
		}

		source.name = resultSF->getName();

		return true;
	}
}

bool PointFeature::computeStat(const CCCoreLib::DgmOctree::NeighboursSet& pointsInNeighbourhood, const IScalarFieldWrapper::Shared& 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::MEDIAN || stat == Feature::MODE || stat == Feature::SKEW);
		double sum = 0.0;
		double sum2 = 0.0;

		CCCoreLib::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:
		{
			CCCoreLib::WeibullDistribution w;
			if (w.computeParameters(values))
				outputValue = w.computeMode();
		}
		break;

		case Feature::MEDIAN:
		{
			size_t medianIndex = values.size() / 2;
			std::nth_element(values.begin(), values.begin() + medianIndex, values.end());
			outputValue = values[medianIndex];
		}
		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:
		{
			CCCoreLib::WeibullDistribution w;
			if (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();

		//update display
		//if (corePoints.cloud->getDisplay())
		{
			int sfIndex1 = corePoints.cloud->getScalarFieldIndexByName(statSF1->getName());
			corePoints.cloud->setCurrentDisplayedScalarField(sfIndex1);
			//corePoints.cloud->getDisplay()->redraw();
			//QCoreApplication::processEvents();
		}
	}

	if (statSF2 && !sf1WasAlreadyExisting)
	{
		//now perform the math operation
		if (op != Feature::NO_OPERATION)
		{
			if (!PerformMathOp(statSF1, statSF2, op))
			{
				error = "Failed to perform the MATH operation";
				success = false;
			}
		}
//		statSF2->computeMinAndMax();

		//DGM: we don't delete it now! As it could be used by other features!
		//if (!keepStatSF2)
		//{
		//	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
	QString description = ToString(type);

	//special case for the 'SF' type
	if (type == SF)
	{
		//'SF' + sf index
		description += QString::number(sourceSFIndex);
	}

	if (scaled())
	{
		description += QString("_SC%1_%2").arg(scale).arg(StatToString(stat));
	}
	else
	{
		description += "_SC0";
	}

	description += "_" + cloud1Label;

	if (cloud2 && !cloud2Label.isEmpty())
	{
		description += "_" + cloud2Label;

		if (op != NO_OPERATION)
		{
			description += "_" + OpToString(op);
		}
	}

	//Point features always have a scale equal to 0 by definition
	return description;
}