613 lines
23 KiB
C++
613 lines
23 KiB
C++
//============================================================================
|
|
// 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 CoordType, typename CoordStorageTag, typename DeviceAdapter>
|
|
void Run(const vtkm::cont::ArrayHandle<vtkm::Vec<CoordType, 3>, CoordStorageTag>& coordi_Handle,
|
|
vtkm::cont::ArrayHandle<vtkm::Id>& pointId_Handle,
|
|
vtkm::cont::ArrayHandle<vtkm::Id>& 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<vtkm::Id> xorder_Handle;
|
|
Algorithm::Copy(counting_Handle, xorder_Handle);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> yorder_Handle;
|
|
Algorithm::Copy(counting_Handle, yorder_Handle);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> zorder_Handle;
|
|
Algorithm::Copy(counting_Handle, zorder_Handle);
|
|
|
|
splitId_Handle.Allocate(nTrainingPoints);
|
|
|
|
vtkm::cont::ArrayHandle<CoordType> xcoordi_Handle;
|
|
vtkm::cont::ArrayHandle<CoordType> ycoordi_Handle;
|
|
vtkm::cont::ArrayHandle<CoordType> 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<vtkm::Id> xrank_Handle =
|
|
ScatterArrayWrapper(pointId_Handle, xorder_Handle, device);
|
|
|
|
Algorithm::SortByKey(ycoordi_Handle, yorder_Handle);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> yrank_Handle =
|
|
ScatterArrayWrapper(pointId_Handle, yorder_Handle, device);
|
|
|
|
Algorithm::SortByKey(zcoordi_Handle, zorder_Handle);
|
|
vtkm::cont::ArrayHandle<vtkm::Id> zrank_Handle =
|
|
ScatterArrayWrapper(pointId_Handle, zorder_Handle, device);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::Id> segId_Handle;
|
|
vtkm::cont::ArrayHandleConstant<vtkm::Id> 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
|