Skip to content

Commit

Permalink
Merge pull request #623 from dkazanc/fixmedian
Browse files Browse the repository at this point in the history 
Fix for remove_outlier3d #622
  • Loading branch information
@carterbox
carterbox committed Sep 28, 2023
2 parents f02ee59 + 2f0ab1c commit 32b6497
Show file tree
Hide file tree
Showing 5 changed files with 48 additions and 38 deletions.
4 changes: 2 additions & 2 deletions .gitmodules
View file
@@ -1,3 +1,3 @@
[submodule "src/PTL"]
[submodule "source/PTL"]
path = source/PTL
url = https://github.com/jrmadsen/PTL.git
url = https://github.com/tomopy/PTL.git
6 changes: 3 additions & 3 deletions azure-pipelines.yml
View file
Expand Up @@ -76,7 +76,7 @@ jobs:
- bash: echo "##vso[task.prependpath]$CONDA/bin"
displayName: Add conda to PATH
- script: |
conda env create -n tomopy --quiet -f envs/linux-$(python.version).yml
conda env create -n tomopy --quiet -f envs/linux-$(python.version).yml --solver=libmamba
displayName: Create build environment
- script: conda list -n tomopy
displayName: List build environment
Expand Down Expand Up @@ -106,7 +106,7 @@ jobs:
- bash: sudo chown -R $USER $CONDA
displayName: Take ownership of conda installation
- script: |
conda env create -n tomopy --quiet -f envs/osx-$(python.version).yml
conda env create -n tomopy --quiet -f envs/osx-$(python.version).yml --solver=libmamba
displayName: Create build environment
- script: conda list -n tomopy
displayName: List build environment
Expand Down Expand Up @@ -134,7 +134,7 @@ jobs:
- powershell: Write-Host "##vso[task.prependpath]$env:CONDA\Scripts"
displayName: Add conda to PATH
- script: |
conda env create -n tomopy --quiet -f envs/win-ci.yml
conda env create -n tomopy --quiet -f envs/win-ci.yml --solver=libmamba
displayName: Create build environment
- script: conda list -n tomopy
displayName: List build environment
Expand Down
67 changes: 36 additions & 31 deletions source/libtomo/misc/median_filt3d.c
View file
Expand Up @@ -64,10 +64,11 @@ 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,
long dimY, long dimZ)
float mu_threshold, long i, long j, long k, size_t index, size_t dimX,
size_t dimY, size_t dimZ)
{
float* ValVec;
long i_m;
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 = dimX * dimY * (size_t)k1 + (size_t)j1 * dimX + (size_t)i1;
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, size_t dimX, size_t 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 + (size_t)(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, size_t dimX, size_t dimY, size_t 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,12 +206,13 @@ 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 = dimX * dimY * (size_t)k1 + (size_t)j1 * dimX + (size_t)i1;
ValVec[counter] = Input[index1];
counter++;
}
}
}

qsort(ValVec, sizefilter_total, sizeof(unsigned short), uint16comp);

if(mu_threshold == 0.0F)
Expand All @@ -224,8 +231,8 @@ 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,
long dimX, long dimY)
int sizefilter_total, float mu_threshold, long i, long j, size_t index,
size_t dimX, size_t dimY)
{
unsigned short* ValVec;
long i_m;
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 + (size_t)(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) * (size_t)(dimY) * (size_t)(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,9 +311,9 @@ 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 + (size_t)(i);
medfilt2D_float(Input, Output, radius, sizefilter_total, mu_threshold, i,
j, index, (long) (dimX), (long) (dimY));
j, index, (size_t) dimX, (size_t) dimY);
}
}
}
Expand All @@ -321,10 +328,10 @@ 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 = dimX * dimY * (size_t)(k) + (size_t)(j) * dimX + (size_t)(i);
medfilt3D_float(Input, Output, radius, sizefilter_total, mu_threshold,
i, j, k, index, (long) (dimX), (long) (dimY),
(long) (dimZ));
i, j, k, index, (size_t) dimX, (size_t) dimY,
(size_t) 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) * (size_t)(dimY) * (size_t)(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,9 +372,9 @@ 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 + (size_t)(i);
medfilt2D_uint16(Input, Output, radius, sizefilter_total, mu_threshold, i,
j, index, (long) (dimX), (long) (dimY));
j, index, (size_t) dimX, (size_t) dimY);
}
}
}
Expand All @@ -384,10 +389,10 @@ 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 = dimX * dimY * (size_t)(k) + (size_t)(j) * dimX + (size_t)(i);
medfilt3D_uint16(Input, Output, radius, sizefilter_total,
mu_threshold, i, j, k, index, (long) (dimX),
(long) (dimY), (long) (dimZ));
mu_threshold, i, j, k, index,(size_t) dimX,
(size_t)dimY, (size_t)dimZ);
}
}
}
Expand Down
4 changes: 2 additions & 2 deletions source/tomopy/misc/corr.py
View file
Expand Up @@ -441,8 +441,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 = np.copy(arr, order='C')
# 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 Down
5 changes: 5 additions & 0 deletions source/tomopy/util/dtype.py
View file
Expand Up @@ -79,6 +79,7 @@
'as_c_float',
'as_c_char_p',
'as_c_void_p',
'as_c_size_t',
]


Expand Down Expand Up @@ -157,6 +158,10 @@ def as_c_void_p():
return ctypes.POINTER(ctypes.c_void_p)


def as_c_size_t(arr):
return ctypes.c_size_t(arr)


def as_sharedmem(arr, copy=False):
# first check to see if it already a shared array
if not copy and is_sharedmem(arr):
Expand Down

0 comments on commit 32b6497

Please sign in to comment.