1 ////////////////////////////////////////////////////////////////
2 //
3 // Green Equilibration via directional average
4 //
5 // copyright (c) 2008-2010 Emil Martinec <ejmartin@uchicago.edu>
6 // optimized for speed 2017 Ingo Weyrich <heckflosse67@gmx.de>
7 //
8 //
9 // code dated: August 25, 2017
10 //
11 // green_equil_RT.cc is free software: you can redistribute it and/or modify
12 // it under the terms of the GNU General Public License as published by
13 // the Free Software Foundation, either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // This program is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU General Public License
22 // along with this program. If not, see <https://www.gnu.org/licenses/>.
23 //
24 ////////////////////////////////////////////////////////////////
25
26 #include <cmath>
27 #include <cstdlib>
28 #include <ctime>
29
30 #include "rt_math.h"
31 #include "rawimagesource.h"
32 #include "opthelper.h"
33
34 namespace rtengine
35 {
36
green_equilibrate_global(array2D<float> & rawData)37 void RawImageSource::green_equilibrate_global(array2D<float> &rawData)
38 {
39 // global correction
40 int ng1 = 0, ng2 = 0;
41 double avgg1 = 0., avgg2 = 0.;
42
43 #ifdef _OPENMP
44 #pragma omp parallel for reduction(+: ng1, ng2, avgg1, avgg2) schedule(dynamic,16)
45 #endif
46
47 for (int i = border; i < H - border; i++) {
48 double avgg = 0.;
49 for (int j = border + ((FC(i, border) & 1) ^ 1); j < W - border; j += 2) {
50 avgg += rawData[i][j];
51 }
52
53 int ng = (W - 2 * border + (FC(i, border) & 1)) / 2;
54
55 if (i & 1) {
56 avgg2 += avgg;
57 ng2 += ng;
58 } else {
59 avgg1 += avgg;
60 ng1 += ng;
61 }
62 }
63
64 // Avoid division by zero
65 if(ng1 == 0 || avgg1 == 0.0) {
66 ng1 = 1;
67 avgg1 = 1.0;
68 }
69 if(ng2 == 0 || avgg2 == 0.0) {
70 ng2 = 1;
71 avgg2 = 1.0;
72 }
73
74 double corrg1 = (avgg1 / ng1 + avgg2 / ng2) / 2.0 / (avgg1 / ng1);
75 double corrg2 = (avgg1 / ng1 + avgg2 / ng2) / 2.0 / (avgg2 / ng2);
76
77 #ifdef _OPENMP
78 #pragma omp parallel for schedule(dynamic,16)
79 #endif
80
81 for (int i = border; i < H - border; i++) {
82 double corrg = (i & 1) ? corrg2 : corrg1;
83
84 for (int j = border + ((FC(i, border) & 1) ^ 1); j < W - border; j += 2) {
85 rawData[i][j] *= corrg;
86 }
87 }
88 }
89
90 //void green_equilibrate()//for dcraw implementation
green_equilibrate(const GreenEqulibrateThreshold & thresh,array2D<float> & rawData)91 void RawImageSource::green_equilibrate(const GreenEqulibrateThreshold &thresh, array2D<float> &rawData)
92 {
93 // thresh = threshold for performing green equilibration; max percentage difference of G1 vs G2
94 // G1-G2 differences larger than this will be assumed to be Nyquist texture, and left untouched
95
96 int height = H, width = W;
97
98 // local variables
99 array2D<float> cfa(width / 2 + (width & 1), height);
100
101 #ifdef _OPENMP
102 #pragma omp parallel for schedule(dynamic,16)
103 #endif
104
105 for (int i = 0; i < height; ++i) {
106 int j = (FC(i, 0) & 1) ^ 1;
107 #ifdef __SSE2__
108
109 for (; j < width - 7; j += 8) {
110 STVFU(cfa[i][j >> 1], LC2VFU(rawData[i][j]));
111 }
112
113 #endif
114
115 for (; j < width; j += 2) {
116 cfa[i][j >> 1] = rawData[i][j];
117 }
118 }
119
120 constexpr float eps = 1.f; //tolerance to avoid dividing by zero
121 // const float thresh6 = 6 * thresh;
122 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
123
124 // Fill G interpolated values with border interpolation and input values
125
126 //int vote1, vote2;
127 //int counter, vtest;
128
129 //The green equilibration algorithm starts here
130 //now smooth the cfa data
131 #ifdef _OPENMP
132 #pragma omp parallel
133 #endif
134 {
135 #ifdef __SSE2__
136 vfloat zd5v = F2V(0.5f);
137 vfloat onev = F2V(1.f);
138 // vfloat threshv = F2V(thresh);
139 // vfloat thresh6v = F2V(thresh6);
140 vfloat epsv = F2V(eps);
141 #endif
142 #ifdef _OPENMP
143 #pragma omp for schedule(dynamic,16)
144 #endif
145
146 for (int rr = 4; rr < height - 4; rr++) {
147 int cc = 5 - (FC(rr, 2) & 1);
148 #ifdef __SSE2__
149
150 for (; cc < width - 12; cc += 8) {
151 //neighbour checking code from Manuel Llorens Garcia
152 vfloat o1_1 = LVFU(cfa[rr - 1][(cc - 1) >> 1]);
153 vfloat o1_2 = LVFU(cfa[rr - 1][(cc + 1) >> 1]);
154 vfloat o1_3 = LVFU(cfa[rr + 1][(cc - 1) >> 1]);
155 vfloat o1_4 = LVFU(cfa[rr + 1][(cc + 1) >> 1]);
156 vfloat o2_1 = LVFU(cfa[rr - 2][cc >> 1]);
157 vfloat o2_2 = LVFU(cfa[rr + 2][cc >> 1]);
158 vfloat o2_3 = LVFU(cfa[rr][(cc >> 1) - 1]);
159 vfloat o2_4 = LVFU(cfa[rr][(cc >> 1) + 1]);
160
161 vfloat d1 = (o1_1 + o1_2 + o1_3 + o1_4);
162 vfloat d2 = (o2_1 + o2_2 + o2_3 + o2_4);
163
164 vfloat c1 = (vabsf(o1_1 - o1_2) + vabsf(o1_1 - o1_3) + vabsf(o1_1 - o1_4) + vabsf(o1_2 - o1_3) + vabsf(o1_3 - o1_4) + vabsf(o1_2 - o1_4));
165 vfloat c2 = (vabsf(o2_1 - o2_2) + vabsf(o2_1 - o2_3) + vabsf(o2_1 - o2_4) + vabsf(o2_2 - o2_3) + vabsf(o2_3 - o2_4) + vabsf(o2_2 - o2_4));
166
167 vfloat tfv;
168 for (int k = 0; k < 4; ++k) {
169 tfv[k] = thresh(rr, cc + 2 * k);
170 }
171 vfloat tf6v = F2V(6.f) * tfv;
172
173 vmask mask1 = vmaskf_lt(c1 + c2, tf6v * vabsf(d1 - d2));
174
175 if (_mm_movemask_ps((vfloat)mask1)) { // if for any of the 4 pixels the condition is true, do the maths for all 4 pixels and mask the unused out at the end
176 //pixel interpolation
177 vfloat gin = LVFU(cfa[rr][cc >> 1]);
178
179 vfloat gmp2p2 = gin - LVFU(cfa[rr + 2][(cc >> 1) + 1]);
180 vfloat gmm2m2 = gin - LVFU(cfa[rr - 2][(cc >> 1) - 1]);
181 vfloat gmm2p2 = gin - LVFU(cfa[rr - 2][(cc >> 1) + 1]);
182 vfloat gmp2m2 = gin - LVFU(cfa[rr + 2][(cc >> 1) - 1]);
183
184 vfloat gse = o1_4 + zd5v * gmp2p2;
185 vfloat gnw = o1_1 + zd5v * gmm2m2;
186 vfloat gne = o1_2 + zd5v * gmm2p2;
187 vfloat gsw = o1_3 + zd5v * gmp2m2;
188
189 vfloat wtse = onev / (epsv + SQRV(gmp2p2) + SQRV(LVFU(cfa[rr + 3][(cc + 3) >> 1]) - o1_4));
190 vfloat wtnw = onev / (epsv + SQRV(gmm2m2) + SQRV(LVFU(cfa[rr - 3][(cc - 3) >> 1]) - o1_1));
191 vfloat wtne = onev / (epsv + SQRV(gmm2p2) + SQRV(LVFU(cfa[rr - 3][(cc + 3) >> 1]) - o1_2));
192 vfloat wtsw = onev / (epsv + SQRV(gmp2m2) + SQRV(LVFU(cfa[rr + 3][(cc - 3) >> 1]) - o1_3));
193
194 vfloat ginterp = (gse * wtse + gnw * wtnw + gne * wtne + gsw * wtsw) / (wtse + wtnw + wtne + wtsw);
195
196 vfloat val = vself(vmaskf_lt(ginterp - gin, tfv * (ginterp + gin)), zd5v * (ginterp + gin), gin);
197 val = vself(mask1, val, gin);
198 STC2VFU(rawData[rr][cc], val);
199 }
200 }
201
202 #endif
203
204 for (; cc < width - 6; cc += 2) {
205 //neighbour checking code from Manuel Llorens Garcia
206 float o1_1 = cfa[rr - 1][(cc - 1) >> 1];
207 float o1_2 = cfa[rr - 1][(cc + 1) >> 1];
208 float o1_3 = cfa[rr + 1][(cc - 1) >> 1];
209 float o1_4 = cfa[rr + 1][(cc + 1) >> 1];
210 float o2_1 = cfa[rr - 2][cc >> 1];
211 float o2_2 = cfa[rr + 2][cc >> 1];
212 float o2_3 = cfa[rr][(cc - 2) >> 1];
213 float o2_4 = cfa[rr][(cc + 2) >> 1];
214
215 float d1 = (o1_1 + o1_2) + (o1_3 + o1_4);
216 float d2 = (o2_1 + o2_2) + (o2_3 + o2_4);
217
218 float c1 = (fabs(o1_1 - o1_2) + fabs(o1_1 - o1_3) + fabs(o1_1 - o1_4) + fabs(o1_2 - o1_3) + fabs(o1_3 - o1_4) + fabs(o1_2 - o1_4));
219 float c2 = (fabs(o2_1 - o2_2) + fabs(o2_1 - o2_3) + fabs(o2_1 - o2_4) + fabs(o2_2 - o2_3) + fabs(o2_3 - o2_4) + fabs(o2_2 - o2_4));
220
221 float tf = thresh(rr, cc);
222
223 if (c1 + c2 < 6 * tf * fabs(d1 - d2)) {
224 //pixel interpolation
225 float gin = cfa[rr][cc >> 1];
226
227 float gmp2p2 = gin - cfa[rr + 2][(cc + 2) >> 1];
228 float gmm2m2 = gin - cfa[rr - 2][(cc - 2) >> 1];
229 float gmm2p2 = gin - cfa[rr - 2][(cc + 2) >> 1];
230 float gmp2m2 = gin - cfa[rr + 2][(cc - 2) >> 1];
231
232 float gse = o1_4 + 0.5f * gmp2p2;
233 float gnw = o1_1 + 0.5f * gmm2m2;
234 float gne = o1_2 + 0.5f * gmm2p2;
235 float gsw = o1_3 + 0.5f * gmp2m2;
236
237 float wtse = 1.f / (eps + SQR(gmp2p2) + SQR(cfa[rr + 3][(cc + 3) >> 1] - o1_4));
238 float wtnw = 1.f / (eps + SQR(gmm2m2) + SQR(cfa[rr - 3][(cc - 3) >> 1] - o1_1));
239 float wtne = 1.f / (eps + SQR(gmm2p2) + SQR(cfa[rr - 3][(cc + 3) >> 1] - o1_2));
240 float wtsw = 1.f / (eps + SQR(gmp2m2) + SQR(cfa[rr + 3][(cc - 3) >> 1] - o1_3));
241
242 float ginterp = (gse * wtse + gnw * wtnw + gne * wtne + gsw * wtsw) / (wtse + wtnw + wtne + wtsw);
243
244 if (ginterp - gin < tf * (ginterp + gin)) {
245 rawData[rr][cc] = 0.5f * (ginterp + gin);
246 }
247 }
248 }
249 }
250 }
251 }
252 }
253