Commit 11faf398 authored by Carsten Richter's avatar Carsten Richter

Merge branch 'median-filter' into 'master'

Make median filter works with mask

Closes #40

See merge request !66
parents d536d417 d31d75a0
Pipeline #4405 passed with stages
in 4 minutes and 7 seconds
......@@ -6,7 +6,7 @@ __pycache__/
*.so
# Cython generated files
xsocs/util/filt_utils/filt_utils.c
xsocs/util/filt_utils/medianfilter.cpp
xsocs/util/filt_utils/histogramnd_lut.c
# Distribution / packaging
......
......@@ -128,6 +128,7 @@ def run_entry_point(entry_point, argv):
target_name = elements[0].strip()
elements = elements[1].split(":")
module_name = elements[0].strip()
# Take care of entry_point optional "extra" requirements declaration
function_name = elements[1].split()[0].strip()
logger.info("Execute target %s (function %s from module %s) using importlib", target_name, function_name, module_name)
......
......@@ -201,11 +201,9 @@ class ConversionParamsWidget(Qt.QWidget):
imgGboxLayout.addRow(self.__medfiltCBox, inputBase)
# Mask and binning/median filer are mutually exclusive
# Mask and binning mutually exclusive
self.__imgMaskCBox.toggled.connect(self.__imgBinCBox.setDisabled)
self.__imgBinCBox.toggled.connect(self.__imgMaskCBox.setDisabled)
self.__imgMaskCBox.toggled.connect(self.__medfiltCBox.setDisabled)
self.__medfiltCBox.toggled.connect(self.__imgMaskCBox.setDisabled)
layout.addWidget(imageGbox)
......
......@@ -41,9 +41,7 @@ import multiprocessing.sharedctypes as mp_sharedctypes
import numpy as np
import xrayutilities as xu
# from scipy.signal import medfilt2d
from ...util.filt_utils import medfilt2D
from ...util.medianfilter import medfilt2d
from ...util.histogramnd_lut import histogramnd_get_lut, histogramnd_from_lut
# from silx.math import histogramnd
from ...io import XsocsH5, QSpaceH5, ShiftH5
......@@ -712,10 +710,6 @@ class QSpaceConverter(object):
not np.all(np.equal(image_binning, (1, 1)))):
raise ValueError(
'Image binning is not implemented with mask')
if (medfilt_dims is not None and
not np.all(np.equal(medfilt_dims, (1, 1)))):
raise ValueError(
'Median filter is not implemented with mask')
shiftH5 = self.__shiftH5
......@@ -1398,6 +1392,9 @@ def _to_qspace(th_idx,
# sum_axis_2 = 1
avg_weight = 1. / (image_binning[0] * image_binning[1])
# Set OpenMP to use a single thread (for median filter)
os.environ["OMP_NUM_THREADS"] = "1"
rc = None
errMsg = None
try:
......@@ -1487,9 +1484,10 @@ def _to_qspace(th_idx,
# intensity = medfilt2d(intensity, 3)
if medfilt_dims[0] != 1 or medfilt_dims[1] != 1:
intensity = medfilt2D(intensity,
kernel=medfilt_dims,
n_threads=None)
intensity = medfilt2d(intensity,
medfilt_dims,
mode='constant',
cval=0)
t_medfilt += time.time() - t0
t0 = time.time()
......
This diff is collapsed.
/*##########################################################################
#
# Copyright (c) 2017-2018 European Synchrotron Radiation Facility
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ###########################################################################*/
// __authors__ = ["H. Payno"]
// __license__ = "MIT"
// __date__ = "10/02/2017"
#ifndef MEDIAN_FILTER
#define MEDIAN_FILTER
#include <vector>
#include <assert.h>
#include <algorithm>
#include <signal.h>
#include <iostream>
#include <cmath>
#include <cfloat>
/* Needed for pytohn2.7 on Windows... */
#ifndef INFINITY
#define INFINITY (DBL_MAX+DBL_MAX)
#endif
#ifndef NAN
#define NAN (INFINITY-INFINITY)
#endif
// Modes for the median filter
enum MODE{
NEAREST=0,
REFLECT=1,
MIRROR=2,
SHRINK=3,
CONSTANT=4,
};
// Simple function browsing a deque and registring the min and max values
// and if those values are unique or not
template<typename T>
void getMinMax(std::vector<T>& v, T& min, T&max,
typename std::vector<T>::const_iterator end){
// init min and max values
typename std::vector<T>::const_iterator it = v.begin();
if (v.size() == 0){
raise(SIGINT);
}else{
min = max = *it;
}
it++;
// Browse all the deque
while(it!=end){
// check if repeated (should always be before min/max setting)
T value = *it;
if(value > max) max = value;
if(value < min) min = value;
it++;
}
}
// apply the median filter only on limited part of the vector
// In case of even number of elements (either due to NaNs in the window
// or for image borders in shrink mode):
// the highest of the 2 central values is returned
template<typename T>
inline T median(std::vector<T>& v, int window_size) {
int pivot = window_size / 2;
std::nth_element(v.begin(), v.begin() + pivot, v.begin()+window_size);
return v[pivot];
}
// return the index into 0, (length_max - 1) in reflect mode
inline int reflect(int index, int length_max){
int res = index;
// if the index is negative get the positive symetrical value
if(res < 0){
res += 1;
res = -res;
}
// then apply the reflect algorithm. Frequence is 2 max length
res = res % (2*length_max);
if(res >= length_max){
res = 2*length_max - res -1;
res = res % length_max;
}
return res;
}
// return the index into 0, (length_max - 1) in mirror mode
inline int mirror(int index, int length_max){
int res = index;
// if the index is negative get the positive symetrical value
if(res < 0){
res = -res;
}
int rightLimit = length_max -1;
// apply the redundancy each two right limit
res = res % (2*rightLimit);
if(res >= length_max){
int distToRedundancy = (2*rightLimit) - res;
res = distToRedundancy;
}
return res;
}
// Browse the column of pixel_x
template<typename T>
void median_filter(
const T* input,
T* output,
int* kernel_dim, // two values : 0:width, 1:height
int* image_dim, // two values : 0:width, 1:height
int y_pixel, // the x pixel to process
int x_pixel_range_min,
int x_pixel_range_max,
bool conditional,
int pMode,
T cval) {
assert(kernel_dim[0] > 0);
assert(kernel_dim[1] > 0);
assert(y_pixel >= 0);
assert(image_dim[0] > 0);
assert(image_dim[1] > 0);
assert(y_pixel >= 0);
assert(y_pixel < image_dim[0]);
assert(x_pixel_range_max < image_dim[1]);
assert(x_pixel_range_min <= x_pixel_range_max);
// kernel odd assertion
assert((kernel_dim[0] - 1)%2 == 0);
assert((kernel_dim[1] - 1)%2 == 0);
// # this should be move up to avoid calculation each time
int halfKernel_x = (kernel_dim[1] - 1) / 2;
int halfKernel_y = (kernel_dim[0] - 1) / 2;
MODE mode = static_cast<MODE>(pMode);
// init buffer
std::vector<T> window_values(kernel_dim[0]*kernel_dim[1]);
bool not_horizontal_border = (y_pixel >= halfKernel_y && y_pixel < image_dim[0] - halfKernel_y);
for(int x_pixel=x_pixel_range_min; x_pixel <= x_pixel_range_max; x_pixel ++ ){
typename std::vector<T>::iterator it = window_values.begin();
// fill the vector
if (not_horizontal_border &&
x_pixel >= halfKernel_x && x_pixel < image_dim[1] - halfKernel_x) {
//This is not a border, just fill it
for(int win_y=y_pixel-halfKernel_y; win_y<= y_pixel+halfKernel_y; win_y++) {
for(int win_x = x_pixel-halfKernel_x; win_x <= x_pixel+halfKernel_x; win_x++){
T value = input[win_y*image_dim[1] + win_x];
if (value == value) { // Ignore NaNs
*it = value;
++it;
}
}
}
} else { // This is a border, handle the special case
for(int win_y=y_pixel-halfKernel_y; win_y<= y_pixel+halfKernel_y; win_y++)
{
for(int win_x = x_pixel-halfKernel_x; win_x <= x_pixel+halfKernel_x; win_x++)
{
T value = 0;
int index_x = win_x;
int index_y = win_y;
switch(mode){
case NEAREST:
index_x = std::min(std::max(win_x, 0), image_dim[1] - 1);
index_y = std::min(std::max(win_y, 0), image_dim[0] - 1);
value = input[index_y*image_dim[1] + index_x];
break;
case REFLECT:
index_x = reflect(win_x, image_dim[1]);
index_y = reflect(win_y, image_dim[0]);
value = input[index_y*image_dim[1] + index_x];
break;
case MIRROR:
index_x = mirror(win_x, image_dim[1]);
index_y = mirror(win_y, image_dim[0]);
value = input[index_y*image_dim[1] + index_x];
break;
case SHRINK:
if ((index_x < 0) || (index_x > image_dim[1] -1) ||
(index_y < 0) || (index_y > image_dim[0] -1)) {
continue;
}
value = input[index_y*image_dim[1] + index_x];
break;
case CONSTANT:
if ((index_x < 0) || (index_x > image_dim[1] -1) ||
(index_y < 0) || (index_y > image_dim[0] -1)) {
value = cval;
} else {
value = input[index_y*image_dim[1] + index_x];
}
break;
}
if (value == value) { // Ignore NaNs
*it = value;
++it;
}
}
}
}
//window_size can be smaller than kernel size in shrink mode or if there is NaNs
int window_size = std::distance(window_values.begin(), it);
if (window_size == 0) {
// Window is empty, this is the case when all values are NaNs
output[image_dim[1]*y_pixel + x_pixel] = NAN;
} else {
// apply the median value if needed for this pixel
const T currentPixelValue = input[image_dim[1]*y_pixel + x_pixel];
if (conditional == true){
typename std::vector<T>::iterator window_end = window_values.begin() + window_size;
T min = 0;
T max = 0;
getMinMax(window_values, min, max, window_end);
// NaNs are propagated through unchanged
if ((currentPixelValue == max) || (currentPixelValue == min)){
output[image_dim[1]*y_pixel + x_pixel] = median<T>(window_values, window_size);
}else{
output[image_dim[1]*y_pixel + x_pixel] = currentPixelValue;
}
}else{
output[image_dim[1]*y_pixel + x_pixel] = median<T>(window_values, window_size);
}
}
}
}
#endif // MEDIAN_FILTER
# coding: utf-8
# /*##########################################################################
#
# Copyright (c) 2015-2018 European Synchrotron Radiation Facility
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#
# ###########################################################################*/
from libcpp cimport bool
# pyx
cdef extern from "median_filter.hpp":
cdef extern void median_filter[T](const T* image,
T* output,
int* kernel_dim,
int* image_dim,
int x_pixel_range_min,
int x_pixel_range_max,
int y_pixel_range_min,
int y_pixel_range_max,
bool conditional,
T cval) nogil;
cdef extern int reflect(int index, int length_max);
cdef extern int mirror(int index, int length_max);
This diff is collapsed.
......@@ -34,30 +34,18 @@ from numpy.distutils.misc_util import Configuration
def configuration(parent_package='', top_path=None):
config = Configuration('util', parent_package, top_path)
# =====================================
# filt_utils
# =====================================
filt_utils_dir = 'filt_utils'
filt_utils_src = [os.path.join(filt_utils_dir, srcf)
for srcf in ['filt_utils.pyx']]
filt_utils_inc = [numpy.get_include()]
config.add_extension('filt_utils',
sources=filt_utils_src,
include_dirs=filt_utils_inc,
language='c')
# =====================================
# histogramnd_lut
# =====================================
filt_utils_dir = 'filt_utils'
histogramnd_src = [os.path.join(filt_utils_dir, srcf)
for srcf in ['histogramnd_lut.pyx']]
histogramnd_inc = [numpy.get_include()]
config.add_extension('medianfilter',
sources=[os.path.join(filt_utils_dir,
'medianfilter.pyx')],
include_dirs=['.', numpy.get_include()],
language='c++')
config.add_extension('histogramnd_lut',
sources=histogramnd_src,
include_dirs=histogramnd_inc,
sources=[os.path.join(filt_utils_dir,
'histogramnd_lut.pyx')],
include_dirs=[numpy.get_include()],
language='c')
return config
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment