Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
19 changes: 18 additions & 1 deletion include/PointCloudTpl.h
Original file line number Diff line number Diff line change
Expand Up @@ -79,8 +79,16 @@ namespace CCCoreLib
return *this;
}

//! Returns the point vector
inline const std::vector<CCVector3>& points() const
{
return m_points;
}

//! Returns the point vector
inline unsigned size() const override { return static_cast<unsigned>(m_points.size()); }

//! Sets all scalar values in the active 'out' scalar field to a given value
void setPointScalarValues(ScalarType value) override
{
ScalarField* currentOutScalarFieldArray = getCurrentOutScalarField();
Expand All @@ -95,6 +103,7 @@ namespace CCCoreLib
currentOutScalarFieldArray->fill(value);
}

//! Returns the bounding box of the cloud
void getBoundingBox(CCVector3& bbMin, CCVector3& bbMax) override
{
if (!m_bbox.isValid())
Expand All @@ -110,10 +119,13 @@ namespace CCCoreLib
bbMax = m_bbox.maxCorner();
}

//! Places the point iterator at the beginning of the cloud
void placeIteratorAtBeginning() override { m_currentPointIndex = 0; }

//! Returns the next point in the cloud (or nullptr if we are at the end)
const CCVector3* getNextPoint() override { return (m_currentPointIndex < m_points.size() ? point(m_currentPointIndex++) : 0); }

//! Enables the (default) scalar field for this cloud
bool enableScalarField() override
{
if (m_points.empty() && m_points.capacity() == 0)
Expand Down Expand Up @@ -164,6 +176,7 @@ namespace CCCoreLib
}
}

//! Returns whether the scalar field is enabled or not
bool isScalarFieldEnabled() const override
{
ScalarField* currentInScalarFieldArray = getCurrentInScalarField();
Expand All @@ -176,6 +189,7 @@ namespace CCCoreLib
return (sfValuesCount != 0 && sfValuesCount >= m_points.size());
}

//! Sets a scalar value to the active 'in' scalar field for a given point
void setPointScalarValue(unsigned pointIndex, ScalarType value) override
{
assert(m_currentInScalarFieldIndex >= 0 && m_currentInScalarFieldIndex < static_cast<int>(m_scalarFields.size()));
Expand All @@ -188,16 +202,19 @@ namespace CCCoreLib
m_scalarFields[m_currentInScalarFieldIndex]->setValue(pointIndex, value);
}

//! Returns the scalar value of the active 'out' scalar field for a given point
ScalarType getPointScalarValue(unsigned pointIndex) const override
{
assert(m_currentOutScalarFieldIndex >= 0 && m_currentOutScalarFieldIndex < static_cast<int>(m_scalarFields.size()));

return m_scalarFields[m_currentOutScalarFieldIndex]->getValue(pointIndex);
}

//! Returns the point at a given index
inline const CCVector3* getPoint(unsigned index) const override { return point(index); }
//! Returns the point at a given index
inline void getPoint(unsigned index, CCVector3& P) const override { P = *point(index); }

//! Returns the point at a given index (as a persitent pointer)
inline const CCVector3* getPointPersistentPtr(unsigned index) const override { return point(index); }

//! Adds a scalar values to the active 'in' scalar field
Expand Down
110 changes: 75 additions & 35 deletions src/StatisticalTestingTools.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,20 +17,20 @@

using namespace CCCoreLib;

//! Max computable Chi2 distance
static double CHI2_MAX = 1e7;
static const double CHI2_MAX = 1e7; //!< Max computable Chi2 distance

//! An element of a double-chained-list structure (used by computeAdaptativeChi2Dist)
struct Chi2Class
{

double pi; /**< Probability Pi **/
int n; /**< Number of elements for the class **/
double pi; //!< Probability Pi
int n; //!< Number of elements for the class

//! Default constructor
Chi2Class() : pi(0.0) , n(0) {}
//! Constructor from parameters
Chi2Class(double _pi, int _n) : pi(_pi) , n(_n) {}
Chi2Class(double _pi = 0.0, int _n = 0)
: pi(_pi)
, n(_n)
{
}

};

Expand All @@ -50,8 +50,10 @@ double StatisticalTestingTools::computeAdaptativeChi2Dist( const GenericDistribu
assert(distrib && cloud);
unsigned n = cloud->size();

if (n==0 || !distrib->isValid())
if (n == 0 || !distrib->isValid())
{
return -1.0;
}

//compute min and max (valid) values
ScalarType minV = 0;
Expand All @@ -72,22 +74,32 @@ double StatisticalTestingTools::computeAdaptativeChi2Dist( const GenericDistribu
else
{
if (V > maxV)
{
maxV = V;
}
else if (V < minV)
{
minV = V;
}
}
++numberOfValidValues;
}
}
}

if (numberOfValidValues == 0)
{
return -1.0;
}

if (histoMin)
{
minV = *histoMin;
}
if (histoMax)
{
maxV = *histoMax;
}

//shall we automatically compute the number of classes?
if (numberOfClasses == 0)
Expand All @@ -96,7 +108,7 @@ double StatisticalTestingTools::computeAdaptativeChi2Dist( const GenericDistribu
}
if (numberOfClasses < 2)
{
return -2.0; //not enough points/classes
return -2.0; // not enough points/classes
}

//try to allocate the histogram values array (if necessary)
Expand All @@ -106,13 +118,13 @@ double StatisticalTestingTools::computeAdaptativeChi2Dist( const GenericDistribu
//not enough memory
return -1.0;
}
memset(histo, 0, sizeof(unsigned)*numberOfClasses);
memset(histo, 0, sizeof(unsigned) * numberOfClasses);

//accumulate histogram
ScalarType dV = maxV - minV;
unsigned histoBefore = 0;
unsigned histoAfter = 0;
if ( GreaterThanEpsilon( dV ) )
if (GreaterThanEpsilon(dV))
{
for (unsigned i = 0; i < n; ++i)
{
Expand All @@ -127,9 +139,13 @@ double StatisticalTestingTools::computeAdaptativeChi2Dist( const GenericDistribu
else if (bin >= static_cast<int>(numberOfClasses))
{
if (V > maxV)
{
histoAfter++;
}
else
{
histo[numberOfClasses - 1]++;
}
}
else
{
Expand Down Expand Up @@ -169,7 +185,9 @@ double StatisticalTestingTools::computeAdaptativeChi2Dist( const GenericDistribu
currentClass.n = histo[k - 1];
currentClass.pi = p2 - p1;
if (npis)
{
npis[k - 1] = currentClass.pi * numberOfValidValues;
}

try
{
Expand Down Expand Up @@ -209,11 +227,17 @@ double StatisticalTestingTools::computeAdaptativeChi2Dist( const GenericDistribu
Chi2ClassList::iterator it = classes.begin();
Chi2ClassList::iterator minIt = it;
for (; it != classes.end(); ++it)
{
if (it->pi < minIt->pi)
{
minIt = it;
}
}

if (minIt->pi >= minPi) //all classes are bigger than the minimum requirement
if (minIt->pi >= minPi) // all classes are bigger than the minimum requirement
{
break;
}

//otherwise we must merge the smallest class with its neighbor (to make the classes repartition more equilibrated)
Chi2ClassList::iterator smallestIt;
Expand Down Expand Up @@ -263,7 +287,9 @@ double StatisticalTestingTools::computeAdaptativeChi2Dist( const GenericDistribu
}

if (!histoValues)
{
delete[] histo;
}

finalNumberOfClasses = static_cast<unsigned>(classes.size());

Expand All @@ -272,12 +298,12 @@ double StatisticalTestingTools::computeAdaptativeChi2Dist( const GenericDistribu

double StatisticalTestingTools::computeChi2Fractile(double p, int d)
{
return Chi2Helper::critchi(p,d);
return Chi2Helper::critchi(p, d);
}

double StatisticalTestingTools::computeChi2Probability(double chi2result, int d)
{
return Chi2Helper::pochisq(chi2result,d);
return Chi2Helper::pochisq(chi2result, d);
}

double StatisticalTestingTools::testCloudWithStatisticalModel(const GenericDistribution* distrib,
Expand All @@ -290,7 +316,9 @@ double StatisticalTestingTools::testCloudWithStatisticalModel(const GenericDistr
assert(theCloud);

if (!distrib->isValid())
{
return -1.0;
}

DgmOctree* theOctree = inputOctree;
if (!theOctree)
Expand All @@ -307,7 +335,9 @@ double StatisticalTestingTools::testCloudWithStatisticalModel(const GenericDistr
if (!theCloud->enableScalarField())
{
if (!inputOctree)
{
delete theOctree;
}
return -3.0;
}

Expand All @@ -325,7 +355,9 @@ double StatisticalTestingTools::testCloudWithStatisticalModel(const GenericDistr
{
//not enough memory
if (!inputOctree)
{
delete theOctree;
}
return -3.0;
}

Expand All @@ -351,12 +383,12 @@ double StatisticalTestingTools::testCloudWithStatisticalModel(const GenericDistr
}

//additional parameters for local process
void* additionalParameters[] = { reinterpret_cast<void*>(const_cast<GenericDistribution*>(distrib)),
reinterpret_cast<void*>(&numberOfNeighbours),
reinterpret_cast<void*>(&numberOfChi2Classes),
reinterpret_cast<void*>(histoValues.data()),
reinterpret_cast<void*>(histoMin),
reinterpret_cast<void*>(histoMax) };
void* additionalParameters[] { reinterpret_cast<void*>(const_cast<GenericDistribution*>(distrib)),
reinterpret_cast<void*>(&numberOfNeighbours),
reinterpret_cast<void*>(&numberOfChi2Classes),
reinterpret_cast<void*>(histoValues.data()),
reinterpret_cast<void*>(histoMin),
reinterpret_cast<void*>(histoMax) };

double maxChi2 = -1.0;

Expand All @@ -373,13 +405,15 @@ double StatisticalTestingTools::testCloudWithStatisticalModel(const GenericDistr
if (!progressCb || !progressCb->isCancelRequested())
{
//theoretical Chi2 fractile
maxChi2 = computeChi2Fractile(pTrust, numberOfChi2Classes-1);
maxChi2 = computeChi2Fractile(pTrust, numberOfChi2Classes - 1);
maxChi2 = sqrt(maxChi2); //on travaille avec les racines carrees des distances du Chi2
}
}

if (!inputOctree)
{
delete theOctree;
}

return maxChi2;
}
Expand All @@ -388,22 +422,22 @@ bool StatisticalTestingTools::computeLocalChi2DistAtLevel( const DgmOctree::octr
void** additionalParameters,
NormalizedProgress* nProgress/*=nullptr*/)
{
//variables additionnelles
GenericDistribution* statModel = reinterpret_cast<GenericDistribution*>(additionalParameters[0]);
unsigned numberOfNeighbours = *reinterpret_cast<unsigned*>(additionalParameters[1]);
unsigned numberOfChi2Classes = *reinterpret_cast<unsigned*>(additionalParameters[2]);
unsigned* histoValues = reinterpret_cast<unsigned*>(additionalParameters[3]);
ScalarType* histoMin = reinterpret_cast<ScalarType*>(additionalParameters[4]);
ScalarType* histoMax = reinterpret_cast<ScalarType*>(additionalParameters[5]);
// variables additionnelles
GenericDistribution* statModel = reinterpret_cast<GenericDistribution*>(additionalParameters[0]);
unsigned numberOfNeighbours = *reinterpret_cast<unsigned*>(additionalParameters[1]);
unsigned numberOfChi2Classes = *reinterpret_cast<unsigned*>(additionalParameters[2]);
unsigned* histoValues = reinterpret_cast<unsigned*>(additionalParameters[3]);
ScalarType* histoMin = reinterpret_cast<ScalarType*>(additionalParameters[4]);
ScalarType* histoMax = reinterpret_cast<ScalarType*>(additionalParameters[5]);

//number of points in the current cell
unsigned n = cell.points->size();

DgmOctree::NearestNeighboursSearchStruct nNSS;
nNSS.level = cell.level;
nNSS.minNumberOfNeighbors = numberOfNeighbours;
cell.parentOctree->getCellPos(cell.truncatedCode,cell.level,nNSS.cellPos,true);
cell.parentOctree->computeCellCenter(nNSS.cellPos,cell.level,nNSS.cellCenter);
nNSS.level = cell.level;
nNSS.minNumberOfNeighbors = numberOfNeighbours;
cell.parentOctree->getCellPos(cell.truncatedCode, cell.level, nNSS.cellPos, true);
cell.parentOctree->computeCellCenter(nNSS.cellPos, cell.level, nNSS.cellCenter);

//we already know the points of the first cell (this is the one we are currently processing!)
{
Expand All @@ -417,9 +451,9 @@ bool StatisticalTestingTools::computeLocalChi2DistAtLevel( const DgmOctree::octr
}

DgmOctree::NeighboursSet::iterator it = nNSS.pointsInNeighbourhood.begin();
for (unsigned j=0;j<n;++j,++it)
for (unsigned j = 0; j < n; ++j, ++it)
{
it->point = cell.points->getPointPersistentPtr(j);
it->point = cell.points->getPointPersistentPtr(j);
it->pointIndex = cell.points->getPointGlobalIndex(j);
}
nNSS.alreadyVisitedNeighbourhoodSize = 1;
Expand All @@ -443,11 +477,15 @@ bool StatisticalTestingTools::computeLocalChi2DistAtLevel( const DgmOctree::octr

unsigned k = cell.parentOctree->findNearestNeighborsStartingFromCell(nNSS, true);
if (k > numberOfNeighbours)
{
k = numberOfNeighbours;
}

neighboursCloud.clear();
for (unsigned j = 0; j < k; ++j)
{
neighboursCloud.addPointIndex(nNSS.pointsInNeighbourhood[j].pointIndex);
}

unsigned finalNumberOfChi2Classes = 0;
//LAZY VERSION (approximate test)
Expand All @@ -462,7 +500,9 @@ bool StatisticalTestingTools::computeLocalChi2DistAtLevel( const DgmOctree::octr
cell.points->setPointScalarValue(i, D);

if (nProgress && !nProgress->oneStep())
{
return false;
}
}

return true;
Expand Down
Loading