From fb906af9494ec916849d3d053541f6c39ec3a761 Mon Sep 17 00:00:00 2001 From: Emil Martinec Date: Thu, 16 Sep 2010 14:07:18 -0500 Subject: [PATCH] Oops! Forgot to add the fast demosaic code itself. Here's the patch... --- rtengine/fast_demo.cc | 272 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 272 insertions(+) create mode 100644 rtengine/fast_demo.cc diff --git a/rtengine/fast_demo.cc b/rtengine/fast_demo.cc new file mode 100644 index 000000000..96e73c844 --- /dev/null +++ b/rtengine/fast_demo.cc @@ -0,0 +1,272 @@ +//////////////////////////////////////////////////////////////// +// +// Fast demosaicing algorythm +// +// copyright (c) 2008-2010 Emil Martinec +// +// +// code dated: August 26, 2010 +// +// fast_demo.cc is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . +// +//////////////////////////////////////////////////////////////// + + + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +void RawImageSource::fast_demo() { + + if (plistener) { + plistener->setProgressStr ("Fast demosaicing..."); + plistener->setProgress (0.0); + } + float progress = 0.0; + + //allocate output arrays + red = new unsigned short*[H]; + for (int i=0; i -1 && i1 < H && j1 > -1) { + c = FC(i1,j1); + sum[c] += ri->data[i1][j1]; + sum[c+3]++; + } + } + c=FC(i,j); + if (c==1) { + red[i][j]=sum[0]/sum[3]; + green[i][j]=ri->data[i][j]; + blue[i][j]=sum[2]/sum[5]; + } else { + green[i][j]=sum[1]/sum[4]; + if (c==0) { + red[i][j]=ri->data[i][j]; + blue[i][j]=sum[2]/sum[5]; + } else { + red[i][j]=sum[0]/sum[3]; + blue[i][j]=ri->data[i][j]; + } + } + }//j + + for (j=W-bord; j -1 && i1 < H && j1 < W) { + c = FC(i1,j1); + sum[c] += ri->data[i1][j1]; + sum[c+3]++; + } + } + c=FC(i,j); + if (c==1) { + red[i][j]=sum[0]/sum[3]; + green[i][j]=ri->data[i][j]; + blue[i][j]=sum[2]/sum[5]; + } else { + green[i][j]=sum[1]/sum[4]; + if (c==0) { + red[i][j]=ri->data[i][j]; + blue[i][j]=sum[2]/sum[5]; + } else { + red[i][j]=sum[0]/sum[3]; + blue[i][j]=ri->data[i][j]; + } + } + }//j + }//i + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + for (j=bord; j -1 && j1 < W && i1 > -1) { + c = FC(i1,j1); + sum[c] += ri->data[i1][j1]; + sum[c+3]++; + } + } + c=FC(i,j); + if (c==1) { + red[i][j]=sum[0]/sum[3]; + green[i][j]=ri->data[i][j]; + blue[i][j]=sum[2]/sum[5]; + } else { + green[i][j]=sum[1]/sum[4]; + if (c==0) { + red[i][j]=ri->data[i][j]; + blue[i][j]=sum[2]/sum[5]; + } else { + red[i][j]=sum[0]/sum[3]; + blue[i][j]=ri->data[i][j]; + } + } + }//i + + for (i=H-bord; i -1 && j1 < W && i1 < H) { + c = FC(i1,j1); + sum[c] += ri->data[i1][j1]; + sum[c+3]++; + } + } + c=FC(i,j); + if (c==1) { + red[i][j]=sum[0]/sum[3]; + green[i][j]=ri->data[i][j]; + blue[i][j]=sum[2]/sum[5]; + } else { + green[i][j]=sum[1]/sum[4]; + if (c==0) { + red[i][j]=ri->data[i][j]; + blue[i][j]=sum[2]/sum[5]; + } else { + red[i][j]=sum[0]/sum[3]; + blue[i][j]=ri->data[i][j]; + } + } + }//i + }//j + + if(plistener) plistener->setProgress(0.05); + + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +#pragma omp parallel private(i,j,c) + { + int rb; + + float wtu, wtd, wtl, wtr; + + float * dirwt = new float [0x20000]; + +#pragma omp for schedule(dynamic) nowait + + //set up directional weight function + for (int i=0; i<0x10000; i++) + dirwt[i] = 1.0/SQR(1.0+i); + + +#pragma omp for schedule(dynamic) nowait + + // interpolate G using gradient weights + for (i=bord; i< H-bord; i++) { + for (j=bord; j < W-bord; j++) { + + if (FC(i,j)==1) { + green[i][j] = ri->data[i][j]; + //red[i][j] = green[i][j]; + //blue[i][j] = green[i][j]; + + } else { + //compute directional weights using image gradients + wtu=dirwt[(abs(ri->data[i+1][j]-ri->data[i-1][j])+abs(ri->data[i][j]-ri->data[i-2][j])+abs(ri->data[i-1][j]-ri->data[i-3][j])) >>8]; + wtd=dirwt[(abs(ri->data[i-1][j]-ri->data[i+1][j])+abs(ri->data[i][j]-ri->data[i+2][j])+abs(ri->data[i+1][j]-ri->data[i+3][j])) >>8]; + wtl=dirwt[(abs(ri->data[i][j+1]-ri->data[i][j-1])+abs(ri->data[i][j]-ri->data[i][j-2])+abs(ri->data[i][j-1]-ri->data[i][j-3])) >>8]; + wtr=dirwt[(abs(ri->data[i][j-1]-ri->data[i][j+1])+abs(ri->data[i][j]-ri->data[i][j+2])+abs(ri->data[i][j+1]-ri->data[i][j+3])) >>8]; + + //wtu=1/SQR(1.0+fabs((int)ri->data[i+1][j]-ri->data[i-1][j])+fabs((int)ri->data[i][j]-ri->data[i-2][j])+fabs((int)ri->data[i-1][j]-ri->data[i-3][j])); + //wtd=1/SQR(1.0+fabs((int)ri->data[i-1][j]-ri->data[i+1][j])+fabs((int)ri->data[i][j]-ri->data[i+2][j])+fabs((int)ri->data[i+1][j]-ri->data[i+3][j])); + //wtl=1/SQR(1.0+fabs((int)ri->data[i][j+1]-ri->data[i][j-1])+fabs((int)ri->data[i][j]-ri->data[i][j-2])+fabs((int)ri->data[i][j-1]-ri->data[i][j-3])); + //wtr=1/SQR(1.0+fabs((int)ri->data[i][j-1]-ri->data[i][j+1])+fabs((int)ri->data[i][j]-ri->data[i][j+2])+fabs((int)ri->data[i][j+1]-ri->data[i][j+3])); + + //store in rgb array the interpolated G value at R/B grid points using directional weighted average + green[i][j]=(int)((wtu*ri->data[i-1][j]+wtd*ri->data[i+1][j]+wtl*ri->data[i][j-1]+wtr*ri->data[i][j+1])/(wtu+wtd+wtl+wtr)); + //red[i][j] = green[i][j]; + //blue[i][j] = green[i][j]; + + } + } + progress+=(double)0.33/(H); + //if(plistener) plistener->setProgress(progress); + } + if(plistener) plistener->setProgress(0.4); + + +#pragma omp for schedule(dynamic) nowait + + for (i=bord; i< H-bord; i++) { + for (j=bord+(FC(i,2)&1); j < W-bord; j+=2) { + + c=FC(i,j); + //interpolate B/R colors at R/B sites + rb = CLIP((int)(green[i][j] - 0.25*((green[i-1][j-1]-ri->data[i-1][j-1])+(green[i-1][j+1]-ri->data[i-1][j+1])+ \ + (green[i+1][j+1]-ri->data[i+1][j+1])+(green[i+1][j-1]-ri->data[i+1][j-1])))); + if (c==0) {//R site + red[i][j] = ri->data[i][j]; + blue[i][j] = rb; + } else {//B site + red[i][j] = rb; + blue[i][j] = ri->data[i][j]; + } + } + progress+=(double)0.33/(H); + //if(plistener) plistener->setProgress(progress); + } + if(plistener) plistener->setProgress(0.7); + + +#pragma omp for schedule(dynamic) nowait + + // interpolate R/B using color differences + for (i=bord; i< H-bord; i++) { + for (j=bord+1-(FC(i,2)&1); j < W-bord; j+=2) { + + //interpolate R and B colors at G sites + red[i][j] = CLIP((int)(green[i][j] - 0.25*((green[i-1][j]-red[i-1][j])+(green[i+1][j]-red[i+1][j])+ \ + (green[i][j-1]-red[i][j-1])+(green[i][j+1]-red[i][j+1])))); + blue[i][j] = CLIP((int)(green[i][j] - 0.25*((green[i-1][j]-blue[i-1][j])+(green[i+1][j]-blue[i+1][j])+ \ + (green[i][j-1]-blue[i][j-1])+(green[i][j+1]-blue[i][j+1])))); + } + progress+=(double)0.33/(H); + //if(plistener) plistener->setProgress(progress); + } + if(plistener) plistener->setProgress(0.99); + + } + + + +#undef bord + + + +}//namespace