Commit b3328534 authored by Andrew Bauer's avatar Andrew Bauer
Browse files

WIP: Improve implementation of Kd-tree

Use double instead float for storing coordinate location and
reduce number of casts betweeen floats and doubles in parts
of Kd tree implementation.
parent aa4975e2
Pipeline #134099 failed with stage
......@@ -122,94 +122,6 @@ void vtkKdNode::SetDataBounds(double x1,double x2,double y1,double y2,
this->MinVal[2] = z1; this->MaxVal[2] = z2;
}
// ----------------------------------------------------------------------------
void vtkKdNode::SetDataBounds(float *v)
{
int x;
double newbounds[6];
vtkIdType numPoints = this->GetNumberOfPoints();
int i;
if (this->Up)
{
double bounds[6];
this->Up->GetDataBounds(bounds);
int dim = this->Up->GetDim();
for (i=0; i<3; i++)
{
if (i == dim)
{
continue;
}
newbounds[i*2] = bounds[i*2];
newbounds[i*2+1] = bounds[i*2+1];
}
newbounds[dim*2] = newbounds[dim*2+1] = static_cast<double>(v[dim]);
for (i = dim+3; i< numPoints*3; i+=3)
{
if (v[i] < newbounds[dim*2])
{
newbounds[dim*2] = static_cast<double>(v[i]);
}
else if (v[i] > newbounds[dim*2+1])
{
newbounds[dim*2+1] = static_cast<double>(v[i]);
}
}
}
else
{
for (i=0; i<3; i++)
{
newbounds[i*2] = newbounds[i*2+1] = static_cast<double>(v[i]);
}
for (x = 3; x< numPoints*3; x+=3)
{
int y=x+1;
int z=x+2;
if (v[x] < newbounds[0])
{
newbounds[0] = static_cast<double>(v[x]);
}
else if (v[x] > newbounds[1])
{
newbounds[1] = static_cast<double>(v[x]);
}
if (v[y] < newbounds[2])
{
newbounds[2] = static_cast<double>(v[y]);
}
else if (v[y] > newbounds[3])
{
newbounds[3] = static_cast<double>(v[y]);
}
if (v[z] < newbounds[4])
{
newbounds[4] = static_cast<double>(v[z]);
}
else if (v[z] > newbounds[5])
{
newbounds[5] = static_cast<double>(v[z]);
}
}
}
this->SetDataBounds(newbounds[0], newbounds[1], newbounds[2],
newbounds[3], newbounds[4], newbounds[5]);
}
// ----------------------------------------------------------------------------
void vtkKdNode::GetDataBounds(double *b) const
{
......
......@@ -98,7 +98,8 @@ public:
* Given a pointer to NumberOfPoints points, set the DataBounds of this
* node to the bounds of these points.
*/
void SetDataBounds(float *v);
template <class TYPE>
void SetDataBounds(TYPE *v);
/**
* Get a pointer to the 3 bound minima (xmin, ymin and zmin) or the
......@@ -310,4 +311,93 @@ private:
void operator=(const vtkKdNode&) = delete;
};
// ----------------------------------------------------------------------------
template <class TYPE>
void vtkKdNode::SetDataBounds(TYPE *v)
{
int x;
double newbounds[6];
vtkIdType numPoints = this->GetNumberOfPoints();
int i;
if (this->Up)
{
double bounds[6];
this->Up->GetDataBounds(bounds);
int dim = this->Up->GetDim();
for (i=0; i<3; i++)
{
if (i == dim)
{
continue;
}
newbounds[i*2] = bounds[i*2];
newbounds[i*2+1] = bounds[i*2+1];
}
newbounds[dim*2] = newbounds[dim*2+1] = static_cast<double>(v[dim]);
for (i = dim+3; i< numPoints*3; i+=3)
{
if (v[i] < newbounds[dim*2])
{
newbounds[dim*2] = static_cast<double>(v[i]);
}
else if (v[i] > newbounds[dim*2+1])
{
newbounds[dim*2+1] = static_cast<double>(v[i]);
}
}
}
else
{
for (i=0; i<3; i++)
{
newbounds[i*2] = newbounds[i*2+1] = static_cast<double>(v[i]);
}
for (x = 3; x< numPoints*3; x+=3)
{
int y=x+1;
int z=x+2;
if (v[x] < newbounds[0])
{
newbounds[0] = static_cast<double>(v[x]);
}
else if (v[x] > newbounds[1])
{
newbounds[1] = static_cast<double>(v[x]);
}
if (v[y] < newbounds[2])
{
newbounds[2] = static_cast<double>(v[y]);
}
else if (v[y] > newbounds[3])
{
newbounds[3] = static_cast<double>(v[y]);
}
if (v[z] < newbounds[4])
{
newbounds[4] = static_cast<double>(v[z]);
}
else if (v[z] > newbounds[5])
{
newbounds[5] = static_cast<double>(v[z]);
}
}
}
this->SetDataBounds(newbounds[0], newbounds[1], newbounds[2],
newbounds[3], newbounds[4], newbounds[5]);
}
#endif
This diff is collapsed.
......@@ -499,7 +499,7 @@ public:
* You must have called BuildLocatorFromPoints() before calling this.
* You are responsible for deleting the returned array.
*/
vtkIdTypeArray *BuildMapForDuplicatePoints(float tolerance);
vtkIdTypeArray *BuildMapForDuplicatePoints(double tolerance);
//@{
/**
......@@ -688,28 +688,24 @@ protected:
* Get back a list of the nodes at a specified level, nodes must
* be preallocated to hold 2^^(level) node structures.
*/
void GetRegionsAtLevel(int level, vtkKdNode **nodes);
/**
* Adds to the vtkIntArray the list of region IDs of all leaf
* nodes in the given node.
*/
static void GetLeafNodeIds(vtkKdNode *node, vtkIntArray *ids);
/**
* Returns the total number of cells in all the data sets
*/
int GetNumberOfCells();
/**
* Returns the total number of cells in data set 1 through
* data set 2.
*/
int GetDataSetsNumberOfCells(int set1, int set2);
/**
......@@ -717,8 +713,6 @@ protected:
* nullptr, the first DataSet is used. This is the point used in
* determining to which spatial region the cell is assigned.
*/
void ComputeCellCenter(vtkDataSet *set, int cellId, float *center);
void ComputeCellCenter(vtkDataSet *set, int cellId, double *center);
/**
......@@ -729,10 +723,9 @@ protected:
* DataSets are returned. The caller should free the list of
* cell centers when done.
*/
float *ComputeCellCenters();
float *ComputeCellCenters(int set);
float *ComputeCellCenters(vtkDataSet *set);
double *ComputeCellCenters();
double *ComputeCellCenters(int set);
double *ComputeCellCenters(vtkDataSet *set);
vtkDataSetCollection *DataSets;
......@@ -783,9 +776,9 @@ protected:
// Recursive helper for public FindPointsInArea
void AddAllPointsInRegion(vtkKdNode* node, vtkIdTypeArray* ids);
int DivideRegion(vtkKdNode *kd, float *c1, int *ids, int nlevels);
int DivideRegion(vtkKdNode *kd, double *c1, int *ids, int nlevels);
void DoMedianFind(vtkKdNode *kd, float *c1, int *ids, int d1, int d2, int d3);
void DoMedianFind(vtkKdNode *kd, double *c1, int *ids, int d1, int d2, int d3);
void SelfRegister(vtkKdNode *kd);
......@@ -815,12 +808,12 @@ protected:
void _printTree(int verbose);
int SearchNeighborsForDuplicate(int regionId, float *point,
int SearchNeighborsForDuplicate(int regionId, double *point,
int **pointsSoFar, int *len,
float tolerance, float tolerance2);
double tolerance, double tolerance2);
int SearchRegionForDuplicate(float *point, int *pointsSoFar,
int len, float tolerance2);
int SearchRegionForDuplicate(double *point, int *pointsSoFar,
int len, double tolerance2);
int _FindClosestPointInRegion(int regionId,
double x, double y, double z, double &dist2);
......@@ -855,21 +848,21 @@ protected:
static void __printTree(vtkKdNode *kd, int depth, int verbose);
static int MidValue(int dim, float *c1, int nvals, double &coord);
static int MidValue(int dim, double *c1, int nvals, double &coord);
static int Select(int dim, float *c1, int *ids, int nvals, double &coord);
static float FindMaxLeftHalf(int dim, float *c1, int K);
static void _Select(int dim, float *X, int *ids, int L, int R, int K);
static int Select(int dim, double *c1, int *ids, int nvals, double &coord);
static double FindMaxLeftHalf(int dim, double *c1, int K);
static void _Select(int dim, double *X, int *ids, int L, int R, int K);
static int ComputeLevel(vtkKdNode *kd);
static int SelfOrder(int id, vtkKdNode *kd);
static int findRegion(vtkKdNode *node, float x, float y, float z);
//static int findRegion(vtkKdNode *node, float x, float y, float z);
static int findRegion(vtkKdNode *node, double x, double y, double z);
static vtkKdNode **_GetRegionsAtLevel(int level, vtkKdNode **nodes,
vtkKdNode *kd);
static void AddNewRegions(vtkKdNode *kd, float *c1,
static void AddNewRegions(vtkKdNode *kd, double *c1,
int midpt, int dim, double coord);
void NewPartitioningRequest(int req);
......@@ -899,11 +892,11 @@ protected:
// to find duplicate points. (BuildLocatorFromPoints)
int NumberOfLocatorPoints;
float *LocatorPoints;
double *LocatorPoints;
int *LocatorIds;
int *LocatorRegionLocation;
float MaxWidth;
double MaxWidth;
// These Last* values are here to save state so we can
// determine later if k-d tree must be rebuilt.
......
......@@ -19,6 +19,7 @@ PRIVATE_DEPENDS
VTK::CommonMath
VTK::CommonSystem
VTK::CommonTransforms
VTK::ParallelCore
TEST_DEPENDS
VTK::FiltersExtraction
VTK::FiltersFlowPaths
......
......@@ -37,6 +37,7 @@
#include <algorithm>
#include <cstdlib>
#include <map>
#include "vtkMultiProcessController.h"
vtkStandardNewMacro(vtkMergeCells);
......@@ -597,7 +598,7 @@ vtkIdType* vtkMergeCells::MapPointsToIdsUsingGlobalIds(vtkDataSet* set)
//-----------------------------------------------------------------------------
// Use a spatial locator to filter out duplicate points and map
// the new ids to their ids in the merged grid.
vtkIdType* vtkMergeCells::MapPointsToIdsUsingLocator(vtkDataSet* set)
vtkIdType* vtkMergeCells::MapPointsToIdsUsingLocator(vtkDataSet* set) // have an unstructured grid
{
vtkUnstructuredGrid* grid = this->UnstructuredGrid;
vtkPoints* points0 = grid->GetPoints();
......@@ -709,8 +710,10 @@ vtkIdType* vtkMergeCells::MapPointsToIdsUsingLocator(vtkDataSet* set)
kd->BuildLocatorFromPoints(ptArrays, numArrays);
int pid = vtkMultiProcessController::GetGlobalController()->GetLocalProcessId();
vtkIdTypeArray* pointToEquivClassMap =
kd->BuildMapForDuplicatePoints(this->PointMergeTolerance);
kd->BuildMapForDuplicatePoints(this->PointMergeTolerance); // called from here...
kd->Delete();
......
......@@ -710,7 +710,7 @@ int vtkPKdTree::DivideRegion(vtkKdNode *kd, int L, int level, int tag)
double bounds[6];
kd->GetBounds(bounds);
float *val = this->GetLocalVal(L);
double *val = this->GetLocalVal(L);
double coord;
if (numpoints > 0)
......@@ -799,7 +799,7 @@ int vtkPKdTree::DivideRegion(vtkKdNode *kd, int L, int level, int tag)
maxdim = newdim;
}
float newDataBounds[12];
double newDataBounds[12];
this->GetDataBounds(L, midpt, R, newDataBounds);
vtkKdNode *left = vtkKdNode::New();
vtkKdNode *right = vtkKdNode::New();
......@@ -853,8 +853,8 @@ void vtkPKdTree::ExchangeVals(int pos1, int pos2)
{
vtkCommunicator *comm = this->Controller->GetCommunicator();
float *myval;
float otherval[3];
double *myval;
double otherval[3];
int player1 = this->WhoHas(pos1);
int player2 = this->WhoHas(pos2);
......@@ -979,9 +979,9 @@ int vtkPKdTree::Select(int dim, int L, int R)
int hasKleft = this->WhoHas(K-1);
int hasKleftrank = this->SubGroup->getLocalRank(hasKleft);
float Kval;
float Kleftval;
float *pt;
double Kval;
double Kleftval;
double *pt;
if (hasK == this->MyId)
{
......@@ -1064,7 +1064,7 @@ int vtkPKdTree::WhoHas(int pos)
}
return _whoHas(0, this->NumProcesses-1, pos);
}
float *vtkPKdTree::GetLocalVal(int pos)
double *vtkPKdTree::GetLocalVal(int pos)
{
if ( (pos < this->StartVal[this->MyId]) || (pos > this->EndVal[this->MyId]))
{
......@@ -1074,7 +1074,7 @@ float *vtkPKdTree::GetLocalVal(int pos)
return this->CurrentPtArray + (3*localPos);
}
float *vtkPKdTree::GetLocalValNext(int pos)
double *vtkPKdTree::GetLocalValNext(int pos)
{
if ( (pos < this->StartVal[this->MyId]) || (pos > this->EndVal[this->MyId]))
{
......@@ -1084,7 +1084,7 @@ float *vtkPKdTree::GetLocalValNext(int pos)
return this->NextPtArray + (3*localPos);
}
void vtkPKdTree::SetLocalVal(int pos, float *val)
void vtkPKdTree::SetLocalVal(int pos, double *val)
{
if ( (pos < this->StartVal[this->MyId]) || (pos > this->EndVal[this->MyId]))
{
......@@ -1100,10 +1100,10 @@ void vtkPKdTree::SetLocalVal(int pos, float *val)
}
void vtkPKdTree::ExchangeLocalVals(int pos1, int pos2)
{
float temp[3];
double temp[3];
float *pt1 = this->GetLocalVal(pos1);
float *pt2 = this->GetLocalVal(pos2);
double *pt1 = this->GetLocalVal(pos1);
double *pt2 = this->GetLocalVal(pos2);
if (!pt1 || !pt2)
{
......@@ -1126,7 +1126,7 @@ void vtkPKdTree::ExchangeLocalVals(int pos1, int pos2)
void vtkPKdTree::DoTransfer(int from, int to, int fromIndex, int toIndex, int count)
{
float *fromPt, *toPt;
double *fromPt, *toPt;
vtkCommunicator *comm = this->Controller->GetCommunicator();
......@@ -1442,9 +1442,9 @@ int *vtkPKdTree::PartitionSubArray(int L, int R, int K, int dim, int p1, int p2)
// that function we know that "T" appears in the array. In this
// function, "T" may or may not appear in the array.
int *vtkPKdTree::PartitionAboutOtherValue(int L, int R, float T, int dim)
int *vtkPKdTree::PartitionAboutOtherValue(int L, int R, double T, int dim)
{
float *Ipt, *Jpt, Lval, Rval;
double *Ipt, *Jpt, Lval, Rval;
int *vals = &this->SelectBuffer[0];
int numTValues = 0;
int numGreater = 0;
......@@ -1653,8 +1653,8 @@ int *vtkPKdTree::PartitionAboutOtherValue(int L, int R, float T, int dim)
int *vtkPKdTree::PartitionAboutMyValue(int L, int R, int K, int dim)
{
float *Ipt, *Jpt;
float T;
double *Ipt, *Jpt;
double T;
int I, J;
int manyTValues = 0;
int *vals = &this->SelectBuffer[0];
......@@ -1665,7 +1665,7 @@ int *vtkPKdTree::PartitionAboutMyValue(int L, int R, int K, int dim)
// X[L] < T and X[R] = T
//
float *pt = this->GetLocalVal(K);
double *pt = this->GetLocalVal(K);
T = pt[dim];
......@@ -1802,7 +1802,7 @@ int *vtkPKdTree::PartitionAboutMyValue(int L, int R, int K, int dim)
//--------------------------------------------------------------------
void vtkPKdTree::GetLocalMinMax(int L, int R, int me,
float *min, float *max)
double *min, double *max)
{
int i, d;
int from = this->StartVal[me];
......@@ -1822,7 +1822,7 @@ void vtkPKdTree::GetLocalMinMax(int L, int R, int me,
from -= this->StartVal[me];
to -= this->StartVal[me];
float *val = this->CurrentPtArray + from*3;
double *val = this->CurrentPtArray + from*3;
for (d=0; d<3; d++)
{
......@@ -1861,16 +1861,16 @@ void vtkPKdTree::GetLocalMinMax(int L, int R, int me,
}
}
}
void vtkPKdTree::GetDataBounds(int L, int K, int R, float globalBounds[12])
void vtkPKdTree::GetDataBounds(int L, int K, int R, double globalBounds[12])
{
float localMinLeft[3]; // Left region is L through K-1
float localMaxLeft[3];
float globalMinLeft[3];
float globalMaxLeft[3];
float localMinRight[3]; // Right region is K through R
float localMaxRight[3];
float globalMinRight[3];
float globalMaxRight[3];
double localMinLeft[3]; // Left region is L through K-1
double localMaxLeft[3];
double globalMinLeft[3];
double globalMaxLeft[3];
double localMinRight[3]; // Right region is K through R
double localMaxRight[3];
double globalMinRight[3];
double globalMaxRight[3];
this->GetLocalMinMax(L, K-1, this->MyId, localMinLeft, localMaxLeft);
......@@ -1888,8 +1888,8 @@ void vtkPKdTree::GetDataBounds(int L, int K, int R, float globalBounds[12])
this->SubGroup->ReduceMax(localMaxRight, globalMaxRight, 3, 0);
this->SubGroup->Broadcast(globalMaxRight, 3, 0);
float *left = globalBounds;
float *right = globalBounds + 6;
double *left = globalBounds;
double *right = globalBounds + 6;
MinMaxToBounds(left, globalMinLeft, globalMaxLeft);
......@@ -2255,7 +2255,7 @@ int vtkPKdTree::AllocateDoubleBuffer()
this->PtArraySize = this->NumCells[this->MyId] * 3;
this->PtArray2 = new float [this->PtArraySize];
this->PtArray2 = new double [this->PtArraySize];
this->CurrentPtArray = this->PtArray;
this->NextPtArray = this->PtArray2;
......@@ -2264,7 +2264,7 @@ int vtkPKdTree::AllocateDoubleBuffer()
}
void vtkPKdTree::SwitchDoubleBuffer()
{
float *temp = this->CurrentPtArray;
double *temp = this->CurrentPtArray;
this->CurrentPtArray = this->NextPtArray;
this->NextPtArray = temp;
......@@ -3155,7 +3155,7 @@ int vtkPKdTree::GetCellArrayGlobalRange(const char *n, double range[2])
}
int vtkPKdTree::GetCellArrayGlobalRange(const char *n, float range[2])
{
double tmp[2] = {0, 0 };
float tmp[2] = {0, 0 };
int fail = this->GetCellArrayGlobalRange(n, tmp);
......@@ -3208,7 +3208,7 @@ int vtkPKdTree::GetPointArrayGlobalRange(const char *n, double range[2])
}
int vtkPKdTree::GetPointArrayGlobalRange(const char *n, float range[2])
{
double tmp[2] = {0, 0};
float tmp[2] = {0, 0};
int fail = this->GetPointArrayGlobalRange(n, tmp);
......@@ -3222,7 +3222,7 @@ int vtkPKdTree::GetPointArrayGlobalRange(const char *n, float range[2])
}
int vtkPKdTree::GetCellArrayGlobalRange(int arrayIndex, float range[2])
{
double tmp[2];
float tmp[2];
int fail = this->GetCellArrayGlobalRange(arrayIndex, tmp);
if (!fail)
{
......@@ -3243,9 +3243,9 @@ int vtkPKdTree::GetCellArrayGlobalRange(int arrayIndex, double range[2])
return 0;
}
int vtkPKdTree::GetPointArrayGlobalRange(int arrayIndex, float range[2])
int vtkPKdTree::GetPointArrayGlobalRange(int arrayIndex, double range[2])
{
double tmp[2];
float tmp[2];
int fail = this->GetPointArrayGlobalRange(arrayIndex, tmp);
if (!fail)
{
......@@ -3254,7 +3254,7 @@ int vtkPKdTree::GetPointArrayGlobalRange(int arrayIndex, float range[2])
}
return fail;
}
int vtkPKdTree::GetPointArrayGlobalRange(int arrayIndex, double range[2])
int vtkPKdTree::GetPointArrayGlobalRange(int arrayIndex, float range[2])
{
if (arrayIndex < 0 || arrayIndex >= this->NumPointArrays || this->PointDataMin.empty())
{
......@@ -3363,8 +3363,8 @@ int vtkPKdTree::GetRegionAssignmentList(int procId, vtkIntArray *list)
return nregions;
}
void vtkPKdTree::GetAllProcessesBorderingOnPoint(float x, float y, float z,
vtkIntArray *list)
void vtkPKdTree::GetAllProcessesBorderingOnPoint(double x, double y, double z,
vtkIntArray *list)