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 20 //#include "rtengine.h" 21 #include "rawimagesource.hh" 22 #include "rt_math.h" 23 //#include "../rtgui/multilangmgr.h" 24 //#include "procparams.h" 25 #include "sleef.c" 26 #include "opthelper.h" 27 28 29 #undef CLIP 30 #define CLIP(x) x 31 32 33 //namespace rtengine { 34 namespace rtengine { 35 igv_demosaic_RT(int winx,int winy,int winw,int winh,int tilex,int tiley,int tilew,int tileh)36 void RawImageSource::igv_demosaic_RT(int winx, int winy, int winw, int winh, 37 int tilex, int tiley, int tilew, int tileh) 38 { 39 //if( tilex>0 || tiley > 0 ) 40 //return; 41 static const float eps=1e-5f, epssq=1e-5f;//mod epssq -10f =>-5f Jacques 3/2013 to prevent artifact (divide by zero) 42 static const int h1=1, h2=2, h3=3, h4=4, h5=5, h6=6; 43 // Image dimensions 44 const int iwidth=winw, iheight=winh; 45 // Input tile dimensions 46 const int width=tilew, height=tileh; 47 const int v1=1*width, v2=2*width, v3=3*width, v4=4*width, v5=5*width, v6=6*width; 48 float* rgb[3]; 49 float* chr[2]; 50 float (*rgbarray), *vdif, *hdif, (*chrarray); 51 52 rgbarray = (float (*)) calloc(width*height*3, sizeof( float)); 53 rgb[0] = rgbarray; 54 rgb[1] = rgbarray + (width*height); 55 rgb[2] = rgbarray + 2*(width*height); 56 57 chrarray = (float (*)) calloc(width*height*2, sizeof( float)); 58 chr[0] = chrarray; 59 chr[1] = chrarray + (width*height); 60 61 vdif = (float (*)) calloc(width*height/2, sizeof *vdif); 62 hdif = (float (*)) calloc(width*height/2, sizeof *hdif); 63 64 //border_interpolate2(winw,winh,7); 65 66 //if (plistener) { 67 // plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::methodstring[RAWParams::BayerSensor::igv])); 68 // plistener->setProgress (0.0); 69 //} 70 /* 71 #ifdef _OPENMP 72 #pragma omp parallel default(none) shared(rgb,vdif,hdif,chr) 73 #endif 74 */ 75 { 76 77 float ng, eg, wg, sg, nv, ev, wv, sv, nwg, neg, swg, seg, nwv, nev, swv, sev; 78 //int row, col, row2, col2; 79 // refcol == 1 if tilex is even 80 // refcol == 0 if tilex is odd 81 int refcol = 1-(tilex%2); 82 //std::cout<<"refcol: "<<refcol<<std::endl; 83 /* 84 #ifdef _OPENMP 85 #pragma omp for 86 #endif 87 */ 88 for (int row=0, row2=tiley; row<height-0; row++, row2++) { 89 for (int col=0, col2=tilex, indx=row*width+col; col<width-0; col++, col2++, indx++) { 90 int c=FC(row,col); 91 float val = rawData[row2][col2]; 92 rgb[c][indx]=CLIP( val ); 93 if( false && row2<16 && col2<16) 94 std::cout<<"row,col="<<row2<<","<<col2<<" rgb["<<c<<"]["<<indx<<"]="<<rgb[c][indx]<<std::endl; 95 } 96 } 97 // border_interpolate2(7, rgb); 98 99 #ifdef _OPENMP 100 #pragma omp single 101 #endif 102 { 103 //if (plistener) plistener->setProgress (0.13); 104 } 105 106 #ifdef _OPENMP 107 #pragma omp for 108 #endif 109 for (int row=5, row2=tiley+row; row<height-5; row++, row2++) 110 for (int col=5+(FC(row,refcol)&1), col2=tilex+col, indx=row*width+col, c=FC(row,col); col<width-5; col+=2, col2+=2, indx+=2) { 111 //N,E,W,S Gradients 112 ng=(eps+(fabsf(rgb[1][indx-v1]-rgb[1][indx-v3])+fabsf(rgb[c][indx]-rgb[c][indx-v2]))/65535.f);; 113 eg=(eps+(fabsf(rgb[1][indx+h1]-rgb[1][indx+h3])+fabsf(rgb[c][indx]-rgb[c][indx+h2]))/65535.f); 114 wg=(eps+(fabsf(rgb[1][indx-h1]-rgb[1][indx-h3])+fabsf(rgb[c][indx]-rgb[c][indx-h2]))/65535.f); 115 sg=(eps+(fabsf(rgb[1][indx+v1]-rgb[1][indx+v3])+fabsf(rgb[c][indx]-rgb[c][indx+v2]))/65535.f); 116 //N,E,W,S High Order Interpolation (Li & Randhawa) 117 //N,E,W,S Hamilton Adams Interpolation 118 // (48.f * 65535.f) = 3145680.f 119 nv=rtengine::LIM(((23.0f*rgb[1][indx-v1]+23.0f*rgb[1][indx-v3]+rgb[1][indx-v5]+rgb[1][indx+v1]+40.0f*rgb[c][indx]-32.0f*rgb[c][indx-v2]-8.0f*rgb[c][indx-v4]))/3145680.f, 0.0f, 1.0f); 120 ev=rtengine::LIM(((23.0f*rgb[1][indx+h1]+23.0f*rgb[1][indx+h3]+rgb[1][indx+h5]+rgb[1][indx-h1]+40.0f*rgb[c][indx]-32.0f*rgb[c][indx+h2]-8.0f*rgb[c][indx+h4]))/3145680.f, 0.0f, 1.0f); 121 wv=rtengine::LIM(((23.0f*rgb[1][indx-h1]+23.0f*rgb[1][indx-h3]+rgb[1][indx-h5]+rgb[1][indx+h1]+40.0f*rgb[c][indx]-32.0f*rgb[c][indx-h2]-8.0f*rgb[c][indx-h4]))/3145680.f, 0.0f, 1.0f); 122 sv=rtengine::LIM(((23.0f*rgb[1][indx+v1]+23.0f*rgb[1][indx+v3]+rgb[1][indx+v5]+rgb[1][indx-v1]+40.0f*rgb[c][indx]-32.0f*rgb[c][indx+v2]-8.0f*rgb[c][indx+v4]))/3145680.f, 0.0f, 1.0f); 123 //Horizontal and vertical color differences 124 vdif[indx>>1]=(sg*nv+ng*sv)/(ng+sg)-(rgb[c][indx])/65535.f; 125 hdif[indx>>1]=(wg*ev+eg*wv)/(eg+wg)-(rgb[c][indx])/65535.f; 126 } 127 128 #ifdef _OPENMP 129 #pragma omp single 130 #endif 131 { 132 //if (plistener) plistener->setProgress (0.26); 133 } 134 135 #ifdef _OPENMP 136 #pragma omp for 137 #endif 138 for (int row=7, row2=tiley+row; row<height-7; row++, row2++) 139 for (int col=7+(FC(row,refcol)&1), col2=tilex+col, indx=row*width+col, c=FC(row,col), d=c/2; col<width-7; col+=2, col2+=2, indx+=2) { 140 //H&V integrated gaussian vector over variance on color differences 141 //Mod Jacques 3/2013 142 ng=rtengine::LIM(epssq+78.0f*SQR(vdif[indx>>1])+69.0f*(SQR(vdif[(indx-v2)>>1])+SQR(vdif[(indx+v2)>>1]))+51.0f*(SQR(vdif[(indx-v4)>>1])+SQR(vdif[(indx+v4)>>1]))+21.0f*(SQR(vdif[(indx-v6)>>1])+SQR(vdif[(indx+v6)>>1]))-6.0f*SQR(vdif[(indx-v2)>>1]+vdif[indx>>1]+vdif[(indx+v2)>>1]) 143 -10.0f*(SQR(vdif[(indx-v4)>>1]+vdif[(indx-v2)>>1]+vdif[indx>>1])+SQR(vdif[indx>>1]+vdif[(indx+v2)>>1]+vdif[(indx+v4)>>1]))-7.0f*(SQR(vdif[(indx-v6)>>1]+vdif[(indx-v4)>>1]+vdif[(indx-v2)>>1])+SQR(vdif[(indx+v2)>>1]+vdif[(indx+v4)>>1]+vdif[(indx+v6)>>1])),0.f,1.f); 144 eg=rtengine::LIM(epssq+78.0f*SQR(hdif[indx>>1])+69.0f*(SQR(hdif[(indx-h2)>>1])+SQR(hdif[(indx+h2)>>1]))+51.0f*(SQR(hdif[(indx-h4)>>1])+SQR(hdif[(indx+h4)>>1]))+21.0f*(SQR(hdif[(indx-h6)>>1])+SQR(hdif[(indx+h6)>>1]))-6.0f*SQR(hdif[(indx-h2)>>1]+hdif[indx>>1]+hdif[(indx+h2)>>1]) 145 -10.0f*(SQR(hdif[(indx-h4)>>1]+hdif[(indx-h2)>>1]+hdif[indx>>1])+SQR(hdif[indx>>1]+hdif[(indx+h2)>>1]+hdif[(indx+h4)>>1]))-7.0f*(SQR(hdif[(indx-h6)>>1]+hdif[(indx-h4)>>1]+hdif[(indx-h2)>>1])+SQR(hdif[(indx+h2)>>1]+hdif[(indx+h4)>>1]+hdif[(indx+h6)>>1])),0.f,1.f); 146 //Limit chrominance using H/V neighbourhood 147 nv=rtengine::ULIM(0.725f*vdif[indx>>1]+0.1375f*vdif[(indx-v2)>>1]+0.1375f*vdif[(indx+v2)>>1],vdif[(indx-v2)>>1],vdif[(indx+v2)>>1]); 148 ev=rtengine::ULIM(0.725f*hdif[indx>>1]+0.1375f*hdif[(indx-h2)>>1]+0.1375f*hdif[(indx+h2)>>1],hdif[(indx-h2)>>1],hdif[(indx+h2)>>1]); 149 //Chrominance estimation 150 chr[d][indx]=(eg*nv+ng*ev)/(ng+eg); 151 //Green channel population 152 rgb[1][indx]=rgb[c][indx]+65535.f*chr[d][indx]; 153 } 154 155 #ifdef _OPENMP 156 #pragma omp single 157 #endif 158 { 159 //if (plistener) plistener->setProgress (0.39); 160 } 161 162 // free(vdif); free(hdif); 163 #ifdef _OPENMP 164 #pragma omp for 165 #endif 166 for (int row=7, row2=tiley+row; row<height-7; row+=2, row2+=2) 167 for (int col=7+(FC(row,refcol)&1), col2=tilex+col, indx=row*width+col, c=1-FC(row,col)/2; col<width-7; col+=2, col2+=2, indx+=2) { 168 //NW,NE,SW,SE Gradients 169 nwg=1.0f/(eps+fabsf(chr[c][indx-v1-h1]-chr[c][indx-v3-h3])+fabsf(chr[c][indx+v1+h1]-chr[c][indx-v3-h3])); 170 neg=1.0f/(eps+fabsf(chr[c][indx-v1+h1]-chr[c][indx-v3+h3])+fabsf(chr[c][indx+v1-h1]-chr[c][indx-v3+h3])); 171 swg=1.0f/(eps+fabsf(chr[c][indx+v1-h1]-chr[c][indx+v3+h3])+fabsf(chr[c][indx-v1+h1]-chr[c][indx+v3-h3])); 172 seg=1.0f/(eps+fabsf(chr[c][indx+v1+h1]-chr[c][indx+v3-h3])+fabsf(chr[c][indx-v1-h1]-chr[c][indx+v3+h3])); 173 //Limit NW,NE,SW,SE Color differences 174 nwv=rtengine::ULIM(chr[c][indx-v1-h1],chr[c][indx-v3-h1],chr[c][indx-v1-h3]); 175 nev=rtengine::ULIM(chr[c][indx-v1+h1],chr[c][indx-v3+h1],chr[c][indx-v1+h3]); 176 swv=rtengine::ULIM(chr[c][indx+v1-h1],chr[c][indx+v3-h1],chr[c][indx+v1-h3]); 177 sev=rtengine::ULIM(chr[c][indx+v1+h1],chr[c][indx+v3+h1],chr[c][indx+v1+h3]); 178 //Interpolate chrominance: R@B and B@R 179 chr[c][indx]=(nwg*nwv+neg*nev+swg*swv+seg*sev)/(nwg+neg+swg+seg); 180 } 181 #ifdef _OPENMP 182 #pragma omp single 183 #endif 184 { 185 //if (plistener) plistener->setProgress (0.52); 186 } 187 #ifdef _OPENMP 188 #pragma omp for 189 #endif 190 for (int row=8, row2=tiley+row; row<height-7; row+=2, row2+=2) 191 for (int col=7+(FC(row,refcol)&1), col2=tilex+col, indx=row*width+col, c=1-FC(row,col)/2; col<width-7; col+=2, col2+=2, indx+=2) { 192 //NW,NE,SW,SE Gradients 193 nwg=1.0f/(eps+fabsf(chr[c][indx-v1-h1]-chr[c][indx-v3-h3])+fabsf(chr[c][indx+v1+h1]-chr[c][indx-v3-h3])); 194 neg=1.0f/(eps+fabsf(chr[c][indx-v1+h1]-chr[c][indx-v3+h3])+fabsf(chr[c][indx+v1-h1]-chr[c][indx-v3+h3])); 195 swg=1.0f/(eps+fabsf(chr[c][indx+v1-h1]-chr[c][indx+v3+h3])+fabsf(chr[c][indx-v1+h1]-chr[c][indx+v3-h3])); 196 seg=1.0f/(eps+fabsf(chr[c][indx+v1+h1]-chr[c][indx+v3-h3])+fabsf(chr[c][indx-v1-h1]-chr[c][indx+v3+h3])); 197 //Limit NW,NE,SW,SE Color differences 198 nwv=rtengine::ULIM(chr[c][indx-v1-h1],chr[c][indx-v3-h1],chr[c][indx-v1-h3]); 199 nev=rtengine::ULIM(chr[c][indx-v1+h1],chr[c][indx-v3+h1],chr[c][indx-v1+h3]); 200 swv=rtengine::ULIM(chr[c][indx+v1-h1],chr[c][indx+v3-h1],chr[c][indx+v1-h3]); 201 sev=rtengine::ULIM(chr[c][indx+v1+h1],chr[c][indx+v3+h1],chr[c][indx+v1+h3]); 202 //Interpolate chrominance: R@B and B@R 203 chr[c][indx]=(nwg*nwv+neg*nev+swg*swv+seg*sev)/(nwg+neg+swg+seg); 204 } 205 #ifdef _OPENMP 206 #pragma omp single 207 #endif 208 { 209 //if (plistener) plistener->setProgress (0.65); 210 } 211 #ifdef _OPENMP 212 #pragma omp for 213 #endif 214 for (int row=7, row2=tiley+row; row<height-7; row++, row2++) 215 for (int col=7+(FC(row,0)&1), col2=tilex+col, indx=row*width+col; col<width-7; col+=2, col2+=2, indx+=2) { 216 //N,E,W,S Gradients 217 ng=1.0f/(eps+fabsf(chr[0][indx-v1]-chr[0][indx-v3])+fabsf(chr[0][indx+v1]-chr[0][indx-v3])); 218 eg=1.0f/(eps+fabsf(chr[0][indx+h1]-chr[0][indx+h3])+fabsf(chr[0][indx-h1]-chr[0][indx+h3])); 219 wg=1.0f/(eps+fabsf(chr[0][indx-h1]-chr[0][indx-h3])+fabsf(chr[0][indx+h1]-chr[0][indx-h3])); 220 sg=1.0f/(eps+fabsf(chr[0][indx+v1]-chr[0][indx+v3])+fabsf(chr[0][indx-v1]-chr[0][indx+v3])); 221 //Interpolate chrominance: R@G and B@G 222 chr[0][indx]=((ng*chr[0][indx-v1]+eg*chr[0][indx+h1]+wg*chr[0][indx-h1]+sg*chr[0][indx+v1])/(ng+eg+wg+sg)); 223 } 224 #ifdef _OPENMP 225 #pragma omp single 226 #endif 227 { 228 //if (plistener) plistener->setProgress (0.78); 229 } 230 #ifdef _OPENMP 231 #pragma omp for 232 #endif 233 for (int row=7, row2=tiley+row; row<height-7; row++, row2++) 234 for (int col=7+(FC(row,0)&1), col2=tilex+col, indx=row*width+col; col<width-7; col+=2, col2+=2, indx+=2) { 235 236 //N,E,W,S Gradients 237 ng=1.0f/(eps+fabsf(chr[1][indx-v1]-chr[1][indx-v3])+fabsf(chr[1][indx+v1]-chr[1][indx-v3])); 238 eg=1.0f/(eps+fabsf(chr[1][indx+h1]-chr[1][indx+h3])+fabsf(chr[1][indx-h1]-chr[1][indx+h3])); 239 wg=1.0f/(eps+fabsf(chr[1][indx-h1]-chr[1][indx-h3])+fabsf(chr[1][indx+h1]-chr[1][indx-h3])); 240 sg=1.0f/(eps+fabsf(chr[1][indx+v1]-chr[1][indx+v3])+fabsf(chr[1][indx-v1]-chr[1][indx+v3])); 241 //Interpolate chrominance: R@G and B@G 242 chr[1][indx]=((ng*chr[1][indx-v1]+eg*chr[1][indx+h1]+wg*chr[1][indx-h1]+sg*chr[1][indx+v1])/(ng+eg+wg+sg)); 243 } 244 #ifdef _OPENMP 245 #pragma omp single 246 #endif 247 { 248 //if (plistener) plistener->setProgress (0.91); 249 250 //Interpolate borders 251 // border_interpolate2(7, rgb); 252 } 253 /* 254 #ifdef _OPENMP 255 #pragma omp for 256 #endif 257 for (int row=0; row < height; row++) //borders 258 for (int col=0; col < width; col++) { 259 if (col==7 && row >= 7 && row < height-7) 260 col = width-7; 261 int indxc=row*width+col; 262 red [row][col] = rgb[indxc][0]; 263 green[row][col] = rgb[indxc][1]; 264 blue [row][col] = rgb[indxc][2]; 265 } 266 */ 267 268 #ifdef _OPENMP 269 #pragma omp for 270 #endif 271 for(int row=7, row2=tiley+row; row<height-7; row++, row2++) 272 for(int col=7, col2=tilex+col, indx=row*width+col; col<width-7; col++, col2++, indx++) { 273 red [row2][col2] = CLIP(rgb[1][indx]-65535.f*chr[0][indx]); 274 green[row2][col2] = CLIP(rgb[1][indx]); 275 blue [row2][col2] = CLIP(rgb[1][indx]-65535.f*chr[1][indx]); 276 //if( row<16 && col<16) 277 //std::cout<<"row,col="<<row<<","<<col<<" red: "<<red [row2][col2]<<std::endl; 278 } 279 }// End of parallelization 280 281 //if (plistener) plistener->setProgress (1.0); 282 283 free(chrarray); free(rgbarray); 284 free(vdif); free(hdif); 285 286 } 287 } 288