Merge branch 'newTestForGradient' into develop

Author: krzysg

src/algorithm/APRConverter.hpp

@@ -104,15 +104,13 @@ bool APRConverter<ImageType>::get_apr_method(APR<ImageType> &aAPR, MeshData<T>&
 
     //assuming uint16, the total memory cost shoudl be approximately (1 + 1 + 1/8 + 2/8 + 2/8) = 2 5/8 original image size in u16bit
     //storage of the particle cell tree for computing the pulling scheme
-    MeshData<ImageType> image_temp; // global image variable useful for passing between methods, or re-using memory (should be the only full sized copy of the image)
-    MeshData<ImageType> grad_temp; // should be a down-sampled image
-    MeshData<float> local_scale_temp; //   Used as down-sampled images for some averaging steps where it is useful to not lose precision, or get over-flow errors
-    MeshData<float> local_scale_temp2;
-
     allocation_timer.start_timer("init and copy image");
-    image_temp.init(input_image);
+    MeshData<ImageType> image_temp(input_image, false /* don't copy */); // global image variable useful for passing between methods, or re-using memory (should be the only full sized copy of the image)
+    MeshData<ImageType> grad_temp; // should be a down-sampled image
     grad_temp.initDownsampled(input_image.y_num, input_image.x_num, input_image.z_num, 0);
+    MeshData<float> local_scale_temp; // Used as down-sampled images for some averaging steps where it is useful to not lose precision, or get over-flow errors
     local_scale_temp.initDownsampled(input_image.y_num, input_image.x_num, input_image.z_num);
+    MeshData<float> local_scale_temp2;
     local_scale_temp2.initDownsampled(input_image.y_num, input_image.x_num, input_image.z_num);
     allocation_timer.stop_timer();
 
@@ -126,7 +124,7 @@ bool APRConverter<ImageType>::get_apr_method(APR<ImageType> &aAPR, MeshData<T>&
     //offset image by factor (this is required if there are zero areas in the background with uint16_t and uint8_t images, as the Bspline co-efficients otherwise may be negative!)
     // Warning both of these could result in over-flow (if your image is non zero, with a 'buffer' and has intensities up to uint16_t maximum value then set image_type = "", i.e. uncomment the following line)
     float bspline_offset = 0;
-    if(std::is_same<uint16_t, ImageType>::value){
+    if (std::is_same<uint16_t, ImageType>::value) {
         bspline_offset = 100;
         image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); });
     } else if (std::is_same<uint8_t, ImageType>::value){
@@ -227,15 +225,10 @@ void APRConverter<ImageType>::get_local_particle_cell_set(MeshData<ImageType> &g
 
 template<typename ImageType>
 void APRConverter<ImageType>::get_gradient(MeshData<ImageType> &image_temp, MeshData<ImageType> &grad_temp, MeshData<float> &local_scale_temp, MeshData<float> &local_scale_temp2, float bspline_offset) {
-    //
     //  Bevan Cheeseman 2018
-    //
     //  Calculate the gradient from the input image. (You could replace this method with your own)
-    //
     //  Input: full sized image.
-    //
     //  Output: down-sampled by 2 gradient magnitude (Note, the gradient is calculated at pixel level then maximum down sampled within the loops below)
-    //
 
     fine_grained_timer.verbose_flag = false;
 

src/algorithm/ComputeGradient.hpp

@@ -6,6 +6,8 @@
 #define PARTPLAY_GRADIENT_HPP
 
 #include "src/data_structures/Mesh/MeshData.hpp"
+#include "src/io/TiffUtils.hpp"
+
 #ifdef HAVE_OPENMP
 	#include "omp.h"
 #endif
@@ -49,9 +51,9 @@ class ComputeGradient {
 
 // Gradient computation
 
-    template<typename T, typename S>
+    template<typename S>
     void
-    calc_bspline_fd_ds_mag(MeshData<T> &input, MeshData<S> &grad, const float hx, const float hy, const float hz);
+    calc_bspline_fd_ds_mag(const MeshData<S> &input, MeshData<S> &grad, const float hx, const float hy, const float hz);
 
     template<typename T,typename S>
     void mask_gradient(MeshData<T>& grad_ds,MeshData<S>& temp_ds,MeshData<T>& temp_full,APRParameters& par);
@@ -719,74 +721,85 @@ void ComputeGradient::get_smooth_bspline_3D(MeshData<T>& input,APRParameters& pa
     spline_timer.stop_timer();
 }
 
-template<typename T,typename S>
-void ComputeGradient::calc_bspline_fd_ds_mag(MeshData<T> &input, MeshData<S> &grad, const float hx, const float hy,const float hz){
-    //
-    //  Bevan Cheeseman 2016
-    //
-    //  Calculate fd filt, for xgrad with bsplines
-
+/**
+ * Calculates downsampled gradient (maximum magnitude) with 'replicate' boundary approach (nearest border value)
+ * @param input - input mesh
+ * @param grad - output gradient (must be initialized)
+ * @param hx - step in x dir
+ * @param hy - step in y dir
+ * @param hz - step in z dir
+ */
+template<typename S>
+void ComputeGradient::calc_bspline_fd_ds_mag(const MeshData<S> &input, MeshData<S> &grad, const float hx, const float hy,const float hz) {
     const size_t z_num = input.z_num;
     const size_t x_num = input.x_num;
     const size_t y_num = input.y_num;
 
     const size_t x_num_ds = grad.x_num;
     const size_t y_num_ds = grad.y_num;
 
-    const float a1 = -1.0/2.0;
-    const float a3 = 1.0/2.0;
-
-    std::vector<S> temp(y_num,0);
+    std::vector<S> temp(y_num, 0);
     const size_t xnumynum = x_num * y_num;
 
     #ifdef HAVE_OPENMP
-	#pragma omp parallel for default(shared)firstprivate(temp)
+    #pragma omp parallel for default(shared) firstprivate(temp)
     #endif
-    for (size_t j = 0;j < z_num; ++j) {
-        S *temp_vec_1 = input.mesh.begin() + j*x_num*y_num + 1*y_num;
-        S *temp_vec_2 = input.mesh.begin() + j*x_num*y_num;
+    for (size_t z = 0; z < z_num; ++z) {
+        // Belows pointers up, down... are forming stencil in X (left <-> right) and Z ( up <-> down) direction and
+        // are pointing to whole Y column. If out of bounds then 'replicate' (nearest array border value) approach is used.
+        //
+        //                 up
+        //   ...   left  center  right ...
+        //                down
+        const S *left = input.mesh.begin() + z * xnumynum + 0 * y_num; // boundary value is chosen
+        const S *center = input.mesh.begin() + z * xnumynum + 0 * y_num;
 
         //LHS boundary condition is accounted for wiht this initialization
-        const size_t j_m = j > 0 ? j - 1 : 0;
-        const size_t j_p = std::min(z_num-1, j+1);
+        const size_t zMinus = z > 0 ? z - 1 : 0 /* boundary */;
+        const size_t zPlus = std::min(z + 1, z_num - 1 /* boundary */);
 
-        for (size_t i = 0; i < x_num-1; ++i) {
-            S *temp_vec_4 = input.mesh.begin() + j_m*xnumynum + i * y_num;
-            S *temp_vec_5 = input.mesh.begin() + j_p*xnumynum + i * y_num;
-            S *temp_vec_3 = input.mesh.begin() + j*x_num*y_num + (i+1)*y_num;
+        for (size_t x = 0; x < x_num; ++x) {
+            const S *up = input.mesh.begin() + zMinus * xnumynum + x * y_num;
+            const S *down = input.mesh.begin() + zPlus * xnumynum + x * y_num;
+            const size_t xPlus = std::min(x + 1, x_num - 1 /* boundary */);
+            const S *right = input.mesh.begin() + z * xnumynum + xPlus * y_num;
 
             //compute the boundary values
-            temp[0] = sqrt(pow((a1*temp_vec_1[0] + a3*temp_vec_3[0])/hx,2.0)  + pow((a1*temp_vec_4[0] + a3*temp_vec_5[0])/hz,2.0));
+            if (y_num >= 2) {
+                temp[0] = sqrt(pow((right[0] - left[0]) / (2 * hx), 2.0) + pow((down[0] - up[0]) / (2 * hz), 2.0) + pow((center[1] - center[0 /* boundary */]) / (2 * hy), 2.0));
+                temp[y_num - 1] = sqrt(pow((right[y_num - 1] - left[y_num - 1]) / (2 * hx), 2.0) + pow((down[y_num - 1] - up[y_num - 1]) / (2 * hz), 2.0) + pow((center[y_num - 1 /* boundary */] - center[y_num - 2]) / (2 * hy), 2.0));
+            }
+            else {
+                temp[0] = 0; // same values minus same values in x/y/z
+            }
 
-            //do the y gradient
+            //do the y gradient in range 1..y_num-2
             #ifdef HAVE_OPENMP
-	        #pragma omp simd
+            #pragma omp simd
             #endif
-            for (size_t k = 1; k < (y_num-1); ++k) {
-                temp[k] =  sqrt(pow((a1*temp_vec_4[k] + a3*temp_vec_5[k])/hz,2.0) +  pow((a1*temp_vec_2[k-1] + a3*temp_vec_2[k+1])/hy,2.0) + pow((a1*temp_vec_1[k] + a3*temp_vec_3[k])/hx,2.0));
+            for (size_t y = 1; y < y_num - 1; ++y) {
+                temp[y] = sqrt(pow((right[y] - left[y]) / (2 * hx), 2.0) + pow((down[y] - up[y]) / (2 * hz), 2.0) + pow((center[y + 1] - center[y - 1]) / (2 * hy), 2.0));
             }
 
-            temp[y_num - 1] = sqrt(pow((a1*temp_vec_1[(y_num-1)] + a3*temp_vec_3[(y_num-1)])/hx,2.0)  + pow((a1*temp_vec_4[(y_num-1)] + a3*temp_vec_5[(y_num-1)])/hz,2.0));
-
-            int64_t j_2 = j/2;
-            int64_t i_2 = i/2;
-
-            for (size_t k = 0; k < (y_num_ds); ++k) {
-                size_t k_s = std::min(2*k+1, y_num-1);
-                const size_t idx = j_2 * x_num_ds * y_num_ds + i_2 * y_num_ds + k;
+            // Set as a downsampled gradient maximum from 2x2x2 gradient cubes
+            int64_t z_2 = z / 2;
+            int64_t x_2 = x / 2;
+            for (size_t k = 0; k < y_num_ds; ++k) {
+                size_t k_s = std::min(2 * k + 1, y_num - 1);
+                const size_t idx = z_2 * x_num_ds * y_num_ds + x_2 * y_num_ds + k;
                 grad.mesh[idx] = std::max(temp[2 * k], std::max(temp[k_s], grad.mesh[idx]));
             }
 
-            std::swap(temp_vec_1, temp_vec_2);
-            std::swap(temp_vec_2, temp_vec_3);
+            // move left, center to current center, right (both +1 to right)
+            std::swap(left, center);
+            std::swap(center, right);
         }
     }
 }
 
-/*
+/**
  * Caclulation of signal value from B-Spline co-efficients
  */
-
 template<typename T>
 void ComputeGradient::calc_inv_bspline_y(MeshData<T>& input){
     //

src/data_structures/Mesh/MeshData.hpp

@@ -17,6 +17,8 @@
 #include <cmath>
 #include <memory>
 #include <sstream>
+#include <iostream>
+#include <iomanip>
 
 #include "src/misc/APRTimer.hpp"
 
@@ -306,6 +308,33 @@ public :
         init(y_num_ds, x_num_ds, z_num_ds, aInitVal);
     }
 
+    /**
+     * Initializes mesh with size of half of provided mesh dimensions (rounding up if not divisible by 2)
+     * @param aMesh - mesh used to get dimensions
+     */
+    template <typename U>
+    void initDownsampled(const MeshData<U> &aMesh) {
+        const int z_num_ds = ceil(1.0*aMesh.z_num/2.0);
+        const int x_num_ds = ceil(1.0*aMesh.x_num/2.0);
+        const int y_num_ds = ceil(1.0*aMesh.y_num/2.0);
+
+        init(y_num_ds, x_num_ds, z_num_ds);
+    }
+
+    /**
+     * Initializes mesh with size of half of provided mesh dimensions (rounding up if not divisible by 2) and initialize values
+     * @param aMesh - mesh used to get dimensions
+     * @param aInitVal
+     */
+    template <typename U>
+    void initDownsampled(const MeshData<U> &aMesh, T aInitVal) {
+        const int z_num_ds = ceil(1.0*aMesh.z_num/2.0);
+        const int x_num_ds = ceil(1.0*aMesh.x_num/2.0);
+        const int y_num_ds = ceil(1.0*aMesh.y_num/2.0);
+
+        init(y_num_ds, x_num_ds, z_num_ds, aInitVal);
+    }
+
     /**
      * Swaps data of meshes this <-> aObj
      * @param aObj
@@ -369,18 +398,32 @@ public :
     }
 
     /**
-     * Prints X-Y planes of mesh (for debug/test purposses - use only on small meshes)
+     * Prints X-Y or X-Z planes of mesh (for debug/test purposses - use only on small meshes)
      */
-    void printMesh() const {
-        for (size_t z = 0; z < z_num; ++z) {
-            std::cout << "z=" << z << "\n";
+    void printMesh(int aColumnWidth, bool aXYplanes = true) const {
+        if (aXYplanes) {
+            for (size_t z = 0; z < z_num; ++z) {
+                std::cout << "z=" << z << "\n";
+                for (size_t y = 0; y < y_num; ++y) {
+                    for (size_t x = 0; x < x_num; ++x) {
+                        std::cout << std::setw(aColumnWidth) << at(y, x, z) << " ";
+                    }
+                    std::cout << "\n";
+                }
+                std::cout << std::endl;
+            }
+        }
+        else { // X-Z planes
             for (size_t y = 0; y < y_num; ++y) {
-                for (size_t x = 0; x < x_num; ++x) {
-                    std::cout << at(y, x, z) << " ";
+                std::cout << "y=" << y << "\n";
+                for (size_t z = 0; z < z_num; ++z) {
+                    for (size_t x = 0; x < x_num; ++x) {
+                        std::cout << std::setw(aColumnWidth) << at(y, x, z) << " ";
+                    }
+                    std::cout << "\n";
                 }
-                std::cout << "\n";
+                std::cout << std::endl;
             }
-            std::cout << std::endl;
         }
     }
 

test/CMakeLists.txt

@@ -1,7 +1,9 @@
 add_executable(testMeshData MeshDataTest.cpp)
 add_executable(testTiff TiffTest.cpp)
 add_executable(testAPR APRTest.cpp)
+add_executable(testComputeGradient ComputeGradientTest.cpp)
 
 target_link_libraries(testMeshData ${GTEST_LIBRARIES} libapr)
 target_link_libraries(testTiff ${GTEST_LIBRARIES} libapr)
 target_link_libraries(testAPR ${GTEST_LIBRARIES} libapr)
+target_link_libraries(testComputeGradient ${GTEST_LIBRARIES} libapr)

test/ComputeGradientTest.cpp

@@ -0,0 +1,112 @@
+/*
+ * Created by Krzysztof Gonciarz 2018
+ */
+
+#include <gtest/gtest.h>
+#include "src/data_structures/Mesh/MeshData.hpp"
+#include "src/algorithm/ComputeGradient.hpp"
+
+namespace {
+    /**
+     * Compares mesh with provided data
+     * @param mesh
+     * @param data - data with [Z][Y][X] structure
+     * @return true if same
+     */
+    template <typename T>
+    bool compare(MeshData<T> &mesh, const float *data, const float epsilon) {
+        size_t dataIdx = 0;
+        for (size_t z = 0; z < mesh.z_num; ++z) {
+            for (size_t y = 0; y < mesh.y_num; ++y) {
+                for (size_t x = 0; x < mesh.x_num; ++x) {
+                    bool v = std::abs(mesh(y, x, z) - data[dataIdx]) < epsilon;
+                    if (v == false) {
+                        std::cerr << "Mesh and expected data differ. First place at (Y, X, Z) = " << y << ", " << x << ", " << z << ") " << mesh(y, x, z) << " vs " << data[dataIdx] << std::endl;
+                        return false;
+                    }
+                    ++dataIdx;
+                }
+            }
+        }
+        return true;
+    }
+
+    TEST(ComputeGradientTest, 2D_XY) {
+        {   // Corner points
+            MeshData<float> m(6, 6, 1, 0);
+            // expect gradient is 3x3 X/Y plane
+            float expect[] = {1.41, 0, 4.24,
+                              0, 0, 0,
+                              2.82, 0, 5.65};
+            // put values in corners
+            m(0, 0, 0) = 2;
+            m(5, 0, 0) = 4;
+            m(0, 5, 0) = 6;
+            m(5, 5, 0) = 8;
+            MeshData<float> grad;
+            grad.initDownsampled(m, 0);
+            ComputeGradient cg;
+            cg.calc_bspline_fd_ds_mag(m, grad, 1, 1, 1);
+            ASSERT_TRUE(compare(grad, expect, 0.01));
+        }
+        {   // In the middle
+            MeshData<float> m(6, 6, 1, 0);
+            // expect gradient is 3x3 X/Y plane
+            float expect[] = {1, 1, 0,
+                              1, 0, 0,
+                              0, 0, 0};
+            // put values in corners
+            m(1, 1, 0) = 2;
+            MeshData<float> grad;
+            grad.initDownsampled(m, 0);
+            ComputeGradient cg;
+            cg.calc_bspline_fd_ds_mag(m, grad, 1, 1, 1);
+            ASSERT_TRUE(compare(grad, expect, 0.01));
+        }
+        {   // One pixel image 1x1x1
+            MeshData<float> m(1, 1, 1, 0);
+            // expect gradient is 3x3 X/Y plane
+            float expect[] = {0};
+            // put values in corners
+            m(0, 0, 0) = 2;
+            MeshData<float> grad;
+            grad.initDownsampled(m, 0);
+            ComputeGradient cg;
+            cg.calc_bspline_fd_ds_mag(m, grad, 1, 1, 1);
+            ASSERT_TRUE(compare(grad, expect, 0.01));
+        }
+
+    }
+
+    TEST(ComputeGradientTest, Corners3D) {
+        MeshData<float> m(6, 6, 4, 0);
+        // expect gradient is 3x3x2 X/Y/Z plane
+        float expect[] = { 1.73, 0, 5.19,
+                           0,    0,    0,
+                           3.46, 0, 6.92,
+
+                           8.66, 0, 12.12,
+                           0,    0,     0,
+                          10.39, 0, 13.85 };
+        // put values in corners
+        m(0, 0, 0) = 2;
+        m(5, 0, 0) = 4;
+        m(0, 5, 0) = 6;
+        m(5, 5, 0) = 8;
+        m(0, 0, 3) = 10;
+        m(5, 0, 3) = 12;
+        m(0, 5, 3) = 14;
+        m(5, 5, 3) = 16;
+
+        MeshData<float> grad;
+        grad.initDownsampled(m, 0);
+        ComputeGradient cg;
+        cg.calc_bspline_fd_ds_mag(m, grad, 1, 1, 1);
+        ASSERT_TRUE(compare(grad, expect, 0.01));
+    }
+}
+
+int main(int argc, char **argv) {
+    testing::InitGoogleTest(&argc, argv);
+    return RUN_ALL_TESTS();
+}