1 /*
2 * This file is part of RawTherapee.
3 *
4 * RawTherapee is free software: you can redistribute it and/or modify
5 * it under the terms of the GNU General Public License as published by
6 * the Free Software Foundation, either version 3 of the License, or
7 * (at your option) any later version.
8 *
9 * RawTherapee is distributed in the hope that it will be useful,
10 * but widthITheightOUT ANY widthARRANTY; without even the implied warranty of
11 * MERCheightANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with RawTherapee. If not, see <http://www.gnu.org/licenses/>.
16 *
17 * 2010 Emil Martinec <ejmartin@uchicago.edu>
18 *
19 */
20
21 #define SQR(x) ((x)*(x))
22
23 #include <cstddef>
24 #include <algorithm>
25
26 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27 // Gabor's implementation of bilateral filtering, without input pixel
28
29 #define NBRWT(a,b) (src[i - a][j - b] * ec[src[i - a][j - b]-src[i][j]+0x10000])
30 #define NORM(a,b) (1 + ec[src[i - a][j - b]-src[i][j]+0x10000])
31
32 #define RB_BEGIN(a,b) double scale = (a); \
33 int* ec = new int [0x20000]; \
34 for (int i=0; i<0x20000; i++) \
35 ec[i] = (int)(exp(-(double)(i-0x10000)*(double)(i-0x10000) / (2.0*rangewidth*rangewidth))*scale); \
36 int rstart = b; \
37 int rend = H-b; \
38 int cstart = b; \
39 int cend = W-b;
40
41 #define RB_END(b) buffer[i][j] = v; }} delete [] ec; \
42 for (int i=0; i<H; i++) \
43 for (int j=0; j<W; j++) \
44 if (i<rstart || j<cstart || i>=rend || j>=cend) \
45 dst[i][j] = src[i][j]; \
46 else \
47 dst[i][j] = buffer[i][j];
48
49 #define RB_OPER5 for (int i=rstart; i<rend; i++) { \
50 for (int j=cstart; j<cend; j++) { \
51 A v = NBRWT(-2,-2) + NBRWT(-2,-1) + NBRWT(-2,0) + NBRWT(-2,1) + NBRWT(-2,2) + \
52 NBRWT(-1,-2) + NBRWT(-1,-1) + NBRWT(-1,0) + NBRWT(-1,1) + NBRWT(-1,2) + \
53 NBRWT(0,-2) + NBRWT(0,-1) /*+ NBRWT(0,0)*/ + NBRWT(0,1) + NBRWT(0,2) + \
54 NBRWT(1,-2) + NBRWT(1,-1) + NBRWT(1,0) + NBRWT(1,1) + NBRWT(1,2) + \
55 NBRWT(2,-2) + NBRWT(2,-1) + NBRWT(2,0) + NBRWT(2,1) + NBRWT(2,2); \
56 v /= NORM(-2,-2) + NORM(-2,-1) + NORM(-2,0) + NORM(-2,1) + NORM(-2,2) + \
57 NORM(-1,-2) + NORM(-1,-1) + NORM(-1,0) + NORM(-1,1) + NORM(-1,2) + \
58 NORM(0,-2) + NORM(0,-1) /*+ NORM(0,0)*/ + NORM(0,1) + NORM(0,2) + \
59 NORM(1,-2) + NORM(1,-1) + NORM(1,0) + NORM(1,1) + NORM(1,2) + \
60 NORM(2,-2) + NORM(2,-1) + NORM(2,0) + NORM(2,1) + NORM(2,2);
61
62
rangeblur(T ** src,T ** dst,T ** buffer,int W,int H,double rangewidth,bool multiThread)63 template<class T, class A> void rangeblur (T** src, T** dst, T** buffer, int W, int H, double rangewidth, bool multiThread) {
64
65 RB_BEGIN(753,2)
66 #pragma omp parallel for if (multiThread)
67 RB_OPER5
68 RB_END(2)
69
70 }
71
72
73
impulse_nr(T ** src,T ** dst,int width,int height,double thresh)74 template<class T> void impulse_nr (T** src, T** dst, int width, int height, double thresh) {
75
76
77 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
78 // impulse noise removal
79 // local variables
80
81 float hpfabs, hfnbrave;
82
83 // buffer for the lowpass image
84 unsigned short ** lpf = new unsigned short *[height];
85 for (int i=0; i<height; i++) {
86 lpf[i] = new unsigned short [width];
87 //memset (lpf[i], 0, width*sizeof(float));
88 }
89
90 // buffer for the highpass image
91 unsigned short ** impish = new unsigned short *[height];
92 for (int i=0; i<height; i++) {
93 impish[i] = new unsigned short [width];
94 //memset (impish[i], 0, width*sizeof(unsigned short));
95 }
96
97 //The cleaning algorithm starts here
98
99 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100 // modified bilateral filter for lowpass image, omitting input pixel
101
102 static float eps = 1.0;
103 float wtdsum, dirwt, norm;
104 int i1, j1;
105
106 //rangeblur<unsigned short, unsigned int> (src, lpf, impish /*used as buffer here*/, width, height, thresh, false);
107
108 AlignedBuffer<double>* buffer = new AlignedBuffer<double> (MAX(width,height));
109
110 gaussHorizontal<unsigned short> (src, lpf, buffer, width, height, 2.0, false /*multiThread*/);
111 gaussVertical<unsigned short> (lpf, lpf, buffer, width, height, 2.0, false);
112
113 delete buffer;
114
115
116 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
117
118 for (int i=0; i < height; i++)
119 for (int j=0; j < width; j++) {
120
121 hpfabs = fabs(src[i][j]-lpf[i][j]);
122 //block average of high pass data
123 for (i1=MAX(0,i-2), hfnbrave=0; i1<=MIN(i+2,height-1); i1++ )
124 for (j1=MAX(0,j-2); j1<=MIN(j+2,width-1); j1++ ) {
125 hfnbrave += fabs(src[i1][j1]-lpf[i1][j1]);
126 }
127 hfnbrave = (hfnbrave-hpfabs)/24;
128 hpfabs>(hfnbrave*(5.5-thresh)) ? impish[i][j]=1 : impish[i][j]=0;
129
130 }//now impulsive values have been corrected
131
132 for (int i=0; i < height; i++)
133 for (int j=0; j < width; j++) {
134 if (!impish[i][j]) continue;
135 norm=0.0;
136 wtdsum=0.0;
137 for (i1=MAX(0,i-2), hfnbrave=0; i1<=MIN(i+2,height-1); i1++ )
138 for (j1=MAX(0,j-2); j1<=MIN(j+2,width-1); j1++ ) {
139 if (i1==i && j1==j) continue;
140 if (impish[i1][j1]) continue;
141 dirwt = 1/(SQR(src[i1][j1]-src[i][j])+eps);//use more sophisticated rangefn???
142 wtdsum += dirwt*src[i1][j1];
143 norm += dirwt;
144 }
145 //wtdsum /= norm;
146 if (norm) {
147 src[i][j]=wtdsum/norm;//low pass filter
148 }
149
150 }//now impulsive values have been corrected
151
152 for (int i=0; i<height; i++)
153 delete [] lpf[i];
154 delete [] lpf;
155
156 for (int i=0; i<height; i++)
157 delete [] impish[i];
158 delete [] impish;
159
160 }
161
162
163
164
165