Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix for remove_outlier3d #622 #623

Merged
merged 10 commits into from Sep 28, 2023
Merged

Fix for remove_outlier3d #622 #623

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
e61dcb8
fixing the memory segfault issue
dkazanc Sep 24, 2023
cb61f9f
fixing the memory allocation issue
dkazanc Sep 24, 2023
1757ef9
removes stdio
dkazanc Sep 24, 2023
4585cd1
corrections
dkazanc Sep 26, 2023
01e5ec5
BLD: Update submodule to tomopy fork
carterbox Sep 26, 2023
1e92caa
REF: Add a ctypes size_t converter
carterbox Sep 26, 2023
3bafa8d
REF: Use size_t in median_flt without breaking ABI
carterbox Sep 26, 2023
fb71b0e
CI: Use libmamba solver to create environment
carterbox Sep 26, 2023
aa8e1cd
Merge branch 'ci-env' into fixmedian
carterbox Sep 26, 2023
2f0ab1c
BUG: Add path key to submodule
carterbox Sep 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
49 changes: 27 additions & 22 deletions source/libtomo/misc/median_filt3d.c
View file
Open in desktop
Expand Up @@ -64,9 +64,10 @@ int uint16comp(const void* elem1, const void* elem2)
return *(const unsigned short*)elem1 > *(const unsigned short*)elem2;
}


void
medfilt3D_float(float* Input, float* Output, int radius, int sizefilter_total,
float mu_threshold, long i, long j, long k, long index, long dimX,
float mu_threshold, long i, long j, long k, size_t index, long dimX,
long dimY, long dimZ)
{
float* ValVec;
Expand All @@ -78,6 +79,7 @@ medfilt3D_float(float* Input, float* Output, int radius, int sizefilter_total,
long k1;
long long counter;
int midval;
size_t index1;
midval = (sizefilter_total / 2);
ValVec = (float*) calloc(sizefilter_total, sizeof(float));

Expand All @@ -98,7 +100,8 @@ medfilt3D_float(float* Input, float* Output, int radius, int sizefilter_total,
k1 = k + k_m;
if((k1 < 0) || (k1 >= dimZ))
k1 = k;
ValVec[counter] = Input[(dimX * dimY) * k1 + j1 * dimX + i1];
index1 = (size_t)(dimX * dimY * k1 ) + (size_t)(j1 * dimX + i1);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What's the pro/con of using size_t for all indexing and size variables instead of a mix of long and size_t with casting? Do some of these indexing variable need to represent negative integers?

My concern with mixing types and casting on lines like this is that it might still overflow before the cast to size_t? Because the order of operations is to perform the computation inside the parenthesis before casting.

Maybe it should be something like this? Where everything that is not cast is already size_t? This way the multiplications are computed using types that are already size_t.

Suggested change
index1 = (size_t)(dimX * dimY * k1 ) + (size_t)(j1 * dimX + i1);
index1 = dimX * dimY * (size_t)k1 + (size_t)j1 * dimX + (size_t)i1;

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks, makes total sense. I guess (size_t)(dimX * dimY * k1 ) could overflow before getting converted to size_t. will try to fix...

ValVec[counter] = Input[index1];
counter++;
}
}
Expand All @@ -122,7 +125,7 @@ medfilt3D_float(float* Input, float* Output, int radius, int sizefilter_total,

void
medfilt2D_float(float* Input, float* Output, int radius, int sizefilter_total,
float mu_threshold, long i, long j, long index, long dimX, long dimY)
float mu_threshold, long i, long j, size_t index, long dimX, long dimY)
{
float* ValVec;
long i_m;
Expand All @@ -131,6 +134,7 @@ medfilt2D_float(float* Input, float* Output, int radius, int sizefilter_total,
long j1;
long long counter;
int midval;
size_t index1;
midval = (sizefilter_total / 2);
ValVec = (float*) calloc(sizefilter_total, sizeof(float));

Expand All @@ -146,7 +150,8 @@ medfilt2D_float(float* Input, float* Output, int radius, int sizefilter_total,
j1 = j + j_m;
if((j1 < 0) || (j1 >= dimY))
j1 = j;
ValVec[counter] = Input[j1 * dimX + i1];
index1 = (size_t)(j1 * dimX + i1);
ValVec[counter] = Input[index1];
counter++;
}
}
Expand All @@ -169,7 +174,7 @@ medfilt2D_float(float* Input, float* Output, int radius, int sizefilter_total,
void
medfilt3D_uint16(unsigned short* Input, unsigned short* Output, int radius,
int sizefilter_total, float mu_threshold, long i, long j, long k,
long index, long dimX, long dimY, long dimZ)
size_t index, long dimX, long dimY, long dimZ)
{
unsigned short* ValVec;
long i_m;
Expand All @@ -180,6 +185,7 @@ medfilt3D_uint16(unsigned short* Input, unsigned short* Output, int radius,
long k1;
long long counter;
int midval;
size_t index1;
midval = (sizefilter_total / 2);
ValVec = (unsigned short*) calloc(sizefilter_total, sizeof(unsigned short));

Expand All @@ -200,7 +206,8 @@ medfilt3D_uint16(unsigned short* Input, unsigned short* Output, int radius,
k1 = k + k_m;
if((k1 < 0) || (k1 >= dimZ))
k1 = k;
ValVec[counter] = Input[(dimX * dimY) * k1 + j1 * dimX + i1];
index1 = (size_t)(dimX * dimY * k1 ) + (size_t)(j1 * dimX + i1);
ValVec[counter] = Input[index1];
counter++;
}
}
Expand All @@ -224,7 +231,7 @@ medfilt3D_uint16(unsigned short* Input, unsigned short* Output, int radius,

void
medfilt2D_uint16(unsigned short* Input, unsigned short* Output, int radius,
int sizefilter_total, float mu_threshold, long i, long j, long index,
int sizefilter_total, float mu_threshold, long i, long j, size_t index,
long dimX, long dimY)
{
unsigned short* ValVec;
Expand All @@ -234,6 +241,7 @@ medfilt2D_uint16(unsigned short* Input, unsigned short* Output, int radius,
long j1;
long long counter;
int midval;
size_t index1;
midval = (sizefilter_total / 2);
ValVec = (unsigned short*) calloc(sizefilter_total, sizeof(unsigned short));

Expand All @@ -249,7 +257,8 @@ medfilt2D_uint16(unsigned short* Input, unsigned short* Output, int radius,
j1 = j + j_m;
if((j1 < 0) || (j1 >= dimY))
j1 = j;
ValVec[counter] = Input[j1 * dimX + i1];
index1 = (size_t)(j1 * dimX + i1);
ValVec[counter] = Input[index1];
counter++;
}
}
Expand Down Expand Up @@ -279,13 +288,11 @@ medianfilter_main_float(float* Input, float* Output, int radius, float mu_thresh
long i;
long j;
long k;
long long index;
size_t index;
size_t totalvoxels;

totalvoxels = (size_t)(dimX * dimY * dimZ);
diameter = (2 * radius + 1); /* diameter of the filter's kernel */
if(mu_threshold != 0.0)
{
memcpy(Output, Input, dimX * dimY * dimZ * sizeof(float));
} /* copy input into output */

/* dealing here with a custom given number of cpu threads */
if(ncores > 0)
Expand All @@ -304,7 +311,7 @@ medianfilter_main_float(float* Input, float* Output, int radius, float mu_thresh
{
for(i = 0; i < dimX; i++)
{
index = (j * dimX + i);
index = (size_t)(j * dimX + i);
medfilt2D_float(Input, Output, radius, sizefilter_total, mu_threshold, i,
j, index, (long) (dimX), (long) (dimY));
}
Expand All @@ -321,7 +328,7 @@ medianfilter_main_float(float* Input, float* Output, int radius, float mu_thresh
{
for(i = 0; i < dimX; i++)
{
index = ((dimX * dimY) * k + j * dimX + i);
index = (size_t)((dimX * dimY) * k) + (size_t)(j * dimX + i);
medfilt3D_float(Input, Output, radius, sizefilter_total, mu_threshold,
i, j, k, index, (long) (dimX), (long) (dimY),
(long) (dimZ));
Expand All @@ -342,13 +349,11 @@ medianfilter_main_uint16(unsigned short* Input, unsigned short* Output, int radi
long i;
long j;
long k;
long long index;
size_t index;
size_t totalvoxels;

totalvoxels = (size_t)(dimX * dimY * dimZ);
diameter = (2 * radius + 1); /* diameter of the filter's kernel */
if(mu_threshold != 0.0)
{
memcpy(Output, Input, dimX * dimY * dimZ * sizeof(unsigned short));
} /* copy input into output */

/* dealing here with a custom given number of cpu threads */
if(ncores > 0)
Expand All @@ -367,7 +372,7 @@ medianfilter_main_uint16(unsigned short* Input, unsigned short* Output, int radi
{
for(i = 0; i < dimX; i++)
{
index = (j * dimX + i);
index = (size_t)(j * dimX + i);
medfilt2D_uint16(Input, Output, radius, sizefilter_total, mu_threshold, i,
j, index, (long) (dimX), (long) (dimY));
}
Expand All @@ -384,7 +389,7 @@ medianfilter_main_uint16(unsigned short* Input, unsigned short* Output, int radi
{
for(i = 0; i < dimX; i++)
{
index = ((dimX * dimY) * k + j * dimX + i);
index = (size_t)(dimX * dimY * k) + (size_t)(j * dimX + i);
medfilt3D_uint16(Input, Output, radius, sizefilter_total,
mu_threshold, i, j, k, index, (long) (dimX),
(long) (dimY), (long) (dimZ));
Expand Down
9 changes: 5 additions & 4 deletions source/tomopy/misc/corr.py
View file
Open in desktop
Expand Up @@ -62,6 +62,7 @@
import numexpr as ne
import concurrent.futures as cf
from scipy.signal import medfilt2d
import copy

logger = logging.getLogger(__name__)

Expand Down Expand Up @@ -441,8 +442,8 @@ def remove_outlier3d(arr, dif, size=3, ncore=None):
input_type = arr.dtype
if (input_type != 'float32') and (input_type != 'uint16'):
arr = dtype.as_float32(arr) # silent convertion to float32 data type
out = np.empty_like(arr)

out = copy.deepcopy(arr)
dkazanc marked this conversation as resolved.
Show resolved Hide resolved
# convert the full kernel size (odd int) to a half size as the C function requires it
kernel_half_size = (max(int(size), 3) - 1) // 2

Expand All @@ -455,10 +456,10 @@ def remove_outlier3d(arr, dif, size=3, ncore=None):

# perform full 3D filtering
if (input_type == 'float32'):
extern.c_median_filt3d_float32(np.ascontiguousarray(arr), out, kernel_half_size, dif, ncore,
extern.c_median_filt3d_float32(np.ascontiguousarray(arr), np.ascontiguousarray(out), kernel_half_size, dif, ncore,
dx, dy, dz)
else:
extern.c_median_filt3d_uint16(np.ascontiguousarray(arr), out, kernel_half_size, dif, ncore,
extern.c_median_filt3d_uint16(np.ascontiguousarray(arr), np.ascontiguousarray(out), kernel_half_size, dif, ncore,
dx, dy, dz)
return out

Expand Down