text: //============================================================================
  //  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 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
  //  Copyright 2014 UT-Battelle, LLC.
  //  Copyright 2014 Los Alamos National Security.
  //
  //  Under the terms of Contract DE-NA0003525 with NTESS,
  //  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      /////////////////////
    /// \brief Segmented split for 3D x, y, z coordinates
    ///
    /// Split \c pointId_Handle, \c X_Handle, \c Y_Handle and \c Z_Handle within each segment
    /// as indicated by \c segId_Handle according to flags in \c flag_Handle.
    ///
    /// \tparam T
    /// \tparam DeviceAdapter
    /// \param pointId_Handle
    /// \param flag_Handle
    /// \param segId_Handle
    /// \param X_Handle
    /// \param Y_Handle
    /// \param Z_Handle
    /// \param device
    template <typename T, typename DeviceAdapter>
    void SegmentedSplitProcedure3D(vtkm::cont::ArrayHandle<T>& pointId_Handle,
                                   vtkm::cont::ArrayHandle<T>& flag_Handle,
                                   vtkm::cont::ArrayHandle<T>& segId_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(segId_Handle, flag_Handle, D_Handle, initValue, vtkm::Add());
  
      vtkm::cont::ArrayHandleCounting<T> Ecouting_Handle(0, 1, pointId_Handle.GetNumberOfValues());
      vtkm::cont::ArrayHandle<T> E_Handle;
      Algorithm::Copy(Ecouting_Handle, E_Handle);
  
      vtkm::cont::ArrayHandle<T> F_Handle;
      Algorithm::ScanInclusiveByKey(segId_Handle, E_Handle, F_Handle, vtkm::Minimum());
  
      vtkm::cont::ArrayHandle<T> InvB_Handle = Inverse01ArrayWrapper(flag_Handle, device);
      vtkm::cont::ArrayHandle<T> G_Handle;
      Algorithm::ScanInclusiveByKey(segId_Handle, InvB_Handle, G_Handle, vtkm::Add());
  
      vtkm::cont::ArrayHandle<T> H_Handle =
        ReverseScanInclusiveByKey(segId_Handle, G_Handle, vtkm::Maximum(), device);
  
      vtkm::cont::ArrayHandle<T> I_Handle;
      SegmentedSplitTransform sstWorklet;
      vtkm::worklet::DispatcherMapField<SegmentedSplitTransform, DeviceAdapter>
        SegmentedSplitTransformDispatcher(sstWorklet);
      SegmentedSplitTransformDispatcher.Invoke(
        flag_Handle, D_Handle, F_Handle, G_Handle, H_Handle, I_Handle);
  
      pointId_Handle = ScatterArrayWrapper(pointId_Handle, I_Handle, device);
  
      flag_Handle = ScatterArrayWrapper(flag_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);
    }
  
    /// \brief Perform one level of KD-Tree construction
    ///
    /// Construct a level of KD-Tree by segemeted splits (partitioning) of \c pointId_Handle,
    /// \c xrank_Handle, \c yrank_Handle and \c zrank_Handle according to the medium element
    /// in each segment as indicated by \c segId_Handle alone the axis determined by \c level.
    /// The split point of each segment will be updated in \c splitId_Handle.
    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);
      }
    }
  
    /// \brief Construct KdTree from x y z coordinate vector.
    ///
    /// This method constructs an array based KD-Tree from x, y, z coordinates of points in \c
    /// coordi_Handle. The method rotates between x, y and z axis and splits input points into
    /// equal halves with respect to the split axis at each level of construction. The indices to
    /// the leaf nodes are returned in \c pointId_Handle and indices to internal nodes (splits)
    /// are returned in splitId_handle.
    ///
    /// \param coordi_Handle (in) x, y, z coordinates of input points
    /// \param pointId_Handle (out) returns indices to leaf nodes of the KD-tree
    /// \param splitId_Handle (out) returns indices to internal nodes of the KD-tree
    /// \param device the device to run the construction on
    // 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
  