Merge branch 'ResliceAndPGSTest' into 'master'

ImageReslice, PGS Solver Tests

See merge request iMSTK/iMSTK!557

Examples/VolumeRendering/VolumeRenderingExample.cpp

@@ -21,6 +21,7 @@
 
 #include "imstkCamera.h"
 #include "imstkImageData.h"
+#include "imstkImageReslice.h"
 #include "imstkKeyboardSceneControl.h"
 #include "imstkLogger.h"
 #include "imstkMeshIO.h"
@@ -50,10 +51,19 @@ main()
     // SDK and Scene
     imstkNew<Scene> scene("VolumeRendering");
 
+    // Read an image
+    auto imageData = MeshIO::read<ImageData>(iMSTK_DATA_ROOT "skullVolume.nrrd");
+
+    // Rotate that image
+    imstkNew<ImageReslice> reslice;
+    reslice->setInputImage(imageData);
+    // Rotate 1radian around y
+    reslice->setTransform(mat4dRotation(Rotd(1.0, Vec3d(0.0, 1.0, 0.0))));
+    reslice->update();
+
     // Create a visual object in the scene for the volume
     imstkNew<SceneObject> volumeObj("VisualVolume");
-    auto                  imageData = MeshIO::read<ImageData>(iMSTK_DATA_ROOT "skullVolume.nrrd");
-    volumeObj->setVisualGeometry(imageData);
+    volumeObj->setVisualGeometry(reslice->getOutputImage());
     scene->addSceneObject(volumeObj);
 
     // Update Camera to position volume close to viewer

Source/Filtering/imstkImageReslice.cpp

+/*=========================================================================
+
+   Library: iMSTK
+
+   Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+   & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+          http://www.apache.org/licenses/LICENSE-2.0.txt
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+
+=========================================================================*/
+
+#include "imstkImageReslice.h"
+#include "imstkGeometryUtilities.h"
+#include "imstkImageData.h"
+#include "imstkLogger.h"
+
+#include <vtkImageData.h>
+#include <vtkImageReslice.h>
+#include <vtkTransform.h>
+
+namespace imstk
+{
+ImageReslice::ImageReslice()
+{
+    setNumberOfInputPorts(1);
+    setNumberOfOutputPorts(1);
+    setOutput(std::make_shared<ImageData>());
+}
+
+std::shared_ptr<ImageData>
+ImageReslice::getOutputImage() const
+{
+    return std::dynamic_pointer_cast<ImageData>(getOutput(0));
+}
+
+void
+ImageReslice::setInputImage(std::shared_ptr<ImageData> inputData)
+{
+    setInput(inputData, 0);
+}
+
+void
+ImageReslice::requestUpdate()
+{
+    std::shared_ptr<ImageData> inputImage = std::dynamic_pointer_cast<ImageData>(getInput(0));
+    if (inputImage == nullptr)
+    {
+        LOG(WARNING) << "No inputImage to resample";
+        return;
+    }
+
+    vtkNew<vtkTransform> transform;
+    Mat4d                test = m_Transform.transpose();
+    transform->SetMatrix(test.data());
+
+    vtkNew<vtkImageReslice> reslice;
+    reslice->SetResliceTransform(transform);
+    if (m_InterpolationType == InterpolateType::NearestNeighbor)
+    {
+        reslice->SetInterpolationModeToNearestNeighbor();
+    }
+    else if (m_InterpolationType == InterpolateType::Linear)
+    {
+        reslice->SetInterpolationModeToLinear();
+    }
+    else if (m_InterpolationType == InterpolateType::Cubic)
+    {
+        reslice->SetInterpolationModeToCubic();
+    }
+    vtkSmartPointer<vtkImageData> vtkInputImage = GeometryUtils::copyToVtkImageData(inputImage);
+    reslice->SetInputData(vtkInputImage);
+    reslice->SetResliceTransform(transform);
+    reslice->SetAutoCropOutput(true);
+    reslice->Update();
+
+    setOutput(GeometryUtils::copyToImageData(reslice->GetOutput()));
+}
+}
\ No newline at end of file

Source/Filtering/imstkImageReslice.h

+/*=========================================================================
+
+   Library: iMSTK
+
+   Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+   & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+          http://www.apache.org/licenses/LICENSE-2.0.txt
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+
+=========================================================================*/
+#pragma once
+
+#include "imstkGeometryAlgorithm.h"
+#include "imstkMath.h"
+
+namespace imstk
+{
+class ImageData;
+
+///
+/// \class ImageReslice
+///
+/// \brief Resamples an image using a transform
+///
+class ImageReslice : public GeometryAlgorithm
+{
+public:
+    enum class InterpolateType
+    {
+        Linear,
+        Cubic,
+        NearestNeighbor
+    };
+
+public:
+    ImageReslice();
+
+    virtual ~ImageReslice() override = default;
+
+public:
+    std::shared_ptr<ImageData> getOutputImage() const;
+
+    void setInputImage(std::shared_ptr<ImageData> inputData);
+
+    imstkGetMacro(Transform, const Mat4d&);
+    imstkGetMacro(InterpolationType, const InterpolateType&);
+
+    ///
+    /// \brief Set the transformation matrix
+    ///
+    imstkSetMacro(Transform, const Mat4d&);
+
+    ///
+    /// \brief Set the interpolation type to use when resampling
+    ///
+    imstkSetMacro(InterpolationType, const InterpolateType&);
+
+protected:
+    void requestUpdate() override;
+
+private:
+    Mat4d m_Transform;
+    InterpolateType m_InterpolationType = InterpolateType::Linear;
+};
+}
\ No newline at end of file

Source/Solvers/CMakeLists.txt

@@ -9,6 +9,7 @@ imstk_add_library( Solvers
 #-----------------------------------------------------------------------------
 # Testing
 #-----------------------------------------------------------------------------
-#if( ${PROJECT_NAME}_BUILD_TESTING )
-#  add_subdirectory( Testing )
-#endif()
+if( ${PROJECT_NAME}_BUILD_TESTING )
+  include(imstkAddTest)
+  imstk_add_test(Solvers)
+endif()

Source/Solvers/Testing/imstkPGSSolverTest.cpp

+/*=========================================================================
+
+   Library: iMSTK
+
+   Copyright (c) Kitware, Inc. & Center for Modeling, Simulation,
+   & Imaging in Medicine, Rensselaer Polytechnic Institute.
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+      http://www.apache.org/licenses/LICENSE-2.0.txt
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+
+=========================================================================*/
+
+#include "imstkProjectedGaussSeidelSolver.h"
+#include "imstkTypes.h"
+
+#include <gtest/gtest.h>
+
+using namespace imstk;
+
+class imstkPGSSolverTest : public ::testing::Test
+{
+protected:
+    ProjectedGaussSeidelSolver<double> m_solver;
+};
+
+///
+/// \brief TODO
+///
+TEST_F(imstkPGSSolverTest, Solve5x5)
+{
+    // Testing Ax=b
+    Eigen::MatrixXd Ad(5, 5);
+    Ad <<
+        1.0, 0.999861, 0.997739, 0.971125, 0.984529,
+        0.999861, 1.0, 0.996607, 0.967639, 0.981667,
+        0.997739, 0.996607, 1.0, 0.984906, 0.994004,
+        0.971125, 0.967639, 0.984906, 1.0, 0.997076,
+        0.984529, 0.981667, 0.994004, 0.997076, 1.0;
+    Eigen::SparseMatrix<double> A = Ad.sparseView();
+
+    Eigen::VectorXd b(5);
+    b(0) = 369.425;
+    b(1) = 370.798;
+    b(2) = 382.972;
+    b(3) = 404.772;
+    b(4) = 393.974;
+    // Not testing projection here (cu clamps the solution)
+    Eigen::MatrixXd cu(5, 2);
+    for (int i = 0; i < 5; i++)
+    {
+        cu(i, 0) = IMSTK_DOUBLE_MIN;
+        cu(i, 1) = IMSTK_DOUBLE_MAX;
+    }
+
+    m_solver.setA(&A);
+    m_solver.setMaxIterations(1000);
+    m_solver.setRelaxation(0.05);
+    m_solver.setEpsilon(1.0e-8);
+
+    Eigen::VectorXd x = m_solver.solve(b, cu);
+
+    // Check that Ax now equals b
+    // Test this way in case multiple solutions exist, here we are only
+    // testing that a solution was found
+    Eigen::VectorXd bPrime = A.toDense() * x;
+    //std::cout << "x: " << x << std::endl << std::endl;
+    //std::cout << "Energy " << m_solver.getEnergy() << std::endl;
+
+    std::cout << "actual b: " << b << std::endl << std::endl;
+    std::cout << "computed b: " << bPrime << std::endl << std::endl;
+    for (int i = 0; i < b.size(); i++)
+    {
+        EXPECT_NEAR(bPrime(i), b(i), 10.0);
+    }
+}
+
+///
+/// \brief TODO
+///
+int
+imstkPGSSolverTest(int argc, char* argv[])
+{
+    // Init Google Test & Mock
+    ::testing::InitGoogleTest(&argc, argv);
+
+    // Run tests with gtest
+    return RUN_ALL_TESTS();
+}

Source/Solvers/imstkProjectedGaussSeidelSolver.h

@@ -24,6 +24,15 @@
 
 namespace imstk
 {
+///
+/// \class ProjectedGaussSeidelSolver
+///
+/// \brief Solves a linear system using the projected gauss seidel method.
+/// Only good for diagonally dominant systems, must have elements on diagonals though.
+/// The initial guess (start) is always zero, convergence value may be specified
+/// with epsilon, relaxation decreases the step size (useful when may rows exist in
+/// A)
+///
 template<typename Scalar>
 class ProjectedGaussSeidelSolver
 {
@@ -31,10 +40,27 @@ public:
     // Sets the vector to be used for the solve
     //void setGuess(Matrix<Scalar, -1, 1>& g) { x = g; }
     void setA(Eigen::SparseMatrix<Scalar>* A) { this->m_A = A; }
+
+    ///
+    /// \brief Set the maximum number of iterations
+    ///
     void setMaxIterations(const unsigned int maxIterations) { this->m_maxIterations = maxIterations; }
+
+    ///
+    /// \brief Similar to step size can be used to avoid overshooting the solution
+    ///
     void setRelaxation(const Scalar relaxation) { this->m_relaxation = relaxation; }
+
+    ///
+    /// \brief Stops when energy=(x_i+1-x_i).norm() < epsilon, when the solution isn't
+    /// changing anymore
+    ///
     void setEpsilon(const Scalar epsilon) { this->m_epsilon = epsilon; }
-    const double getEnergy() { return m_conv; } // Get final convergence
+
+    ///
+    /// \brief Energy is defined as energy=(x_i+1-x_i).norm()
+    ///
+    const double getEnergy() const { return m_conv; }
 
     Eigen::Matrix<Scalar, -1, 1>& solve(const Eigen::Matrix<Scalar, -1, 1>& b, const Eigen::Matrix<Scalar, -1, 2>& cu)
     {
@@ -66,6 +92,8 @@ public:
                     delta += A.coeff(r, c) * m_x[c];
                 }
 
+                // PGS can't converge for non-diagonal elements so its assumed
+                // we have these
                 delta = (b[r] - delta) / A.coeff(r, r);
                 // Apply relaxation factor
                 m_x(r) += m_relaxation * (delta - m_x(r));