PageRenderTime 101ms CodeModel.GetById 24ms app.highlight 66ms RepoModel.GetById 1ms app.codeStats 0ms

/src/FreeImage/Source/LibRawLite/internal/dcb_demosaicing.c

https://bitbucket.org/cabalistic/ogredeps/
C | 710 lines | 433 code | 211 blank | 66 comment | 74 complexity | 5f561b30ea41dd4f99f371f6f27d6c61 MD5 | raw file
  1/*
  2 *    Copyright (C) 2010,  Jacek Gozdz (cuniek@kft.umcs.lublin.pl)
  3 * 
  4 *    This code is licensed under a (3-clause) BSD license as follows :
  5 *
  6 *    Redistribution and use in source and binary forms, with or without
  7 *    modification, are permitted provided that the following 
  8 *	  conditions are met:
  9 *     
 10 *    * Redistributions of source code must retain the above copyright
 11 *      notice, this list of conditions and the following disclaimer.
 12 *    * Redistributions in binary form must reproduce the above
 13 *      copyright notice, this list of conditions and the following 
 14 *		disclaimer in the documentation and/or other materials provided 
 15 * 	    with the distribution.
 16 *    * Neither the name of the author nor the names of its
 17 *      contributors may be used to endorse or promote products 
 18 * 		derived from this software without specific prior written permission.
 19 *     
 20 *    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 21 *    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 22 *    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
 23 * 	  FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL 
 24 * 	  THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, 
 25 * 	  INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES 
 26 * 	  (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR 
 27 * 	  SERVICES; LOSS OF USE,
 28 *    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND 
 29 * 	  ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, 
 30 *	  OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
 31 * 	  OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY 
 32 *    OF SUCH DAMAGE.
 33 */
 34
 35// DCB demosaicing by Jacek Gozdz (cuniek@kft.umcs.lublin.pl)
 36
 37// FBDD denoising by Jacek Gozdz (cuniek@kft.umcs.lublin.pl) and 
 38// Luis Sanz Rodríguez (luis.sanz.rodriguez@gmail.com)
 39 
 40// last modification: 11.07.2010 
 41 
 42
 43// interpolates green vertically and saves it to image3
 44void CLASS dcb_ver(float (*image3)[3])
 45{
 46	int row, col, u=width, v=2*u, indx;
 47
 48	for (row=2; row < height-2; row++)
 49		for (col=2+(FC(row,2)&1),indx=row*width+col; col < u-2; col+=2,indx+=2) {
 50	
 51			image3[indx][1] = CLIP((image[indx+u][1] + image[indx-u][1])/2.0);
 52
 53	}	
 54}
 55
 56
 57// interpolates green horizontally and saves it to image2
 58void CLASS dcb_hor(float (*image2)[3])
 59{
 60	int row, col, u=width, v=2*u, indx;
 61	
 62	for (row=2; row < height-2; row++)
 63		for (col=2+(FC(row,2)&1),indx=row*width+col; col < u-2; col+=2,indx+=2) {
 64			
 65			image2[indx][1] = CLIP((image[indx+1][1] + image[indx-1][1])/2.0);
 66
 67	}	
 68}
 69	
 70
 71
 72// missing colors are interpolated
 73void CLASS dcb_color()
 74{
 75	int row, col, c, d, u=width, indx;
 76
 77
 78	for (row=1; row < height-1; row++)
 79		for (col=1+(FC(row,1) & 1), indx=row*width+col, c=2-FC(row,col); col < u-1; col+=2, indx+=2) {
 80
 81			
 82			image[indx][c] = CLIP(( 
 83			4*image[indx][1] 
 84			- image[indx+u+1][1] - image[indx+u-1][1] - image[indx-u+1][1] - image[indx-u-1][1] 
 85			+ image[indx+u+1][c] + image[indx+u-1][c] + image[indx-u+1][c] + image[indx-u-1][c] )/4.0);
 86		}
 87
 88	for (row=1; row<height-1; row++)
 89		for (col=1+(FC(row,2) & 1), indx=row*width+col,c=FC(row,col+1),d=2-c; col<width-1; col+=2, indx+=2) {
 90			
 91			image[indx][c] = CLIP((2*image[indx][1] - image[indx+1][1] - image[indx-1][1] + image[indx+1][c] + image[indx-1][c])/2.0);
 92			image[indx][d] = CLIP((2*image[indx][1] - image[indx+u][1] - image[indx-u][1] + image[indx+u][d] + image[indx-u][d])/2.0);
 93		}	
 94}
 95
 96
 97// missing R and B are interpolated horizontally and saved in image2
 98void CLASS dcb_color2(float (*image2)[3])
 99{
100	int row, col, c, d, u=width, indx;
101
102
103	for (row=1; row < height-1; row++)
104		for (col=1+(FC(row,1) & 1), indx=row*width+col, c=2-FC(row,col); col < u-1; col+=2, indx+=2) {
105
106			
107			image2[indx][c] = CLIP(( 
108			4*image2[indx][1] 
109			- image2[indx+u+1][1] - image2[indx+u-1][1] - image2[indx-u+1][1] - image2[indx-u-1][1] 
110			+ image[indx+u+1][c] + image[indx+u-1][c] + image[indx-u+1][c] + image[indx-u-1][c] )/4.0);
111		}
112
113	for (row=1; row<height-1; row++)
114		for (col=1+(FC(row,2) & 1), indx=row*width+col,c=FC(row,col+1),d=2-c; col<width-1; col+=2, indx+=2) {
115			
116			image2[indx][c] = CLIP((image[indx+1][c] + image[indx-1][c])/2.0);
117			image2[indx][d] = CLIP((2*image2[indx][1] - image2[indx+u][1] - image2[indx-u][1] + image[indx+u][d] + image[indx-u][d])/2.0);
118		}	
119}
120
121
122// missing R and B are interpolated vertically and saved in image3
123void CLASS dcb_color3(float (*image3)[3])
124{
125	int row, col, c, d, u=width, indx;
126
127
128	for (row=1; row < height-1; row++)
129		for (col=1+(FC(row,1) & 1), indx=row*width+col, c=2-FC(row,col); col < u-1; col+=2, indx+=2) {
130
131			
132			image3[indx][c] = CLIP(( 
133			4*image3[indx][1] 
134			- image3[indx+u+1][1] - image3[indx+u-1][1] - image3[indx-u+1][1] - image3[indx-u-1][1] 
135			+ image[indx+u+1][c] + image[indx+u-1][c] + image[indx-u+1][c] + image[indx-u-1][c] )/4.0);
136		}
137
138	for (row=1; row<height-1; row++)
139		for (col=1+(FC(row,2) & 1), indx=row*width+col,c=FC(row,col+1),d=2-c; col<width-1; col+=2, indx+=2) {
140			
141			image3[indx][c] = CLIP((2*image3[indx][1] - image3[indx+1][1] - image3[indx-1][1] + image[indx+1][c] + image[indx-1][c])/2.0);
142			image3[indx][d] = CLIP((image[indx+u][d] + image[indx-u][d])/2.0);
143		}	
144}
145
146
147// decides the primary green interpolation direction
148void CLASS dcb_decide(float (*image2)[3], float (*image3)[3])
149{
150	int row, col, c, d, u=width, v=2*u, indx;
151	float current, current2, current3;
152
153	for (row=2; row < height-2; row++)
154		for (col=2+(FC(row,2)&1),indx=row*width+col, c=FC(row,col); col < u-2; col+=2,indx+=2) {
155	
156		d=ABS(c-2);
157
158   current = MAX(image[indx+v][c], MAX(image[indx-v][c], MAX(image[indx-2][c], image[indx+2][c]))) -
159			 MIN(image[indx+v][c], MIN(image[indx-v][c], MIN(image[indx-2][c], image[indx+2][c]))) +
160			 MAX(image[indx+1+u][d], MAX(image[indx+1-u][d], MAX(image[indx-1+u][d], image[indx-1-u][d]))) -
161			 MIN(image[indx+1+u][d], MIN(image[indx+1-u][d], MIN(image[indx-1+u][d], image[indx-1-u][d])));
162		
163  current2 = MAX(image2[indx+v][d], MAX(image2[indx-v][d], MAX(image2[indx-2][d], image2[indx+2][d]))) -
164			 MIN(image2[indx+v][d], MIN(image2[indx-v][d], MIN(image2[indx-2][d], image2[indx+2][d]))) +
165			 MAX(image2[indx+1+u][c], MAX(image2[indx+1-u][c], MAX(image2[indx-1+u][c], image2[indx-1-u][c]))) -
166			 MIN(image2[indx+1+u][c], MIN(image2[indx+1-u][c], MIN(image2[indx-1+u][c], image2[indx-1-u][c])));
167			  
168  current3 = MAX(image3[indx+v][d], MAX(image3[indx-v][d], MAX(image3[indx-2][d], image3[indx+2][d]))) -
169			 MIN(image3[indx+v][d], MIN(image3[indx-v][d], MIN(image3[indx-2][d], image3[indx+2][d]))) +
170			 MAX(image3[indx+1+u][c], MAX(image3[indx+1-u][c], MAX(image3[indx-1+u][c], image3[indx-1-u][c]))) -
171			 MIN(image3[indx+1+u][c], MIN(image3[indx+1-u][c], MIN(image3[indx-1+u][c], image3[indx-1-u][c])));
172
173		
174			if (ABS(current-current2) < ABS(current-current3))
175				image[indx][1] = image2[indx][1];
176			else
177				image[indx][1] = image3[indx][1];				
178			
179		
180	}
181}
182
183
184// saves red and blue in image2
185void CLASS dcb_copy_to_buffer(float (*image2)[3])
186{
187	int indx;
188
189	for (indx=0; indx < height*width; indx++) {
190		image2[indx][0]=image[indx][0]; //R
191		image2[indx][2]=image[indx][2]; //B
192	}
193}
194
195
196
197// restores red and blue from image2
198void CLASS dcb_restore_from_buffer(float (*image2)[3])
199{
200	int indx;
201
202	for (indx=0; indx < height*width; indx++) {
203		image[indx][0]=image2[indx][0]; //R
204		image[indx][2]=image2[indx][2]; //B
205	}
206}
207
208
209// R and B smoothing using green contrast, all pixels except 2 pixel wide border
210void CLASS dcb_pp()
211{
212	int g1, r1, b1, u=width, indx, row, col;
213
214	
215	for (row=2; row < height-2; row++) 
216	for (col=2, indx=row*u+col; col < width-2; col++, indx++) { 
217
218		r1 = ( image[indx-1][0] + image[indx+1][0] + image[indx-u][0] + image[indx+u][0] + image[indx-u-1][0] + image[indx+u+1][0] + image[indx-u+1][0] + image[indx+u-1][0])/8.0;
219		g1 = ( image[indx-1][1] + image[indx+1][1] + image[indx-u][1] + image[indx+u][1] + image[indx-u-1][1] + image[indx+u+1][1] + image[indx-u+1][1] + image[indx+u-1][1])/8.0;
220		b1 = ( image[indx-1][2] + image[indx+1][2] + image[indx-u][2] + image[indx+u][2] + image[indx-u-1][2] + image[indx+u+1][2] + image[indx-u+1][2] + image[indx+u-1][2])/8.0;
221		 
222		image[indx][0] = CLIP(r1 + ( image[indx][1] - g1 ));
223		image[indx][2] = CLIP(b1 + ( image[indx][1] - g1 ));
224	
225	}
226}
227
228
229
230// green blurring correction, helps to get the nyquist right
231void CLASS dcb_nyquist()
232{
233	int row, col, c, u=width, v=2*u, indx;
234
235	
236	for (row=2; row < height-2; row++)
237		for (col=2+(FC(row,2)&1),indx=row*width+col, c=FC(row,col); col < u-2; col+=2,indx+=2) {
238
239			image[indx][1] = CLIP((image[indx+v][1] + image[indx-v][1] + image[indx-2][1] + image[indx+2][1])/4.0 + 
240							  image[indx][c] - ( image[indx+v][c] + image[indx-v][c] + image[indx-2][c] + image[indx+2][c])/4.0);
241	
242	}	
243
244}
245
246
247
248
249
250// missing colors are interpolated using high quality algorithm by Luis Sanz Rodríguez
251void CLASS dcb_color_full()
252{
253	int row,col,c,d,i,j,u=width,v=2*u,w=3*u,indx, g1, g2;
254	float f[4],g[4],(*chroma)[2];
255
256	chroma = (float (*)[2]) calloc(width*height,sizeof *chroma); merror (chroma, "dcb_color_full()");
257
258	for (row=1; row < height-1; row++)
259		for (col=1+(FC(row,1)&1),indx=row*width+col,c=FC(row,col),d=c/2; col < u-1; col+=2,indx+=2)
260			chroma[indx][d]=image[indx][c]-image[indx][1];
261
262	for (row=3; row<height-3; row++)
263		for (col=3+(FC(row,1)&1),indx=row*width+col,c=1-FC(row,col)/2,d=1-c; col<u-3; col+=2,indx+=2) {
264			f[0]=1.0/(float)(1.0+fabs(chroma[indx-u-1][c]-chroma[indx+u+1][c])+fabs(chroma[indx-u-1][c]-chroma[indx-w-3][c])+fabs(chroma[indx+u+1][c]-chroma[indx-w-3][c]));
265			f[1]=1.0/(float)(1.0+fabs(chroma[indx-u+1][c]-chroma[indx+u-1][c])+fabs(chroma[indx-u+1][c]-chroma[indx-w+3][c])+fabs(chroma[indx+u-1][c]-chroma[indx-w+3][c]));
266			f[2]=1.0/(float)(1.0+fabs(chroma[indx+u-1][c]-chroma[indx-u+1][c])+fabs(chroma[indx+u-1][c]-chroma[indx+w+3][c])+fabs(chroma[indx-u+1][c]-chroma[indx+w-3][c]));
267			f[3]=1.0/(float)(1.0+fabs(chroma[indx+u+1][c]-chroma[indx-u-1][c])+fabs(chroma[indx+u+1][c]-chroma[indx+w-3][c])+fabs(chroma[indx-u-1][c]-chroma[indx+w+3][c]));
268			g[0]=1.325*chroma[indx-u-1][c]-0.175*chroma[indx-w-3][c]-0.075*chroma[indx-w-1][c]-0.075*chroma[indx-u-3][c];
269			g[1]=1.325*chroma[indx-u+1][c]-0.175*chroma[indx-w+3][c]-0.075*chroma[indx-w+1][c]-0.075*chroma[indx-u+3][c];
270			g[2]=1.325*chroma[indx+u-1][c]-0.175*chroma[indx+w-3][c]-0.075*chroma[indx+w-1][c]-0.075*chroma[indx+u-3][c];
271			g[3]=1.325*chroma[indx+u+1][c]-0.175*chroma[indx+w+3][c]-0.075*chroma[indx+w+1][c]-0.075*chroma[indx+u+3][c];
272			chroma[indx][c]=(f[0]*g[0]+f[1]*g[1]+f[2]*g[2]+f[3]*g[3])/(f[0]+f[1]+f[2]+f[3]);
273		}
274	for (row=3; row<height-3; row++)
275		for (col=3+(FC(row,2)&1),indx=row*width+col,c=FC(row,col+1)/2; col<u-3; col+=2,indx+=2)
276			for(d=0;d<=1;c=1-c,d++){
277				f[0]=1.0/(float)(1.0+fabs(chroma[indx-u][c]-chroma[indx+u][c])+fabs(chroma[indx-u][c]-chroma[indx-w][c])+fabs(chroma[indx+u][c]-chroma[indx-w][c]));
278				f[1]=1.0/(float)(1.0+fabs(chroma[indx+1][c]-chroma[indx-1][c])+fabs(chroma[indx+1][c]-chroma[indx+3][c])+fabs(chroma[indx-1][c]-chroma[indx+3][c]));
279				f[2]=1.0/(float)(1.0+fabs(chroma[indx-1][c]-chroma[indx+1][c])+fabs(chroma[indx-1][c]-chroma[indx-3][c])+fabs(chroma[indx+1][c]-chroma[indx-3][c]));
280				f[3]=1.0/(float)(1.0+fabs(chroma[indx+u][c]-chroma[indx-u][c])+fabs(chroma[indx+u][c]-chroma[indx+w][c])+fabs(chroma[indx-u][c]-chroma[indx+w][c]));
281			
282				g[0]=0.875*chroma[indx-u][c]+0.125*chroma[indx-w][c];
283				g[1]=0.875*chroma[indx+1][c]+0.125*chroma[indx+3][c];
284				g[2]=0.875*chroma[indx-1][c]+0.125*chroma[indx-3][c];
285				g[3]=0.875*chroma[indx+u][c]+0.125*chroma[indx+w][c];				
286
287				chroma[indx][c]=(f[0]*g[0]+f[1]*g[1]+f[2]*g[2]+f[3]*g[3])/(f[0]+f[1]+f[2]+f[3]);
288			}
289
290	for(row=6; row<height-6; row++)
291		for(col=6,indx=row*width+col; col<width-6; col++,indx++){
292			image[indx][0]=CLIP(chroma[indx][0]+image[indx][1]);
293			image[indx][2]=CLIP(chroma[indx][1]+image[indx][1]);
294			
295	g1 = MIN(image[indx+1+u][0], MIN(image[indx+1-u][0], MIN(image[indx-1+u][0], MIN(image[indx-1-u][0], MIN(image[indx-1][0], MIN(image[indx+1][0], MIN(image[indx-u][0], image[indx+u][0])))))));
296
297	g2 = MAX(image[indx+1+u][0], MAX(image[indx+1-u][0], MAX(image[indx-1+u][0], MAX(image[indx-1-u][0], MAX(image[indx-1][0], MAX(image[indx+1][0], MAX(image[indx-u][0], image[indx+u][0])))))));
298
299
300	image[indx][0] =  ULIM(image[indx][0], g2, g1);
301	
302
303	
304	g1 = MIN(image[indx+1+u][2], MIN(image[indx+1-u][2], MIN(image[indx-1+u][2], MIN(image[indx-1-u][2], MIN(image[indx-1][2], MIN(image[indx+1][2], MIN(image[indx-u][2], image[indx+u][2])))))));
305
306	g2 = MAX(image[indx+1+u][2], MAX(image[indx+1-u][2], MAX(image[indx-1+u][2], MAX(image[indx-1-u][2], MAX(image[indx-1][2], MAX(image[indx+1][2], MAX(image[indx-u][2], image[indx+u][2])))))));
307
308	image[indx][2] =  ULIM(image[indx][2], g2, g1);
309				
310			
311			
312		}
313
314	free(chroma);
315}
316
317
318
319
320// green is used to create an interpolation direction map saved in image[][3]
321// 1 = vertical
322// 0 = horizontal
323void CLASS dcb_map()
324{	
325	int current, row, col, c, u=width, v=2*u, indx;
326
327	for (row=1; row < height-1; row++) {
328	for (col=1, indx=row*width+col; col < width-1; col++, indx++) { 
329
330		if (image[indx][1] > ( image[indx-1][1] + image[indx+1][1] + image[indx-u][1] + image[indx+u][1])/4.0)
331			image[indx][3] = ((MIN( image[indx-1][1], image[indx+1][1]) + image[indx-1][1] + image[indx+1][1]) < 
332							  (MIN( image[indx-u][1], image[indx+u][1]) + image[indx-u][1] + image[indx+u][1]));   
333		else
334			image[indx][3] = ((MAX( image[indx-1][1], image[indx+1][1]) + image[indx-1][1] + image[indx+1][1]) > 
335							  (MAX( image[indx-u][1], image[indx+u][1]) + image[indx-u][1] + image[indx+u][1])); 
336	}
337	}
338}
339
340
341// interpolated green pixels are corrected using the map
342void CLASS dcb_correction()
343{
344	int current, row, col, u=width, v=2*u, indx;
345
346	for (row=2; row < height-2; row++)
347		for (col=2+(FC(row,2)&1),indx=row*width+col; col < u-2; col+=2,indx+=2) {
348
349			current = 4*image[indx][3] + 
350				      2*(image[indx+u][3] + image[indx-u][3] + image[indx+1][3] + image[indx-1][3]) + 
351					    image[indx+v][3] + image[indx-v][3] + image[indx+2][3] + image[indx-2][3];
352						
353			image[indx][1] = ((16-current)*(image[indx-1][1] + image[indx+1][1])/2.0 + current*(image[indx-u][1] + image[indx+u][1])/2.0)/16.0;
354
355	}
356
357}
358
359
360// interpolated green pixels are corrected using the map
361// with contrast correction
362void CLASS dcb_correction2()
363{
364	int current, row, col, c, u=width, v=2*u, indx;
365	ushort (*pix)[4];
366
367	for (row=4; row < height-4; row++)
368		for (col=4+(FC(row,2)&1),indx=row*width+col, c=FC(row,col); col < u-4; col+=2,indx+=2) {
369			
370			current = 4*image[indx][3] + 
371				      2*(image[indx+u][3] + image[indx-u][3] + image[indx+1][3] + image[indx-1][3]) + 
372					    image[indx+v][3] + image[indx-v][3] + image[indx+2][3] + image[indx-2][3];
373						
374			image[indx][1] = CLIP(((16-current)*((image[indx-1][1] + image[indx+1][1])/2.0 + image[indx][c] - (image[indx+2][c] + image[indx-2][c])/2.0) + current*((image[indx-u][1] + image[indx+u][1])/2.0 + image[indx][c] - (image[indx+v][c] + image[indx-v][c])/2.0))/16.0);			   
375		
376	}
377
378}
379
380
381void CLASS dcb_refinement()
382{
383	int row, col, c, u=width, v=2*u, w=3*u, indx, current;
384	float f[5], g1, g2, tmp, tmp2=0, tmp3=0;
385	
386	for (row=4; row < height-4; row++)
387		for (col=4+(FC(row,2)&1),indx=row*width+col, c=FC(row,col); col < u-4; col+=2,indx+=2) {
388
389			current = 4*image[indx][3] + 
390				      2*(image[indx+u][3] + image[indx-u][3] + image[indx+1][3] + image[indx-1][3]) 
391				      +image[indx+v][3] + image[indx-v][3] + image[indx-2][3] + image[indx+2][3];
392
393if (image[indx][c] > 1)
394{
395
396f[0] = (float)(image[indx-u][1] + image[indx+u][1])/(2*image[indx][c]);
397
398
399if (image[indx-v][c] > 0)
400f[1] = 2*(float)image[indx-u][1]/(image[indx-v][c] + image[indx][c]);
401else
402f[1] = f[0];
403
404if (image[indx-v][c] > 0)
405f[2] = (float)(image[indx-u][1] + image[indx-w][1])/(2*image[indx-v][c]);
406else
407f[2] = f[0];
408
409if (image[indx+v][c] > 0)
410f[3] = 2*(float)image[indx+u][1]/(image[indx+v][c] + image[indx][c]);
411else
412f[3] = f[0];
413
414if (image[indx+v][c] > 0)
415f[4] = (float)(image[indx+u][1] + image[indx+w][1])/(2*image[indx+v][c]);
416else
417f[4] = f[0];
418
419g1 = (5*f[0] + 3*f[1] + f[2] + 3*f[3] + f[4])/13.0;
420
421
422
423f[0] = (float)(image[indx-1][1] + image[indx+1][1])/(2*image[indx][c]);
424
425if (image[indx-2][c] > 0)
426f[1] = 2*(float)image[indx-1][1]/(image[indx-2][c] + image[indx][c]);
427else
428f[1] = f[0];
429
430if (image[indx-2][c] > 0)
431f[2] = (float)(image[indx-1][1] + image[indx-3][1])/(2*image[indx-2][c]);
432else
433f[2] = f[0];
434
435if (image[indx+2][c] > 0)
436f[3] = 2*(float)image[indx+1][1]/(image[indx+2][c] + image[indx][c]);
437else
438f[3] = f[0];
439
440if (image[indx+2][c] > 0)
441f[4] = (float)(image[indx+1][1] + image[indx+3][1])/(2*image[indx+2][c]);
442else
443f[4] = f[0];
444
445g2 = (5*f[0] + 3*f[1] + f[2] + 3*f[3] + f[4])/13.0;
446
447
448image[indx][1] = CLIP((image[indx][c])*(current*g1 + (16-current)*g2)/16.0);
449}
450else
451image[indx][1] = image[indx][c];
452			
453	// get rid of overshooted pixels
454
455	g1 = MIN(image[indx+1+u][1], MIN(image[indx+1-u][1], MIN(image[indx-1+u][1], MIN(image[indx-1-u][1], MIN(image[indx-1][1], MIN(image[indx+1][1], MIN(image[indx-u][1], image[indx+u][1])))))));
456
457	g2 = MAX(image[indx+1+u][1], MAX(image[indx+1-u][1], MAX(image[indx-1+u][1], MAX(image[indx-1-u][1], MAX(image[indx-1][1], MAX(image[indx+1][1], MAX(image[indx-u][1], image[indx+u][1])))))));
458
459	image[indx][1] =  ULIM(image[indx][1], g2, g1);
460	
461	}	
462
463
464
465}
466
467// converts RGB to LCH colorspace and saves it to image3
468void CLASS rgb_to_lch(double (*image2)[3])
469{
470	int indx;
471	for (indx=0; indx < height*width; indx++) {
472
473            image2[indx][0] = image[indx][0] + image[indx][1] + image[indx][2]; 		// L
474            image2[indx][1] = 1.732050808 *(image[indx][0] - image[indx][1]);			// C
475            image2[indx][2] = 2.0*image[indx][2] - image[indx][0] - image[indx][1];		// H
476	}
477}
478
479// converts LCH to RGB colorspace and saves it back to image
480void CLASS lch_to_rgb(double (*image2)[3])
481{
482	int indx;
483	for (indx=0; indx < height*width; indx++) {
484
485            image[indx][0] = CLIP(image2[indx][0] / 3.0 - image2[indx][2] / 6.0 + image2[indx][1] / 3.464101615);
486            image[indx][1] = CLIP(image2[indx][0] / 3.0 - image2[indx][2] / 6.0 - image2[indx][1] / 3.464101615);
487            image[indx][2] = CLIP(image2[indx][0] / 3.0 + image2[indx][2] / 3.0);
488	}
489}
490
491
492
493
494// denoising using interpolated neighbours
495void CLASS fbdd_correction()
496{
497	int row, col, c, u=width, indx;
498	ushort (*pix)[4];
499
500	for (row=2; row < height-2; row++) {
501	for (col=2, indx=row*width+col; col < width-2; col++, indx++) { 	
502
503		c =  fc(row,col);
504
505		image[indx][c] = ULIM(image[indx][c], 
506					MAX(image[indx-1][c], MAX(image[indx+1][c], MAX(image[indx-u][c], image[indx+u][c]))), 
507					MIN(image[indx-1][c], MIN(image[indx+1][c], MIN(image[indx-u][c], image[indx+u][c]))));
508
509	} 
510	}	
511}
512
513
514// corrects chroma noise
515void CLASS fbdd_correction2(double (*image2)[3])
516{
517    int indx, u=width, v=2*width; 
518    int col, row;
519    double Co, Ho, ratio;
520
521    for (row=6; row < height-6; row++)
522    {
523       for (col=6; col < width-6; col++)
524       {
525          indx = row*width+col;
526
527          if ( image2[indx][1]*image2[indx][2] != 0 ) 
528          {
529             Co = (image2[indx+v][1] + image2[indx-v][1] + image2[indx-2][1] + image2[indx+2][1] - 
530                   MAX(image2[indx-2][1], MAX(image2[indx+2][1], MAX(image2[indx-v][1], image2[indx+v][1]))) -
531                   MIN(image2[indx-2][1], MIN(image2[indx+2][1], MIN(image2[indx-v][1], image2[indx+v][1]))))/2.0;
532             Ho = (image2[indx+v][2] + image2[indx-v][2] + image2[indx-2][2] + image2[indx+2][2] - 
533                   MAX(image2[indx-2][2], MAX(image2[indx+2][2], MAX(image2[indx-v][2], image2[indx+v][2]))) -
534                   MIN(image2[indx-2][2], MIN(image2[indx+2][2], MIN(image2[indx-v][2], image2[indx+v][2]))))/2.0;
535             ratio = sqrt ((Co*Co+Ho*Ho) / (image2[indx][1]*image2[indx][1] + image2[indx][2]*image2[indx][2]));
536            
537             if (ratio < 0.85)
538             {
539                 image2[indx][0] = -(image2[indx][1] + image2[indx][2] - Co - Ho) + image2[indx][0];
540                 image2[indx][1] = Co;
541                 image2[indx][2] = Ho;
542	     }
543          }
544       }
545    }
546}
547
548
549// Cubic Spline Interpolation by Li and Randhawa, modified by Jacek Gozdz and Luis Sanz Rodríguez
550void CLASS fbdd_green()
551{
552	int row, col, c, u=width, v=2*u, w=3*u, x=4*u, y=5*u, indx, min, max, current;
553	float f[4], g[4];
554
555	for (row=5; row < height-5; row++)
556		for (col=5+(FC(row,1)&1),indx=row*width+col,c=FC(row,col); col < u-5; col+=2,indx+=2) {
557			
558		
559f[0]=1.0/(1.0+abs(image[indx-u][1]-image[indx-w][1])+abs(image[indx-w][1]-image[indx+y][1]));
560f[1]=1.0/(1.0+abs(image[indx+1][1]-image[indx+3][1])+abs(image[indx+3][1]-image[indx-5][1]));
561f[2]=1.0/(1.0+abs(image[indx-1][1]-image[indx-3][1])+abs(image[indx-3][1]-image[indx+5][1]));
562f[3]=1.0/(1.0+abs(image[indx+u][1]-image[indx+w][1])+abs(image[indx+w][1]-image[indx-y][1]));
563
564g[0]=CLIP((23*image[indx-u][1]+23*image[indx-w][1]+2*image[indx-y][1]+8*(image[indx-v][c]-image[indx-x][c])+40*(image[indx][c]-image[indx-v][c]))/48.0);
565g[1]=CLIP((23*image[indx+1][1]+23*image[indx+3][1]+2*image[indx+5][1]+8*(image[indx+2][c]-image[indx+4][c])+40*(image[indx][c]-image[indx+2][c]))/48.0);
566g[2]=CLIP((23*image[indx-1][1]+23*image[indx-3][1]+2*image[indx-5][1]+8*(image[indx-2][c]-image[indx-4][c])+40*(image[indx][c]-image[indx-2][c]))/48.0);
567g[3]=CLIP((23*image[indx+u][1]+23*image[indx+w][1]+2*image[indx+y][1]+8*(image[indx+v][c]-image[indx+x][c])+40*(image[indx][c]-image[indx+v][c]))/48.0);
568
569	image[indx][1]=CLIP((f[0]*g[0]+f[1]*g[1]+f[2]*g[2]+f[3]*g[3])/(f[0]+f[1]+f[2]+f[3]));
570	
571	min = MIN(image[indx+1+u][1], MIN(image[indx+1-u][1], MIN(image[indx-1+u][1], MIN(image[indx-1-u][1], MIN(image[indx-1][1], MIN(image[indx+1][1], MIN(image[indx-u][1], image[indx+u][1])))))));
572
573	max = MAX(image[indx+1+u][1], MAX(image[indx+1-u][1], MAX(image[indx-1+u][1], MAX(image[indx-1-u][1], MAX(image[indx-1][1], MAX(image[indx+1][1], MAX(image[indx-u][1], image[indx+u][1])))))));
574
575	image[indx][1] = ULIM(image[indx][1], max, min);				
576		}
577}
578
579
580
581
582// FBDD (Fake Before Demosaicing Denoising)
583void CLASS fbdd(int noiserd)
584{ 
585	double (*image2)[3];
586        // safety net: disable for 4-color bayer or full-color images
587        if(colors!=3 || !filters)
588            return;
589	image2 = (double (*)[3]) calloc(width*height, sizeof *image2);
590
591	border_interpolate(4);
592
593if (noiserd>1)
594{
595#ifdef DCRAW_VERBOSE
596	if (verbose) fprintf (stderr,_("FBDD full noise reduction...\n"));
597#endif
598	fbdd_green();
599	//dcb_color_full(image2);
600	dcb_color_full();
601	fbdd_correction();	
602	
603
604	dcb_color();
605	rgb_to_lch(image2);
606	fbdd_correction2(image2);
607	fbdd_correction2(image2);
608	lch_to_rgb(image2); 
609
610}
611else
612{
613#ifdef DCRAW_VERBOSE
614	if (verbose) fprintf (stderr,_("FBDD noise reduction...\n"));
615#endif
616
617	fbdd_green();
618	//dcb_color_full(image2);
619	dcb_color_full();
620	fbdd_correction();	
621}
622
623free(image2);
624
625}
626
627
628
629
630// DCB demosaicing main routine
631void CLASS dcb(int iterations, int dcb_enhance)
632{
633
634
635	int i=1;
636	
637	float (*image2)[3];
638	image2 = (float (*)[3]) calloc(width*height, sizeof *image2);
639
640	float (*image3)[3];
641	image3 = (float (*)[3]) calloc(width*height, sizeof *image3);
642	
643#ifdef DCRAW_VERBOSE
644	if (verbose) fprintf (stderr,_("DCB demosaicing...\n"));
645#endif
646
647 		border_interpolate(6);
648
649		dcb_hor(image2);
650 		dcb_color2(image2);
651  		
652 		dcb_ver(image3);
653 		dcb_color3(image3);
654  		
655 		dcb_decide(image2, image3);
656
657		free(image3);
658
659		dcb_copy_to_buffer(image2);
660
661		while (i<=iterations)
662		{
663#ifdef DCRAW_VERBOSE
664			if (verbose) fprintf (stderr,_("DCB correction pass %d...\n"), i);
665#endif                        
666			dcb_nyquist();
667			dcb_nyquist();
668			dcb_nyquist();
669			dcb_map();
670			dcb_correction();	
671			i++;
672		}
673		
674		dcb_color();
675		dcb_pp();	
676		
677#ifdef DCRAW_VERBOSE
678	if (verbose) fprintf (stderr,_("finishing DCB...\n"));
679#endif		
680		dcb_map();
681		dcb_correction2();
682
683		dcb_map();
684		dcb_correction();
685
686		dcb_map();
687		dcb_correction();
688
689		dcb_map();
690		dcb_correction();	
691
692		dcb_map();
693		dcb_restore_from_buffer(image2);
694		dcb_color();
695
696	if (dcb_enhance)
697	{	
698#ifdef DCRAW_VERBOSE
699		if (verbose) fprintf (stderr,_("optional DCB refinement...\n"));
700#endif
701		dcb_refinement();
702 		//dcb_color_full(image2);		
703 		dcb_color_full();		
704	}
705
706
707		free(image2);
708
709}
710