/src/FreeImage/Source/LibRawLite/internal/dcb_demosaicing.c
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