/usr/share/ufo/med-mad-reject-2d.cl is in ufo-filters-data 0.14.1+dfsg1-1.1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
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 = 0; num_px != (1+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 = 0; num_px != (1+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;
}
|