diff --git a/.github/workflows/cmake.yml b/.github/workflows/cmake.yml index 43e86178..6c7c3f83 100644 --- a/.github/workflows/cmake.yml +++ b/.github/workflows/cmake.yml @@ -37,7 +37,7 @@ jobs: VCPKG_ROOT: ${{ github.workspace }}/vcpkg steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 with: submodules: true @@ -45,7 +45,7 @@ jobs: - 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,7 +75,7 @@ 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 @@ -85,7 +85,7 @@ jobs: # 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 diff --git a/.gitmodules b/.gitmodules index e08c9430..09d41da6 100644 --- a/.gitmodules +++ b/.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 diff --git a/CMakeLists.txt b/CMakeLists.txt index c4912458..70cf734b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ ############################################################################### # 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}") @@ -9,6 +9,10 @@ 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) @@ -17,7 +21,7 @@ option(APR_BUILD_EXAMPLES "Build APR examples" OFF) option(APR_USE_LIBTIFF "Use LibTIFF" ON) option(APR_TESTS "Build APR tests" OFF) option(APR_PREFER_EXTERNAL_GTEST "When found, use the installed GTEST libs instead of included sources" ON) -option(APR_PREFER_EXTERNAL_BLOSC "When found, use the installed BLOSC libs instead of included sources" ON) +option(APR_PREFER_EXTERNAL_BLOSC "When found, use the installed BLOSC libs instead of included sources" OFF) option(APR_USE_CUDA "should APR use CUDA? (experimental - under development)" OFF) option(APR_USE_OPENMP "should APR use OpenMP?" ON) option(APR_BENCHMARK "add benchmarking code" OFF) @@ -359,13 +363,24 @@ if(APR_TESTS) endif() if(GTEST_FOUND) include_directories(${GTEST_INCLUDE_DIRS}) + message(STATUS "Gtest found: ${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() diff --git a/examples/Example_get_apr.cpp b/examples/Example_get_apr.cpp index bfbaf31d..6ac8e9c7 100644 --- a/examples/Example_get_apr.cpp +++ b/examples/Example_get_apr.cpp @@ -48,7 +48,7 @@ int runAPR(cmdLineOptions options) { //the apr datastructure APR apr; - APRConverter aprConverter; + APRConverter aprConverter; //read in the command line options into the parameters file aprConverter.par.Ip_th = options.Ip_th; diff --git a/examples/Example_get_apr_by_block.cpp b/examples/Example_get_apr_by_block.cpp index 70805a4c..10217a56 100644 --- a/examples/Example_get_apr_by_block.cpp +++ b/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 aprConverter; + APRConverterBatch 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; } \ No newline at end of file diff --git a/examples/Example_get_apr_by_block.hpp b/examples/Example_get_apr_by_block.hpp index 019cddc2..1902f00b 100644 --- a/examples/Example_get_apr_by_block.hpp +++ b/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 diff --git a/examples/Example_reconstruct_image.cpp b/examples/Example_reconstruct_image.cpp index dcebde0e..f821d3c4 100644 --- a/examples/Example_reconstruct_image.cpp +++ b/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 -void add_random_to_img(PixelData& img,float sd){ +void add_random_to_img(PixelData& img, float sd){ std::default_random_engine generator; std::normal_distribution 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 recon_pc; + if(options.output_pc_recon) { + //create mesh data structure for reconstruction + PixelData 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); } ////////////////////////// diff --git a/external/c-blosc b/external/c-blosc index 5352508a..77be44a5 160000 --- a/external/c-blosc +++ b/external/c-blosc @@ -1 +1 @@ -Subproject commit 5352508a420bec865c57cda116a5e5303df898d2 +Subproject commit 77be44a58da10b2261c346b514878ba4d5c020e4 diff --git a/external/gtest b/external/gtest deleted file mode 160000 index 6b74da47..00000000 --- a/external/gtest +++ /dev/null @@ -1 +0,0 @@ -Subproject commit 6b74da4757a549563d7c37c8fae3e704662a043b diff --git a/src/algorithm/APRConverter.hpp b/src/algorithm/APRConverter.hpp index d6728f5f..b34e1b74 100644 --- a/src/algorithm/APRConverter.hpp +++ b/src/algorithm/APRConverter.hpp @@ -99,6 +99,9 @@ class APRConverter { return true; } + float bspline_offset = 0; + + protected: template @@ -108,8 +111,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 grad_temp; // should be a down-sampled image @@ -133,26 +134,6 @@ class APRConverter { }; -template -struct MinMax{T min; T max; }; - -template -static MinMax getMinMax(const PixelData& input_image) { - T minVal = std::numeric_limits::max(); - T maxVal = std::numeric_limits::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{minVal, maxVal}; -} - template void APRConverter::initPipelineMemory(int y_num,int x_num,int z_num){ //initializes the internal memory to be used in the pipeline. @@ -224,7 +205,7 @@ void APRConverter::computeL(APR& aAPR,PixelData& 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 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 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 +214,12 @@ void APRConverter::computeL(APR& aAPR,PixelData& 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::value) { - bspline_offset = 100; - image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); }); - } else if (std::is_same::value){ - bspline_offset = 5; - image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); }); - } else { + if (std::is_floating_point::value) { image_temp.copyFromMesh(input_image); + } else { + bspline_offset = compute_bspline_offset(input_image, par.lambda); + image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); }); } fine_grained_timer.stop_timer(); @@ -276,6 +252,8 @@ void APRConverter::applyParameters(APR& aAPR,APRParameters& aprParame // Apply the main parameters // + aprParameters.validate_parameters(); + fine_grained_timer.start_timer("load_and_apply_mask"); // Apply mask if given if(par.mask_file != ""){ diff --git a/src/algorithm/APRConverterBatch.hpp b/src/algorithm/APRConverterBatch.hpp index 3d611a9b..eab034d3 100644 --- a/src/algorithm/APRConverterBatch.hpp +++ b/src/algorithm/APRConverterBatch.hpp @@ -214,17 +214,14 @@ bool APRConverterBatch::get_apr_method_patch(APR &aAPR, PixelData& 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::value) { - bspline_offset = 100; - image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); }); - } else if (std::is_same::value){ - bspline_offset = 5; - image_temp.copyFromMeshWithUnaryOp(input_image, [=](const auto &a) { return (a + bspline_offset); }); - } else { + + if (std::is_floating_point::value) { image_temp.copyFromMesh(input_image); + } else { + bspline_offset = compute_bspline_offset(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"); @@ -376,6 +373,8 @@ template void APRConverterBatch::applyParameters(PixelData &grad_temp, PixelData &local_scale_temp, PixelData &local_scale_temp2, APRParameters& aprParameters) { + aprParameters.validate_parameters(); + fine_grained_timer.start_timer("load_and_apply_mask"); // Apply mask if given if(par.mask_file != ""){ diff --git a/src/algorithm/APRParameters.hpp b/src/algorithm/APRParameters.hpp index 9ab53eb3..f99a151b 100644 --- a/src/algorithm/APRParameters.hpp +++ b/src/algorithm/APRParameters.hpp @@ -61,6 +61,12 @@ class APRParameters { return os; } + + void validate_parameters(){ + if (sigma_th == 0){ + std::cerr << "Warning: sigma_th is set to 0, this may result in unexpected results due to divide by zero errors. Consider setting this to a non-zero small value, if it is not needed." << std::endl; + } + } }; diff --git a/src/algorithm/AutoParameters.hpp b/src/algorithm/AutoParameters.hpp index 948d0636..14d87cd1 100644 --- a/src/algorithm/AutoParameters.hpp +++ b/src/algorithm/AutoParameters.hpp @@ -157,4 +157,43 @@ void autoParametersLiEntropy(APRParameters& par, } } + + +template +struct MinMax{T min; T max; }; + +template +static MinMax getMinMax(const PixelData& input_image) { + T minVal = std::numeric_limits::max(); + T maxVal = std::numeric_limits::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{minVal, maxVal}; +} + +/** + * Compute bspline offset for APRConverter of integer type ImageType + */ +template +float compute_bspline_offset(PixelData& 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::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 diff --git a/src/data_structures/APR/APR.hpp b/src/data_structures/APR/APR.hpp index 83ae95b3..b2b34fd5 100644 --- a/src/data_structures/APR/APR.hpp +++ b/src/data_structures/APR/APR.hpp @@ -23,7 +23,7 @@ class APR { friend class APRConverterBatch; friend class APRBenchHelper; -protected: +public: // initialize tree RandomAccess void initialize_tree_random_sparse(); @@ -60,8 +60,6 @@ class APR { APRParameters parameters; // this is here to keep a record of what parameters were used, to then be written if needed. -public: - #ifdef APR_USE_CUDA @@ -86,6 +84,10 @@ class APR { * @param with_tree include the tree access */ void init_cuda(bool with_tree=true) { + gpuAccess.genInfo = &aprInfo; + gpuTreeAccess.genInfo = &treeInfo; + linearAccess.genInfo = &aprInfo; + linearAccessTree.genInfo = &treeInfo; auto apr_helper = gpuAPRHelper(); if(with_tree) { auto tree_helper = gpuTreeHelper(); @@ -190,6 +192,7 @@ class APR { tree_initialized = apr2copy.tree_initialized; apr_initialized = apr2copy.apr_initialized; name = apr2copy.name; + parameters = apr2copy.parameters; //old data structures apr_access = apr2copy.apr_access; diff --git a/src/data_structures/APR/access/GPUAccess.hpp b/src/data_structures/APR/access/GPUAccess.hpp index 9d57b09d..4d5f676a 100644 --- a/src/data_structures/APR/access/GPUAccess.hpp +++ b/src/data_structures/APR/access/GPUAccess.hpp @@ -67,7 +67,6 @@ class GPUAccessHelper { gpuAccess->init_y_vec(linearAccess->y_vec); gpuAccess->init_level_xz_vec(linearAccess->level_xz_vec); gpuAccess->init_xz_end_vec(linearAccess->xz_end_vec); - gpuAccess->genInfo = linearAccess->genInfo; gpuAccess->copy2Device(); gpuAccess->initialized = true; } @@ -78,7 +77,6 @@ class GPUAccessHelper { gpuAccess->init_y_vec(linearAccess->y_vec); gpuAccess->init_level_xz_vec(linearAccess->level_xz_vec); gpuAccess->init_xz_end_vec(linearAccess->xz_end_vec); - gpuAccess->genInfo = linearAccess->genInfo; gpuAccess->copy2Device(total_number_particles(tree_access.level_max()), tree_access.gpuAccess); gpuAccess->initialized = true; } @@ -88,7 +86,7 @@ class GPUAccessHelper { gpuAccess->copy2Host(); } - uint64_t total_number_particles() { return gpuAccess->total_number_particles(); } + uint64_t total_number_particles() { return gpuAccess->genInfo->total_number_particles; } uint64_t total_number_particles(const int level) { uint64_t index = linearAccess->level_xz_vec[level] + linearAccess->x_num(level) - 1 + (linearAccess->z_num(level)-1)*linearAccess->x_num(level); diff --git a/src/io/APRWriter.hpp b/src/io/APRWriter.hpp index 4ad6231a..2fbe4dee 100644 --- a/src/io/APRWriter.hpp +++ b/src/io/APRWriter.hpp @@ -702,7 +702,8 @@ class APRWriter { break; case Operation::WRITE: - fileId = hdf5_create_file_blosc(aFileName); +// fileId = hdf5_create_file_blosc(aFileName); + fileId = H5Fcreate(aFileName.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); if (fileId == -1) { std::cerr << "Could not create file [" << aFileName << "]" << std::endl; diff --git a/src/numerics/APRDownsampleGPU.cu b/src/numerics/APRDownsampleGPU.cu index 909728d4..9da3410e 100644 --- a/src/numerics/APRDownsampleGPU.cu +++ b/src/numerics/APRDownsampleGPU.cu @@ -900,6 +900,7 @@ template void compute_ne_rows_tree_cuda(GPUAccessHelper& tree_access, VectorData& ne_count, ScopedCudaMemHandler& ne_rows_gpu) { ne_count.resize(tree_access.level_max() + 3); + ne_count[0] = 0; int z_blocks_max = (tree_access.z_num(tree_access.level_max()) + blockSize_z - 1) / blockSize_z; int num_levels = tree_access.level_max() - tree_access.level_min() + 1; @@ -973,12 +974,13 @@ void compute_ne_rows_tree_cuda(GPUAccessHelper& tree_access, VectorData& ne ne_rows_gpu.get()); } - error_check(cudaFree(block_sums_device) ) + error_check(cudaFree(block_sums_device)) } void compute_ne_rows_tree(GPUAccessHelper& tree_access, VectorData& ne_counter, VectorData& ne_rows) { ne_counter.resize(tree_access.level_max() + 3); + ne_counter[0] = 0; int z = 0; int x = 0; diff --git a/src/numerics/miscCuda.cu b/src/numerics/miscCuda.cu index 4dc043b6..93b5b94e 100644 --- a/src/numerics/miscCuda.cu +++ b/src/numerics/miscCuda.cu @@ -237,6 +237,7 @@ __global__ void fill_ne_rows_cuda(const uint64_t* level_xz_vec, } + template void compute_ne_rows_cuda(GPUAccessHelper& access, VectorData& ne_count, ScopedCudaMemHandler& ne_rows_gpu, int blockSize) { @@ -264,12 +265,12 @@ void compute_ne_rows_cuda(GPUAccessHelper& access, VectorData& ne_count, Sc count_ne_rows_cuda << < grid_dim, block_dim >> > (access.get_level_xz_vec_ptr(), - access.get_xz_end_vec_ptr(), - access.z_num(level), - access.x_num(level), - level, - blockSize, - block_sums_device + offset); + access.get_xz_end_vec_ptr(), + access.z_num(level), + access.x_num(level), + level, + blockSize, + block_sums_device + offset); offset += z_blocks_max; } @@ -305,14 +306,14 @@ void compute_ne_rows_cuda(GPUAccessHelper& access, VectorData& ne_count, Sc fill_ne_rows_cuda<<< grid_dim, block_dim >>> (access.get_level_xz_vec_ptr(), - access.get_xz_end_vec_ptr(), - access.z_num(level), - access.x_num(level), - level, - blockSize, - ne_sz, - ne_count[level], - ne_rows_gpu.get()); + access.get_xz_end_vec_ptr(), + access.z_num(level), + access.x_num(level), + level, + blockSize, + ne_sz, + ne_count[level], + ne_rows_gpu.get()); } error_check( cudaFree(block_sums_device) ) diff --git a/test/APRParametersTest.cpp b/test/APRParametersTest.cpp new file mode 100644 index 00000000..6fc37a04 --- /dev/null +++ b/test/APRParametersTest.cpp @@ -0,0 +1,40 @@ +#include +#include +#include +#include +#include "algorithm/APRParameters.hpp" + +class APRParametersTest : public ::testing::Test { +protected: + // This can be used to install a custom std::cerr buffer + std::streambuf* original_cerr; + + std::stringstream cerr_content; + + void SetUp() override { + // Before each test, backup the stream buffer for std::cerr + original_cerr = std::cerr.rdbuf(); + // Redirect std::cerr to our stringstream buffer or any other stream + std::cerr.rdbuf(cerr_content.rdbuf()); + } + + void TearDown() override { + // Restore the original buffer before exiting the test + std::cerr.rdbuf(original_cerr); + } + + // Helper function to check if the warning message is in the buffer + bool hasWarningMessage(const std::string& message) { + return cerr_content.str().find(message) != std::string::npos; + } +}; + +TEST_F(APRParametersTest, ValidateParametersWarnsOnZeroSigmaTh) { + APRParameters params; + params.sigma_th = 0; // Set sigma_th to zero to trigger the warning + + params.validate_parameters(); // Call the validation method + + // Check if the correct warning message was logged to std::cerr + EXPECT_TRUE(hasWarningMessage("Warning: sigma_th is set to 0")); +} \ No newline at end of file diff --git a/test/APRTest.cpp b/test/APRTest.cpp index 9f2d990e..33ea37d6 100644 --- a/test/APRTest.cpp +++ b/test/APRTest.cpp @@ -123,106 +123,62 @@ class CreateGTSmall1DTest : public CreateAPRTest }; template -bool compare_two_iterators(Iterator1& it1, Iterator2& it2,bool success = true){ +bool compare_two_iterators(Iterator1& it1, Iterator2& it2, int maxNumOfErrPrinted = 10){ - if(it1.total_number_particles() != it2.total_number_particles()){ - success = false; - std::cout << "Number of particles mismatch" << std::endl; + if(it1.total_number_particles() != it2.total_number_particles()) { + std::cout << "Number of particles mismatch: " << it1.total_number_particles() << " vs " << it2.total_number_particles() << std::endl; + return false; } uint64_t counter_1 = 0; uint64_t counter_2 = 0; - for (int level = it1.level_min(); level <= it1.level_max(); ++level) { - int z = 0; - int x = 0; + uint64_t errors = 0; - for (z = 0; z < it1.z_num(level); z++) { - for (x = 0; x < it1.x_num(level); ++x) { + for (int level = it1.level_min(); level <= it1.level_max(); ++level) { + for (int z = 0; z < it1.z_num(level); z++) { + for (int x = 0; x < it1.x_num(level); ++x) { it2.begin(level, z, x); - for (it1.begin(level, z, x); it1 < it1.end(); - it1++) { + for (it1.begin(level, z, x); it1 < it1.end(); it1++) { counter_1++; - if(it1 != it1){ - - uint64_t new_index = it1; - uint64_t org_index = it2; - - (void) new_index; - (void) org_index; - - - success = false; -// std::cout << "1" << std::endl; - } - - if(it1.y() != it2.y()){ - - auto y_new = it1.y(); - auto y_org = it2.y(); - - (void) y_new; - (void) y_org; - - success = false; -// std::cout << "y_new" << y_new << std::endl; -// std::cout << "y_org" << y_org << std::endl; + if( (it1 != it2) || (it1.y() != it2.y()) ){ + if(errors < maxNumOfErrPrinted || maxNumOfErrPrinted == -1) { + std::cout << "iterator mismatch, idx " << (uint64_t) it1 << " vs " << (uint64_t) it2 << ", y value " << it2.y() << " vs " << it1.y() << std::endl; + } + errors++; } if(it2 < it2.end()){ it2++; } - } } } } - - for (int level = it2.level_min(); level <= it2.level_max(); ++level) { - int z = 0; - int x = 0; - - for (z = 0; z < it2.z_num(level); z++) { - for (x = 0; x < it2.x_num(level); ++x) { + for (int z = 0; z < it2.z_num(level); z++) { + for (int x = 0; x < it2.x_num(level); ++x) { it1.begin(level, z, x); - for (it2.begin(level, z, x); it2 < it2.end(); - it2++) { + for (it2.begin(level, z, x); it2 < it2.end(); it2++) { counter_2++; - if(it1 != it1){ - - uint64_t new_index = it1; - uint64_t org_index = it2; - - (void) new_index; - (void) org_index; - - - success = false; - } - - if(it1.y() != it2.y()){ - - auto y_new = it1.y(); - auto y_org = it2.y(); - - (void) y_new; - (void) y_org; - - success = false; + if( (it1 != it2) || (it1.y() != it2.y()) ){ + if(errors < maxNumOfErrPrinted || maxNumOfErrPrinted == -1) { + std::cout << "iterator mismatch, idx " << (uint64_t) it1 << " vs " << (uint64_t) it2 << " y value " << it1.y() << " vs " << it2.y() << std::endl; + } + errors++; } - if(it1 < it1.end()){ + if(it1 < it1.end()) { it1++; } @@ -231,13 +187,12 @@ bool compare_two_iterators(Iterator1& it1, Iterator2& it2,bool success = true){ } } - if(counter_1 != counter_2){ - success = false; + if((counter_1 != counter_2) || (errors > 0)){ std::cout << "Iteration mismatch" << std::endl; + return false; } - - return success; + return true; } bool check_neighbours(APR& apr,APRIterator ¤t, APRIterator &neigh){ @@ -870,9 +825,9 @@ bool test_pulling_scheme_sparse(TestData& test_data){ auto sparse_lin_it = apr_lin_sparse.iterator(); - success = compare_two_iterators(org_it,sparse_it,success); - success = compare_two_iterators(sparse_lin_it,sparse_it,success); - success = compare_two_iterators(org_it,sparse_lin_it,success); + success = success && compare_two_iterators(org_it, sparse_it); + success = success && compare_two_iterators(sparse_lin_it, sparse_it); + success = success && compare_two_iterators(org_it, sparse_lin_it); return success; } @@ -973,12 +928,12 @@ bool test_linear_access_create(TestData& test_data) { //apr_lin.init_linear(); auto it_new = apr_lin.iterator(); - success = compare_two_iterators(it_lin_old,it_new,success); + success = success && compare_two_iterators(it_lin_old, it_new); //Test Linear -> Random generation auto it_new_random = apr_lin.random_iterator(); - success = compare_two_iterators(it_new_random,it_new,success); + success = success && compare_two_iterators(it_new_random, it_new); //Test the APR Tree construction. @@ -992,7 +947,7 @@ bool test_linear_access_create(TestData& test_data) { std::cout << "PARTS: " << total_number_parts << " " << total_number_parts_lin << std::endl; - success = compare_two_iterators(tree_it_org,tree_it_lin,success); + success = success && compare_two_iterators(tree_it_org, tree_it_lin); return success; @@ -1298,16 +1253,16 @@ bool test_linear_access_io(TestData& test_data) { auto it_new = apr_lin.iterator(); auto it_read = test_data.apr.iterator(); - success = compare_two_iterators(it_org,it_new,success); - success = compare_two_iterators(it_org,it_read,success); - success = compare_two_iterators(it_new,it_read,success); + success = success && compare_two_iterators(it_org, it_new); + success = success && compare_two_iterators(it_org, it_read); + success = success && compare_two_iterators(it_new, it_read); // Test the tree IO auto tree_it_org = apr_random.random_tree_iterator(); auto tree_it_lin = apr_lin.tree_iterator(); - success = compare_two_iterators(tree_it_org,tree_it_lin,success); + success = success && compare_two_iterators(tree_it_org, tree_it_lin); return success; } @@ -1325,11 +1280,11 @@ bool test_apr_tree(TestData& test_data) { auto it_lin = test_data.apr.iterator(); auto it_random = test_data.apr.random_iterator(); - success = compare_two_iterators(it_lin,it_random,success); + success = success && compare_two_iterators(it_lin, it_random); auto it_tree_t = test_data.apr.random_tree_iterator(); - success = compare_two_iterators(it_lin,it_random,success); + success = success && compare_two_iterators(it_lin, it_random); ParticleData tree_data; @@ -1338,7 +1293,7 @@ bool test_apr_tree(TestData& test_data) { auto apr_tree_iterator = test_data.apr.random_tree_iterator(); auto apr_tree_iterator_lin = test_data.apr.tree_iterator(); - success = compare_two_iterators(apr_tree_iterator,apr_tree_iterator_lin,success); + success = success && compare_two_iterators(apr_tree_iterator, apr_tree_iterator_lin); for (int level = (apr_tree_iterator.level_max()); level >= apr_tree_iterator.level_min(); --level) { @@ -1620,7 +1575,7 @@ bool test_apr_file(TestData& test_data){ auto apr_iterator = test_data.apr.iterator(); auto apr_iterator_read = aprRead.iterator(); - success = compare_two_iterators(apr_iterator,apr_iterator_read,success); + success = success && compare_two_iterators(apr_iterator, apr_iterator_read); //test apr iterator with channel writeFile.open(file_name,"WRITE"); @@ -1662,11 +1617,11 @@ bool test_apr_file(TestData& test_data){ auto it1 = apr_channel_0.iterator(); auto it2 = test_data.apr.iterator(); - success = compare_two_iterators(it1,it2,success); + success = success && compare_two_iterators(it1, it2); it1 = apr_channel_0_55.iterator(); it2 = test_data.apr.iterator(); - success = compare_two_iterators(it1,it2,success); + success = success && compare_two_iterators(it1, it2); @@ -1710,7 +1665,7 @@ bool test_apr_file(TestData& test_data){ auto tree_it = aprRead2.random_tree_iterator(); auto tree_it_org = test_data.apr.random_tree_iterator(); - success = compare_two_iterators(tree_it,tree_it_org,success); + success = success && compare_two_iterators(tree_it, tree_it_org); //Test file list std::vector correct_names = {"tree_parts","tree_parts1","tree_parts2"}; @@ -2819,8 +2774,6 @@ bool test_pipeline_u16(TestData& test_data){ // // - bool success = true; - //the apr datastructure APR apr; APRConverter aprConverter; @@ -2831,12 +2784,9 @@ bool test_pipeline_u16(TestData& test_data){ aprConverter.par.Ip_th = 0; aprConverter.par.rel_error = readPars.rel_error; aprConverter.par.lambda = readPars.lambda; - - aprConverter.par.sigma_th_max = readPars.sigma_th_max; - aprConverter.par.sigma_th = readPars.sigma_th; - + aprConverter.par.sigma_th_max = std::fmax(readPars.sigma_th_max,0.1); + aprConverter.par.sigma_th = std::fmax(readPars.sigma_th,0.1); aprConverter.par.grad_th = readPars.grad_th; - aprConverter.par.auto_parameters = false; //where things are @@ -2855,34 +2805,24 @@ bool test_pipeline_u16(TestData& test_data){ PixelData scale_saved = TiffUtils::getMesh(test_data.output_dir + "scale_saved.tif"); PixelData gradient_saved = TiffUtils::getMesh(test_data.output_dir + "gradient_saved.tif"); - for (size_t i = 0; i < scale_computed.mesh.size(); ++i) { - float computed_val = scale_computed.mesh[i]; - float saved_val = scale_saved.mesh[i]; + std::cout << "Comparing local intensity scales" << std::endl; + int scale_errors = compareMeshes(scale_saved, scale_computed, 1); - if(std::abs(computed_val - saved_val) > 1){ - success = false; - } - } - - for (size_t i = 0; i < gradient_computed.mesh.size(); ++i) { - float computed_val = gradient_computed.mesh[i]; - float saved_val = gradient_saved.mesh[i]; + std::cout << "Comparing gradients" << std::endl; + int grad_errors = compareMeshes(gradient_saved, gradient_computed, 1); - if(std::abs(computed_val - saved_val) > 1){ - success = false; - } - } + bool success = (scale_errors == 0) && (grad_errors == 0); APR apr_c; + aprConverter.bspline_offset = 0; aprConverter.initPipelineAPR(apr_c, test_data.img_original); - - aprConverter.get_apr_custom_grad_scale(apr_c,gradient_saved,scale_saved); + aprConverter.get_apr_custom_grad_scale(apr_c, gradient_saved, scale_saved); auto it_org = test_data.apr.iterator(); auto it_gen = apr_c.iterator(); //test the access - success = compare_two_iterators(it_org,it_gen,success); + success = success && compare_two_iterators(it_org, it_gen); ParticleData particles_intensities; @@ -3423,16 +3363,12 @@ bool test_iterator_methods(TestData &test_data){ bool test_apr_copy(TestData &test_data){ - bool success = true; - APR aprCopy(test_data.apr); auto it_org = test_data.apr.iterator(); auto it_copy = aprCopy.iterator(); - success = compare_two_iterators(it_org,it_copy,success); - - return success; + return compare_two_iterators(it_org, it_copy); } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7bc8f6cc..dc3e5a11 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -10,6 +10,7 @@ buildTarget(testAPR APRTest.cpp) buildTarget(testComputeGradient ComputeGradientTest.cpp) buildTarget(testLocalIntensityScale LocalIntensityScaleTest.cpp) buildTarget(testPullingScheme PullingSchemeTest.cpp) +buildTarget(testAPRParameters APRParametersTest.cpp) #APR GPU Tests if(APR_USE_CUDA) diff --git a/vcpkg b/vcpkg index d72783cb..486a4640 160000 --- a/vcpkg +++ b/vcpkg @@ -1 +1 @@ -Subproject commit d72783cb3aeddfd667861caef1060e54ca6fa7a9 +Subproject commit 486a4640db740f5994e492eb60748111dfc48de7 diff --git a/vcpkg.json b/vcpkg.json index f7b8b7f0..193ac9b0 100644 --- a/vcpkg.json +++ b/vcpkg.json @@ -1,11 +1,11 @@ { "name": "libapr", - "version-string": "0.0.1", + "version-string": "2.1.0", "dependencies": [ "blosc", "hdf5", "szip", - "gtest", "tiff" - ] -} \ No newline at end of file + ], + "builtin-baseline": "486a4640db740f5994e492eb60748111dfc48de7" +}