Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge topic 'kdtree'

Browse Source 
37ded1cf Resolve code review issues, 1. cudaDeviceSetLimit issue, 2. Namespace issue
bf11195d Resolve the code review issues
75c794f6 Add 3D kd tree reconstruction and nearest neighbor searh

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Robert Maynard <robert.maynard@kitware.com>
Merge-request: !797
This commit is contained in:
Li-Ta Lo 2017-07-11 15:40:15 +00:00 committed by Kitware Robot
commit 5bf0bfe11a
7 changed files with 1106 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
vtkm/worklet/CMakeLists.txt
View File

@ -36,6 +36,7 @@ set(headers
FieldHistogram.h
FieldStatistics.h
Gradient.h
KdTree3D.h
KernelSplatter.h
Keys.h
Magnitude.h
@ -68,6 +69,7 @@ add_subdirectory(internal)
add_subdirectory(contourtree)
add_subdirectory(gradient)
add_subdirectory(splatkernels)
add_subdirectory(spatialstructure)
add_subdirectory(tetrahedralize)
add_subdirectory(triangulate)
add_subdirectory(wavelets)

69
vtkm/worklet/KdTree3D.h Normal file
View File

@ -0,0 +1,69 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtkm_m_worklet_KdTree3D_h
#define vtkm_m_worklet_KdTree3D_h
#include <vtkm/worklet/spatialstructure/KdTree3DConstruction.h>
#include <vtkm/worklet/spatialstructure/KdTree3DNNSearch.h>
namespace vtkm
{
namespace worklet
{
class KdTree3D
{
public:
KdTree3D() {}
// Execute the 3d kd tree construction given x y z coordinate vectors
// Returns:
// Leaf Node vector and internal node (split) vectpr
template <typename CoordiType, typename TreeIdType, typename DeviceAdapter>
void Run(vtkm::cont::ArrayHandle<vtkm::Vec<CoordiType, 3>>& coordi_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& pointId_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& splitId_Handle,
DeviceAdapter device)
{
vtkm::worklet::spatialstructure::KdTree3DConstruction kdtree3DConstruction;
kdtree3DConstruction.Run(coordi_Handle, pointId_Handle, splitId_Handle, device);
}
// Execute the Neaseat Neighbor Search given kdtree and search points
// Returns:
// Vectors of NN point index and NNpoint distance
template <typename CoordiType, typename TreeIdType, typename DeviceAdapter>
void Run(vtkm::cont::ArrayHandle<vtkm::Vec<CoordiType, 3>>& coordi_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& pointId_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& splitId_Handle,
vtkm::cont::ArrayHandle<vtkm::Vec<CoordiType, 3>>& qc_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& nnId_Handle,
vtkm::cont::ArrayHandle<CoordiType>& nnDis_Handle,
DeviceAdapter device)
{
vtkm::worklet::spatialstructure::KdTree3DNNSearch kdtree3DNNS;
kdtree3DNNS.Run(
coordi_Handle, pointId_Handle, splitId_Handle, qc_Handle, nnId_Handle, nnDis_Handle, device);
}
};
}
} // namespace vtkm::worklet
#endif // vtkm_m_worklet_Kdtree3D_h

26
vtkm/worklet/spatialstructure/CMakeLists.txt Normal file
View File

@ -0,0 +1,26 @@
##============================================================================
## Copyright (c) Kitware, Inc.
## All rights reserved.
## See LICENSE.txt for details.
## This software is distributed WITHOUT ANY WARRANTY; without even
## the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
## PURPOSE. See the above copyright notice for more information.
##
## Copyright 2014 Sandia Corporation.
## Copyright 2014 UT-Battelle, LLC.
## Copyright 2014 Los Alamos National Security.
##
## Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
## the U.S. Government retains certain rights in this software.
##
## Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
## Laboratory (LANL), the U.S. Government retains certain rights in
## this software.
##============================================================================
set(headers
KdTree3DConstruction.h
KdTree3DNNSearch.h
)
vtkm_declare_headers(${headers})

612
vtkm/worklet/spatialstructure/KdTree3DConstruction.h Normal file
View File

@ -0,0 +1,612 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_KdTree3DConstruction_h
#define vtk_m_worklet_KdTree3DConstruction_h
#include <vtkm/Math.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/ArrayHandleReverse.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/arg/ControlSignatureTagBase.h>
#include <vtkm/cont/serial/DeviceAdapterSerial.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/internal/DispatcherBase.h>
#include <vtkm/worklet/internal/WorkletBase.h>
namespace vtkm
{
namespace worklet
{
namespace spatialstructure
{
class KdTree3DConstruction
{
public:
////////// General WORKLET for Kd-tree //////
class ComputeFlag : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> rank, FieldIn<> pointCountInSeg, FieldOut<> flag);
typedef void ExecutionSignature(_1, _2, _3);
VTKM_CONT
ComputeFlag() {}
template <typename T>
VTKM_EXEC void operator()(const T& rank, const T& pointCountInSeg, T& flag) const
{
if (rank >= pointCountInSeg / 2.0f)
flag = 1; //right subtree
else
flag = 0; //left subtree
}
};
class InverseArray : public vtkm::worklet::WorkletMapField
{ //only for 0/1 array
public:
typedef void ControlSignature(FieldIn<> in, FieldOut<> out);
typedef void ExecutionSignature(_1, _2);
VTKM_CONT
InverseArray() {}
template <typename T>
VTKM_EXEC void operator()(const T& in, T& out) const
{
if (in == 0)
out = 1;
else
out = 0;
}
};
class SegmentedSplitTransform : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> B,
FieldIn<> D,
FieldIn<> F,
FieldIn<> G,
FieldIn<> H,
FieldOut<> I);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6);
VTKM_CONT
SegmentedSplitTransform() {}
template <typename T>
VTKM_EXEC void operator()(const T& B, const T& D, const T& F, const T& G, const T& H, T& I)
const
{
if (B == 1)
{
I = F + H + D;
}
else
{
I = F + G - 1;
}
}
};
class ScatterArray : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> in, FieldIn<> index, WholeArrayOut<> out);
typedef void ExecutionSignature(_1, _2, _3);
VTKM_CONT
ScatterArray() {}
template <typename T, typename OutputArrayPortalType>
VTKM_EXEC void operator()(const T& in,
const T& index,
const OutputArrayPortalType& outputPortal) const
{
outputPortal.Set(index, in);
}
};
class NewSegmentId : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> inSegmentId, FieldIn<> flag, FieldOut<> outSegmentId);
typedef void ExecutionSignature(_1, _2, _3);
VTKM_CONT
NewSegmentId() {}
template <typename T>
VTKM_EXEC void operator()(const T& oldSegId, const T& flag, T& newSegId) const
{
if (flag == 0)
newSegId = oldSegId * 2;
else
newSegId = oldSegId * 2 + 1;
}
};
class SaveSplitPointId : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> pointId,
FieldIn<> flag,
FieldIn<> oldSplitPointId,
FieldOut<> newSplitPointId);
typedef void ExecutionSignature(_1, _2, _3, _4);
VTKM_CONT
SaveSplitPointId() {}
template <typename T>
VTKM_EXEC void operator()(const T& pointId,
const T& flag,
const T& oldSplitPointId,
T& newSplitPointId) const
{
if (flag == 0) //do not change
newSplitPointId = oldSplitPointId;
else //split point id
newSplitPointId = pointId;
}
};
class FindSplitPointId : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> pointId, FieldIn<> rank, FieldOut<> splitIdInsegment);
typedef void ExecutionSignature(_1, _2, _3);
VTKM_CONT
FindSplitPointId() {}
template <typename T>
VTKM_EXEC void operator()(const T& pointId, const T& rank, T& splitIdInsegment) const
{
if (rank == 0) //do not change
splitIdInsegment = pointId;
else //split point id
splitIdInsegment = -1; //indicate this is not split point
}
};
class ArrayAdd : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> inArray0, FieldIn<> inArray1, FieldOut<> outArray);
typedef void ExecutionSignature(_1, _2, _3);
VTKM_CONT
ArrayAdd() {}
template <typename T>
VTKM_EXEC void operator()(const T& in0, const T& in1, T& out) const
{
out = in0 + in1;
}
};
class SeprateVec3AryHandle : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> inVec3,
FieldOut<> out0,
FieldOut<> out1,
FieldOut<> out2);
typedef void ExecutionSignature(_1, _2, _3, _4);
VTKM_CONT
SeprateVec3AryHandle() {}
template <typename T>
VTKM_EXEC void operator()(const Vec<T, 3>& inVec3, T& out0, T& out1, T& out2) const
{
out0 = inVec3[0];
out1 = inVec3[1];
out2 = inVec3[2];
}
};
////////// General worklet WRAPPER for Kd-tree //////
template <typename T, class BinaryFunctor, typename DeviceAdapter>
vtkm::cont::ArrayHandle<T> ReverseScanInclusiveByKey(vtkm::cont::ArrayHandle<T>& keyHandle,
vtkm::cont::ArrayHandle<T>& dataHandle,
BinaryFunctor binary_functor,
DeviceAdapter vtkmNotUsed(device))
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
vtkm::cont::ArrayHandle<T> resultHandle;
auto reversedResultHandle = vtkm::cont::make_ArrayHandleReverse(resultHandle);
Algorithm::ScanInclusiveByKey(vtkm::cont::make_ArrayHandleReverse(keyHandle),
vtkm::cont::make_ArrayHandleReverse(dataHandle),
reversedResultHandle,
binary_functor);
return resultHandle;
}
template <typename T, typename DeviceAdapter>
vtkm::cont::ArrayHandle<T> Inverse01ArrayWrapper(vtkm::cont::ArrayHandle<T>& inputHandle,
DeviceAdapter vtkmNotUsed(device))
{
vtkm::cont::ArrayHandle<T> InverseHandle;
InverseArray invWorklet;
vtkm::worklet::DispatcherMapField<InverseArray, DeviceAdapter> InverseArrayDispatcher(
invWorklet);
InverseArrayDispatcher.Invoke(inputHandle, InverseHandle);
return InverseHandle;
}
template <typename T, typename DeviceAdapter>
vtkm::cont::ArrayHandle<T> ScatterArrayWrapper(vtkm::cont::ArrayHandle<T>& inputHandle,
vtkm::cont::ArrayHandle<T>& indexHandle,
DeviceAdapter vtkmNotUsed(device))
{
vtkm::cont::ArrayHandle<T> outputHandle;
outputHandle.Allocate(inputHandle.GetNumberOfValues());
ScatterArray scatterWorklet;
vtkm::worklet::DispatcherMapField<ScatterArray, DeviceAdapter> ScatterArrayDispatcher(
scatterWorklet);
ScatterArrayDispatcher.Invoke(inputHandle, indexHandle, outputHandle);
return outputHandle;
}
template <typename T, typename DeviceAdapter>
vtkm::cont::ArrayHandle<T> NewKeyWrapper(vtkm::cont::ArrayHandle<T>& oldSegIdHandle,
vtkm::cont::ArrayHandle<T>& flagHandle,
DeviceAdapter vtkmNotUsed(device))
{
vtkm::cont::ArrayHandle<T> newSegIdHandle;
NewSegmentId newsegidWorklet;
vtkm::worklet::DispatcherMapField<NewSegmentId, DeviceAdapter> newSegIdDispatcher(
newsegidWorklet);
newSegIdDispatcher.Invoke(oldSegIdHandle, flagHandle, newSegIdHandle);
return newSegIdHandle;
}
template <typename T, typename DeviceAdapter>
vtkm::cont::ArrayHandle<T> SaveSplitPointIdWrapper(vtkm::cont::ArrayHandle<T>& pointIdHandle,
vtkm::cont::ArrayHandle<T>& flagHandle,
vtkm::cont::ArrayHandle<T>& rankHandle,
vtkm::cont::ArrayHandle<T>& oldSplitIdHandle,
DeviceAdapter device)
{
vtkm::cont::ArrayHandle<T> splitIdInSegmentHandle;
FindSplitPointId findSplitPointIdWorklet;
vtkm::worklet::DispatcherMapField<FindSplitPointId, DeviceAdapter>
findSplitPointIdWorkletDispatcher(findSplitPointIdWorklet);
findSplitPointIdWorkletDispatcher.Invoke(pointIdHandle, rankHandle, splitIdInSegmentHandle);
vtkm::cont::ArrayHandle<T> splitIdInSegmentByScanHandle =
ReverseScanInclusiveByKey(flagHandle, splitIdInSegmentHandle, vtkm::Maximum(), device);
vtkm::cont::ArrayHandle<T> splitIdHandle;
SaveSplitPointId saveSplitPointIdWorklet;
vtkm::worklet::DispatcherMapField<SaveSplitPointId, DeviceAdapter>
saveSplitPointIdWorkletDispatcher(saveSplitPointIdWorklet);
saveSplitPointIdWorkletDispatcher.Invoke(
splitIdInSegmentByScanHandle, flagHandle, oldSplitIdHandle, splitIdHandle);
return splitIdHandle;
}
template <typename T, typename DeviceAdapter>
vtkm::cont::ArrayHandle<T> ArrayAddWrapper(vtkm::cont::ArrayHandle<T>& array0Handle,
vtkm::cont::ArrayHandle<T>& array1Handle,
DeviceAdapter vtkmNotUsed(device))
{
vtkm::cont::ArrayHandle<T> resultHandle;
ArrayAdd arrayAddWorklet;
vtkm::worklet::DispatcherMapField<ArrayAdd, DeviceAdapter> arrayAddDispatcher(arrayAddWorklet);
arrayAddDispatcher.Invoke(array0Handle, array1Handle, resultHandle);
return resultHandle;
}
///////////////////////////////////////////////////
////////General Kd tree function //////////////////
///////////////////////////////////////////////////
template <typename T, typename DeviceAdapter>
vtkm::cont::ArrayHandle<T> ComputeFlagProcedure(vtkm::cont::ArrayHandle<T>& rankHandle,
vtkm::cont::ArrayHandle<T>& segIdHandle,
DeviceAdapter device)
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
vtkm::cont::ArrayHandle<T> segCountAryHandle;
{
vtkm::cont::ArrayHandle<T> tmpAryHandle;
vtkm::cont::ArrayHandleConstant<T> constHandle(1, rankHandle.GetNumberOfValues());
Algorithm::ScanInclusiveByKey(
segIdHandle, constHandle, tmpAryHandle, vtkm::Add()); //compute ttl segs in segment
segCountAryHandle =
ReverseScanInclusiveByKey(segIdHandle, tmpAryHandle, vtkm::Maximum(), device);
}
vtkm::cont::ArrayHandle<T> flagHandle;
vtkm::worklet::DispatcherMapField<ComputeFlag, DeviceAdapter> ComputeFlagDispatcher;
ComputeFlagDispatcher.Invoke(rankHandle, segCountAryHandle, flagHandle);
return flagHandle;
}
template <typename T, typename DeviceAdapter>
vtkm::cont::ArrayHandle<T> SegmentedSplitProcedure(vtkm::cont::ArrayHandle<T>& A_Handle,
vtkm::cont::ArrayHandle<T>& B_Handle,
vtkm::cont::ArrayHandle<T>& C_Handle,
DeviceAdapter device)
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
vtkm::cont::ArrayHandle<T> D_Handle;
T initValue = 0;
Algorithm::ScanExclusiveByKey(C_Handle, B_Handle, D_Handle, initValue, vtkm::Add());
vtkm::cont::ArrayHandleCounting<T> Ecouting_Handle(0, 1, A_Handle.GetNumberOfValues());
vtkm::cont::ArrayHandle<T> E_Handle;
Algorithm::Copy(Ecouting_Handle, E_Handle);
vtkm::cont::ArrayHandle<T> F_Handle;
Algorithm::ScanInclusiveByKey(C_Handle, E_Handle, F_Handle, vtkm::Minimum());
vtkm::cont::ArrayHandle<T> InvB_Handle = Inverse01ArrayWrapper(B_Handle, device);
vtkm::cont::ArrayHandle<T> G_Handle;
Algorithm::ScanInclusiveByKey(C_Handle, InvB_Handle, G_Handle, vtkm::Add());
vtkm::cont::ArrayHandle<T> H_Handle =
ReverseScanInclusiveByKey(C_Handle, G_Handle, vtkm::Maximum(), device);
vtkm::cont::ArrayHandle<T> I_Handle;
SegmentedSplitTransform sstWorklet;
vtkm::worklet::DispatcherMapField<SegmentedSplitTransform, DeviceAdapter>
SegmentedSplitTransformDispatcher(sstWorklet);
SegmentedSplitTransformDispatcher.Invoke(
B_Handle, D_Handle, F_Handle, G_Handle, H_Handle, I_Handle);
return ScatterArrayWrapper(A_Handle, I_Handle, device);
}
template <typename T, typename DeviceAdapter>
void RenumberRanksProcedure(vtkm::cont::ArrayHandle<T>& A_Handle,
vtkm::cont::ArrayHandle<T>& B_Handle,
vtkm::cont::ArrayHandle<T>& C_Handle,
vtkm::cont::ArrayHandle<T>& D_Handle,
DeviceAdapter device)
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
vtkm::Id nPoints = A_Handle.GetNumberOfValues();
vtkm::cont::ArrayHandleCounting<T> Ecouting_Handle(0, 1, nPoints);
vtkm::cont::ArrayHandle<T> E_Handle;
Algorithm::Copy(Ecouting_Handle, E_Handle);
vtkm::cont::ArrayHandle<T> F_Handle;
Algorithm::ScanInclusiveByKey(D_Handle, E_Handle, F_Handle, vtkm::Minimum());
vtkm::cont::ArrayHandle<T> G_Handle;
G_Handle = ArrayAddWrapper(A_Handle, F_Handle, device);
vtkm::cont::ArrayHandleConstant<T> HConstant_Handle(1, nPoints);
vtkm::cont::ArrayHandle<T> H_Handle;
Algorithm::Copy(HConstant_Handle, H_Handle);
vtkm::cont::ArrayHandle<T> I_Handle;
T initValue = 0;
Algorithm::ScanExclusiveByKey(C_Handle, H_Handle, I_Handle, initValue, vtkm::Add());
vtkm::cont::ArrayHandle<T> J_Handle;
J_Handle = ScatterArrayWrapper(I_Handle, G_Handle, device);
vtkm::cont::ArrayHandle<T> K_Handle;
K_Handle = ScatterArrayWrapper(B_Handle, G_Handle, device);
vtkm::cont::ArrayHandle<T> L_Handle;
L_Handle = SegmentedSplitProcedure(J_Handle, K_Handle, D_Handle, device);
vtkm::cont::ArrayHandle<T> M_Handle;
Algorithm::ScanInclusiveByKey(C_Handle, E_Handle, M_Handle, vtkm::Minimum());
vtkm::cont::ArrayHandle<T> N_Handle;
N_Handle = ArrayAddWrapper(L_Handle, M_Handle, device);
A_Handle = ScatterArrayWrapper(I_Handle, N_Handle, device);
}
/////////////3D construction /////////////////////
template <typename T, typename DeviceAdapter>
void SegmentedSplitProcedure3D(vtkm::cont::ArrayHandle<T>& A_Handle,
vtkm::cont::ArrayHandle<T>& B_Handle,
vtkm::cont::ArrayHandle<T>& C_Handle,
vtkm::cont::ArrayHandle<T>& X_Handle,
vtkm::cont::ArrayHandle<T>& Y_Handle,
vtkm::cont::ArrayHandle<T>& Z_Handle,
DeviceAdapter device)
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
vtkm::cont::ArrayHandle<T> D_Handle;
T initValue = 0;
Algorithm::ScanExclusiveByKey(C_Handle, B_Handle, D_Handle, initValue, vtkm::Add());
vtkm::cont::ArrayHandleCounting<T> Ecouting_Handle(0, 1, A_Handle.GetNumberOfValues());
vtkm::cont::ArrayHandle<T> E_Handle;
Algorithm::Copy(Ecouting_Handle, E_Handle);
vtkm::cont::ArrayHandle<T> F_Handle;
Algorithm::ScanInclusiveByKey(C_Handle, E_Handle, F_Handle, vtkm::Minimum());
vtkm::cont::ArrayHandle<T> InvB_Handle = Inverse01ArrayWrapper(B_Handle, device);
vtkm::cont::ArrayHandle<T> G_Handle;
Algorithm::ScanInclusiveByKey(C_Handle, InvB_Handle, G_Handle, vtkm::Add());
vtkm::cont::ArrayHandle<T> H_Handle =
ReverseScanInclusiveByKey(C_Handle, G_Handle, vtkm::Maximum(), device);
vtkm::cont::ArrayHandle<T> I_Handle;
SegmentedSplitTransform sstWorklet;
vtkm::worklet::DispatcherMapField<SegmentedSplitTransform, DeviceAdapter>
SegmentedSplitTransformDispatcher(sstWorklet);
SegmentedSplitTransformDispatcher.Invoke(
B_Handle, D_Handle, F_Handle, G_Handle, H_Handle, I_Handle);
A_Handle = ScatterArrayWrapper(A_Handle, I_Handle, device);
B_Handle = ScatterArrayWrapper(B_Handle, I_Handle, device);
X_Handle = ScatterArrayWrapper(X_Handle, I_Handle, device);
Y_Handle = ScatterArrayWrapper(Y_Handle, I_Handle, device);
Z_Handle = ScatterArrayWrapper(Z_Handle, I_Handle, device);
}
template <typename T, typename DeviceAdapter>
void OneLevelSplit3D(vtkm::cont::ArrayHandle<T>& pointId_Handle,
vtkm::cont::ArrayHandle<T>& xrank_Handle,
vtkm::cont::ArrayHandle<T>& yrank_Handle,
vtkm::cont::ArrayHandle<T>& zrank_Handle,
vtkm::cont::ArrayHandle<T>& segId_Handle,
vtkm::cont::ArrayHandle<T>& splitId_Handle,
vtkm::Int32 level,
DeviceAdapter device)
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
vtkm::cont::ArrayHandle<T> flag_Handle;
if (level % 3 == 0)
{
flag_Handle = ComputeFlagProcedure(xrank_Handle, segId_Handle, device);
}
else if (level % 3 == 1)
{
flag_Handle = ComputeFlagProcedure(yrank_Handle, segId_Handle, device);
}
else
{
flag_Handle = ComputeFlagProcedure(zrank_Handle, segId_Handle, device);
}
SegmentedSplitProcedure3D(
pointId_Handle, flag_Handle, segId_Handle, xrank_Handle, yrank_Handle, zrank_Handle, device);
vtkm::cont::ArrayHandle<T> segIdOld_Handle;
Algorithm::Copy(segId_Handle, segIdOld_Handle);
segId_Handle = NewKeyWrapper(segIdOld_Handle, flag_Handle, device);
RenumberRanksProcedure(xrank_Handle, flag_Handle, segId_Handle, segIdOld_Handle, device);
RenumberRanksProcedure(yrank_Handle, flag_Handle, segId_Handle, segIdOld_Handle, device);
RenumberRanksProcedure(zrank_Handle, flag_Handle, segId_Handle, segIdOld_Handle, device);
if (level % 3 == 0)
{
splitId_Handle =
SaveSplitPointIdWrapper(pointId_Handle, flag_Handle, xrank_Handle, splitId_Handle, device);
}
else if (level % 3 == 1)
{
splitId_Handle =
SaveSplitPointIdWrapper(pointId_Handle, flag_Handle, yrank_Handle, splitId_Handle, device);
}
else
{
splitId_Handle =
SaveSplitPointIdWrapper(pointId_Handle, flag_Handle, zrank_Handle, splitId_Handle, device);
}
}
// Execute the 3d kd tree construction given x y z coordinate vectors
// Returns:
// Leaf Node vector and internal node (split) vectpr
template <typename CoordiType, typename TreeIdType, typename DeviceAdapter>
void Run(vtkm::cont::ArrayHandle<vtkm::Vec<CoordiType, 3>>& coordi_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& pointId_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& splitId_Handle,
DeviceAdapter device)
{
typedef typename vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter> Algorithm;
vtkm::Id nTrainingPoints = coordi_Handle.GetNumberOfValues();
vtkm::cont::ArrayHandleCounting<vtkm::Id> counting_Handle(0, 1, nTrainingPoints);
Algorithm::Copy(counting_Handle, pointId_Handle);
vtkm::cont::ArrayHandle<TreeIdType> xorder_Handle;
Algorithm::Copy(counting_Handle, xorder_Handle);
vtkm::cont::ArrayHandle<TreeIdType> yorder_Handle;
Algorithm::Copy(counting_Handle, yorder_Handle);
vtkm::cont::ArrayHandle<TreeIdType> zorder_Handle;
Algorithm::Copy(counting_Handle, zorder_Handle);
splitId_Handle.Allocate(nTrainingPoints);
vtkm::cont::ArrayHandle<CoordiType> xcoordi_Handle;
vtkm::cont::ArrayHandle<CoordiType> ycoordi_Handle;
vtkm::cont::ArrayHandle<CoordiType> zcoordi_Handle;
SeprateVec3AryHandle sepVec3Worklet;
vtkm::worklet::DispatcherMapField<SeprateVec3AryHandle, DeviceAdapter> sepVec3Dispatcher(
sepVec3Worklet);
sepVec3Dispatcher.Invoke(coordi_Handle, xcoordi_Handle, ycoordi_Handle, zcoordi_Handle);
Algorithm::SortByKey(xcoordi_Handle, xorder_Handle);
vtkm::cont::ArrayHandle<TreeIdType> xrank_Handle =
ScatterArrayWrapper(pointId_Handle, xorder_Handle, device);
Algorithm::SortByKey(ycoordi_Handle, yorder_Handle);
vtkm::cont::ArrayHandle<TreeIdType> yrank_Handle =
ScatterArrayWrapper(pointId_Handle, yorder_Handle, device);
Algorithm::SortByKey(zcoordi_Handle, zorder_Handle);
vtkm::cont::ArrayHandle<TreeIdType> zrank_Handle =
ScatterArrayWrapper(pointId_Handle, zorder_Handle, device);
vtkm::cont::ArrayHandle<TreeIdType> segId_Handle;
vtkm::cont::ArrayHandleConstant<TreeIdType> constHandle(0, nTrainingPoints);
Algorithm::Copy(constHandle, segId_Handle);
///// build kd tree /////
vtkm::Int32 maxLevel = static_cast<vtkm::Int32>(ceil(vtkm::Log2(nTrainingPoints) + 1));
for (vtkm::Int32 i = 0; i < maxLevel - 1; i++)
{
OneLevelSplit3D(pointId_Handle,
xrank_Handle,
yrank_Handle,
zrank_Handle,
segId_Handle,
splitId_Handle,
i,
device);
}
}
};
}
}
} // namespace vtkm::worklet
#endif // vtk_m_worklet_KdTree3DConstruction_h

225
vtkm/worklet/spatialstructure/KdTree3DNNSearch.h Normal file
View File

@ -0,0 +1,225 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_KdTree3DNNSearch_h
#define vtk_m_worklet_KdTree3DNNSearch_h
#include <vtkm/Math.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCounting.h>
#include <vtkm/cont/ArrayHandleReverse.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/arg/ControlSignatureTagBase.h>
#include <vtkm/cont/serial/DeviceAdapterSerial.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/internal/DispatcherBase.h>
#include <vtkm/worklet/internal/WorkletBase.h>
namespace vtkm
{
namespace worklet
{
namespace spatialstructure
{
class KdTree3DNNSearch
{
public:
class NearestNeighborSearch3DWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> qcIn,
WholeArrayIn<> treeIdIn,
WholeArrayIn<> treeSplitIdIn,
WholeArrayIn<> treeCoordiIn,
FieldOut<> nnIdOut,
FieldOut<> nnDisOut);
typedef void ExecutionSignature(_1, _2, _3, _4, _5, _6);
VTKM_CONT
NearestNeighborSearch3DWorklet() {}
template <typename CooriVecT, typename CooriT, typename IdPortalT, typename CoordiPortalT>
VTKM_EXEC_CONT void NearestNeighborSearch3D(CooriVecT qc,
CooriT& dis,
vtkm::Id& nnpIdx,
vtkm::Int32 level,
vtkm::Id sIdx,
vtkm::Id tIdx,
IdPortalT treePortal,
IdPortalT splitIdPortal,
CoordiPortalT coordiPortal) const
{
CooriT qx = qc[0];
CooriT qy = qc[1];
CooriT qz = qc[2];
if (tIdx - sIdx == 1)
{ ///// leaf node
vtkm::Id leafNodeIdx = treePortal.Get(sIdx);
CooriT leafX = coordiPortal.Get(leafNodeIdx)[0];
CooriT leafY = coordiPortal.Get(leafNodeIdx)[1];
CooriT leafZ = coordiPortal.Get(leafNodeIdx)[2];
CooriT _dis = vtkm::Sqrt((leafX - qx) * (leafX - qx) + (leafY - qy) * (leafY - qy) +
(leafZ - qz) * (leafZ - qz));
if (_dis < dis)
{
dis = _dis;
nnpIdx = leafNodeIdx;
}
}
else
{ //normal Node
vtkm::Id splitNodeLoc = static_cast<vtkm::Id>(vtkm::Ceil((sIdx + tIdx) / 2.0));
CooriT splitX = coordiPortal.Get(splitIdPortal.Get(splitNodeLoc))[0];
CooriT splitY = coordiPortal.Get(splitIdPortal.Get(splitNodeLoc))[1];
CooriT splitZ = coordiPortal.Get(splitIdPortal.Get(splitNodeLoc))[2];
CooriT splitAxis;
CooriT queryCoordi;
if (level % 3 == 0)
{ //x axis level
splitAxis = splitX;
queryCoordi = qx;
}
else if (level % 3 == 1)
{
splitAxis = splitY;
queryCoordi = qy;
}
else
{
splitAxis = splitZ;
queryCoordi = qz;
}
if (queryCoordi <= splitAxis)
{ //left tree first
if (queryCoordi - dis <= splitAxis)
NearestNeighborSearch3D(qc,
dis,
nnpIdx,
level + 1,
sIdx,
static_cast<vtkm::Id>(vtkm::Ceil((sIdx + tIdx) / 2.0)),
treePortal,
splitIdPortal,
coordiPortal);
if (queryCoordi + dis > splitAxis)
NearestNeighborSearch3D(qc,
dis,
nnpIdx,
level + 1,
static_cast<vtkm::Id>(vtkm::Ceil((sIdx + tIdx) / 2.0)),
tIdx,
treePortal,
splitIdPortal,
coordiPortal);
}
else
{ //right tree first
if (queryCoordi + dis > splitAxis)
NearestNeighborSearch3D(qc,
dis,
nnpIdx,
level + 1,
static_cast<vtkm::Id>(vtkm::Ceil((sIdx + tIdx) / 2.0)),
tIdx,
treePortal,
splitIdPortal,
coordiPortal);
if (queryCoordi - dis <= splitAxis)
NearestNeighborSearch3D(qc,
dis,
nnpIdx,
level + 1,
sIdx,
static_cast<vtkm::Id>(vtkm::Ceil((sIdx + tIdx) / 2.0)),
treePortal,
splitIdPortal,
coordiPortal);
}
}
}
template <typename CoordiVecType,
typename IdPortalType,
typename CoordiPortalType,
typename IdType,
typename CoordiType>
VTKM_EXEC void operator()(const CoordiVecType& qc,
const IdPortalType& treeIdPortal,
const IdPortalType& treeSplitIdPortal,
const CoordiPortalType& treeCoordiPortal,
IdType& nnId,
CoordiType& nnDis) const
{
nnDis = std::numeric_limits<CoordiType>::max();
NearestNeighborSearch3D(qc,
nnDis,
nnId,
0,
0,
treeIdPortal.GetNumberOfValues(),
treeIdPortal,
treeSplitIdPortal,
treeCoordiPortal);
}
};
// Execute the Neaseat Neighbor Search given kdtree and search points
// Returns:
// Vectors of NN point index and NNpoint distance
template <typename CoordiType, typename TreeIdType, typename DeviceAdapter>
void Run(vtkm::cont::ArrayHandle<vtkm::Vec<CoordiType, 3>>& coordi_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& pointId_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& splitId_Handle,
vtkm::cont::ArrayHandle<vtkm::Vec<CoordiType, 3>>& qc_Handle,
vtkm::cont::ArrayHandle<TreeIdType>& nnId_Handle,
vtkm::cont::ArrayHandle<CoordiType>& nnDis_Handle,
DeviceAdapter vtkmNotUsed(device))
{
#if VTKM_DEVICE_ADAPTER == VTKM_DEVICE_ADAPTER_CUDA
//set up stack size for cuda envinroment
size_t stackSizeBackup;
cudaDeviceGetLimit(&stackSizeBackup, cudaLimitStackSize);
cudaDeviceSetLimit(cudaLimitStackSize, 1024 * 16);
#endif
NearestNeighborSearch3DWorklet nns3dWorklet;
vtkm::worklet::DispatcherMapField<NearestNeighborSearch3DWorklet, DeviceAdapter>
nns3DDispatcher(nns3dWorklet);
nns3DDispatcher.Invoke(
qc_Handle, pointId_Handle, splitId_Handle, coordi_Handle, nnId_Handle, nnDis_Handle);
#if VTKM_DEVICE_ADAPTER == VTKM_DEVICE_ADAPTER_CUDA
cudaDeviceSetLimit(cudaLimitStackSize, stackSizeBackup);
#endif
}
};
}
}
} // namespace vtkm::worklet
#endif // vtk_m_worklet_KdTree3DNNSearch_h

1
vtkm/worklet/testing/CMakeLists.txt
View File

@ -30,6 +30,7 @@ set(unit_tests
UnitTestExtractStructured.cxx
UnitTestFieldHistogram.cxx
UnitTestFieldStatistics.cxx
UnitTestKdTreeBuildNNS.cxx
UnitTestKeys.cxx
UnitTestMagnitude.cxx
UnitTestMarchingCubes.cxx

171
vtkm/worklet/testing/UnitTestKdTreeBuildNNS.cxx Normal file
View File

@ -0,0 +1,171 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 Sandia Corporation.
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#include <random>
#include <vtkm/worklet/KdTree3D.h>
namespace
{
typedef vtkm::cont::DeviceAdapterAlgorithm<VTKM_DEFAULT_DEVICE_ADAPTER_TAG> Algorithm;
////brute force method /////
template <typename CoordiVecT, typename CoordiPortalT, typename CoordiT>
VTKM_EXEC_CONT vtkm::Id NNSVerify3D(CoordiVecT qc, CoordiPortalT coordiPortal, CoordiT& dis)
{
dis = std::numeric_limits<CoordiT>::max();
vtkm::Id nnpIdx;
for (vtkm::Int32 i = 0; i < coordiPortal.GetNumberOfValues(); i++)
{
CoordiT splitX = coordiPortal.Get(i)[0];
CoordiT splitY = coordiPortal.Get(i)[1];
CoordiT splitZ = coordiPortal.Get(i)[2];
CoordiT _dis =
vtkm::Sqrt((splitX - qc[0]) * (splitX - qc[0]) + (splitY - qc[1]) * (splitY - qc[1]) +
(splitZ - qc[2]) * (splitZ - qc[2]));
if (_dis < dis)
{
dis = _dis;
nnpIdx = i;
}
}
return nnpIdx;
}
class NearestNeighborSearchBruteForce3DWorklet : public vtkm::worklet::WorkletMapField
{
public:
typedef void ControlSignature(FieldIn<> qcIn,
WholeArrayIn<> treeCoordiIn,
FieldOut<> nnIdOut,
FieldOut<> nnDisOut);
typedef void ExecutionSignature(_1, _2, _3, _4);
VTKM_CONT
NearestNeighborSearchBruteForce3DWorklet() {}
template <typename CoordiVecType, typename CoordiPortalType, typename IdType, typename CoordiType>
VTKM_EXEC void operator()(const CoordiVecType& qc,
const CoordiPortalType& coordiPortal,
IdType& nnId,
CoordiType& nnDis) const
{
nnDis = std::numeric_limits<CoordiType>::max();
nnId = NNSVerify3D(qc, coordiPortal, nnDis);
}
};
void TestKdTreeBuildNNS()
{
vtkm::Int32 nTrainingPoints = 1000;
vtkm::Int32 nTestingPoint = 1000;
std::vector<vtkm::Float32> xcoordi;
std::vector<vtkm::Float32> ycoordi;
std::vector<vtkm::Float32> zcoordi;
std::vector<vtkm::Vec<vtkm::Float32, 3>> coordi;
///// randomly genarate training points/////
std::default_random_engine dre;
std::uniform_real_distribution<vtkm::Float32> dr(0.0f, 10.0f);
for (vtkm::Int32 i = 0; i < nTrainingPoints; i++)
{
xcoordi.push_back(dr(dre));
ycoordi.push_back(dr(dre));
zcoordi.push_back(dr(dre));
vtkm::Vec<vtkm::Float32, 3> c;
c[0] = xcoordi[xcoordi.size() - 1];
c[1] = ycoordi[ycoordi.size() - 1];
c[2] = zcoordi[zcoordi.size() - 1];
coordi.push_back(c);
}
///// preprare data to build 3D kd tree /////
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> coordi_Handle;
Algorithm::Copy(vtkm::cont::make_ArrayHandle(coordi), coordi_Handle);
vtkm::cont::ArrayHandle<vtkm::Id> pointId_Handle;
vtkm::cont::ArrayHandle<vtkm::Id> splitId_Handle;
// Run data
vtkm::worklet::KdTree3D kdtree3d;
kdtree3d.Run(coordi_Handle, pointId_Handle, splitId_Handle, VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
//Nearest Neighbor worklet Testing
/// randomly generate testing points /////
std::vector<vtkm::Vec<vtkm::Float32, 3>> qcVec;
for (vtkm::Int32 i = 0; i < nTestingPoint; i++)
{
vtkm::Vec<vtkm::Float32, 3> qc;
qc[0] = dr(dre);
qc[1] = dr(dre);
qc[2] = dr(dre);
qcVec.push_back(qc);
}
///// preprare testing data /////
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32, 3>> qc_Handle;
Algorithm::Copy(vtkm::cont::make_ArrayHandle(qcVec), qc_Handle);
vtkm::cont::ArrayHandle<vtkm::Id> nnId_Handle;
vtkm::cont::ArrayHandle<vtkm::Float32> nnDis_Handle;
kdtree3d.Run(coordi_Handle,
pointId_Handle,
splitId_Handle,
qc_Handle,
nnId_Handle,
nnDis_Handle,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
vtkm::cont::ArrayHandle<vtkm::Id> bfnnId_Handle;
vtkm::cont::ArrayHandle<vtkm::Float32> bfnnDis_Handle;
NearestNeighborSearchBruteForce3DWorklet nnsbf3dWorklet;
vtkm::worklet::DispatcherMapField<NearestNeighborSearchBruteForce3DWorklet> nnsbf3DDispatcher(
nnsbf3dWorklet);
nnsbf3DDispatcher.Invoke(
qc_Handle, vtkm::cont::make_ArrayHandle(coordi), bfnnId_Handle, bfnnDis_Handle);
///// verfity search result /////
bool passTest = true;
for (vtkm::Int32 i = 0; i < nTestingPoint; i++)
{
vtkm::Id workletIdx = nnId_Handle.GetPortalControl().Get(i);
vtkm::Id bfworkletIdx = bfnnId_Handle.GetPortalControl().Get(i);
if (workletIdx != bfworkletIdx)
{
passTest = false;
}
}
VTKM_TEST_ASSERT(passTest, "Kd tree NN search result incorrect.");
}
} // anonymous namespace
int UnitTestKdTreeBuildNNS(int, char* [])
{
return vtkm::cont::testing::Testing::Run(TestKdTreeBuildNNS);
}