diff --git a/include/libtomo/stripe.h b/include/libtomo/stripe.h index 73c52b02d..4642388c6 100644 --- a/include/libtomo/stripe.h +++ b/include/libtomo/stripe.h @@ -51,5 +51,17 @@ # define DLL #endif +#include + DLL void remove_stripe_sf(float* data, int dx, int dy, int dz, int size, int istart, int iend); + +DLL int +stripesdetect3d_main_float(float* Input, float* Output, int window_halflength_vertical, + int ratio_radius, int ncores, int dimX, int dimY, int dimZ); + +DLL int +stripesmask3d_main_float(float* Input, bool* Output, float threshold_val, + int stripe_length_min, int stripe_depth_min, + int stripe_width_min, float sensitivity, int ncores, int dimX, + int dimY, int dimZ); diff --git a/source/libtomo/misc/CMakeLists.txt b/source/libtomo/misc/CMakeLists.txt index 18c2ea516..a3c6431dd 100755 --- a/source/libtomo/misc/CMakeLists.txt +++ b/source/libtomo/misc/CMakeLists.txt @@ -9,8 +9,6 @@ tomopy_add_library( morph.c remove_ring.c median_filt3d.c - utils.c - utils.h ${HEADERS}) find_package(OpenMP REQUIRED COMPONENTS C) diff --git a/source/libtomo/misc/median_filt3d.c b/source/libtomo/misc/median_filt3d.c index 024d632d7..0b1aadfeb 100644 --- a/source/libtomo/misc/median_filt3d.c +++ b/source/libtomo/misc/median_filt3d.c @@ -47,9 +47,22 @@ #include #include #include +#include #include "libtomo/median_filt3d.h" -#include "utils.h" + +int floatcomp(const void* elem1, const void* elem2) +{ + if(*(const float*)elem1 < *(const float*)elem2) + return -1; + return *(const float*)elem1 > *(const float*)elem2; +} +int uint16comp(const void* elem1, const void* elem2) +{ + if(*(const unsigned short*)elem1 < *(const unsigned short*)elem2) + return -1; + return *(const unsigned short*)elem1 > *(const unsigned short*)elem2; +} void medfilt3D_float(float* Input, float* Output, int radius, int sizefilter_total, @@ -90,7 +103,8 @@ medfilt3D_float(float* Input, float* Output, int radius, int sizefilter_total, } } } - quicksort_float(ValVec, 0, sizefilter_total - 1); /* perform sorting */ + + qsort(ValVec, sizefilter_total, sizeof(float), floatcomp); if(mu_threshold == 0.0F) { @@ -136,7 +150,7 @@ medfilt2D_float(float* Input, float* Output, int radius, int sizefilter_total, counter++; } } - quicksort_float(ValVec, 0, sizefilter_total - 1); /* perform sorting */ + qsort(ValVec, sizefilter_total, sizeof(float), floatcomp); if(mu_threshold == 0.0F) { @@ -191,7 +205,8 @@ medfilt3D_uint16(unsigned short* Input, unsigned short* Output, int radius, } } } - quicksort_uint16(ValVec, 0, sizefilter_total - 1); /* perform sorting */ + + qsort(ValVec, sizefilter_total, sizeof(unsigned short), uint16comp); if(mu_threshold == 0.0F) { @@ -238,7 +253,8 @@ medfilt2D_uint16(unsigned short* Input, unsigned short* Output, int radius, counter++; } } - quicksort_uint16(ValVec, 0, sizefilter_total - 1); /* perform sorting */ + + qsort(ValVec, sizefilter_total, sizeof(unsigned short), uint16comp); if(mu_threshold == 0.0F) { @@ -268,7 +284,7 @@ medianfilter_main_float(float* Input, float* Output, int radius, float mu_thresh diameter = (2 * radius + 1); /* diameter of the filter's kernel */ if(mu_threshold != 0.0) { - copyIm(Input, Output, (long) (dimX), (long) (dimY), (long) (dimZ)); + memcpy(Output, Input, dimX * dimY * dimZ * sizeof(float)); } /* copy input into output */ /* dealing here with a custom given number of cpu threads */ @@ -331,7 +347,7 @@ medianfilter_main_uint16(unsigned short* Input, unsigned short* Output, int radi diameter = (2 * radius + 1); /* diameter of the filter's kernel */ if(mu_threshold != 0.0) { - copyIm_unshort(Input, Output, (long) (dimX), (long) (dimY), (long) (dimZ)); + memcpy(Output, Input, dimX * dimY * dimZ * sizeof(unsigned short)); } /* copy input into output */ /* dealing here with a custom given number of cpu threads */ diff --git a/source/libtomo/misc/utils.c b/source/libtomo/misc/utils.c deleted file mode 100644 index a598163d0..000000000 --- a/source/libtomo/misc/utils.c +++ /dev/null @@ -1,202 +0,0 @@ -// Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved. - -// Copyright 2015. UChicago Argonne, LLC. This software was produced -// under U.S. Government contract DE-AC02-06CH11357 for Argonne National -// Laboratory (ANL), which is operated by UChicago Argonne, LLC for the -// U.S. Department of Energy. The U.S. Government has rights to use, -// reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR -// UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR -// ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is -// modified to produce derivative works, such modified software should -// be clearly marked, so as not to confuse it with the version available -// from ANL. - -// Additionally, redistribution and use in source and binary forms, with -// or without modification, are permitted provided that the following -// conditions are met: - -// * Redistributions of source code must retain the above copyright -// notice, this list of conditions and the following disclaimer. - -// * Redistributions in binary form must reproduce the above copyright -// notice, this list of conditions and the following disclaimer in -// the documentation and/or other materials provided with the -// distribution. - -// * Neither the name of UChicago Argonne, LLC, Argonne National -// Laboratory, ANL, the U.S. Government, nor the names of its -// contributors may be used to endorse or promote products derived -// from this software without specific prior written permission. - -// THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS -// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago -// Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, -// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, -// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT -// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN -// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -// POSSIBILITY OF SUCH DAMAGE. - -#include "utils.h" - -// for windows build -#ifdef WIN32 -# ifdef PY3K -void -PyInit_libtomo(void) -{ -} -# else -void -initlibtomo(void) -{ -} -# endif -#endif - -/* Copy Image (float) */ -void -copyIm(const float* A, float* U, long dimX, long dimY, long dimZ) -{ - long j; - for(j = 0; j < dimX * dimY * dimZ; j++) - U[j] = A[j]; -} - -/* Copy Image -unsigned char (8bit)*/ -void -copyIm_unchar(const unsigned char* A, unsigned char* U, int dimX, int dimY, int dimZ) -{ - int j; - for(j = 0; j < dimX * dimY * dimZ; j++) - U[j] = A[j]; -} - -/* Copy Image - unsigned short (16bit)*/ -void -copyIm_unshort(const unsigned short* A, unsigned short* U, int dimX, int dimY, int dimZ) -{ - int j; - for(j = 0; j < dimX * dimY * dimZ; j++) - U[j] = A[j]; -} - -/* sorting using bubble method (float)*/ -void -sort_bubble_float(float* x, int n_size) -{ - int i; - int j; - float temp; - - for(i = 0; i < n_size - 1; i++) - { - for(j = 0; j < n_size - i - 1; j++) - { - if(x[j] > x[j + 1]) - { - temp = x[j]; - x[j] = x[j + 1]; - x[j + 1] = temp; - } - } - } -} - -/* sorting using bubble method (uint16)*/ -void -sort_bubble_uint16(unsigned short* x, int n_size) -{ - int i; - int j; - unsigned short temp; - - for(i = 0; i < n_size - 1; i++) - { - for(j = 0; j < n_size - i - 1; j++) - { - if(x[j] > x[j + 1]) - { - temp = x[j]; - x[j] = x[j + 1]; - x[j + 1] = temp; - } - } - } -} - -void -quicksort_float(float* x, int first, int last) -{ - int i; - int j; - int pivot; - float temp; - - if(first < last) - { - pivot = first; - i = first; - j = last; - - while(i < j) - { - while(x[i] <= x[pivot] && i < last) - i++; - while(x[j] > x[pivot]) - j--; - if(i < j) - { - temp = x[i]; - x[i] = x[j]; - x[j] = temp; - } - } - - temp = x[pivot]; - x[pivot] = x[j]; - x[j] = temp; - quicksort_float(x, first, j - 1); - quicksort_float(x, j + 1, last); - } -} - -void -quicksort_uint16(unsigned short* x, int first, int last) -{ - int i; - int j; - int pivot; - unsigned short temp; - - if(first < last) - { - pivot = first; - i = first; - j = last; - - while(i < j) - { - while(x[i] <= x[pivot] && i < last) - i++; - while(x[j] > x[pivot]) - j--; - if(i < j) - { - temp = x[i]; - x[i] = x[j]; - x[j] = temp; - } - } - - temp = x[pivot]; - x[pivot] = x[j]; - x[j] = temp; - quicksort_uint16(x, first, j - 1); - quicksort_uint16(x, j + 1, last); - } -} diff --git a/source/libtomo/misc/utils.h b/source/libtomo/misc/utils.h deleted file mode 100644 index 80ea12fb7..000000000 --- a/source/libtomo/misc/utils.h +++ /dev/null @@ -1,65 +0,0 @@ -// Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved. - -// Copyright 2015. UChicago Argonne, LLC. This software was produced -// under U.S. Government contract DE-AC02-06CH11357 for Argonne National -// Laboratory (ANL), which is operated by UChicago Argonne, LLC for the -// U.S. Department of Energy. The U.S. Government has rights to use, -// reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR -// UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR -// ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is -// modified to produce derivative works, such modified software should -// be clearly marked, so as not to confuse it with the version available -// from ANL. - -// Additionally, redistribution and use in source and binary forms, with -// or without modification, are permitted provided that the following -// conditions are met: - -// * Redistributions of source code must retain the above copyright -// notice, this list of conditions and the following disclaimer. - -// * Redistributions in binary form must reproduce the above copyright -// notice, this list of conditions and the following disclaimer in -// the documentation and/or other materials provided with the -// distribution. - -// * Neither the name of UChicago Argonne, LLC, Argonne National -// Laboratory, ANL, the U.S. Government, nor the names of its -// contributors may be used to endorse or promote products derived -// from this software without specific prior written permission. - -// THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS -// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT -// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS -// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago -// Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, -// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, -// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER -// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT -// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN -// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE -// POSSIBILITY OF SUCH DAMAGE. - -#pragma once - -#ifdef WIN32 -# define DLL __declspec(dllexport) -#else -# define DLL -#endif - -void DLL -copyIm(const float* A, float* U, long dimX, long dimY, long dimZ); -void DLL -copyIm_unchar(const unsigned char* A, unsigned char* U, int dimX, int dimY, int dimZ); -void DLL -copyIm_unshort(const unsigned short* A, unsigned short* U, int dimX, int dimY, int dimZ); -void DLL -sort_bubble_float(float* x, int n_size); -void DLL -sort_bubble_uint16(unsigned short* x, int n_size); -void DLL -quicksort_float(float* x, int first, int last); -void DLL -quicksort_uint16(unsigned short* x, int first, int last); diff --git a/source/libtomo/prep/CMakeLists.txt b/source/libtomo/prep/CMakeLists.txt index b946fbc2d..2148ffd6b 100755 --- a/source/libtomo/prep/CMakeLists.txt +++ b/source/libtomo/prep/CMakeLists.txt @@ -1,7 +1,20 @@ set(HEADERS "${tomopy_SOURCE_DIR}/include/libtomo/prep.h" "${tomopy_SOURCE_DIR}/include/libtomo/stripe.h") -tomopy_add_library(tomo-prep SHARED prep.c stripe.c ${HEADERS}) +tomopy_add_library( + tomo-prep + SHARED + prep.c + stripe.c + stripes_detect3d.c + ${HEADERS}) + +find_package(OpenMP REQUIRED COMPONENTS C) +target_link_libraries(tomo-prep PRIVATE OpenMP::OpenMP_C) +if (WIN32) +target_compile_options( + tomo-prep PRIVATE $<$:/openmp:experimental>) +endif() target_include_directories( tomo-prep @@ -14,6 +27,10 @@ target_compile_definitions(tomo-prep PRIVATE ${${PROJECT_NAME}_DEFINITIONS}) target_compile_options( tomo-prep PRIVATE $<$:${${PROJECT_NAME}_C_FLAGS}>) +set_target_properties(tomo-prep PROPERTIES + C_STANDARD 99 # C99 or later requried for boolean type support + C_STANDARD_REQUIRED ON) + install(TARGETS tomo-prep EXPORT libtomoTargets) install( diff --git a/source/libtomo/prep/stripes_detect3d.c b/source/libtomo/prep/stripes_detect3d.c new file mode 100644 index 000000000..b51ec8380 --- /dev/null +++ b/source/libtomo/prep/stripes_detect3d.c @@ -0,0 +1,723 @@ +// Copyright (c) 2015, UChicago Argonne, LLC. All rights reserved. + +// Copyright 2015. UChicago Argonne, LLC. This software was produced +// under U.S. Government contract DE-AC02-06CH11357 for Argonne National +// Laboratory (ANL), which is operated by UChicago Argonne, LLC for the +// U.S. Department of Energy. The U.S. Government has rights to use, +// reproduce, and distribute this software. NEITHER THE GOVERNMENT NOR +// UChicago Argonne, LLC MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR +// ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE. If software is +// modified to produce derivative works, such modified software should +// be clearly marked, so as not to confuse it with the version available +// from ANL. + +// Additionally, redistribution and use in source and binary forms, with +// or without modification, are permitted provided that the following +// conditions are met: + +// * Redistributions of source code must retain the above copyright +// notice, this list of conditions and the following disclaimer. + +// * Redistributions in binary form must reproduce the above copyright +// notice, this list of conditions and the following disclaimer in +// the documentation and/or other materials provided with the +// distribution. + +// * Neither the name of UChicago Argonne, LLC, Argonne National +// Laboratory, ANL, the U.S. Government, nor the names of its +// contributors may be used to endorse or promote products derived +// from this software without specific prior written permission. + +// THIS SOFTWARE IS PROVIDED BY UChicago Argonne, LLC AND CONTRIBUTORS +// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS +// FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL UChicago +// Argonne, LLC OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, +// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, +// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT +// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN +// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +// POSSIBILITY OF SUCH DAMAGE. + +// C-module for detecting and emphasising stripes present in the data (3D case) +// Original author: Daniil Kazantsev, Diamond Light Source Ltd. + +#include +#include +#include +#include +#include +#include + +#include "libtomo/stripe.h" + +/********************************************************************/ +/**********************Supporting Functions**************************/ +/********************************************************************/ +/* Calculate the forward difference derrivative of the 3D input in the +direction of the "axis" parameter using the step_size in pixels to skip pixels +(i.e. step_size = 1 is the classical gradient) +axis = 0: horizontal direction +axis = 1: depth direction +axis = 2: vertical direction +*/ +void +gradient3D_local(const float* input, float* output, long dimX, long dimY, long dimZ, + int axis, int step_size) +{ + long i; + long j; + long k; + long i1; + long j1; + long k1; + size_t index; + +#pragma omp parallel for shared(input, output) private(i, j, k, i1, j1, k1, index) + for(j = 0; j < dimY; j++) + { + for(i = 0; i < dimX; i++) + { + for(k = 0; k < dimZ; k++) + { + index = (size_t) (dimX * dimY * k) + (size_t) (j * dimX + i); + /* Forward differences */ + if(axis == 0) + { + i1 = i + step_size; + if(i1 >= dimX) + i1 = i - step_size; + output[index] = + input[(size_t) (dimX * dimY * k) + (size_t) (j * dimX + i1)] - + input[index]; + } + else if(axis == 1) + { + j1 = j + step_size; + if(j1 >= dimY) + j1 = j - step_size; + output[index] = + input[(size_t) (dimX * dimY * k) + (size_t) (j1 * dimX + i)] - + input[index]; + } + else + { + k1 = k + step_size; + if(k1 >= dimZ) + k1 = k - step_size; + output[index] = + input[(size_t) (dimX * dimY * k1) + (size_t) (j * dimX + i)] - + input[index]; + } + } + } + } +} + +void +ratio_mean_stride3d(float* input, float* output, int radius, long i, long j, long k, + long dimX, long dimY, long dimZ) +{ + float mean_plate; + float mean_horiz; + float mean_horiz2; + float min_val; + int diameter = 2 * radius + 1; + int all_pixels_window = diameter * diameter; + long i_m; + long j_m; + long k_m; + long i1; + long j1; + long k1; + size_t index; + size_t newindex; + + index = (size_t) (dimX * dimY * k) + (size_t) (j * dimX + i); + + min_val = 0.0F; + /* calculate mean of gradientX in a 2D plate parallel to stripes direction */ + mean_plate = 0.0F; + for(j_m = -radius; j_m <= radius; j_m++) + { + j1 = j + j_m; + if((j1 < 0) || (j1 >= dimY)) + j1 = j - j_m; + for(k_m = -radius; k_m <= radius; k_m++) + { + k1 = k + k_m; + if((k1 < 0) || (k1 >= dimZ)) + k1 = k - k_m; + newindex = (size_t) (dimX * dimY * k1) + (size_t) (j1 * dimX + i); + mean_plate += fabsf(input[newindex]); + } + } + mean_plate /= (float) (all_pixels_window); + + /* calculate mean of gradientX in a 2D plate orthogonal to stripes direction */ + mean_horiz = 0.0F; + for(j_m = -1; j_m <= 1; j_m++) + { + j1 = j + j_m; + if((j1 < 0) || (j1 >= dimY)) + j1 = j - j_m; + for(i_m = 1; i_m <= radius; i_m++) + { + i1 = i + i_m; + if(i1 >= dimX) + i1 = i - i_m; + newindex = (size_t) (dimX * dimY * k) + (size_t) (j1 * dimX + i1); + mean_horiz += fabsf(input[newindex]); + } + } + mean_horiz /= (float) (radius * 3); + + /* Calculate another mean symmetrically */ + mean_horiz2 = 0.0F; + for(j_m = -1; j_m <= 1; j_m++) + { + j1 = j + j_m; + if((j1 < 0) || (j1 >= dimY)) + j1 = j - j_m; + for(i_m = -radius; i_m <= -1; i_m++) + { + i1 = i + i_m; + if(i1 < 0) + i1 = i - i_m; + newindex = (size_t) (dimX * dimY * k) + (size_t) (j1 * dimX + i1); + mean_horiz2 += fabsf(input[newindex]); + } + } + mean_horiz2 /= (float) (radius * 3); + + /* calculate the ratio between two means assuming that the mean + orthogonal to stripes direction should be larger than the mean + parallel to it */ + if((mean_horiz >= mean_plate) && (mean_horiz != 0.0F)) + output[index] = mean_plate / mean_horiz; + if((mean_horiz < mean_plate) && (mean_plate != 0.0F)) + output[index] = mean_horiz / mean_plate; + if((mean_horiz2 >= mean_plate) && (mean_horiz2 != 0.0F)) + min_val = mean_plate / mean_horiz2; + if((mean_horiz2 < mean_plate) && (mean_plate != 0.0F)) + min_val = mean_horiz2 / mean_plate; + + /* accepting the smallest value */ + if(output[index] > min_val) + output[index] = min_val; +} + +int +floatcomp(const void* elem1, const void* elem2) +{ + if(*(const float*) elem1 < *(const float*) elem2) + return -1; + return *(const float*) elem1 > *(const float*) elem2; +} + +void +vertical_median_stride3d(const float* input, float* output, + int window_halflength_vertical, int window_fulllength, + int midval_window_index, long i, long j, long k, long dimX, + long dimY, long dimZ) +{ + int counter; + long k_m; + long k1; + size_t index; + + index = (size_t) (dimX * dimY * k) + (size_t) (j * dimX + i); + + float* _values; + _values = (float*) calloc(window_fulllength, sizeof(float)); + + counter = 0; + for(k_m = -window_halflength_vertical; k_m <= window_halflength_vertical; k_m++) + { + k1 = k + k_m; + if((k1 < 0) || (k1 >= dimZ)) + k1 = k - k_m; + _values[counter] = input[(size_t) (dimX * dimY * k1) + (size_t) (j * dimX + i)]; + counter++; + } + qsort(_values, window_fulllength, sizeof(float), floatcomp); + output[index] = _values[midval_window_index]; + + free(_values); +} + +void +mean_stride3d(const float* input, float* output, long i, long j, long k, long dimX, + long dimY, long dimZ) +{ + /* a 3d mean to enusre a more stable gradient */ + long i1; + long i2; + long j1; + long j2; + long k1; + long k2; + float val1; + float val2; + float val3; + float val4; + float val5; + float val6; + size_t index; + + index = (size_t) (dimX * dimY * k) + (size_t) (j * dimX + i); + + i1 = i - 1; + i2 = i + 1; + j1 = j - 1; + j2 = j + 1; + k1 = k - 1; + k2 = k + 1; + + if(i1 < 0) + i1 = i2; + if(i2 >= dimX) + i2 = i1; + if(j1 < 0) + j1 = j2; + if(j2 >= dimY) + j2 = j1; + if(k1 < 0) + k1 = k2; + if(k2 >= dimZ) + k2 = k1; + + val1 = input[(size_t) (dimX * dimY * k) + (size_t) (j * dimX + i1)]; + val2 = input[(size_t) (dimX * dimY * k) + (size_t) (j * dimX + i2)]; + val3 = input[(size_t) (dimX * dimY * k) + (size_t) (j1 * dimX + i)]; + val4 = input[(size_t) (dimX * dimY * k) + (size_t) (j2 * dimX + i)]; + val5 = input[(size_t) (dimX * dimY * k1) + (size_t) (j * dimX + i)]; + val6 = input[(size_t) (dimX * dimY * k2) + (size_t) (j * dimX + i)]; + + output[index] = 0.1428F * (input[index] + val1 + val2 + val3 + val4 + val5 + val6); +} + +void +remove_inconsistent_stripes(const bool* mask, bool* out, int stripe_length_min, + int stripe_depth_min, float sensitivity, int switch_dim, + long i, long j, long k, long dimX, long dimY, long dimZ) +{ + int counter_vert_voxels; + int counter_depth_voxels; + int halfstripe_length = stripe_length_min / 2; + int halfstripe_depth = stripe_depth_min / 2; + long k_m; + long k1; + long y_m; + long y1; + size_t index; + index = (size_t) (dimX * dimY * k) + (size_t) (j * dimX + i); + + int threshold_vertical = (int) ((0.01F * sensitivity) * stripe_length_min); + int threshold_depth = (int) ((0.01F * sensitivity) * stripe_depth_min); + + /* start by considering vertical features */ + if(switch_dim == 0) + { + if(mask[index] == true) + { + counter_vert_voxels = 0; + for(k_m = -halfstripe_length; k_m <= halfstripe_length; k_m++) + { + k1 = k + k_m; + if((k1 < 0) || (k1 >= dimZ)) + k1 = k - k_m; + if(mask[(size_t) (dimX * dimY * k1) + (size_t) (j * dimX + i)] == true) + counter_vert_voxels++; + } + if(counter_vert_voxels < threshold_vertical) + out[index] = false; + } + } + else + { + /* + Considering the depth of features an removing the deep ones + Here we assume that the stripes do not normally extend far + in the depth dimension compared to the features that belong to a + sample. + */ + if(mask[index] == true) + { + if(stripe_depth_min != 0) + { + counter_depth_voxels = 0; + for(y_m = -halfstripe_depth; y_m <= halfstripe_depth; y_m++) + { + y1 = j + y_m; + if((y1 < 0) || (y1 >= dimY)) + y1 = j - y_m; + if(mask[(size_t) (dimX * dimY * k) + (size_t) (y1 * dimX + i)] == + true) + counter_depth_voxels++; + } + if(counter_depth_voxels > threshold_depth) + out[index] = false; + } + } + } +} + +void +remove_short_stripes(const bool* mask, bool* out, int stripe_length_min, long i, long j, + long k, long dimX, long dimY, long dimZ) +{ + int counter_vert_voxels; + int halfstripe_length = stripe_length_min / 2; + long k_m; + long k1; + size_t index; + index = (size_t) (dimX * dimY * k) + (size_t) (j * dimX + i); + + if(mask[index] == true) + { + counter_vert_voxels = 0; + for(k_m = -halfstripe_length; k_m <= halfstripe_length; k_m++) + { + k1 = k + k_m; + if((k1 < 0) || (k1 >= dimZ)) + k1 = k - k_m; + if(mask[(size_t) (dimX * dimY * k1) + (size_t) (j * dimX + i)] == true) + counter_vert_voxels++; + } + if(counter_vert_voxels < halfstripe_length) + out[index] = false; + } +} + +void +merge_stripes(const bool* mask, bool* out, int stripe_length_min, int stripe_width_min, + long i, long j, long k, long dimX, long dimY, long dimZ) +{ + int halfstripe_width = stripe_width_min / 2; + int vertical_length = 2 * stripe_width_min; + + long x; + long x_l; + long x_r; + long k_u; + long k_d; + int mask_left; + int mask_right; + int mask_up; + int mask_down; + size_t index; + index = (size_t) (dimX * dimY * k) + (size_t) (j * dimX + i); + + if(mask[index] == false) + { + /* checking if there is a mask to the left of False */ + mask_left = 0; + for(x = -halfstripe_width; x <= 0; x++) + { + x_l = i + x; + if(x_l < 0) + x_l = i - x; + if(mask[(size_t) (dimX * dimY * k) + (size_t) (j * dimX + x_l)] == true) + { + mask_left = 1; + break; + } + } + /* checking if there is a mask to the right of False */ + mask_right = 0; + for(x = 0; x <= halfstripe_width; x++) + { + x_r = i + x; + if(x_r >= dimX) + x_r = i - x; + if(mask[(size_t) (dimX * dimY * k) + (size_t) (j * dimX + x_r)] == true) + { + mask_right = 1; + break; + } + } + /* now if there is a mask from the left and from the right side of True value make + * it True */ + if((mask_left == 1) && (mask_right == 1)) + out[index] = true; + + /* perform vertical merging */ + if(out[index] == false) + { + /* checking if there is a mask up of True */ + mask_up = 0; + for(x = -vertical_length; x <= 0; x++) + { + k_u = k + x; + if(k_u < 0) + k_u = k - x; + if(mask[(size_t) (dimX * dimY * k_u) + (size_t) (j * dimX + i)] == true) + { + mask_up = 1; + break; + } + } + /* checking if there is a mask down of False */ + mask_down = 0; + for(x = 0; x <= vertical_length; x++) + { + k_d = k + x; + if(k_d >= dimZ) + k_d = k - x; + if(mask[(size_t) (dimX * dimY * k_d) + (size_t) (j * dimX + i)] == true) + { + mask_down = 1; + break; + } + } + /* now if there is a mask above and bellow of the False make it True */ + if((mask_up == 1) && (mask_down == 1)) + out[index] = true; + } + } +} + +/********************************************************************/ +/*************************stripesdetect3d****************************/ +/********************************************************************/ +DLL int +stripesdetect3d_main_float(float* Input, float* Output, int window_halflength_vertical, + int ratio_radius, int ncores, int dimX, int dimY, int dimZ) +{ + long i; + long j; + long k; + long long totalvoxels; + totalvoxels = (long long) ((long) (dimX) * (long) (dimY) * (long) (dimZ)); + + int window_fulllength = (2 * window_halflength_vertical + 1); + int midval_window_index = (int) (0.5F * window_fulllength) - 1; + + float* temp3d_arr; + temp3d_arr = malloc(totalvoxels * sizeof(float)); + if(temp3d_arr == NULL) + { + printf("Memory allocation of the 'temp3d_arr' array in " + "'stripesdetect3d_main_float' failed"); + } + + /* dealing here with a custom given number of cpu threads */ + if(ncores > 0) + { + // Explicitly disable dynamic teams + omp_set_dynamic(0); + // Use a number of threads for all consecutive parallel regions + omp_set_num_threads(ncores); + } + +/* Perform a gentle (6-stencil) 3d mean smoothing of the data to ensure more stability in + * the gradient calculation */ +#pragma omp parallel for shared(temp3d_arr) private(i, j, k) + for(k = 0; k < dimZ; k++) + { + for(j = 0; j < dimY; j++) + { + for(i = 0; i < dimX; i++) + { + mean_stride3d(Input, temp3d_arr, i, j, k, dimX, dimY, dimZ); + } + } + } + + /* Take the gradient in the horizontal direction, axis = 0, step = 2*/ + gradient3D_local(Input, Output, dimX, dimY, dimZ, 0, 2); + + /* + Here we calculate a ratio between the mean in a small 2D neighbourhood parallel to the + stripe and the mean orthogonal to the stripe. The gradient variation in the direction + orthogonal to the stripe is expected to be large (a jump), while in parallel direction + small. Therefore at the edges of a stripe we should get a ratio small/large or + large/small. + */ +#pragma omp parallel for shared(Output, temp3d_arr) private(i, j, k) + for(k = 0; k < dimZ; k++) + { + for(j = 0; j < dimY; j++) + { + for(i = 0; i < dimX; i++) + { + ratio_mean_stride3d(Output, temp3d_arr, ratio_radius, i, j, k, dimX, dimY, + dimZ); + } + } + } + + /* + We process the resulting ratio map with a vertical median filter which removes + inconsistent from longer stripes features + */ +#pragma omp parallel for shared(temp3d_arr, Output) private(i, j, k) + for(k = 0; k < dimZ; k++) + { + for(j = 0; j < dimY; j++) + { + for(i = 0; i < dimX; i++) + { + vertical_median_stride3d(temp3d_arr, Output, window_halflength_vertical, + window_fulllength, midval_window_index, i, j, k, + dimX, dimY, dimZ); + } + } + } + + free(temp3d_arr); + return 0; +} + +/********************************************************************/ +/*************************stripesmask3d******************************/ +/********************************************************************/ +DLL int +stripesmask3d_main_float(float* Input, bool* Output, float threshold_val, + int stripe_length_min, int stripe_depth_min, + int stripe_width_min, float sensitivity, int ncores, int dimX, + int dimY, int dimZ) +{ + long i; + long j; + long k; + int iter_merge; + int switch_dim; + size_t index; + size_t totalvoxels = (long) (dimX) * (long) (dimY) * (long) (dimZ); + + bool* mask; + mask = malloc(totalvoxels * sizeof(bool)); + if(mask == NULL) + { + printf( + "Memory allocation of the 'mask' array in 'stripesmask3d_main_float' failed"); + } + + /* dealing here with a custom given number of cpu threads */ + if(ncores > 0) + { + // Explicitly disable dynamic teams + omp_set_dynamic(0); + // Use a number of threads for all consecutive parallel regions + omp_set_num_threads(ncores); + } + + /* + First step is to mask all the values in the given weights input image + that are bellow a given "threshold_val" parameter + */ +#pragma omp parallel for shared(Input, mask) private(i, j, k, index) + for(k = 0; k < dimZ; k++) + { + for(j = 0; j < dimY; j++) + { + for(i = 0; i < dimX; i++) + { + index = (size_t) (dimX * dimY * k) + (size_t) (j * dimX + i); + if(Input[index] <= threshold_val) + { + mask[index] = true; + } + else + { + mask[index] = false; + } + } + } + } + + /* Copy mask to output */ + memcpy(Output, mask, totalvoxels * sizeof(bool)); + + /* the depth consistency for features */ + switch_dim = 1; +#pragma omp parallel for shared(mask, Output) private(i, j, k) + for(k = 0; k < dimZ; k++) + { + for(j = 0; j < dimY; j++) + { + for(i = 0; i < dimX; i++) + { + remove_inconsistent_stripes(mask, Output, stripe_length_min, + stripe_depth_min, sensitivity, switch_dim, i, + j, k, dimX, dimY, dimZ); + } + } + } + /* Copy output to mask */ + memcpy(mask, Output, totalvoxels * sizeof(bool)); + + /* + Now we need to remove stripes that are shorter than "stripe_length_min" parameter + or inconsistent otherwise. For every pixel we will run a 1D vertical window to count + nonzero values in the mask. We also check for the depth of the mask's value, + assuming that the stripes are normally shorter in depth compare to the features that + belong to true data. + */ + + /*continue by including longer vertical features and discarding the shorter ones */ + switch_dim = 0; +#pragma omp parallel for shared(mask, Output) private(i, j, k) + for(k = 0; k < dimZ; k++) + { + for(j = 0; j < dimY; j++) + { + for(i = 0; i < dimX; i++) + { + remove_inconsistent_stripes(mask, Output, stripe_length_min, + stripe_depth_min, sensitivity, switch_dim, i, + j, k, dimX, dimY, dimZ); + } + } + } + /* Copy output to mask */ + memcpy(mask, Output, totalvoxels * sizeof(bool)); + + /* now we clean the obtained mask if the features do not hold our assumptions about + * the lengths */ + +#pragma omp parallel for shared(mask, Output) private(i, j, k) + for(k = 0; k < dimZ; k++) + { + for(j = 0; j < dimY; j++) + { + for(i = 0; i < dimX; i++) + { + remove_short_stripes(mask, Output, stripe_length_min, i, j, k, dimX, dimY, + dimZ); + } + } + } + + /* Copy output to mask */ + memcpy(mask, Output, totalvoxels * sizeof(bool)); + + /* + We can merge stripes together if they are relatively close to each other + horizontally and vertically. We do that iteratively. + */ + for(iter_merge = 0; iter_merge < stripe_width_min; iter_merge++) + { +#pragma omp parallel for shared(mask, Output) private(i, j, k) + for(k = 0; k < dimZ; k++) + { + for(j = 0; j < dimY; j++) + { + for(i = 0; i < dimX; i++) + { + merge_stripes(mask, Output, stripe_length_min, stripe_width_min, i, j, + k, dimX, dimY, dimZ); + } + } + } + /* Copy output to mask */ + memcpy(mask, Output, totalvoxels * sizeof(bool)); + } + + free(mask); + return 0; +} diff --git a/source/tomopy/misc/corr.py b/source/tomopy/misc/corr.py index 53d80fb9d..0ebd2ef70 100644 --- a/source/tomopy/misc/corr.py +++ b/source/tomopy/misc/corr.py @@ -396,10 +396,10 @@ def median_filter3d(arr, size=3, ncore=None): # perform full 3D filtering if (input_type == 'float32'): - extern.c_median_filt3d_float32(arr, out, kernel_half_size, dif, ncore, + extern.c_median_filt3d_float32(np.ascontiguousarray(arr), out, kernel_half_size, dif, ncore, dx, dy, dz) else: - extern.c_median_filt3d_uint16(arr, out, kernel_half_size, dif, ncore, + extern.c_median_filt3d_uint16(np.ascontiguousarray(arr), out, kernel_half_size, dif, ncore, dx, dy, dz) return out @@ -455,10 +455,10 @@ def remove_outlier3d(arr, dif, size=3, ncore=None): # perform full 3D filtering if (input_type == 'float32'): - extern.c_median_filt3d_float32(arr, out, kernel_half_size, dif, ncore, + extern.c_median_filt3d_float32(np.ascontiguousarray(arr), out, kernel_half_size, dif, ncore, dx, dy, dz) else: - extern.c_median_filt3d_uint16(arr, out, kernel_half_size, dif, ncore, + extern.c_median_filt3d_uint16(np.ascontiguousarray(arr), out, kernel_half_size, dif, ncore, dx, dy, dz) return out diff --git a/source/tomopy/prep/stripe.py b/source/tomopy/prep/stripe.py index cd80d47b4..b5e81e468 100644 --- a/source/tomopy/prep/stripe.py +++ b/source/tomopy/prep/stripe.py @@ -45,7 +45,6 @@ # ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE # # POSSIBILITY OF SUCH DAMAGE. # # ######################################################################### - """ Module for pre-processing tasks. """ @@ -66,28 +65,36 @@ from scipy.ndimage import uniform_filter1d from scipy import interpolate import logging -logger = logging.getLogger(__name__) +logger = logging.getLogger(__name__) __author__ = "Doga Gursoy, Eduardo X. Miqueles, Nghia Vo" __credits__ = "Juan V. Bermudez, Hugo H. Slepicka" __copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' -__all__ = ['remove_stripe_fw', - 'remove_stripe_ti', - 'remove_stripe_sf', - 'remove_stripe_based_sorting', - 'remove_stripe_based_filtering', - 'remove_stripe_based_fitting', - 'remove_large_stripe', - 'remove_dead_stripe', - 'remove_all_stripe', - 'remove_stripe_based_interpolation'] - - -def remove_stripe_fw( - tomo, level=None, wname='db5', sigma=2, - pad=True, ncore=None, nchunk=None): +__all__ = [ + 'remove_stripe_fw', + 'remove_stripe_ti', + 'remove_stripe_sf', + 'remove_stripe_based_sorting', + 'remove_stripe_based_filtering', + 'remove_stripe_based_fitting', + 'remove_large_stripe', + 'remove_dead_stripe', + 'remove_all_stripe', + 'remove_stripe_based_interpolation', + 'stripes_detect3d', + 'stripes_mask3d', +] + + +def remove_stripe_fw(tomo, + level=None, + wname='db5', + sigma=2, + pad=True, + ncore=None, + nchunk=None): """ Remove horizontal stripes from sinogram using the Fourier-Wavelet (FW) based method :cite:`Munch:09`. @@ -118,13 +125,12 @@ def remove_stripe_fw( size = np.max(tomo.shape) level = int(np.ceil(np.log2(size))) - arr = mproc.distribute_jobs( - tomo, - func=_remove_stripe_fw, - args=(level, wname, sigma, pad), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=_remove_stripe_fw, + args=(level, wname, sigma, pad), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr @@ -163,8 +169,8 @@ def _remove_stripe_fw(tomo, level, wname, sigma, pad): fcV *= np.transpose(np.tile(damp, (mx, 1))) # Inverse FFT. - cV[n] = np.real(ifft(np.fft.ifftshift( - fcV), axis=0, extra_info=num_jobs)) + cV[n] = np.real( + ifft(np.fft.ifftshift(fcV), axis=0, extra_info=num_jobs)) # Wavelet reconstruction. for n in range(level)[::-1]: @@ -196,13 +202,12 @@ def remove_stripe_ti(tomo, nblock=0, alpha=1.5, ncore=None, nchunk=None): ndarray Corrected 3D tomographic data. """ - arr = mproc.distribute_jobs( - tomo, - func=_remove_stripe_ti, - args=(nblock, alpha), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=_remove_stripe_ti, + args=(nblock, alpha), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr @@ -224,12 +229,20 @@ def _remove_stripe_ti(tomo, nblock, alpha): def _kernel(m, n): - v = [[np.array([1, -1]), - np.array([-3 / 2, 2, -1 / 2]), - np.array([-11 / 6, 3, -3 / 2, 1 / 3])], - [np.array([-1, 2, -1]), - np.array([2, -5, 4, -1])], - [np.array([-1, 3, -3, 1])]] + v = [ + [ + np.array([1, -1]), + np.array([-3 / 2, 2, -1 / 2]), + np.array([-11 / 6, 3, -3 / 2, 1 / 3]), + ], + [ + np.array([-1, 2, -1]), + np.array([2, -5, 4, -1]), + ], + [ + np.array([-1, 3, -3, 1]), + ], + ] return v[m - 1][n - 1] @@ -342,17 +355,20 @@ def remove_stripe_sf(tomo, size=5, ncore=None, nchunk=None): Corrected 3D tomographic data. """ tomo = dtype.as_float32(tomo) - arr = mproc.distribute_jobs( - tomo, - func=extern.c_remove_stripe_sf, - args=(size,), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=extern.c_remove_stripe_sf, + args=(size,), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr -def remove_stripe_based_sorting(tomo, size=None, dim=1, ncore=None, nchunk=None): +def remove_stripe_based_sorting(tomo, + size=None, + dim=1, + ncore=None, + nchunk=None): """ Remove full and partial stripe artifacts from sinogram using Nghia Vo's approach :cite:`Vo:18` (algorithm 3). @@ -376,13 +392,12 @@ def remove_stripe_based_sorting(tomo, size=None, dim=1, ncore=None, nchunk=None) ndarray Corrected 3D tomographic data. """ - arr = mproc.distribute_jobs( - tomo, - func=_remove_stripe_based_sorting, - args=(size, dim), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=_remove_stripe_based_sorting, + args=(size, dim), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr @@ -401,14 +416,12 @@ def _rs_sort(sinogram, size, matindex, dim): """ sinogram = np.transpose(sinogram) matcomb = np.asarray(np.dstack((matindex, sinogram))) - matsort = np.asarray( - [row[row[:, 1].argsort()] for row in matcomb]) + matsort = np.asarray([row[row[:, 1].argsort()] for row in matcomb]) if dim == 1: matsort[:, :, 1] = median_filter(matsort[:, :, 1], (size, 1)) else: matsort[:, :, 1] = median_filter(matsort[:, :, 1], (size, size)) - matsortback = np.asarray( - [row[row[:, 0].argsort()] for row in matsort]) + matsortback = np.asarray([row[row[:, 0].argsort()] for row in matsort]) sino_corrected = matsortback[:, :, 1] return np.transpose(sino_corrected) @@ -425,8 +438,12 @@ def _remove_stripe_based_sorting(tomo, size, dim): tomo[:, m, :] = _rs_sort(sino, size, matindex, dim) -def remove_stripe_based_filtering( - tomo, sigma=3, size=None, dim=1, ncore=None, nchunk=None): +def remove_stripe_based_filtering(tomo, + sigma=3, + size=None, + dim=1, + ncore=None, + nchunk=None): """ Remove stripe artifacts from sinogram using Nghia Vo's approach :cite:`Vo:18` (algorithm 2). @@ -453,13 +470,12 @@ def remove_stripe_based_filtering( ndarray Corrected 3D tomographic data. """ - arr = mproc.distribute_jobs( - tomo, - func=_remove_stripe_based_filtering, - args=(sigma, size, dim), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=_remove_stripe_based_filtering, + args=(sigma, size, dim), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr @@ -486,8 +502,8 @@ def _rs_filter(sinogram, window, listsign, size, dim, pad): ifft(fft(sinolist * listsign) * window) * listsign)[pad:ncol - pad] sinosharp = sinogram - sinosmooth matindex = _create_matindex(nrow, ncol - 2 * pad) - sinosmooth_cor = np.transpose(_rs_sort( - np.transpose(sinosmooth), size, matindex, dim)) + sinosmooth_cor = np.transpose( + _rs_sort(np.transpose(sinosmooth), size, matindex, dim)) return np.transpose(sinosmooth_cor + sinosharp) @@ -502,12 +518,14 @@ def _remove_stripe_based_filtering(tomo, sigma, size, dim): size = max(5, int(0.01 * tomo.shape[2])) for m in range(tomo.shape[1]): sino = tomo[:, m, :] - tomo[:, m, :] = _rs_filter( - sino, window, listsign, size, dim, pad) + tomo[:, m, :] = _rs_filter(sino, window, listsign, size, dim, pad) -def remove_stripe_based_fitting( - tomo, order=3, sigma=(5, 20), ncore=None, nchunk=None): +def remove_stripe_based_fitting(tomo, + order=3, + sigma=(5, 20), + ncore=None, + nchunk=None): """ Remove stripe artifacts from sinogram using Nghia Vo's approach :cite:`Vo:18` (algorithm 1). @@ -532,13 +550,12 @@ def remove_stripe_based_fitting( ndarray Corrected 3D tomographic data. """ - arr = mproc.distribute_jobs( - tomo, - func=_remove_stripe_based_fitting, - args=(order, sigma), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=_remove_stripe_based_fitting, + args=(order, sigma), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr @@ -637,8 +654,13 @@ def _remove_stripe_based_fitting(tomo, order, sigma): tomo[:, m, :] = _rs_fit(sino, order, win2d, matsign, pad) -def remove_large_stripe(tomo, snr=3, size=51, drop_ratio=0.1, norm=True, - ncore=None, nchunk=None): +def remove_large_stripe(tomo, + snr=3, + size=51, + drop_ratio=0.1, + norm=True, + ncore=None, + nchunk=None): """ Remove large stripe artifacts from sinogram using Nghia Vo's approach :cite:`Vo:18` (algorithm 5). @@ -667,13 +689,12 @@ def remove_large_stripe(tomo, snr=3, size=51, drop_ratio=0.1, norm=True, ndarray Corrected 3D tomographic data. """ - arr = mproc.distribute_jobs( - tomo, - func=_remove_large_stripe, - args=(snr, size, drop_ratio, norm), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=_remove_large_stripe, + args=(snr, size, drop_ratio, norm), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr @@ -685,8 +706,8 @@ def _detect_stripe(listdata, snr): listsorted = np.sort(listdata)[::-1] xlist = np.arange(0, numdata, 1.0) ndrop = np.int16(0.25 * numdata) - (_slope, _intercept) = np.polyfit( - xlist[ndrop:-ndrop - 1], listsorted[ndrop:-ndrop - 1], 1) + (_slope, _intercept) = np.polyfit(xlist[ndrop:-ndrop - 1], + listsorted[ndrop:-ndrop - 1], 1) numt1 = _intercept + _slope * xlist[-1] noiselevel = np.abs(numt1 - _intercept) noiselevel = np.clip(noiselevel, 1e-6, None) @@ -713,8 +734,10 @@ def _rs_large(sinogram, snr, size, matindex, drop_ratio=0.1, norm=True): sinosmooth = median_filter(sinosort, (1, size)) list1 = np.mean(sinosort[ndrop:nrow - ndrop], axis=0) list2 = np.mean(sinosmooth[ndrop:nrow - ndrop], axis=0) - listfact = np.divide(list1, list2, - out=np.ones_like(list1), where=list2 != 0) + listfact = np.divide(list1, + list2, + out=np.ones_like(list1), + where=list2 != 0) # Locate stripes listmask = _detect_stripe(listfact, snr) listmask = binary_dilation(listmask, iterations=1).astype(listmask.dtype) @@ -724,11 +747,9 @@ def _rs_large(sinogram, snr, size, matindex, drop_ratio=0.1, norm=True): sinogram = sinogram / matfact sinogram1 = np.transpose(sinogram) matcombine = np.asarray(np.dstack((matindex, sinogram1))) - matsort = np.asarray( - [row[row[:, 1].argsort()] for row in matcombine]) + matsort = np.asarray([row[row[:, 1].argsort()] for row in matcombine]) matsort[:, :, 1] = np.transpose(sinosmooth) - matsortback = np.asarray( - [row[row[:, 0].argsort()] for row in matsort]) + matsortback = np.asarray([row[row[:, 0].argsort()] for row in matsort]) sino_corrected = np.transpose(matsortback[:, :, 1]) listxmiss = np.where(listmask > 0.0)[0] sinogram[:, listxmiss] = sino_corrected[:, listxmiss] @@ -742,8 +763,12 @@ def _remove_large_stripe(tomo, snr, size, drop_ratio, norm): tomo[:, m, :] = _rs_large(sino, snr, size, matindex, drop_ratio, norm) -def remove_dead_stripe(tomo, snr=3, size=51, norm=True, - ncore=None, nchunk=None): +def remove_dead_stripe(tomo, + snr=3, + size=51, + norm=True, + ncore=None, + nchunk=None): """ Remove unresponsive and fluctuating stripe artifacts from sinogram using Nghia Vo's approach :cite:`Vo:18` (algorithm 6). @@ -769,13 +794,12 @@ def remove_dead_stripe(tomo, snr=3, size=51, norm=True, ndarray Corrected 3D tomographic data. """ - arr = mproc.distribute_jobs( - tomo, - func=_remove_dead_stripe, - args=(snr, size, norm), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=_remove_dead_stripe, + args=(snr, size, norm), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr @@ -788,8 +812,10 @@ def _rs_dead(sinogram, snr, size, matindex, norm=True): sinosmooth = np.apply_along_axis(uniform_filter1d, 0, sinogram, 10) listdiff = np.sum(np.abs(sinogram - sinosmooth), axis=0) listdiffbck = median_filter(listdiff, size) - listfact = np.divide(listdiff, listdiffbck, - out=np.ones_like(listdiff), where=listdiffbck != 0) + listfact = np.divide(listdiff, + listdiffbck, + out=np.ones_like(listdiff), + where=listdiffbck != 0) listmask = _detect_stripe(listfact, snr) listmask = binary_dilation(listmask, iterations=1).astype(listmask.dtype) listmask[0:2] = 0.0 @@ -814,8 +840,13 @@ def _remove_dead_stripe(tomo, snr, size, norm): tomo[:, m, :] = _rs_dead(sino, snr, size, matindex, norm) -def remove_all_stripe(tomo, snr=3, la_size=61, sm_size=21, dim=1, - ncore=None, nchunk=None): +def remove_all_stripe(tomo, + snr=3, + la_size=61, + sm_size=21, + dim=1, + ncore=None, + nchunk=None): """ Remove all types of stripe artifacts from sinogram using Nghia Vo's approach :cite:`Vo:18` (combination of algorithm 3,4,5, and 6). @@ -843,13 +874,12 @@ def remove_all_stripe(tomo, snr=3, la_size=61, sm_size=21, dim=1, ndarray Corrected 3D tomographic data. """ - arr = mproc.distribute_jobs( - tomo, - func=_remove_all_stripe, - args=(snr, la_size, sm_size, dim), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=_remove_all_stripe, + args=(snr, la_size, sm_size, dim), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr @@ -862,8 +892,13 @@ def _remove_all_stripe(tomo, snr, la_size, sm_size, dim): tomo[:, m, :] = sino -def remove_stripe_based_interpolation(tomo, snr=3, size=31, drop_ratio=0.1, - norm=True, ncore=None, nchunk=None): +def remove_stripe_based_interpolation(tomo, + snr=3, + size=31, + drop_ratio=0.1, + norm=True, + ncore=None, + nchunk=None): """ Remove most types of stripe artifacts from sinograms based on interpolation. Derived from algorithm 4, 5, and 6 in :cite:`Vo:18`. @@ -893,13 +928,12 @@ def remove_stripe_based_interpolation(tomo, snr=3, size=31, drop_ratio=0.1, ndarray Corrected 3D tomographic data. """ - arr = mproc.distribute_jobs( - tomo, - func=_remove_stripe_based_interpolation, - args=(snr, size, drop_ratio, norm), - axis=1, - ncore=ncore, - nchunk=nchunk) + arr = mproc.distribute_jobs(tomo, + func=_remove_stripe_based_interpolation, + args=(snr, size, drop_ratio, norm), + axis=1, + ncore=ncore, + nchunk=nchunk) return arr @@ -915,8 +949,10 @@ def _rs_interpolation(sinogram, snr, size, drop_ratio=0.1, norm=True): sinosmooth = median_filter(sinosort, (1, size)) list1 = np.mean(sinosort[ndrop:nrow - ndrop], axis=0) list2 = np.mean(sinosmooth[ndrop:nrow - ndrop], axis=0) - listfact = np.divide(list1, list2, - out=np.ones_like(list1), where=list2 != 0) + listfact = np.divide(list1, + list2, + out=np.ones_like(list1), + where=list2 != 0) listmask = _detect_stripe(listfact, snr) listmask = np.float32(binary_dilation(listmask, iterations=1)) matfact = np.tile(listfact, (nrow, 1)) @@ -939,3 +975,205 @@ def _remove_stripe_based_interpolation(tomo, snr, size, drop_ratio, norm): sino = tomo[:, m, :] sino = _rs_interpolation(sino, snr, size, drop_ratio, norm) tomo[:, m, :] = sino + + +def stripes_detect3d(tomo, size=10, radius=3, ncore=None): + """ + Detect stripes in a 3D array. Usually applied to normalized projection + data. + + The method works with full and partial stripes of constant ot varying + intensity. + + .. versionadded:: 1.14 + + Parameters + ---------- + tomo : ndarray + 3D tomographic data of float32 data type, preferably in the [0, 1] + range, although reasonable deviations accepted (e.g. the result of the + normalization and the negative log taken of the raw data). The + projection data should be given with [angle, detY(depth), + detX(horizontal)] axis orientation. With this orientation, the stripes + are features along the angle axis. + size : int, optional + The pixel size of the 1D median filter orthogonal to stripes + orientation to minimise false detections. Increase it if you have + longer or full stripes in the data. + radius : int, optional + The pixel size of the 3D stencil to calculate the mean ratio between + the angular and detX orientations of the detX gradient. The larger + values can affect the width of the detected stripe, use 1,2,3 values. + ncore : int, optional + Number of cores that will be assigned to jobs. All cores will be used + if unspecified. + + Returns + ------- + ndarray + Weights in the range of [0, 1] of float32 data type where stripe's + edges are highlighted with the smaller (e.g. < 0.5) values. The weights + can be manually thresholded or passed to stripes_mask3d function for + further processing and a binary mask generation. + + Raises + ------ + ValueError + If the `tomo` is not three dimensional. + + If the `size` is invalid. + """ + if ncore is None: + ncore = mproc.mp.cpu_count() + + input_type = tomo.dtype + if (input_type != 'float32'): + tomo = dtype.as_float32(tomo) # silent convertion to float32 data type + out = np.empty_like(tomo, order='C') + + if tomo.ndim == 3: + dz, dy, dx = tomo.shape + if (dz == 0) or (dy == 0) or (dx == 0): + msg = "The length of one of dimensions is equal to zero" + raise ValueError(msg) + else: + msg = "The input array must be a 3D array" + raise ValueError(msg) + + if size <= 0 or size > dz // 2: + msg = ("The size of the filter should be larger than zero " + "and smaller than the half of the vertical dimension") + raise ValueError(msg) + + # perform stripes detection + extern.c_stripes_detect3d(np.ascontiguousarray(tomo), out, size, radius, + ncore, dx, dy, dz) + return out + + +def stripes_mask3d(weights, + threshold=0.6, + min_stripe_length=20, + min_stripe_depth=10, + min_stripe_width=5, + sensitivity_perc=85.0, + ncore=None): + """ + Takes the result of the stripes_detect3d module as an input and generates a + binary 3D mask with ones where stripes present. + + The method tries to eliminate non-stripe features in data by checking the + consistency of weights in three directions. + + .. versionadded:: 1.14 + + Parameters + ---------- + weights : ndarray + 3D weights array, a result of stripes_detect3d module given in [angles, + detY(depth), detX] axis orientation. + threshold : float, optional + Threshold for the given weights. This parameter defines what weights + will be considered as potential candidates for stripes. The lower value + (< 0.5) will result in only the most prominent stripes in the data. + Increase the threshold cautiously because increasing the threshold + increases the probability of false detections. The good range to try is + between 0.5 and 0.7. + min_stripe_length : int, optional + Minimum length of a stripe in pixels with respect to the "angles" axis. + If there are full stripes in the data, then this could be >50% of the + size of the the "angles" axis. + min_stripe_depth : int, optional + Minimum depth of a stripe in pixels with respect to the "detY" axis. + The stripes do not extend very deep normally in the data. By setting + this parameter to the approximate depth of the stripe more false alarms + can be removed. + min_stripe_width : int, optional + Minimum width of a stripe in pixels with respect to the "detX" axis. + The stripes that close to each other can be merged together with this + parameter. + sensitivity_perc : float, optional + The value in the range [0, 100] that controls the strictness of the + minimum length, depth and width parameters of a stripe. 0 is + less-strict. 100 is more-strict. + ncore : int, optional + Number of cores that will be assigned to jobs. All cores will be used + if unspecified. + + Returns + ------- + ndarray + A binary mask of bool data type with stripes highlighted as True + values. + + Raises + ------ + ValueError + If the input array is not three dimensional. + + If a min_stripe_length parameter is negative, zero or longer than its + corresponding dimension ("angle") + + If a min_stripe_depth parameter is negative or longer than its + corresponding dimension ("detY") + + If a min_stripe_width parameter is negative, zero or longer than its + corresponding dimension ("detX") + + If a sensitivity_perc parameter doesn't lie in the (0,100] range + + """ + if ncore is None: + ncore = mproc.mp.cpu_count() + + input_type = weights.dtype + if (input_type != 'float32'): + weights = dtype.as_float32( + weights) # silent convertion to float32 data type + out = np.zeros(np.shape(weights), dtype=bool, order='C') + + if weights.ndim == 3: + dz, dy, dx = weights.shape + if (dz == 0) or (dy == 0) or (dx == 0): + msg = "The length of one of dimensions is equal to zero" + raise ValueError(msg) + else: + msg = "The input array must be a 3D array" + raise ValueError(msg) + + if min_stripe_length <= 0 or min_stripe_length >= dz: + msg = ("The minimum length of a stripe cannot be zero " + "or exceed the size of the angular dimension") + raise ValueError(msg) + + if min_stripe_depth < 0 or min_stripe_depth >= dy: + msg = ("The minimum depth of a stripe cannot exceed " + "the size of the depth dimension") + raise ValueError(msg) + + if min_stripe_width <= 0 or min_stripe_width >= dx: + msg = ("The minimum width of a stripe cannot be zero " + "or exceed the size of the horizontal dimension") + raise ValueError(msg) + + if 0.0 < sensitivity_perc <= 100.0: + pass + else: + msg = "sensitivity_perc value must be in (0, 100] percentage range" + raise ValueError(msg) + + # perform mask creation based on the input provided by stripes_detect3d module + extern.c_stripesmask3d( + np.ascontiguousarray(weights), + out, + threshold, + min_stripe_length, + min_stripe_depth, + min_stripe_width, + sensitivity_perc, + ncore, + dx, + dy, + dz, + ) + return out diff --git a/source/tomopy/util/dtype.py b/source/tomopy/util/dtype.py index a17ec535d..f4de12ade 100644 --- a/source/tomopy/util/dtype.py +++ b/source/tomopy/util/dtype.py @@ -71,6 +71,8 @@ 'as_uint8', 'as_uint16', 'as_c_float_p', + 'as_c_bool_p', + 'as_c_uint8_p', 'as_c_uint16_p', 'as_c_int', 'as_c_int_p', @@ -117,6 +119,15 @@ def as_c_float_p(arr): return arr.ctypes.data_as(c_float_p) +def as_c_bool_p(arr): + c_bool_p = ctypes.POINTER(ctypes.c_bool) + return arr.ctypes.data_as(c_bool_p) + +def as_c_uint8_p(arr): + c_uint8_p = ctypes.POINTER(ctypes.c_uint8) + return arr.ctypes.data_as(c_uint8_p) + + def as_c_uint16_p(arr): c_uint16_p = ctypes.POINTER(ctypes.c_uint16) return arr.ctypes.data_as(c_uint16_p) @@ -125,6 +136,8 @@ def as_c_uint16_p(arr): def as_c_int(arr): return ctypes.c_int(arr) +def as_c_long(arr): + return ctypes.c_long(arr) def as_c_int_p(arr): arr = arr.astype(np.intc, copy=False) diff --git a/source/tomopy/util/extern/prep.py b/source/tomopy/util/extern/prep.py index 886cb534e..82f3f13cb 100644 --- a/source/tomopy/util/extern/prep.py +++ b/source/tomopy/util/extern/prep.py @@ -60,7 +60,9 @@ __copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC." __docformat__ = 'restructuredtext en' __all__ = ['c_normalize_bg', - 'c_remove_stripe_sf'] + 'c_remove_stripe_sf', + 'c_stripes_detect3d', + 'c_stripesmask3d'] LIB_TOMOPY_PREP = c_shared_lib("tomo-prep") @@ -96,3 +98,55 @@ def c_remove_stripe_sf(tomo, size): dtype.as_c_int(istart), dtype.as_c_int(iend)) tomo[:] = contiguous_tomo[:] + +def c_stripes_detect3d( + input, + output, + size, + radius, + ncore, + dx, + dy, + dz, +): + LIB_TOMOPY_PREP.stripesdetect3d_main_float.restype = dtype.as_c_void_p() + LIB_TOMOPY_PREP.stripesdetect3d_main_float( + dtype.as_c_float_p(input), + dtype.as_c_float_p(output), + dtype.as_c_int(size), + dtype.as_c_int(radius), + dtype.as_c_int(ncore), + dtype.as_c_int(dx), + dtype.as_c_int(dy), + dtype.as_c_int(dz), + ) + return output + +def c_stripesmask3d( + input, + output, + threshold_val, + min_stripe_length, + min_stripe_depth, + min_stripe_width, + sensitivity_perc, + ncore, + dx, + dy, + dz, +): + LIB_TOMOPY_PREP.stripesmask3d_main_float.restype = dtype.as_c_void_p() + LIB_TOMOPY_PREP.stripesmask3d_main_float( + dtype.as_c_float_p(input), + dtype.as_c_bool_p(output), + dtype.as_c_float(threshold_val), + dtype.as_c_int(min_stripe_length), + dtype.as_c_int(min_stripe_depth), + dtype.as_c_int(min_stripe_width), + dtype.as_c_float(sensitivity_perc), + dtype.as_c_int(ncore), + dtype.as_c_int(dx), + dtype.as_c_int(dy), + dtype.as_c_int(dz), + ) + return output diff --git a/test/test_tomopy/test_data/stripes_detect3d.npy b/test/test_tomopy/test_data/stripes_detect3d.npy new file mode 100644 index 000000000..eb6bb13bd Binary files /dev/null and b/test/test_tomopy/test_data/stripes_detect3d.npy differ diff --git a/test/test_tomopy/test_data/stripes_mask3d.npy b/test/test_tomopy/test_data/stripes_mask3d.npy new file mode 100644 index 000000000..9f1882c66 Binary files /dev/null and b/test/test_tomopy/test_data/stripes_mask3d.npy differ diff --git a/test/test_tomopy/test_data/test_stripe_data.npy b/test/test_tomopy/test_data/test_stripe_data.npy new file mode 100644 index 000000000..13ccf7fa1 Binary files /dev/null and b/test/test_tomopy/test_data/test_stripe_data.npy differ diff --git a/test/test_tomopy/test_prep/test_stripe.py b/test/test_tomopy/test_prep/test_stripe.py index 0419fc634..73256ffaf 100644 --- a/test/test_tomopy/test_prep/test_stripe.py +++ b/test/test_tomopy/test_prep/test_stripe.py @@ -145,3 +145,21 @@ def test_remove_stripe_based_interpolation(self): np.expand_dims(mat, 1), 1.5, 5)[:, 0, :] num = np.abs(np.mean(mat_corr[:, self.b:self.e]) - 6.0) self.assertTrue(num > self.eps) + + def test_stripe_detection(self): + assert_allclose( + srm.stripes_detect3d(read_file('test_stripe_data.npy'), + size=10, + radius=1), + read_file('stripes_detect3d.npy'), rtol=1e-6) + + def test_stripe_mask(self): + assert_allclose( + srm.stripes_mask3d(read_file('stripes_detect3d.npy'), + threshold=0.6, + min_stripe_length = 10, + min_stripe_depth = 0, + min_stripe_width = 5, + sensitivity_perc=85.0), + read_file('stripes_mask3d.npy'), rtol=1e-6) +