1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
|
/*
* Replacing pixels by local median value when detected as outliers based on MAD (2D box, variable size)
* This file is part of ufo-serge filter set.
* Copyright (C) 2016 Serge Cohen
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
* Serge Cohen <serge.cohen@synchrotron-soleil.fr>
*/
// The size of the box in which to perform the filtering
// This one should be provided at compile time through a -D directive.
// #define BOXSIZE 15
// The threshold to apply on the MAD to reject (or not) the pixel's value
typedef int ssize_t;
constant const ssize_t half_box = (BOXSIZE-1)>>1;
__kernel void
med_mad_rej_2D (
__global float *iImage,
__global float *oFilteredImage,
const float iRejThresh
)
{
size_t sizeX = get_global_size(0);
size_t sizeY = get_global_size(1);
size_t globalX = get_global_id(0);
size_t globalY = get_global_id(1);
// Getting the element of an image (both lvalue and rvalue) from index :
#define FromImage(x,y,image) image[(x) + (y)*sizeX]
/* Getting all the value within the box and in the image : */
float v[BOXSIZE*BOXSIZE];
ssize_t cX, cY;
size_t index=0;
for ( cX=-half_box + globalX; (half_box + globalX + 1) != cX; ++cX) {
for ( cY=-half_box + globalY; (half_box + globalY + 1) != cY; ++cY) {
if ( (0 <= cX) && (sizeX > cX) && (0 <= cY) && (sizeY > cY) ) {
v[index] = FromImage(cX, cY, iImage);
++index;
}
}
}
size_t num_px = index;
float swapper;
/* Computing the median : */
for (index = 1; num_px != index; ++index) {
for (size_t j = 0; j!= (num_px-index); ++j ) {
swapper = v[j];
v[j] = min(v[j], v[j+1]);
v[j+1] = max(swapper, v[j+1]);
}
}
float med = v[num_px>>1];
/* Computing the MAD */
for ( index = 0; (num_px>>1) != index; ++index) {
v[index] = med - v[index];
}
for ( index = (num_px>>1) ; num_px != index; ++index) {
v[index] = v[index] - med;
}
for (index = 1; num_px != index; ++index) {
for (size_t j = 0; j!= (num_px-index); ++j ) {
swapper = v[j];
v[j] = min(v[j], v[j+1]);
v[j+1] = max(swapper, v[j+1]);
}
}
float mad = v[num_px>>1];
/* Testing the current value : */
if ( fabs(FromImage(globalX, globalY, iImage) - med) < (mad*iRejThresh) )
FromImage(globalX, globalY, oFilteredImage) = FromImage(globalX, globalY, iImage);
else
FromImage(globalX, globalY, oFilteredImage) = med;
}
|
