Interactive debug

.github/workflows/cmake.yml

@@ -37,15 +37,15 @@ jobs:
       VCPKG_ROOT: ${{ github.workspace }}/vcpkg
 
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
         with:
           submodules: true
 
       # Setup the build machine with the most recent versions of CMake and Ninja. Both are cached if not already: on subsequent runs both will be quickly restored from GitHub cache service.
       - uses: lukka/get-cmake@latest
       # Restore both vcpkg and its artifacts from the GitHub cache service.
       - name: Restore vcpkg and its artifacts.
-        uses: actions/cache@v2
+        uses: actions/cache@v3
         with:
           # The first path is where vcpkg generates artifacts while consuming the vcpkg.json manifest file.
           # The second path is the location of vcpkg (it contains the vcpkg executable and data files).
@@ -64,7 +64,7 @@ jobs:
 
       - name: Install OpenMP dependencies with brew for OSX
         if: contains(matrix.os,'macos')
-        run: brew install libomp #todo caching
+        run: brew install libomp pkg-config #todo caching
 
       - name: Show content of workspace after cache has been restored
         run: find $RUNNER_WORKSPACE
@@ -75,17 +75,21 @@ jobs:
       # Run CMake to generate Ninja project files, using the vcpkg's toolchain file to resolve and install the dependencies as specified in vcpkg.json.
       - name: Install dependencies and generate project files
         run: |
-          cmake -S "${{ github.workspace }}" -B "${{ env.CMAKE_BUILD_DIR }}" -DCMAKE_TOOLCHAIN_FILE="${{ env.VCPKG_ROOT }}/scripts/buildsystems/vcpkg.cmake" -DAPR_BUILD_STATIC_LIB=ON -DAPR_BUILD_SHARED_LIB=OFF -DAPR_BUILD_EXAMPLES=ON -DAPR_TESTS=ON -DAPR_USE_CUDA=OFF -DAPR_PREFER_EXTERNAL_BLOSC=${{ matrix.extblosc }} -DAPR_PREFER_EXTERNAL_BLOSC=${{ matrix.extblosc }} -DAPR_USE_OPENMP=${{ matrix.openmp }} ${{ matrix.buildargs }}
+          cmake -S "${{ github.workspace }}" -B "${{ env.CMAKE_BUILD_DIR }}" -DCMAKE_TOOLCHAIN_FILE="${{ env.VCPKG_ROOT }}/scripts/buildsystems/vcpkg.cmake" -DCMAKE_BUILD_TYPE=Release -DAPR_BUILD_STATIC_LIB=ON -DAPR_BUILD_SHARED_LIB=OFF -DAPR_BUILD_EXAMPLES=ON -DAPR_TESTS=ON -DAPR_USE_CUDA=OFF -DAPR_PREFER_EXTERNAL_BLOSC=${{ matrix.extblosc }} -DAPR_USE_OPENMP=${{ matrix.openmp }} ${{ matrix.buildargs }}
 
       # Build the whole project with Ninja (which is spawn by CMake).
       - name: Build
         run: |
           cmake --build "${{ env.CMAKE_BUILD_DIR }}" --parallel 4 --config Release
 
+
+      - name: Setup tmate session
+        uses: mxschmitt/action-tmate@v3
+
       # Run tests
       - name: tests
         run: |
-          ctest --test-dir "${{ env.CMAKE_BUILD_DIR }}"
+          ctest --test-dir "${{ env.CMAKE_BUILD_DIR }}" --output-on-failure
 
       - name: Show content of workspace at its completion
         run: find $RUNNER_WORKSPACE

.gitmodules

@@ -1,9 +1,6 @@
 [submodule "external/glm"]
 	path = external/glm
 	url = https://github.com/g-truc/glm.git
-[submodule "external/gtest"]
-	path = external/gtest
-	url = https://github.com/google/googletest
 [submodule "external/c-blosc"]
 	path = external/c-blosc
 	url = https://github.com/Blosc/c-blosc

CMakeLists.txt

@@ -1,14 +1,18 @@
 ###############################################################################
 # APR - Adaptive Particle Representation
 ###############################################################################
-cmake_minimum_required(VERSION 3.10)
+cmake_minimum_required(VERSION 3.14)
 project(APR DESCRIPTION "Adaptive Particle Representation library")
 
 message(STATUS "CMAKE VERSION ${CMAKE_VERSION}")
 
 set(CMAKE_CXX_STANDARD 14)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
 
+if(POLICY CMP0135)
+    cmake_policy(SET CMP0135 NEW)
+endif(POLICY CMP0135)
+
 # APR build options:
 option(APR_INSTALL "Install APR library" OFF)
 option(APR_BUILD_SHARED_LIB "Builds shared library" OFF)
@@ -360,12 +364,22 @@ if(APR_TESTS)
     if(GTEST_FOUND)
         include_directories(${GTEST_INCLUDE_DIRS})
     else(GTEST_FOUND)
+        message(STATUS "APR: GTest not found, it will be downloaded and built")
+
+        include(FetchContent)
+        FetchContent_Declare(
+            googletest
+            # Specify the commit you depend on and update it regularly.
+            URL https://github.com/google/googletest/archive/refs/tags/v1.13.0.zip
+        )
+
+        # For Windows: Prevent overriding the parent project's compiler/linker settings
+        set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
         set(BUILD_GMOCK OFF CACHE BOOL "" FORCE)
-        set(BUILD_GTEST ON CACHE BOOL "" FORCE)
         set(INSTALL_GTEST OFF CACHE BOOL "" FORCE)
-        message(STATUS "APR: GTest not found, using internal gtest")
-        add_subdirectory("external/gtest")
-        set(GTEST_LIBRARIES gtest)
+        FetchContent_MakeAvailable(googletest)
+
+        set(GTEST_LIBRARIES GTest::gtest_main)
     endif(GTEST_FOUND)
 
     enable_testing()

examples/Example_get_apr.cpp

@@ -48,7 +48,7 @@ int runAPR(cmdLineOptions options) {
     //the apr datastructure
     APR apr;
 
-    APRConverter<uint16_t> aprConverter;
+    APRConverter<float> aprConverter;
 
     //read in the command line options into the parameters file
     aprConverter.par.Ip_th = options.Ip_th;

examples/Example_get_apr_by_block.cpp

@@ -44,7 +44,7 @@ must be used. The exact influence of this has not yet been studied.
 int runAPR(cmdLineOptions options) {
 
     APR apr;
-    APRConverterBatch<uint16_t> aprConverter;
+    APRConverterBatch<float> aprConverter;
 
     //read in the command line options into the parameters file
     aprConverter.par.Ip_th = options.Ip_th;
@@ -55,6 +55,7 @@ int runAPR(cmdLineOptions options) {
     aprConverter.par.neighborhood_optimization = options.neighborhood_optimization;
     aprConverter.par.output_steps = options.output_steps;
     aprConverter.par.grad_th = options.grad_th;
+    aprConverter.par.auto_parameters = options.auto_parameters;
 
     //where things are
     aprConverter.par.input_image_name = options.input;
@@ -234,7 +235,6 @@ cmdLineOptions read_command_line_options(int argc, char **argv){
     if(command_option_exists(argv, argv + argc, "-neighborhood_optimization_off"))
     {
         result.neighborhood_optimization = false;
-
     }
 
     if(command_option_exists(argv, argv + argc, "-output_steps"))
@@ -247,5 +247,10 @@ cmdLineOptions read_command_line_options(int argc, char **argv){
         result.store_tree = true;
     }
 
+    if(command_option_exists(argv, argv + argc, "-auto_parameters"))
+    {
+        result.auto_parameters = true;
+    }
+
     return result;
 }

examples/Example_get_apr_by_block.hpp

@@ -30,6 +30,8 @@ struct cmdLineOptions{
     float rel_error = 0.1;
     float grad_th = 1;
 
+    bool auto_parameters = false;
+
     int z_block_size = 128;
     int z_ghost = 16; // number of "ghost slices" to use in the APR pipeline
     int z_ghost_sampling = 64; // number of "ghost slices" to use when sampling intensities

examples/Example_reconstruct_image.cpp

@@ -47,7 +47,7 @@ struct cmdLineOptions{
     bool output_spatial_properties = false;
     bool output_pc_recon = false;
     bool output_smooth_recon = false;
-
+    float gaussian_noise_sigma = 0.0f;
 };
 
 static bool command_option_exists(char **begin, char **end, const std::string &option) {
@@ -99,6 +99,10 @@ static cmdLineOptions read_command_line_options(int argc, char **argv) {
         result.output_spatial_properties = true;
     }
 
+    if(command_option_exists(argv, argv + argc, "-noise")) {
+        result.gaussian_noise_sigma = std::stof(std::string(get_command_option(argv, argv + argc, "-noise")));
+    }
+
     if(!(result.output_pc_recon || result.output_smooth_recon || result.output_spatial_properties)){
         //default is pc recon
         result.output_pc_recon = true;
@@ -107,7 +111,7 @@ static cmdLineOptions read_command_line_options(int argc, char **argv) {
     return result;
 }
 template<typename T>
-void add_random_to_img(PixelData<T>& img,float sd){
+void add_random_to_img(PixelData<T>& img, float sd){
 
     std::default_random_engine generator;
     std::normal_distribution<float> distribution(0.0,sd);
@@ -149,32 +153,26 @@ int main(int argc, char **argv) {
     apr.name = options.output;
     timer.stop_timer();
 
-    // Intentionaly block-scoped since local recon_pc will be destructed when block ends and release memory.
-    {
-
-        if(options.output_pc_recon) {
-            //create mesh data structure for reconstruction
-            bool add_random_gitter = true;
-
-            PixelData<uint16_t> recon_pc;
+    if(options.output_pc_recon) {
+        //create mesh data structure for reconstruction
+        PixelData<uint16_t> recon_pc;
 
-            timer.start_timer("pc interp");
-            //perform piece-wise constant interpolation
-            APRReconstruction::reconstruct_constant(apr,recon_pc, parts);
-            timer.stop_timer();
+        timer.start_timer("pc interp");
+        //perform piece-wise constant interpolation
+        APRReconstruction::reconstruct_constant(apr, recon_pc, parts);
+        timer.stop_timer();
 
-            if(add_random_gitter){
-                add_random_to_img(recon_pc,1.0f);
-            }
+        if(options.gaussian_noise_sigma > 0.0f) {
+            add_random_to_img(recon_pc, options.gaussian_noise_sigma);
+        }
 
-            float elapsed_seconds = timer.t2 - timer.t1;
-            std::cout << "PC recon "
-                      << (recon_pc.x_num * recon_pc.y_num * recon_pc.z_num * 2) / (elapsed_seconds * 1000000.0f)
-                      << " MB per second" << std::endl;
+        float elapsed_seconds = timer.t2 - timer.t1;
+        std::cout << "PC recon "
+                    << (recon_pc.x_num * recon_pc.y_num * recon_pc.z_num * 2) / (elapsed_seconds * 1000000.0f)
+                    << " MB per second" << std::endl;
 
-            // write output as tiff
-            TiffUtils::saveMeshAsTiff(options.directory + apr.name + "_pc.tif", recon_pc);
-        }
+        // write output as tiff
+        TiffUtils::saveMeshAsTiff(options.directory + apr.name + "_pc.tif", recon_pc);
     }
 
     //////////////////////////

src/algorithm/APRConverter.hpp

@@ -99,6 +99,8 @@ class APRConverter {
         return true;
     }
 
+    float bspline_offset = 0;
+
 protected:
 
     template<typename T>
@@ -108,8 +110,6 @@ class APRConverter {
 
     //get apr without setting parameters, and with an already loaded image.
 
-    float bspline_offset = 0;
-
     //DATA (so it can be re-used)
 
     PixelData<ImageType> grad_temp; // should be a down-sampled image
@@ -133,26 +133,6 @@ class APRConverter {
 };
 
 
-template <typename T>
-struct MinMax{T min; T max; };
-
-template <typename T>
-static MinMax<T> getMinMax(const PixelData<T>& input_image) {
-    T minVal = std::numeric_limits<T>::max();
-    T maxVal = std::numeric_limits<T>::min();
-
-#ifdef HAVE_OPENMP
-#pragma omp parallel for default(shared) reduction(max:maxVal) reduction(min:minVal)
-#endif
-    for (size_t i = 0; i < input_image.mesh.size(); ++i) {
-        T val = input_image.mesh[i];
-        if (val > maxVal) maxVal = val;
-        if (val < minVal) minVal = val;
-    }
-
-    return MinMax<T>{minVal, maxVal};
-}
-
 template<typename ImageType>
 void APRConverter<ImageType>::initPipelineMemory(int y_num,int x_num,int z_num){
     //initializes the internal memory to be used in the pipeline.
@@ -224,7 +204,7 @@ void APRConverter<ImageType>::computeL(APR& aAPR,PixelData<T>& input_image){
     //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
     allocation_timer.start_timer("init and copy image");
-    PixelData<ImageType> image_temp(input_image, false /* don't copy */, true /* pinned memory */); // global image variable useful for passing between methods, or re-using memory (should be the only full sized copy of the image)
+    PixelData<ImageType> image_temp(input_image, false /* don't copy */, false /* pinned memory */); // global image variable useful for passing between methods, or re-using memory (should be the only full sized copy of the image)
 
     allocation_timer.stop_timer();
 
@@ -233,17 +213,12 @@ void APRConverter<ImageType>::computeL(APR& aAPR,PixelData<T>& input_image){
     ////////////////////////
 
     fine_grained_timer.start_timer("offset image");
-    //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)
 
-    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){
-        bspline_offset = 5;
-        image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); });
-    } else {
+    if (std::is_floating_point<ImageType>::value) {
         image_temp.copyFromMesh(input_image);
+    } else {
+        bspline_offset = compute_bspline_offset<ImageType>(input_image, par.lambda);
+        image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); });
     }
 
     fine_grained_timer.stop_timer();

src/algorithm/APRConverterBatch.hpp

@@ -214,17 +214,14 @@ bool APRConverterBatch<ImageType>::get_apr_method_patch(APR &aAPR, PixelData<T>&
     computation_timer.start_timer("Calculations");
 
     fine_grained_timer.start_timer("offset image");
-    //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)
-    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){
-        bspline_offset = 5;
-        image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); });
-    } else {
+
+    if (std::is_floating_point<ImageType>::value) {
         image_temp.copyFromMesh(input_image);
+    } else {
+        bspline_offset = compute_bspline_offset<ImageType>(input_image, par.lambda);
+        image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); });
     }
+
     fine_grained_timer.stop_timer();
 
     method_timer.start_timer("compute_gradient_magnitude_using_bsplines");

src/algorithm/AutoParameters.hpp

@@ -157,4 +157,43 @@ void autoParametersLiEntropy(APRParameters& par,
     }
 }
 
+
+
+template <typename T>
+struct MinMax{T min; T max; };
+
+template <typename T>
+static MinMax<T> getMinMax(const PixelData<T>& input_image) {
+    T minVal = std::numeric_limits<T>::max();
+    T maxVal = std::numeric_limits<T>::min();
+
+#ifdef HAVE_OPENMP
+#pragma omp parallel for default(shared) reduction(max:maxVal) reduction(min:minVal)
+#endif
+    for (size_t i = 0; i < input_image.mesh.size(); ++i) {
+        const T val = input_image.mesh[i];
+        if (val > maxVal) maxVal = val;
+        if (val < minVal) minVal = val;
+    }
+
+    return MinMax<T>{minVal, maxVal};
+}
+
+/**
+ * Compute bspline offset for APRConverter of integer type ImageType
+ */
+template<typename ImageType, typename T>
+float compute_bspline_offset(PixelData<T>& input_image, float lambda) {
+    // if bspline smoothing is disabled, there is no need for an offset
+    if(lambda <= 0) return 0;
+
+    // compute offset to center the intensities in the ImageType range (can be negative)
+    auto img_range = getMinMax(input_image);
+    float offset = (std::numeric_limits<ImageType>::max() - (img_range.max - img_range.min)) / 2 - img_range.min;
+
+    // clamp the offset to [-100, 100]
+    return std::max(std::min(offset, 100.f), -100.f);
+}
+
+
 #endif //APR_AUTOPARAMETERS_HPP