summaryrefslogtreecommitdiff
path: root/silx/math/medianfilter/include/median_filter.hpp
diff options
context:
space:
mode:
Diffstat (limited to 'silx/math/medianfilter/include/median_filter.hpp')
-rw-r--r--silx/math/medianfilter/include/median_filter.hpp217
1 files changed, 129 insertions, 88 deletions
diff --git a/silx/math/medianfilter/include/median_filter.hpp b/silx/math/medianfilter/include/median_filter.hpp
index 78d6c5a..b4d953a 100644
--- a/silx/math/medianfilter/include/median_filter.hpp
+++ b/silx/math/medianfilter/include/median_filter.hpp
@@ -1,6 +1,6 @@
/*##########################################################################
#
-# Copyright (c) 2017 European Synchrotron Radiation Facility
+# 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
@@ -33,70 +33,74 @@
#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
+ SHRINK=3,
+ CONSTANT=4,
};
-// Simple function browsing a deque and registring the min and max values
+// Simple function browsing a deque and registering the min and max values
// and if those values are unique or not
template<typename T>
-void getMinMax(std::vector<const T*>& v, T& min, T&max,
- typename std::vector<const T*>::const_iterator end){
+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<const T*>::const_iterator it = v.begin();
+ typename std::vector<T>::const_iterator it = v.begin();
if (v.size() == 0){
raise(SIGINT);
}else{
- min = max = *(*it);
+ min = max = *it;
}
it++;
// Browse all the deque
while(it!=end){
// check if repeated (should always be before min/max setting)
- if(*(*it) > max) max = *(*it);
- if(*(*it) < min) min = *(*it);
+ T value = *it;
+ if(value > max) max = value;
+ if(value < min) min = value;
it++;
}
}
-template<typename T>
-bool cmp(const T* a, const T* b){
- return *a < *b;
-}
// 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>
-const T* median(std::vector<const T*>& v, int window_size) {
- std::nth_element(v.begin(), v.begin() + window_size/2, v.begin()+window_size, cmp<T>);
- return v[window_size/2];
+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];
}
-template<typename T>
-void print_window(std::vector<const T*>& v,
- typename std::vector<const T*>::const_iterator end){
- typename std::vector<const T*>::const_iterator it;
- for(it = v.begin(); it != end; ++it){
- std::cout << *(*it) << " ";
- }
- std::cout << std::endl;
-}
// return the index into 0, (length_max - 1) in reflect mode
-int reflect(int index, int length_max){
+inline int reflect(int index, int length_max){
int res = index;
- // if the index is negative get the positive symetrical value
+ // if the index is negative get the positive symmetrical value
if(res < 0){
res += 1;
res = -res;
}
- // then apply the reflect algorithm. Frequence is 2 max length
+ // then apply the reflect algorithm. Frequency is 2 max length
res = res % (2*length_max);
if(res >= length_max){
res = 2*length_max - res -1;
@@ -106,13 +110,12 @@ int reflect(int index, int length_max){
}
// return the index into 0, (length_max - 1) in mirror mode
-int mirror(int index, int length_max){
+inline int mirror(int index, int length_max){
int res = index;
- // if the index is negative get the positive symetrical value
+ // if the index is negative get the positive symmetrical value
if(res < 0){
res = -res;
}
-
int rightLimit = length_max -1;
// apply the redundancy each two right limit
res = res % (2*rightLimit);
@@ -126,7 +129,7 @@ int mirror(int index, int length_max){
// Browse the column of pixel_x
template<typename T>
void median_filter(
- const T* input,
+ const T* input,
T* output,
int* kernel_dim, // two values : 0:width, 1:height
int* image_dim, // two values : 0:width, 1:height
@@ -134,7 +137,8 @@ void median_filter(
int x_pixel_range_min,
int x_pixel_range_max,
bool conditional,
- int pMode){
+ int pMode,
+ T cval) {
assert(kernel_dim[0] > 0);
assert(kernel_dim[1] > 0);
@@ -156,72 +160,109 @@ void median_filter(
MODE mode = static_cast<MODE>(pMode);
// init buffer
- std::vector<const T*> window_values(kernel_dim[0]*kernel_dim[1]);
+ 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<const T*>::iterator it = window_values.begin();
+ typename std::vector<T>::iterator it = window_values.begin();
// fill the vector
- 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++)
+
+ 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++)
{
- 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);
- break;
-
- case REFLECT:
- index_x = reflect(win_x, image_dim[1]);
- index_y = reflect(win_y, image_dim[0]);
- break;
-
- case MIRROR:
- index_x = mirror(win_x, image_dim[1]);
- index_y = mirror(win_y, image_dim[0]);
- break;
- case SHRINK:
- if((index_x < 0) || (index_x > image_dim[1] -1)){
- continue;
- }
- if((index_y < 0) || (index_y > image_dim[0] -1)){
- continue;
- }
- break;
+ 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]);
+ // deal with 1d case
+ if(win_y == 0 && image_dim[0] == 1){
+ index_y = 0;
+ }else{
+ 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;
+ }
}
- *it = (&input[index_y*image_dim[1] + index_x]);
- ++it;
}
}
- // get end of the windows. This is needed since in shrink mode we are
- // not sure to fill the entire window.
- typename std::vector<const T*>::iterator window_end;
- int window_size = kernel_dim[0]*kernel_dim[1];
- if(mode == SHRINK){
- int x_shrink_ker_dim = std::min(x_pixel+halfKernel_x, image_dim[1]-1) - std::max(0, x_pixel-halfKernel_x)+1;
- int y_shrink_ker_dim = std::min(y_pixel+halfKernel_y, image_dim[0]-1) - std::max(0, y_pixel-halfKernel_y)+1;
- window_size = x_shrink_ker_dim*y_shrink_ker_dim;
- window_end = window_values.begin() + window_size;
- }else{
- window_end = window_values.end();
- }
+ //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);
- // apply the median value if needed for this pixel
- const T* currentPixelValue = &input[image_dim[1]*y_pixel + x_pixel];
- if (conditional == true){
- T min = 0;
- T max = 0;
- getMinMax(window_values, min, max, window_end);
- if ((*currentPixelValue == max) || (*currentPixelValue == min)){
- output[image_dim[1]*y_pixel + x_pixel] = *(median<T>(window_values, window_size));
+ 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] = *currentPixelValue;
+ output[image_dim[1]*y_pixel + x_pixel] = median<T>(window_values, window_size);
}
- }else{
- output[image_dim[1]*y_pixel + x_pixel] = *(median<T>(window_values, window_size));
}
}
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-04-25 15:23:39 +0000