New test file added

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
@@ -720,7 +722,7 @@ void ComputeGradient::get_smooth_bspline_3D(MeshData<T>& input,APRParameters& pa
 }
 
 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){
+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
     //
@@ -733,52 +735,62 @@ void ComputeGradient::calc_bspline_fd_ds_mag(MeshData<T> &input, MeshData<S> &gr
     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;
 
+    // 4                 4
+    // 2  1,3   ...   1  2  3 ...
+    // 5                 5
     #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 j = 0; j < z_num; ++j) {
+        S *left = input.mesh.begin() + j * xnumynum + 1 * y_num;
+        S *center = input.mesh.begin() + j * xnumynum;
 
         //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);
-
-        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;
+        const size_t j_p = std::min(z_num - 1, j + 1);
+        float g[x_num][y_num];
+        for (size_t i = 0; i < x_num - 1; ++i) {
+            S *up = input.mesh.begin() + j_m * xnumynum + i * y_num;
+            S *down = input.mesh.begin() + j_p * xnumynum + i * y_num;
+            S *right = input.mesh.begin() + j * xnumynum + (i + 1) * 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));
-
+            temp[0] = sqrt(pow((right[0] - left[0]) / (2 * hx), 2.0) + pow((down[0] - up[0]) / (2 * hz), 2.0));
+            g[0][i] = temp[0];
             //do the y gradient
-            #ifdef HAVE_OPENMP
-	        #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));
+#ifdef HAVE_OPENMP
+#pragma omp simd
+#endif
+            for (size_t k = 1; k < y_num - 1; ++k) {
+                temp[k] = sqrt(pow((right[k] - left[k]) / (2 * hx), 2.0) + pow((down[k] - up[k]) / (2 * hz), 2.0) +
+                               pow((center[k + 1] - center[k - 1]) / (2 * hy), 2.0));
+                g[k][i] = temp[k];
             }
 
-            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));
+            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));
+            g[y_num - 1][i] = temp[y_num - 1];
 
-            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);
+            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;
                 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);
+        }
+        for (int y = 0; y < y_num; ++y) {
+            for (int x = 0; x < x_num; ++x) {
+                std::cout << g[y][x] << " ";
+            }
+            std::cout << "\n";
         }
     }
 }

src/data_structures/Mesh/MeshData.hpp

@@ -306,6 +306,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

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/MeshDataTest.cpp

@@ -170,7 +170,7 @@ namespace {
         }
     }
 
-    TEST_F(MeshDataTest, PreallocateTest) {
+    TEST_F(MeshDataTest, InitDownsampledTest) {
         {
             MeshData<int> md;
             md.initDownsampled(3, 5, 7);
@@ -189,6 +189,25 @@ namespace {
             int size = 2 * 3 * 4;
             ASSERT_EQ(md.mesh.size(), size);
         }
+        {
+            MeshData<int> md;
+            md.initDownsampled(MeshData<char>(4, 6, 8));
+            ASSERT_EQ(md.y_num, 2);
+            ASSERT_EQ(md.x_num, 3);
+            ASSERT_EQ(md.z_num, 4);
+            int size = 2 * 3 * 4;
+            ASSERT_EQ(md.mesh.size(), size);
+        }
+        {
+            MeshData<int> md;
+            md.initDownsampled(MeshData<float>(3, 5, 7), 2);
+            ASSERT_EQ(md.y_num, 2);
+            ASSERT_EQ(md.x_num, 3);
+            ASSERT_EQ(md.z_num, 4);
+            int size = 2 * 3 * 4;
+            ASSERT_EQ(md.mesh.size(), size);
+            for (size_t i = 0; i < md.mesh.size(); ++i) ASSERT_EQ(md.mesh[i], 2);
+        }
     }
 
     TEST(MeshDataSimpleTest, UnaryOpTest) {