1 /***
2 *
3 * Bayer CFA Demosaicing using Integrated Gaussian Vector on Color Differences
4 * Revision 1.0 - 2013/02/28
5 *
6 * Copyright (c) 2007-2013 Luis Sanz Rodriguez
7 * Using High Order Interpolation technique by Jim S, Jimmy Li, and Sharmil Randhawa
8 *
9 * Contact info: luis.sanz.rodriguez@gmail.com
10 *
11 * This code is distributed under a GNU General Public License, version 3.
12 * Visit <http://www.gnu.org/licenses/> for more information.
13 *
14 ***/
15 // Adapted to RT by Jacques Desmis 3/2013
16 // SSE version by Ingo Weyrich 5/2013
17 // Adapted to PhotoFlow 08/2014
18 #include <string.h>
19 #include <vector>
20 #include <array>
21
22 //#include "rtengine.h"
23 #include "rawimagesource.hh"
24 #include "rt_math.h"
25 //#include "../rtgui/multilangmgr.h"
26 //#include "procparams.h"
27 #include "sleef.c"
28 #include "opthelper.h"
29
30
31 //#undef CLIP
32 //#define CLIP(x) x
33
34
35
36
37 /**
38 * RATIO CORRECTED DEMOSAICING
39 * Luis Sanz Rodriguez (luis.sanz.rodriguez(at)gmail(dot)com)
40 *
41 * Release 2.2 @ 171117
42 *
43 * Original code from https://github.com/LuisSR/RCD-Demosaicing
44 * Adapted from RawTherapee: https://raw.githubusercontent.com/Beep6581/RawTherapee/rcd-demosaic/rtengine/demosaic_algos.cc
45 * Licensed under the GNU GPL version 3
46 */
47
48 #undef _OPENMP
49
50
51 //namespace rtengine {
52 namespace rtengine {
53
rcd_demosaic_RT(int winx,int winy,int winw,int winh,int tilex,int tiley,int tilew,int tileh)54 void RawImageSource::rcd_demosaic_RT(int winx, int winy, int winw, int winh,
55 int tilex, int tiley, int tilew, int tileh)
56 {
57 // Image dimensions
58 const int iwidth=winw, iheight=winh;
59 // Input tile dimensions
60 const int width=tilew, height=tileh;
61
62 std::vector<float> cfa(width * height);
63 std::vector< std::array<float, 3> > rgb(width * height);
64
65 int border = 8;
66 bool verbose = false;
67
68 #ifdef _OPENMP
69 #pragma omp parallel for
70 #endif
71 for (int row = 0, row2 = tiley; row < height; row++, row2++) {
72 for (int col = 0, col2 = tilex, indx = row * width + col;
73 col < width; col++, col2++, indx++) {
74 int c = FC(row, col);
75 cfa[indx] = rgb[indx][c] = CLIP(rawData[row2][col2]) / 65535.f;
76 }
77 }
78
79 if (plistener) {
80 plistener->setProgress(0.05);
81 }
82 // ------------------------------------------------------------------------
83 /* RT
84 int row, col, indx, c;
85 */
86 int w1 = width, w2 = 2 * width, w3 = 3 * width, w4 = 4 * width;
87
88 //Tolerance to avoid dividing by zero
89 static const float eps = 1e-5, epssq = 1e-10;
90
91
92 /* RT
93 //Gradients
94 float N_Grad, E_Grad, W_Grad, S_Grad, NW_Grad, NE_Grad, SW_Grad, SE_Grad;
95
96 //Pixel estimation
97 float N_Est, E_Est, W_Est, S_Est, NW_Est, NE_Est, SW_Est, SE_Est, V_Est, H_Est, P_Est, Q_Est;
98
99 //Directional discrimination
100 //float V_Stat, H_Stat, P_Stat, Q_Stat;
101 float VH_Central_Value, VH_Neighbour_Value, PQ_Central_Value, PQ_Neighbour_Value;
102 */
103 float ( *VH_Dir ), ( *VH_Disc ), ( *PQ_Dir ), ( *PQ_Disc );
104
105 //Low pass filter
106 float ( *lpf );
107
108
109 /**
110 * STEP 1: Find cardinal and diagonal interpolation directions
111 */
112
113 int x0 = 221, y0 = 509;
114
115 VH_Dir = ( float ( * ) ) calloc( width * height, sizeof *VH_Dir ); //merror ( VH_Dir, "rcd_demosaicing_171117()" );
116 PQ_Dir = ( float ( * ) ) calloc( width * height, sizeof *PQ_Dir ); //merror ( PQ_Dir, "rcd_demosaicing_171117()" );
117
118 #ifdef _OPENMP
119 #pragma omp parallel for
120 #endif
121 for (int row = 4; row < height - 4; row++ ) {
122 for (int col = 4, indx = row * width + col; col < width - 4; col++, indx++ ) {
123
124 //Calculate h/v local discrimination
125 float V_Stat = 0.f - 18.0f * cfa[indx] * cfa[indx - w1] - 18.0f * cfa[indx] * cfa[indx + w1] - 36.0f * cfa[indx] * cfa[indx - w2] - 36.0f * cfa[indx] * cfa[indx + w2] + 18.0f * cfa[indx] * cfa[indx - w3] + 18.0f * cfa[indx] * cfa[indx + w3] - 2.0f * cfa[indx] * cfa[indx - w4] - 2.0f * cfa[indx] * cfa[indx + w4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx - w1] * cfa[indx + w1] - 12.0f * cfa[indx - w1] * cfa[indx - w2] + 24.0f * cfa[indx - w1] * cfa[indx + w2] - 38.0f * cfa[indx - w1] * cfa[indx - w3] + 16.0f * cfa[indx - w1] * cfa[indx + w3] + 12.0f * cfa[indx - w1] * cfa[indx - w4] - 6.0f * cfa[indx - w1] * cfa[indx + w4] + 46.0f * cfa[indx - w1] * cfa[indx - w1] + 24.0f * cfa[indx + w1] * cfa[indx - w2] - 12.0f * cfa[indx + w1] * cfa[indx + w2] + 16.0f * cfa[indx + w1] * cfa[indx - w3] - 38.0f * cfa[indx + w1] * cfa[indx + w3] - 6.0f * cfa[indx + w1] * cfa[indx - w4] + 12.0f * cfa[indx + w1] * cfa[indx + w4] + 46.0f * cfa[indx + w1] * cfa[indx + w1] + 14.0f * cfa[indx - w2] * cfa[indx + w2] - 12.0f * cfa[indx - w2] * cfa[indx + w3] - 2.0f * cfa[indx - w2] * cfa[indx - w4] + 2.0f * cfa[indx - w2] * cfa[indx + w4] + 11.0f * cfa[indx - w2] * cfa[indx - w2] - 12.0f * cfa[indx + w2] * cfa[indx - w3] + 2.0f * cfa[indx + w2] * cfa[indx - w4] - 2.0f * cfa[indx + w2] * cfa[indx + w4] + 11.0f * cfa[indx + w2] * cfa[indx + w2] + 2.0f * cfa[indx - w3] * cfa[indx + w3] - 6.0f * cfa[indx - w3] * cfa[indx - w4] + 10.0f * cfa[indx - w3] * cfa[indx - w3] - 6.0f * cfa[indx + w3] * cfa[indx + w4] + 10.0f * cfa[indx + w3] * cfa[indx + w3] + 1.0f * cfa[indx - w4] * cfa[indx - w4] + 1.0f * cfa[indx + w4] * cfa[indx + w4];
126 float H_Stat = 0.f - 18.0f * cfa[indx] * cfa[indx - 1] - 18.0f * cfa[indx] * cfa[indx + 1] - 36.0f * cfa[indx] * cfa[indx - 2] - 36.0f * cfa[indx] * cfa[indx + 2] + 18.0f * cfa[indx] * cfa[indx - 3] + 18.0f * cfa[indx] * cfa[indx + 3] - 2.0f * cfa[indx] * cfa[indx - 4] - 2.0f * cfa[indx] * cfa[indx + 4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx - 1] * cfa[indx + 1] - 12.0f * cfa[indx - 1] * cfa[indx - 2] + 24.0f * cfa[indx - 1] * cfa[indx + 2] - 38.0f * cfa[indx - 1] * cfa[indx - 3] + 16.0f * cfa[indx - 1] * cfa[indx + 3] + 12.0f * cfa[indx - 1] * cfa[indx - 4] - 6.0f * cfa[indx - 1] * cfa[indx + 4] + 46.0f * cfa[indx - 1] * cfa[indx - 1] + 24.0f * cfa[indx + 1] * cfa[indx - 2] - 12.0f * cfa[indx + 1] * cfa[indx + 2] + 16.0f * cfa[indx + 1] * cfa[indx - 3] - 38.0f * cfa[indx + 1] * cfa[indx + 3] - 6.0f * cfa[indx + 1] * cfa[indx - 4] + 12.0f * cfa[indx + 1] * cfa[indx + 4] + 46.0f * cfa[indx + 1] * cfa[indx + 1] + 14.0f * cfa[indx - 2] * cfa[indx + 2] - 12.0f * cfa[indx - 2] * cfa[indx + 3] - 2.0f * cfa[indx - 2] * cfa[indx - 4] + 2.0f * cfa[indx - 2] * cfa[indx + 4] + 11.0f * cfa[indx - 2] * cfa[indx - 2] - 12.0f * cfa[indx + 2] * cfa[indx - 3] + 2.0f * cfa[indx + 2] * cfa[indx - 4] - 2.0f * cfa[indx + 2] * cfa[indx + 4] + 11.0f * cfa[indx + 2] * cfa[indx + 2] + 2.0f * cfa[indx - 3] * cfa[indx + 3] - 6.0f * cfa[indx - 3] * cfa[indx - 4] + 10.0f * cfa[indx - 3] * cfa[indx - 3] - 6.0f * cfa[indx + 3] * cfa[indx + 4] + 10.0f * cfa[indx + 3] * cfa[indx + 3] + 1.0f * cfa[indx - 4] * cfa[indx - 4] + 1.0f * cfa[indx + 4] * cfa[indx + 4];
127 V_Stat += epssq;
128 H_Stat += epssq;
129
130 VH_Dir[indx] = V_Stat / (fabs(V_Stat) + fabs(H_Stat));
131
132 //Calculate m/p local discrimination
133 float P_Stat = 0.f - 18.0f * cfa[indx] * cfa[indx - w1 - 1] - 18.0f * cfa[indx] * cfa[indx + w1 + 1] - 36.0f * cfa[indx] * cfa[indx - w2 - 2] - 36.0f * cfa[indx] * cfa[indx + w2 + 2] + 18.0f * cfa[indx] * cfa[indx - w3 - 3] + 18.0f * cfa[indx] * cfa[indx + w3 + 3] - 2.0f * cfa[indx] * cfa[indx - w4 - 4] - 2.0f * cfa[indx] * cfa[indx + w4 + 4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx - w1 - 1] * cfa[indx + w1 + 1] - 12.0f * cfa[indx - w1 - 1] * cfa[indx - w2 - 2] + 24.0f * cfa[indx - w1 - 1] * cfa[indx + w2 + 2] - 38.0f * cfa[indx - w1 - 1] * cfa[indx - w3 - 3] + 16.0f * cfa[indx - w1 - 1] * cfa[indx + w3 + 3] + 12.0f * cfa[indx - w1 - 1] * cfa[indx - w4 - 4] - 6.0f * cfa[indx - w1 - 1] * cfa[indx + w4 + 4] + 46.0f * cfa[indx - w1 - 1] * cfa[indx - w1 - 1] + 24.0f * cfa[indx + w1 + 1] * cfa[indx - w2 - 2] - 12.0f * cfa[indx + w1 + 1] * cfa[indx + w2 + 2] + 16.0f * cfa[indx + w1 + 1] * cfa[indx - w3 - 3] - 38.0f * cfa[indx + w1 + 1] * cfa[indx + w3 + 3] - 6.0f * cfa[indx + w1 + 1] * cfa[indx - w4 - 4] + 12.0f * cfa[indx + w1 + 1] * cfa[indx + w4 + 4] + 46.0f * cfa[indx + w1 + 1] * cfa[indx + w1 + 1] + 14.0f * cfa[indx - w2 - 2] * cfa[indx + w2 + 2] - 12.0f * cfa[indx - w2 - 2] * cfa[indx + w3 + 3] - 2.0f * cfa[indx - w2 - 2] * cfa[indx - w4 - 4] + 2.0f * cfa[indx - w2 - 2] * cfa[indx + w4 + 4] + 11.0f * cfa[indx - w2 - 2] * cfa[indx - w2 - 2] - 12.0f * cfa[indx + w2 + 2] * cfa[indx - w3 - 3] + 2 * cfa[indx + w2 + 2] * cfa[indx - w4 - 4] - 2.0f * cfa[indx + w2 + 2] * cfa[indx + w4 + 4] + 11.0f * cfa[indx + w2 + 2] * cfa[indx + w2 + 2] + 2.0f * cfa[indx - w3 - 3] * cfa[indx + w3 + 3] - 6.0f * cfa[indx - w3 - 3] * cfa[indx - w4 - 4] + 10.0f * cfa[indx - w3 - 3] * cfa[indx - w3 - 3] - 6.0f * cfa[indx + w3 + 3] * cfa[indx + w4 + 4] + 10.0f * cfa[indx + w3 + 3] * cfa[indx + w3 + 3] + 1.0f * cfa[indx - w4 - 4] * cfa[indx - w4 - 4] + 1.0f * cfa[indx + w4 + 4] * cfa[indx + w4 + 4];
134 float Q_Stat = 0.f - 18.0f * cfa[indx] * cfa[indx + w1 - 1] - 18.0f * cfa[indx] * cfa[indx - w1 + 1] - 36.0f * cfa[indx] * cfa[indx + w2 - 2] - 36.0f * cfa[indx] * cfa[indx - w2 + 2] + 18.0f * cfa[indx] * cfa[indx + w3 - 3] + 18.0f * cfa[indx] * cfa[indx - w3 + 3] - 2.0f * cfa[indx] * cfa[indx + w4 - 4] - 2.0f * cfa[indx] * cfa[indx - w4 + 4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx + w1 - 1] * cfa[indx - w1 + 1] - 12.0f * cfa[indx + w1 - 1] * cfa[indx + w2 - 2] + 24.0f * cfa[indx + w1 - 1] * cfa[indx - w2 + 2] - 38.0f * cfa[indx + w1 - 1] * cfa[indx + w3 - 3] + 16.0f * cfa[indx + w1 - 1] * cfa[indx - w3 + 3] + 12.0f * cfa[indx + w1 - 1] * cfa[indx + w4 - 4] - 6.0f * cfa[indx + w1 - 1] * cfa[indx - w4 + 4] + 46.0f * cfa[indx + w1 - 1] * cfa[indx + w1 - 1] + 24.0f * cfa[indx - w1 + 1] * cfa[indx + w2 - 2] - 12.0f * cfa[indx - w1 + 1] * cfa[indx - w2 + 2] + 16.0f * cfa[indx - w1 + 1] * cfa[indx + w3 - 3] - 38.0f * cfa[indx - w1 + 1] * cfa[indx - w3 + 3] - 6.0f * cfa[indx - w1 + 1] * cfa[indx + w4 - 4] + 12.0f * cfa[indx - w1 + 1] * cfa[indx - w4 + 4] + 46.0f * cfa[indx - w1 + 1] * cfa[indx - w1 + 1] + 14.0f * cfa[indx + w2 - 2] * cfa[indx - w2 + 2] - 12.0f * cfa[indx + w2 - 2] * cfa[indx - w3 + 3] - 2.0f * cfa[indx + w2 - 2] * cfa[indx + w4 - 4] + 2.0f * cfa[indx + w2 - 2] * cfa[indx - w4 + 4] + 11.0f * cfa[indx + w2 - 2] * cfa[indx + w2 - 2] - 12.0f * cfa[indx - w2 + 2] * cfa[indx + w3 - 3] + 2 * cfa[indx - w2 + 2] * cfa[indx + w4 - 4] - 2.0f * cfa[indx - w2 + 2] * cfa[indx - w4 + 4] + 11.0f * cfa[indx - w2 + 2] * cfa[indx - w2 + 2] + 2.0f * cfa[indx + w3 - 3] * cfa[indx - w3 + 3] - 6.0f * cfa[indx + w3 - 3] * cfa[indx + w4 - 4] + 10.0f * cfa[indx + w3 - 3] * cfa[indx + w3 - 3] - 6.0f * cfa[indx - w3 + 3] * cfa[indx - w4 + 4] + 10.0f * cfa[indx - w3 + 3] * cfa[indx - w3 + 3] + 1.0f * cfa[indx + w4 - 4] * cfa[indx + w4 - 4] + 1.0f * cfa[indx - w4 + 4] * cfa[indx - w4 + 4];
135 P_Stat += epssq;
136 Q_Stat += epssq;
137
138
139 PQ_Dir[indx] = P_Stat / ( fabs(P_Stat) + fabs(Q_Stat) );
140
141 //if(tilex<10 && tiley<10 && row<6 && col<6)
142 // printf("indx=%d VH_Dir=%f V_Stat=%e H_Stat=%e PQ_Dir=%f\n", indx, VH_Dir[indx], V_Stat, H_Stat, PQ_Dir[indx]);
143
144 }
145 }
146
147 // RT ---------------------------------------------------------------------
148 if (plistener) {
149 plistener->setProgress(0.2);
150 }
151 // -------------------------------------------------------------------------
152
153
154 VH_Disc = ( float ( * ) ) calloc( width * height, sizeof *VH_Disc ); //merror ( VH_Disc, "rcd_demosaicing_171117()" );
155 PQ_Disc = ( float ( * ) ) calloc( width * height, sizeof *PQ_Disc ); //merror ( PQ_Disc, "rcd_demosaicing_171117()" );
156
157 #ifdef _OPENMP
158 #pragma omp parallel for
159 #endif
160 for ( int row = 4; row < height - 4; row++ ) {
161 for ( int col = 4, indx = row * width + col; col < width - 4; col++, indx++ ) {
162
163 //Refined h/v local discrimination
164 float VH_Central_Value = VH_Dir[indx];
165 float VH_Neighbour_Value = 0.25f * (VH_Dir[indx - w1 - 1] + VH_Dir[indx - w1 + 1] + VH_Dir[indx + w1 - 1] + VH_Dir[indx + w1 + 1]);
166
167 VH_Disc[indx] = ( fabs( 0.5f - VH_Central_Value ) < fabs( 0.5f - VH_Neighbour_Value ) ) ? VH_Neighbour_Value : VH_Central_Value;
168
169 //Refined m/p local discrimination
170 float PQ_Central_Value = PQ_Dir[indx];
171 float PQ_Neighbour_Value = 0.25f * (PQ_Dir[indx - w1 - 1] + PQ_Dir[indx - w1 + 1] + PQ_Dir[indx + w1 - 1] + PQ_Dir[indx + w1 + 1]);
172
173 PQ_Disc[indx] = ( fabs( 0.5f - PQ_Central_Value ) < fabs( 0.5f - PQ_Neighbour_Value ) ) ? PQ_Neighbour_Value : PQ_Central_Value;
174
175 }
176 }
177
178 free( VH_Dir );
179 free( PQ_Dir );
180
181 // RT ---------------------------------------------------------------------
182 if (plistener) {
183 plistener->setProgress(0.4);
184 }
185 // ------------------------------------------------------------------------
186
187 /**
188 * STEP 2: Calculate the low pass filter
189 */
190
191 lpf = ( float ( * ) ) calloc( width * height, sizeof *lpf ); //merror ( lpf, "rcd_demosaicing_171117()" );
192
193 #ifdef _OPENMP
194 #pragma omp parallel for
195 #endif
196 for ( int row = 1; row < height - 1; row++ ) {
197 for ( int col = 1, indx = row * width + col; col < width - 1; col++, indx++ ) {
198
199 //Low pass filter incorporating red and blue local samples
200 lpf[indx] = 0.25f * cfa[indx] + 0.125f * ( cfa[indx - w1] + cfa[indx + w1] + cfa[indx - 1] + cfa[indx + 1] ) + 0.0625f * ( cfa[indx - w1 - 1] + cfa[indx - w1 + 1] + cfa[indx + w1 - 1] + cfa[indx + w1 + 1] );
201
202 }
203 }
204
205 // RT ---------------------------------------------------------------------
206 if (plistener) {
207 plistener->setProgress(0.5);
208 }
209 // ------------------------------------------------------------------------
210
211 /**
212 * STEP 3: Populate the green channel
213 */
214 #ifdef _OPENMP
215 #pragma omp parallel for
216 #endif
217 for ( int row = 4; row < height - 4; row++ ) {
218 for ( int col = 4 + ( FC( row, 0 )&1 ), indx = row * width + col; col < width - 4; col += 2, indx += 2 ) {
219
220 //Cardinal gradients
221 float N_Grad = eps + fabs( cfa[indx - w1] - cfa[indx + w1] ) + fabs( cfa[indx] - cfa[indx - w2] ) + fabs( cfa[indx - w1] - cfa[indx - w3] ) + fabs( cfa[indx - w2] - cfa[indx - w4] );
222 float S_Grad = eps + fabs( cfa[indx + w1] - cfa[indx - w1] ) + fabs( cfa[indx] - cfa[indx + w2] ) + fabs( cfa[indx + w1] - cfa[indx + w3] ) + fabs( cfa[indx + w2] - cfa[indx + w4] );
223 float W_Grad = eps + fabs( cfa[indx - 1] - cfa[indx + 1] ) + fabs( cfa[indx] - cfa[indx - 2] ) + fabs( cfa[indx - 1] - cfa[indx - 3] ) + fabs( cfa[indx - 2] - cfa[indx - 4] );
224 float E_Grad = eps + fabs( cfa[indx + 1] - cfa[indx - 1] ) + fabs( cfa[indx] - cfa[indx + 2] ) + fabs( cfa[indx + 1] - cfa[indx + 3] ) + fabs( cfa[indx + 2] - cfa[indx + 4] );
225
226 //Cardinal pixel estimations
227 float N_Est = cfa[indx - w1] * ( 1.f + ( lpf[indx] - lpf[indx - w2] ) / ( eps + lpf[indx] + lpf[indx - w2] ) );
228 float S_Est = cfa[indx + w1] * ( 1.f + ( lpf[indx] - lpf[indx + w2] ) / ( eps + lpf[indx] + lpf[indx + w2] ) );
229 float W_Est = cfa[indx - 1] * ( 1.f + ( lpf[indx] - lpf[indx - 2] ) / ( eps + lpf[indx] + lpf[indx - 2] ) );
230 float E_Est = cfa[indx + 1] * ( 1.f + ( lpf[indx] - lpf[indx + 2] ) / ( eps + lpf[indx] + lpf[indx + 2] ) );
231
232 //Interpolate G@R & G@B
233 float V_Est = ( S_Grad * N_Est + N_Grad * S_Est ) / ( N_Grad + S_Grad );
234 float H_Est = ( W_Grad * E_Est + E_Grad * W_Est ) / ( E_Grad + W_Grad );
235
236 rgb[indx][1] = LIM( VH_Disc[indx] * H_Est + ( 1.0f - VH_Disc[indx] ) * V_Est, 0.f, 1.f );
237
238 }
239 }
240
241 free( lpf );
242
243 // RT ---------------------------------------------------------------------
244 if (plistener) {
245 plistener->setProgress(0.7);
246 }
247 // -------------------------------------------------------------------------
248
249 /**
250 * STEP 4: Populate the red and blue channel
251 */
252 for ( int row = 4; row < height - 4; row++ ) {
253 for ( int col = 4 + ( FC( row, 0 )&1 ), indx = row * width + col, c = 2 - FC( row, col ); col < width - 4; col += 2, indx += 2 ) {
254
255 //Diagonal gradients
256 float NW_Grad = eps + fabs( rgb[indx - w1 - 1][c] - rgb[indx + w1 + 1][c] ) + fabs( rgb[indx - w1 - 1][c] - rgb[indx - w3 - 3][c] ) + fabs( rgb[indx][1] - rgb[indx - w2 - 2][1] );
257 float NE_Grad = eps + fabs( rgb[indx - w1 + 1][c] - rgb[indx + w1 - 1][c] ) + fabs( rgb[indx - w1 + 1][c] - rgb[indx - w3 + 3][c] ) + fabs( rgb[indx][1] - rgb[indx - w2 + 2][1] );
258 float SW_Grad = eps + fabs( rgb[indx + w1 - 1][c] - rgb[indx - w1 + 1][c] ) + fabs( rgb[indx + w1 - 1][c] - rgb[indx + w3 - 3][c] ) + fabs( rgb[indx][1] - rgb[indx + w2 - 2][1] );
259 float SE_Grad = eps + fabs( rgb[indx + w1 + 1][c] - rgb[indx - w1 - 1][c] ) + fabs( rgb[indx + w1 + 1][c] - rgb[indx + w3 + 3][c] ) + fabs( rgb[indx][1] - rgb[indx + w2 + 2][1] );
260
261 //Diagonal colour differences
262 float NW_Est = rgb[indx - w1 - 1][c] - rgb[indx - w1 - 1][1];
263 float NE_Est = rgb[indx - w1 + 1][c] - rgb[indx - w1 + 1][1];
264 float SW_Est = rgb[indx + w1 - 1][c] - rgb[indx + w1 - 1][1];
265 float SE_Est = rgb[indx + w1 + 1][c] - rgb[indx + w1 + 1][1];
266
267 //Interpolate R@B and B@R
268 float P_Est = ( NW_Grad * SE_Est + SE_Grad * NW_Est ) / ( NW_Grad + SE_Grad );
269 float Q_Est = ( NE_Grad * SW_Est + SW_Grad * NE_Est ) / ( NE_Grad + SW_Grad );
270
271 rgb[indx][c] = LIM( rgb[indx][1] + ( 1.0f - PQ_Disc[indx] ) * P_Est + PQ_Disc[indx] * Q_Est, 0.f, 1.f );
272
273 }
274 }
275
276 // RT ---------------------------------------------------------------------
277 if (plistener) {
278 plistener->setProgress(0.825);
279 }
280 // -------------------------------------------------------------------------
281
282 for ( int row = 4; row < height - 4; row++ ) {
283 for ( int col = 4 + ( FC( row, 1 )&1 ), indx = row * width + col; col < width - 4; col += 2, indx += 2 ) {
284
285 for ( int c = 0; c <= 2; c += 2 ) {
286
287 //Cardinal gradients
288 float N_Grad = eps + fabs( rgb[indx][1] - rgb[indx - w2][1] ) + fabs( rgb[indx - w1][c] - rgb[indx + w1][c] ) + fabs( rgb[indx - w1][c] - rgb[indx - w3][c] );
289 float S_Grad = eps + fabs( rgb[indx][1] - rgb[indx + w2][1] ) + fabs( rgb[indx + w1][c] - rgb[indx - w1][c] ) + fabs( rgb[indx + w1][c] - rgb[indx + w3][c] );
290 float W_Grad = eps + fabs( rgb[indx][1] - rgb[indx - 2][1] ) + fabs( rgb[indx - 1][c] - rgb[indx + 1][c] ) + fabs( rgb[indx - 1][c] - rgb[indx - 3][c] );
291 float E_Grad = eps + fabs( rgb[indx][1] - rgb[indx + 2][1] ) + fabs( rgb[indx + 1][c] - rgb[indx - 1][c] ) + fabs( rgb[indx + 1][c] - rgb[indx + 3][c] );
292
293 //Cardinal colour differences
294 float N_Est = rgb[indx - w1][c] - rgb[indx - w1][1];
295 float S_Est = rgb[indx + w1][c] - rgb[indx + w1][1];
296 float W_Est = rgb[indx - 1][c] - rgb[indx - 1][1];
297 float E_Est = rgb[indx + 1][c] - rgb[indx + 1][1];
298
299 //Interpolate R@G and B@G
300 float V_Est = ( N_Grad * S_Est + S_Grad * N_Est ) / ( N_Grad + S_Grad );
301 float H_Est = ( E_Grad * W_Est + W_Grad * E_Est ) / ( E_Grad + W_Grad );
302
303 rgb[indx][c] = LIM( rgb[indx][1] + ( 1.0f - VH_Disc[indx] ) * V_Est + VH_Disc[indx] * H_Est, 0.f, 1.f );
304
305 }
306 }
307 }
308
309 free( PQ_Disc );
310 free( VH_Disc );
311
312 // RT ---------------------------------------------------------------------
313 if (plistener) {
314 plistener->setProgress(0.95);
315 }
316
317 #ifdef _OPENMP
318 #pragma omp parallel for
319 #endif
320 for (int row = border, row2 = tiley+row; row < height-border; ++row, ++row2) {
321 for (int col = border, col2 = tilex+col, idx = row * width + col;
322 col < width-border; ++col, ++col2, ++idx) {
323 red[row2][col2] = CLIP(rgb[idx][0] * 65535.f);
324 green[row2][col2] = CLIP(rgb[idx][1] * 65535.f);
325 blue[row2][col2] = CLIP(rgb[idx][2] * 65535.f);
326 }
327 }
328
329 if (plistener) {
330 plistener->setProgress(1);
331 }
332 // -------------------------------------------------------------------------
333 }
334
335 }
336