diff --git a/vtkm/filter/CMakeLists.txt b/vtkm/filter/CMakeLists.txt index 9abeb8390..75c6967d8 100644 --- a/vtkm/filter/CMakeLists.txt +++ b/vtkm/filter/CMakeLists.txt @@ -40,6 +40,7 @@ set(headers Histogram.h ImageConnectivity.h Lagrangian.h + LagrangianStructures.h MarchingCubes.h Mask.h MaskPoints.h @@ -103,6 +104,7 @@ set(header_template_sources Histogram.hxx ImageConnectivity.hxx Lagrangian.hxx + LagrangianStructures.hxx MarchingCubes.hxx Mask.hxx MaskPoints.hxx diff --git a/vtkm/filter/LagrangianStructures.h b/vtkm/filter/LagrangianStructures.h index 893b0c7e7..115122903 100644 --- a/vtkm/filter/LagrangianStructures.h +++ b/vtkm/filter/LagrangianStructures.h @@ -27,39 +27,29 @@ public: using Scalar = vtkm::worklet::particleadvection::ScalarType; using Vector = vtkm::Vec; - VTKM_CONT LagrangianStructures(); - VTKM_CONT void SetStepSize(ScalarType s) { this->StepSize = s; } - VTKM_CONT - ScalarType GetStepSize() { return this->StepSize; } + Scalar GetStepSize() { return this->StepSize; } - VTKM_CONT void SetNumberOfSteps(vtkm::Id n) { this->NumberOfSteps = n; } - VTKM_CONT vtkm::Id GetNumberOfSteps() { return this->NumberOfSteps; } - VTKM_CONT + void SetAdvectionTime(Scalar advectionTime) { this->AdvectionTime = advectionTime; } + Scalar GetAdvectionTime() { return this->AdvectionTime; } + void SetUseAuxiliaryGrid(bool useAuxiliaryGrid){ this->UseAuxiliaryGrid = useAuxiliaryGrid }; - VTKM_CONT bool GetUseAuxiliaryGrid() { return this->UseAuxiliaryGrid; } - VTKM_CONT void SetAuxiliaryGridDimensions(vtkm::Id3 auxiliaryDims){ this->AuxiliaryDims = auxiliaryDims }; - VTKM_CONT vtkm::Id3 GetAuxiliaryGridDimensions(){ return this->AuxiliaryDims }; - VTKM_CONT void SetUseFlowMapOutput(bool useFLowMapOutput){ this->UseFlowMapOutput = useFlowMapOutput }; - VTKM_CONT bool GetUseFlowMapOutput() { return this->UseFlowMapOutput; } - VTKM_CONT inline void SetFlowMapOutput(vtkm::cont::ArrayHandle& flowMap){ this->FlowMapOutput = flowMap }; - VTKM_CONT inline vtkm::cont::ArrayHandle GetFlowMapOutput(){ return this->FlowMapOutput }; template @@ -81,6 +71,7 @@ private: vtkm::worklet::LagrangianStructures Worklet; Scalar StepSize; vtkm::Id NumberOfSteps; + Scalar AdvectionTime; bool UseAuxiliaryGrid = false; vtkm::Id3 AuxiliaryDims; bool UseFlowMapOutput; diff --git a/vtkm/filter/LagrangianStructures.hxx b/vtkm/filter/LagrangianStructures.hxx index 7a079c7c4..c19ef06dc 100644 --- a/vtkm/filter/LagrangianStructures.hxx +++ b/vtkm/filter/LagrangianStructures.hxx @@ -15,6 +15,8 @@ #include #include +#include + namespace vtkm { namespace filter @@ -23,17 +25,9 @@ namespace filter //----------------------------------------------------------------------------- inline VTKM_CONT LagrangianStructures::LagrangianStructures() : vtkm::filter::FilterDataSetWithField() - , Worklet() { } -//----------------------------------------------------------------------------- -inline VTKM_CONT void LagrangianStructures::SetSeeds( - vtkm::cont::ArrayHandle>& seeds) -{ - this->Seeds = seeds; -} - //----------------------------------------------------------------------------- template inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute( @@ -55,53 +49,84 @@ inline VTKM_CONT vtkm::cont::DataSet LagrangianStructures::DoExecute( using GridEvaluator = vtkm::worklet::particleadvection::GridEvaluator; using Integrator = vtkm::worklet::particleadvection::RK4Integrator; - Scalar stepLength = static_cast(0.025); - vtkm::Id numberOfSteps = 500; + Scalar stepLength = this->GetStepLength(); + vtkm::Id numberOfSteps = this->GetNumberOfStep(); vtkm::cont::CoordinateSystem coordinates = input.GetCoordinateSystem(); vtkm::cont::DynamicCellSet cellset = input.GetCellSet(); - GridEvaluator evaluator(input.GetCoordinateSystem(), input.GetCellSet(), vectors); - Integrator integrator(evaluator, stepLength); - - vtkm::worklet::ParticleAdvection particles; - vtkm::worklet::ParticleAdvectionResult res; - - std::cout << "Advecting particles" << std::endl; - res = particles.Run(integrator, points, numberOfSteps); - std::cout << "Advected particles" << std::endl; - - vtkm::cont::ArrayHandle outputPoints; - vtkm::cont::ArrayCopy(res.positions, outputPoints); - vtkm::cont::ArrayCopy(input.GetCoordinateSystem().GetData(), points); - + vtkm::cont::Dataset lcsInput; + if (this->GetUseAuxiliaryGrid()) + { + vtkm::Id3 lcsGridDims = this->GetAuxiliaryGridDimensions(); + vtkm::Bounds inputBounds = this->input.GetBounds(); + vtkm::Vec origin(inputBounds.X.Min, inputBounds.Y.Min, inputBounds.Z.Min); + vtkm::Vec spacing; + spacing[0] = inputBounds.X.Length() / static_cast(lcsGridDims[0] - 1); + spacing[1] = inputBounds.Y.Length() / static_cast(lcsGridDims[1] - 1); + spacing[2] = inputBounds.Z.Length() / static_cast(lcsGridDims[2] - 1); + vtkm::cont::DataSetBuilderUniform uniformDatasetBuilder; + lcsInput = uniformDatasetBuilder.Create(lcsGridDims, origin, spacing); + } + else + { + // Check if input dataset is structured. + // If not, we cannot proceed. + if (!(cellset.IsType() || cellset.IsType())) + throw vtkm::cont::ErrorFilterExecution("Provided data is not structured, + provide parameters for an auxiliary grid."); + lcsInput = input; + } + vtkm::cont::ArrayHandle lcsInputPoints, lcsOutputPoints; + vtkm::cont::ArrayCopy(lcsInput.GetCoordinateSystem().GetData(), lcsInputPoints); + if (this->GetUseFlowMapOutput()) + { + // Check if there is a 1:1 correspondense between the flow map + // and the input points + lcsOutputPoints = this->GetFlowMapOutput(); + if (lcsInputPoints.GetNumberOfValues() != lcsOutputPoints.GetNumberOfValues()) + throw vtkm::cont::ErrorFilterExecution("Provided flow map does not correspond + to the input points for LCS filter."); + } + else + { + std::cout << "Advecting particles" << std::endl; + GridEvaluator evaluator(input.GetCoordinateSystem(), input.GetCellSet(), field); + Integrator integrator(evaluator, stepLength); + vtkm::worklet::ParticleAdvection particles; + vtkm::worklet::ParticleAdvectionResult advectionResult; + vtkm::cont::ArrayHandle advectionPoints; + vtkm::cont::ArrayCopy(lcsInputPoints, advectionPoints); + advectionResult = particles.Run(integrator, advectionPoints, numberOfSteps); + std::cout << "Advected particles" << std::endl; + lcsOutputPoints = advectionResult.positions; + } // FTLE output field vtkm::cont::ArrayHandle outputField; - //std::cout << "Calculating FTLE field" << std::endl; - vtkm::FloatDefault advectionTime = stepLength * numberOfSteps; - std::cout << "Calculating FTLE" << std::endl; - detail::CalculateLCSField(cellset, points, outputPoints, advectionTime, outputField); - std::cout << "Calculated FTLE" << std::endl; + std::cout << "Calculating FTLE field" << std::endl; + vtkm::FloatDefault advectionTime = this->GetAdvectionTime(); - if (cellSet.IsType()) + vtkm::cont::DynamicCellSet lcsCellSet = lcsInput.GetCellSet(); + if (lcsCellSet.IsType()) { - using AnalysisType = detail::FTLECalc<2>; - AnalysisType ftleCalculator(advectionTime, cellSet); + using AnalysisType = vtkm::worklet::LagrangianStructures<2>; + AnalysisType ftleCalculator(advectionTime, lcsCellSet); vtkm::worklet::DispatcherMapField dispatcher(ftleCalculator); - dispatcher.Invoke(inputPoints, outputPoints, outputField); + dispatcher.Invoke(lcsInputPoints, lcsOutputPoints, outputField); } - else if (cellSet.IsType()) + else if (lcsCellSet.IsType()) { - using AnalysisType = detail::FTLECalc<3>; - AnalysisType ftleCalculator(advectionTime, cellSet); + using AnalysisType = vtkm::worklet::LagrangianStructures<3>; + AnalysisType ftleCalculator(advectionTime, lcsCellSet); vtkm::worklet::DispatcherMapField dispatcher(ftleCalculator); - dispatcher.Invoke(inputPoints, outputPoints, outputField); + dispatcher.Invoke(lcsInputPoints, lcsOutputPoints, outputField); } + std::cout << "Calculated FTLE field" << std::endl; vtkm::cont::DataSet output; vtkm::cont::DataSetFieldAdd fieldAdder; - output.AddCoordinateSystem(input.GetCoordinateSystem()); - output.AddCellSet(input.GetCellSet()); + output.AddCoordinateSystem(lcsInput.GetCoordinateSystem()); + output.AddCellSet(lcsInput.GetCellSet()); fieldAdder.AddPointField(output, "FTLE", outputField); return output; } diff --git a/vtkm/worklet/CMakeLists.txt b/vtkm/worklet/CMakeLists.txt index 500252a03..8ad81fa99 100644 --- a/vtkm/worklet/CMakeLists.txt +++ b/vtkm/worklet/CMakeLists.txt @@ -34,6 +34,7 @@ set(headers FieldStatistics.h Gradient.h Invoker.h + LagrangianStructures.h KdTree3D.h KernelSplatter.h Keys.h @@ -121,6 +122,7 @@ add_subdirectory(contourtree_augmented) add_subdirectory(cosmotools) add_subdirectory(gradient) add_subdirectory(histogram) +add_subdirectory(lcs) add_subdirectory(splatkernels) add_subdirectory(spatialstructure) add_subdirectory(tetrahedralize) diff --git a/vtkm/worklet/lcs/CMakeLists.txt b/vtkm/worklet/lcs/CMakeLists.txt new file mode 100644 index 000000000..82aad13ae --- /dev/null +++ b/vtkm/worklet/lcs/CMakeLists.txt @@ -0,0 +1,15 @@ +##============================================================================ +## 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. +##============================================================================ + +set(headers + LagrangianStrucutreHelpers.h + ) + +vtkm_declare_headers(${headers})