From 08885a7d7713053612ed8593716c778dbf029044 Mon Sep 17 00:00:00 2001 From: Emil Martinec Date: Wed, 15 Sep 2010 21:39:55 -0500 Subject: [PATCH] Minor improvements to denoise tools. Fixed a minor annoyance with zoom panel. --- rtengine/dirpyrLab_denoise.cc | 23 +- rtengine/rawimagesource.cc | 6452 ++++++++++++++++----------------- rtengine/simpleprocess.cc | 4 +- rtgui/cropwindow.cc | 10 +- rtgui/zoompanel.cc | 6 +- 5 files changed, 3246 insertions(+), 3249 deletions(-) diff --git a/rtengine/dirpyrLab_denoise.cc b/rtengine/dirpyrLab_denoise.cc index ee949c68b..d63d35d24 100644 --- a/rtengine/dirpyrLab_denoise.cc +++ b/rtengine/dirpyrLab_denoise.cc @@ -39,6 +39,7 @@ #define DIRWT_L(i1,j1,i,j) (/*domker[(i1-i)/scale+halfwin][(j1-j)/scale+halfwin] */ rangefn_L[(int)(data_fine->L[i1][j1]-data_fine->L[i][j]+0x10000)] ) #define DIRWT_AB(i1,j1,i,j) ( /*domker[(i1-i)/scale+halfwin][(j1-j)/scale+halfwin]*/ rangefn_ab[(int)(data_fine->a[i1][j1]-data_fine->a[i][j]+0x10000)] * \ +rangefn_ab[(int)(data_fine->L[i1][j1]-data_fine->L[i][j]+0x10000)] * \ rangefn_ab[(int)(data_fine->b[i1][j1]-data_fine->b[i][j]+0x10000)] ) #define NRWT_L(a) (nrwt_l[a] ) @@ -106,15 +107,13 @@ namespace rtengine { int * rangefn_L = new int [0x20000]; - int * irangefn_L = new int [0x20000]; float * nrwt_l = new float [0x10000]; int * rangefn_ab = new int [0x20000]; - int * irangefn_ab = new int [0x20000]; float * nrwt_ab = new float [0x20000]; - int intfactor = 16384; + int intfactor = 1024;//16384; //set up NR weight functions @@ -143,10 +142,6 @@ namespace rtengine { for (int i=0; i<0x20000; i++) rangefn_ab[i] = (int)(( exp(-(double)fabs(i-0x10000) * tonefactor / (1+3*noise_ab)) * noisevar_ab/((double)(i-0x10000)*(double)(i-0x10000) + noisevar_ab))*intfactor); - for (int i=0; i<0x20000; i++) - irangefn_L[i] = 1+(int)( exp(-(double)fabs(i-0x10000) / (1+16*noise_L) )*intfactor); - for (int i=0; i<0x20000; i++) - irangefn_ab[i] = 1+(int)( exp(-(double)fabs(i-0x10000)/ (1+64*noise_ab) )*intfactor); for (int i=0; i<0x20000; i++) nrwt_ab[i] = ((1+abs(i-0x10000)/(1+8*noise_ab)) * exp(-(double)fabs(i-0x10000)/ (1+8*noise_ab) ) ); @@ -251,8 +246,6 @@ namespace rtengine { delete [] rangefn_L; delete [] rangefn_ab; - delete [] irangefn_L; - delete [] irangefn_ab; delete [] nrwt_l; delete [] nrwt_ab; @@ -478,10 +471,11 @@ namespace rtengine { //Wiener filter //luma - if (level<3) { + if (level<2) { hipass[0] = data_fine->L[i][j]-smooth->L[i][j]; hpffluct[0]=SQR(hipass[0])+0.001; hipass[0] *= hpffluct[0]/(hpffluct[0]+noisevar_L); + data_fine->L[i][j] = CLIP(hipass[0]+smooth->L[i][j]); } //chroma @@ -497,12 +491,9 @@ namespace rtengine { }*/ hipass[1] *= nrfactor; hipass[2] *= nrfactor; - - //hipass[0] = hipass[1] = hipass[2] = 0.0;//for testing - - wtdsum[0]=data_fine->L[i][j] = CLIP(hipass[0]+smooth->L[i][j]); - wtdsum[1]=data_fine->a[i][j] = hipass[1]+smooth->a[i][j]; - wtdsum[2]=data_fine->b[i][j] = hipass[2]+smooth->b[i][j]; + + data_fine->a[i][j] = hipass[1]+smooth->a[i][j]; + data_fine->b[i][j] = hipass[2]+smooth->b[i][j]; } delete smooth; diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc index a96a834e5..e5d30cb53 100644 --- a/rtengine/rawimagesource.cc +++ b/rtengine/rawimagesource.cc @@ -1,3226 +1,3226 @@ -/* - * This file is part of RawTherapee. - * - * Copyright (c) 2004-2010 Gabor Horvath - * - * RawTherapee 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. - * - * RawTherapee 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 RawTherapee. If not, see . - */ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#ifdef _OPENMP -#include -#endif - -namespace rtengine { - -int loadRaw (const char* fname, struct RawImage* ri); - -extern const Settings* settings; - -#undef ABS -#undef MAX -#undef MIN -#undef DIST - -#define ABS(a) ((a)<0?-(a):(a)) -#define MAX(a,b) ((a)<(b)?(b):(a)) -#define MIN(a,b) ((a)>(b)?(b):(a)) -#define DIST(a,b) (ABS(a-b)) - -RawImageSource::RawImageSource () : ImageSource(), plistener(NULL), border(4), cache(NULL), green(NULL), red(NULL), blue(NULL) { - - hrmap[0] = NULL; - hrmap[1] = NULL; - hrmap[2] = NULL; - needhr = NULL; - hpmap = NULL; - oldmethod = "None"; - camProfile = NULL; - embProfile = NULL; -} - -RawImageSource::~RawImageSource () { - - delete idata; - if (ri) { - if (ri->allocation) - free (ri->allocation); - if (ri->data) - free (ri->data); - if (ri->profile_data) - free (ri->profile_data); - delete ri; - } - if (green) - freeArray(green, H); - if (red) - freeArray(red, H); - if (blue) - freeArray(blue, H); - - delete [] cache; - if (hrmap[0]!=NULL) { - int dh = H/HR_SCALE; - freeArray(hrmap[0], dh); - freeArray(hrmap[1], dh); - freeArray(hrmap[2], dh); - } - if (needhr) - freeArray(needhr, H); - if (hpmap) - freeArray(hpmap, H); - if (camProfile) - cmsCloseProfile (camProfile); - if (embProfile) - cmsCloseProfile (embProfile); -} - -void RawImageSource::transformRect (PreviewProps pp, int tran, int &ssx1, int &ssy1, int &width, int &height, int &fw) { - - pp.x += border; - pp.y += border; - - if (d1x) { - if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) { - pp.x /= 2; - pp.w = pp.w/2+1; - } - else { - pp.y /= 2; - pp.h = pp.h/2+1; - } - } - - int w = W, h = H; - if (fuji) { - w = ri->fuji_width * 2 + 1; - h = (H - ri->fuji_width)*2 + 1; - } - - int sw = w, sh = h; - if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) { - sw = h; - sh = w; - } - if( pp.w > sw-2*border) pp.w = sw-2*border; - if( pp.h > sh-2*border) pp.h = sh-2*border; - - int ppx = pp.x, ppy = pp.y; - if (tran & TR_HFLIP) - ppx = sw - pp.x - pp.w; - if (tran & TR_VFLIP) - ppy = sh - pp.y - pp.h; - - int sx1 = ppx; - int sy1 = ppy; - int sx2 = ppx + pp.w; - int sy2 = ppy + pp.h; - - if ((tran & TR_ROT) == TR_R180) { - sx1 = w - ppx - pp.w; - sy1 = h - ppy - pp.h; - sx2 = sx1 + pp.w; - sy2 = sy1 + pp.h; - } - else if ((tran & TR_ROT) == TR_R90) { - sx1 = ppy; - sy1 = h - ppx - pp.w; - sx2 = sx1 + pp.h; - sy2 = sy1 + pp.w; - } - else if ((tran & TR_ROT) == TR_R270) { - sx1 = w - ppy - pp.h; - sy1 = ppx; - sx2 = sx1 + pp.h; - sy2 = sy1 + pp.w; - } - - if (fuji) { - // atszamoljuk a koordinatakat fuji-ra: - ssx1 = (sx1+sy1) / 2; - ssy1 = (sy1 - sx2 ) / 2 + ri->fuji_width; - int ssx2 = (sx2+sy2) / 2 + 1; - int ssy2 = (sy2 - sx1) / 2 + ri->fuji_width; - fw = (sx2 - sx1) / 2 / pp.skip; - width = (ssx2 - ssx1) / pp.skip + ((ssx2 - ssx1) % pp.skip > 0); - height = (ssy2 - ssy1) / pp.skip + ((ssy2 - ssy1) % pp.skip > 0); - } - else { - ssx1 = sx1; - ssy1 = sy1; - width = (sx2 - sx1) / pp.skip + ((sx2 - sx1) % pp.skip > 0); - height = (sy2 - sy1) / pp.skip + ((sy2 - sy1) % pp.skip > 0); - } -} - -void RawImageSource::interpolate_image(Image16* image, HRecParams hrp, double rm, double gm, double bm, int skip, int tran, int fw, int imwidth, int imheight, int sx1, int sy1, int start, int end) -{ - - unsigned short* red = new unsigned short[imwidth]; - unsigned short* grn = new unsigned short[imwidth]; - unsigned short* blue = new unsigned short[imwidth]; - - for (int i=sy1 + start*skip ,ix=start; ixfilters && this->red && this->blue) { - for (int j=0,jx=sx1; jred[i][jx]); - grn[j] = CLIP(gm*green[i][jx]); - blue[j] = CLIP(bm*this->blue[i][jx]); - } - } else if(ri->filters) { - if (i==0) - interpolate_row_rb_mul_pp (red, blue, NULL, green[i], green[i+1], i, rm, gm, bm, sx1, imwidth, skip); - else if (i==H-1) - interpolate_row_rb_mul_pp (red, blue, green[i-1], green[i], NULL, i, rm, gm, bm, sx1, imwidth, skip); - else - interpolate_row_rb_mul_pp (red, blue, green[i-1], green[i], green[i+1], i, rm, gm, bm, sx1, imwidth, skip); - - for (int j=0,jx=sx1; jdata[i][jx*3+0]); - grn[j] = CLIP(gm*ri->data[i][jx*3+1]); - blue[j] = CLIP(bm*ri->data[i][jx*3+2]); - } - } - - if (hrp.enabled) - hlRecovery (hrp.method, red, grn, blue, i, sx1, imwidth, skip); - - transLine (red, grn, blue, ix, image, tran, imwidth, imheight, fw); - - } - delete [] red; - delete [] grn; - delete [] blue; -} - -void RawImageSource::getImage (ColorTemp ctemp, int tran, Image16* image, PreviewProps pp, HRecParams hrp, ColorManagementParams cmp) { - - isrcMutex.lock (); - - tran = defTransform (tran); - - // compute channel multipliers - double r, g, b, rm, gm, bm; - ctemp.getMultipliers (r, g, b); - rm = icoeff[0][0]*r + icoeff[0][1]*g + icoeff[0][2]*b; - gm = icoeff[1][0]*r + icoeff[1][1]*g + icoeff[1][2]*b; - bm = icoeff[2][0]*r + icoeff[2][1]*g + icoeff[2][2]*b; - rm = ri->camwb_red / rm; - gm = ri->camwb_green / gm; - bm = ri->camwb_blue / bm; - double mul_lum = 0.299*rm + 0.587*gm + 0.114*bm; - rm /= mul_lum; - gm /= mul_lum; - bm /= mul_lum; - - if (hrp.enabled) - defGain = log(ri->defgain) / log(2.0); - else { - defGain = 0.0; - rm *= ri->defgain; - gm *= ri->defgain; - bm *= ri->defgain; - } - - if (hrp.enabled==true && hrp.method=="Color" && hrmap[0]==NULL) - updateHLRecoveryMap_ColorPropagation (); - - // compute image area to render in order to provide the requested part of the image - int sx1, sy1, imwidth, imheight, fw; - transformRect (pp, tran, sx1, sy1, imwidth, imheight, fw); - - // check possible overflows - int maximwidth, maximheight; - if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) { - maximwidth = image->height; - maximheight = image->width; - } - else { - maximwidth = image->width; - maximheight = image->height; - } - if (d1x) - maximheight /= 2; - - // correct if overflow (very rare), but not fuji because it is corrected in transline - if (!fuji && imwidth>maximwidth) - imwidth = maximwidth; - if (!fuji && imheight>maximheight) - imheight = maximheight; - -#ifdef _OPENMP -#pragma omp parallel -{ - int tid = omp_get_thread_num(); - int nthreads = omp_get_num_threads(); - int blk = imheight/nthreads; - - - if (tidwidth%2==0) || ((tran & TR_ROT) == TR_R180 && image->height%2+image->width%2==1) || ((tran & TR_ROT) == TR_R270 && image->height%2==0); - // first row - for (int j=1+a; jwidth-1; j+=2) { - image->r[0][j] = (image->r[1][j] + image->r[0][j+1] + image->r[0][j-1]) / 3; - image->g[0][j] = (image->g[1][j] + image->g[0][j+1] + image->g[0][j-1]) / 3; - image->b[0][j] = (image->b[1][j] + image->b[0][j+1] + image->b[0][j-1]) / 3; - } - // other rows - for (int i=1; iheight-1; i++) { - for (int j=2-(a+i+1)%2; jwidth-1; j+=2) { - // edge-adaptive interpolation - double dh = (ABS(image->r[i][j+1] - image->r[i][j-1]) + ABS(image->g[i][j+1] - image->g[i][j-1]) + ABS(image->b[i][j+1] - image->b[i][j-1])) / 1.0; - double dv = (ABS(image->r[i+1][j] - image->r[i-1][j]) + ABS(image->g[i+1][j] - image->g[i-1][j]) + ABS(image->b[i+1][j] - image->b[i-1][j])) / 1.0; - double eh = 1.0 / (1.0 + dh); - double ev = 1.0 / (1.0 + dv); - image->r[i][j] = (eh * (image->r[i][j+1] + image->r[i][j-1]) + ev * (image->r[i+1][j] + image->r[i-1][j])) / (2.0 * (eh + ev)); - image->g[i][j] = (eh * (image->g[i][j+1] + image->g[i][j-1]) + ev * (image->g[i+1][j] + image->g[i-1][j])) / (2.0 * (eh + ev)); - image->b[i][j] = (eh * (image->b[i][j+1] + image->b[i][j-1]) + ev * (image->b[i+1][j] + image->b[i-1][j])) / (2.0 * (eh + ev)); - } - // first pixel - if (2-(a+i+1)%2==2) { - image->r[i][0] = (image->r[i+1][0] + image->r[i-1][0] + image->r[i][1]) / 3; - image->g[i][0] = (image->g[i+1][0] + image->g[i-1][0] + image->g[i][1]) / 3; - image->b[i][0] = (image->b[i+1][0] + image->b[i-1][0] + image->b[i][1]) / 3; - } - // last pixel - if (2-(a+i+image->width)%2==2) { - image->r[i][image->width-1] = (image->r[i+1][image->width-1] + image->r[i-1][image->width-1] + image->r[i][image->width-2]) / 3; - image->g[i][image->width-1] = (image->g[i+1][image->width-1] + image->g[i-1][image->width-1] + image->g[i][image->width-2]) / 3; - image->b[i][image->width-1] = (image->b[i+1][image->width-1] + image->b[i-1][image->width-1] + image->b[i][image->width-2]) / 3; - } - } - // last row - int b = (a==1 && image->height%2) || (a==0 && image->height%2==0); - for (int j=1+b; jwidth-1; j+=2) { - image->r[image->height-1][j] = (image->r[image->height-2][j] + image->r[image->height-1][j+1] + image->r[image->height-1][j-1]) / 3; - image->g[image->height-1][j] = (image->g[image->height-2][j] + image->g[image->height-1][j+1] + image->g[image->height-1][j-1]) / 3; - image->b[image->height-1][j] = (image->b[image->height-2][j] + image->b[image->height-1][j+1] + image->b[image->height-1][j-1]) / 3; - } - } - - // Flip if needed - if (tran & TR_HFLIP) - hflip (image); - if (tran & TR_VFLIP) - vflip (image); - - // Color correction - if (ri->filters && pp.skip==1) - correction_YIQ_LQ (image, settings->colorCorrectionSteps); - - // Applying postmul - colorSpaceConversion (image, cmp, embProfile, camProfile, cam, defGain); - - isrcMutex.unlock (); -} - - -void RawImageSource::cfa_clean(float thresh) //Emil's hot/dead pixel removal -- only filters egregiously impulsive pixels - { - // local variables - int rr, cc; - int gin, g[8]; - - float eps=1e-10;//tolerance to avoid dividing by zero - float p[8]; - float pixave, pixratio; - - // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //The cleaning algorithm starts here - - for (rr=4; rr < H-4; rr++) - for (cc=4; cc < W-4; cc++) { - - //pixel neighbor average - gin=ri->data[rr][cc]; - g[0]=ri->data[rr-2][cc-2]; - g[1]=ri->data[rr-2][cc]; - g[2]=ri->data[rr-2][cc+2]; - g[3]=ri->data[rr][cc-2]; - g[4]=ri->data[rr][cc+2]; - g[5]=ri->data[rr+2][cc-2]; - g[6]=ri->data[rr+2][cc]; - g[7]=ri->data[rr+2][cc+2]; - - pixave=(float)(g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6]+g[7])/8; - pixratio=MIN(gin,pixave)/(eps+MAX(gin,pixave)); - - if (pixratio > thresh) continue; - - p[0]=1/(eps+fabs(g[0]-gin)+fabs(g[0]-ri->data[rr-4][cc-4])+fabs(ri->data[rr-1][cc-1]-ri->data[rr-3][cc-3])); - p[1]=1/(eps+fabs(g[1]-gin)+fabs(g[1]-ri->data[rr-4][cc])+fabs(ri->data[rr-1][cc]-ri->data[rr-3][cc])); - p[2]=1/(eps+fabs(g[2]-gin)+fabs(g[2]-ri->data[rr-4][cc+4])+fabs(ri->data[rr-1][cc+1]-ri->data[rr-3][cc+3])); - p[3]=1/(eps+fabs(g[3]-gin)+fabs(g[3]-ri->data[rr][cc-4])+fabs(ri->data[rr][cc-1]-ri->data[rr][cc-3])); - p[4]=1/(eps+fabs(g[4]-gin)+fabs(g[4]-ri->data[rr][cc+4])+fabs(ri->data[rr][cc+1]-ri->data[rr][cc+3])); - p[5]=1/(eps+fabs(g[5]-gin)+fabs(g[5]-ri->data[rr+4][cc-4])+fabs(ri->data[rr+1][cc-1]-ri->data[rr+3][cc-3])); - p[6]=1/(eps+fabs(g[6]-gin)+fabs(g[6]-ri->data[rr+4][cc])+fabs(ri->data[rr+1][cc]-ri->data[rr+3][cc])); - p[7]=1/(eps+fabs(g[7]-gin)+fabs(g[7]-ri->data[rr+4][cc+4])+fabs(ri->data[rr+1][cc+1]-ri->data[rr+3][cc+3])); - - ri->data[rr][cc] = (int)((g[0]*p[0]+g[1]*p[1]+g[2]*p[2]+g[3]*p[3]+g[4]*p[4]+ \ - g[5]*p[5]+g[6]*p[6]+g[7]*p[7])/(p[0]+p[1]+p[2]+p[3]+p[4]+p[5]+p[6]+p[7])); - - }//now impulsive values have been corrected - - - // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - } - - -void RawImageSource::rotateLine (unsigned short* line, unsigned short** channel, int tran, int i, int w, int h) { - - if ((tran & TR_ROT) == TR_R180) - for (int j=0; j=0 && yheight && y>=0 && xwidth) { - image->r[image->height-1-y][image->width-1-x] = red[j]; - image->g[image->height-1-y][image->width-1-x] = green[j]; - image->b[image->height-1-y][image->width-1-x] = blue[j]; - } - } - } - else if ((tran & TR_ROT) == TR_R270) { - int end = MIN(h+fw-i, w-fw+i); - for (int j=start; j=0 && xheight && y>=0 && ywidth) { - image->r[image->height-1-x][y] = red[j]; - image->g[image->height-1-x][y] = green[j]; - image->b[image->height-1-x][y] = blue[j]; - } - } - } - else if ((tran & TR_ROT) == TR_R90) { - int end = MIN(h+fw-i, w-fw+i); - for (int j=start; j=0 && ywidth && y>=0 && xheight) { - image->r[x][image->width-1-y] = red[j]; - image->g[x][image->width-1-y] = green[j]; - image->b[x][image->width-1-y] = blue[j]; - } - } - } - else { - int end = MIN(h+fw-i, w-fw+i); - for (int j=start; j=0 && yheight && y>=0 && xwidth) { - image->r[y][x] = red[j]; - image->g[y][x] = green[j]; - image->b[y][x] = blue[j]; - } - } - } - } - // Nikon D1X vertical interpolation + coarse rotation - else if (d1x) { - // copy new pixels - if ((tran & TR_ROT) == TR_R180) { - for (int j=0; jr[2*imheight-2-2*i][imwidth-1-j] = red[j]; - image->g[2*imheight-2-2*i][imwidth-1-j] = green[j]; - image->b[2*imheight-2-2*i][imwidth-1-j] = blue[j]; - } - - if (i==1 || i==2) { // linear interpolation - int row = 2*imheight-1-2*i; - for (int j=0; jr[row][col] = (red[j] + image->r[row+1][col]) >> 1; - image->g[row][col] = (green[j] + image->g[row+1][col]) >> 1; - image->b[row][col] = (blue[j] + image->b[row+1][col]) >> 1; - } - } - else if (i==imheight-1) { - int row = 2*imheight-1-2*i; - for (int j=0; jr[row][col] = (red[j] + image->r[row+1][col]) >> 1; - image->g[row][col] = (green[j] + image->g[row+1][col]) >> 1; - image->b[row][col] = (blue[j] + image->b[row+1][col]) >> 1; - } - row = 2*imheight-1-2*i+2; - for (int j=0; jr[row][col] = (red[j] + image->r[row+1][col]) >> 1; - image->g[row][col] = (green[j] + image->g[row+1][col]) >> 1; - image->b[row][col] = (blue[j] + image->b[row+1][col]) >> 1; - } - } - else if (i>2 && ir[row][col] = CLIP((int)(-0.0625*red[j] + 0.5625*image->r[row-1][col] + 0.5625*image->r[row+1][col] - 0.0625*image->r[row+3][col])); - image->g[row][col] = CLIP((int)(-0.0625*green[j] + 0.5625*image->g[row-1][col] + 0.5625*image->g[row+1][col] - 0.0625*image->g[row+3][col])); - image->b[row][col] = CLIP((int)(-0.0625*blue[j] + 0.5625*image->b[row-1][col] + 0.5625*image->b[row+1][col] - 0.0625*image->b[row+3][col])); - } - } - } - else if ((tran & TR_ROT) == TR_R90) { - for (int j=0; jr[j][2*imheight-2-2*i] = red[j]; - image->g[j][2*imheight-2-2*i] = green[j]; - image->b[j][2*imheight-2-2*i] = blue[j]; - } - if (i==1 || i==2) { // linear interpolation - int col = 2*imheight-1-2*i; - for (int j=0; jr[j][col] = (red[j] + image->r[j][col+1]) >> 1; - image->g[j][col] = (green[j] + image->g[j][col+1]) >> 1; - image->b[j][col] = (blue[j] + image->b[j][col+1]) >> 1; - } - } - else if (i==imheight-1) { - int col = 2*imheight-1-2*i; - for (int j=0; jr[j][col] = (red[j] + image->r[j][col+1]) >> 1; - image->g[j][col] = (green[j] + image->g[j][col+1]) >> 1; - image->b[j][col] = (blue[j] + image->b[j][col+1]) >> 1; - } - col = 2*imheight-1-2*i+2; - for (int j=0; jr[j][col] = (red[j] + image->r[j][col+1]) >> 1; - image->g[j][col] = (green[j] + image->g[j][col+1]) >> 1; - image->b[j][col] = (blue[j] + image->b[j][col+1]) >> 1; - } - } - else if (i>2 && ir[j][col] = CLIP((int)(-0.0625*red[j] + 0.5625*image->r[j][col-1] + 0.5625*image->r[j][col+1] - 0.0625*image->r[j][col+3])); - image->g[j][col] = CLIP((int)(-0.0625*green[j] + 0.5625*image->g[j][col-1] + 0.5625*image->g[j][col+1] - 0.0625*image->g[j][col+3])); - image->b[j][col] = CLIP((int)(-0.0625*blue[j] + 0.5625*image->b[j][col-1] + 0.5625*image->b[j][col+1] - 0.0625*image->b[j][col+3])); - } - } - } - else if ((tran & TR_ROT) == TR_R270) { - for (int j=0; jr[imwidth-1-j][2*i] = red[j]; - image->g[imwidth-1-j][2*i] = green[j]; - image->b[imwidth-1-j][2*i] = blue[j]; - } - if (i==1 || i==2) { // linear interpolation - for (int j=0; jr[row][2*i-1] = (red[j] + image->r[row][2*i-2]) >> 1; - image->g[row][2*i-1] = (green[j] + image->g[row][2*i-2]) >> 1; - image->b[row][2*i-1] = (blue[j] + image->b[row][2*i-2]) >> 1; - } - } - else if (i==imheight-1) { - for (int j=0; jr[row][2*i-1] = (red[j] + image->r[row][2*i-2]) >> 1; - image->g[row][2*i-1] = (green[j] + image->g[row][2*i-2]) >> 1; - image->b[row][2*i-1] = (blue[j] + image->b[row][2*i-2]) >> 1; - image->r[row][2*i-3] = (image->r[row][2*i-2] + image->r[row][2*i-4]) >> 1; - image->g[row][2*i-3] = (image->g[row][2*i-2] + image->g[row][2*i-4]) >> 1; - image->b[row][2*i-3] = (image->b[row][2*i-2] + image->b[row][2*i-4]) >> 1; - } - } - else if (i>0 && ir[row][2*i-3] = CLIP((int)(-0.0625*red[j] + 0.5625*image->r[row][2*i-2] + 0.5625*image->r[row][2*i-4] - 0.0625*image->r[row][2*i-6])); - image->g[row][2*i-3] = CLIP((int)(-0.0625*green[j] + 0.5625*image->g[row][2*i-2] + 0.5625*image->g[row][2*i-4] - 0.0625*image->g[row][2*i-6])); - image->b[row][2*i-3] = CLIP((int)(-0.0625*blue[j] + 0.5625*image->b[row][2*i-2] + 0.5625*image->b[row][2*i-4] - 0.0625*image->b[row][2*i-6])); - } - } - } - else { - rotateLine (red, image->r, tran, 2*i, imwidth, imheight); - rotateLine (green, image->g, tran, 2*i, imwidth, imheight); - rotateLine (blue, image->b, tran, 2*i, imwidth, imheight); - - if (i==1 || i==2) { // linear interpolation - for (int j=0; jr[2*i-1][j] = (red[j] + image->r[2*i-2][j]) >> 1; - image->g[2*i-1][j] = (green[j] + image->g[2*i-2][j]) >> 1; - image->b[2*i-1][j] = (blue[j] + image->b[2*i-2][j]) >> 1; - } - } - else if (i==imheight-1) { - for (int j=0; jr[2*i-3][j] = (image->r[2*i-4][j] + image->r[2*i-2][j]) >> 1; - image->g[2*i-3][j] = (image->g[2*i-4][j] + image->g[2*i-2][j]) >> 1; - image->b[2*i-3][j] = (image->b[2*i-4][j] + image->b[2*i-2][j]) >> 1; - image->r[2*i-1][j] = (red[j] + image->r[2*i-2][j]) >> 1; - image->g[2*i-1][j] = (green[j] + image->g[2*i-2][j]) >> 1; - image->b[2*i-1][j] = (blue[j] + image->b[2*i-2][j]) >> 1; - } - } - else if (i>2 && ir[2*i-3][j] = CLIP((int)(-0.0625*red[j] + 0.5625*image->r[2*i-2][j] + 0.5625*image->r[2*i-4][j] - 0.0625*image->r[2*i-6][j])); - image->g[2*i-3][j] = CLIP((int)(-0.0625*green[j] + 0.5625*image->g[2*i-2][j] + 0.5625*image->g[2*i-4][j] - 0.0625*image->g[2*i-6][j])); - image->b[2*i-3][j] = CLIP((int)(-0.0625*blue[j] + 0.5625*image->b[2*i-2][j] + 0.5625*image->b[2*i-4][j] - 0.0625*image->b[2*i-6][j])); - } - } - } - } - // other (conventional) CCD coarse rotation - else { - rotateLine (red, image->r, tran, i, imwidth, imheight); - rotateLine (green, image->g, tran, i, imwidth, imheight); - rotateLine (blue, image->b, tran, i, imwidth, imheight); - } -} - -void RawImageSource::getFullSize (int& w, int& h, int tr) { - - tr = defTransform (tr); - - if (fuji) { - w = ri->fuji_width * 2 + 1; - h = (H - ri->fuji_width)*2 + 1; - } - else if (d1x) { - w = W; - h = 2*H-1; - } - else { - w = W; - h = H; - } - - if ((tr & TR_ROT) == TR_R90 || (tr & TR_ROT) == TR_R270) { - int tmp = w; - w = h; - h = tmp; - } - w -= 2 * border; - h -= 2 * border; -} - -void RawImageSource::getSize (int tran, PreviewProps pp, int& w, int& h) { - - tran = defTransform (tran); - -// if (fuji) { -// return; -// } -// else if (d1x) { -// return; -// } -// else { - w = pp.w / pp.skip + (pp.w % pp.skip > 0); - h = pp.h / pp.skip + (pp.h % pp.skip > 0); -// } -} - -void RawImageSource::hflip (Image16* image) { - int width = image->width; - int height = image->height; - - unsigned short* rowr = new unsigned short[width]; - unsigned short* rowg = new unsigned short[width]; - unsigned short* rowb = new unsigned short[width]; - for (int i=0; ir[i][width-1-j]; - rowg[j] = image->g[i][width-1-j]; - rowb[j] = image->b[i][width-1-j]; - } - memcpy (image->r[i], rowr, width*sizeof(unsigned short)); - memcpy (image->g[i], rowg, width*sizeof(unsigned short)); - memcpy (image->b[i], rowb, width*sizeof(unsigned short)); - } - delete [] rowr; - delete [] rowg; - delete [] rowb; -} - -void RawImageSource::vflip (Image16* image) { - int width = image->width; - int height = image->height; - - register unsigned short tmp; - for (int i=0; ir[i][j]; - image->r[i][j] = image->r[height-1-i][j]; - image->r[height-1-i][j] = tmp; - tmp = image->g[i][j]; - image->g[i][j] = image->g[height-1-i][j]; - image->g[height-1-i][j] = tmp; - tmp = image->b[i][j]; - image->b[i][j] = image->b[height-1-i][j]; - image->b[height-1-i][j] = tmp; - } -} - -void RawImageSource::inverse33 (double (*coeff)[3], double (*icoeff)[3]) { - double nom = coeff[0][2]*coeff[1][1]*coeff[2][0] - coeff[0][1]*coeff[1][2]*coeff[2][0] - coeff[0][2]*coeff[1][0]*coeff[2][1] + coeff[0][0]*coeff[1][2]*coeff[2][1] + coeff[0][1]*coeff[1][0]*coeff[2][2] - coeff[0][0]*coeff[1][1]*coeff[2][2]; - icoeff[0][0] = (coeff[1][2]*coeff[2][1]-coeff[1][1]*coeff[2][2]) / nom; - icoeff[0][1] = -(coeff[0][2]*coeff[2][1]-coeff[0][1]*coeff[2][2]) / nom; - icoeff[0][2] = (coeff[0][2]*coeff[1][1]-coeff[0][1]*coeff[1][2]) / nom; - icoeff[1][0] = -(coeff[1][2]*coeff[2][0]-coeff[1][0]*coeff[2][2]) / nom; - icoeff[1][1] = (coeff[0][2]*coeff[2][0]-coeff[0][0]*coeff[2][2]) / nom; - icoeff[1][2] = -(coeff[0][2]*coeff[1][0]-coeff[0][0]*coeff[1][2]) / nom; - icoeff[2][0] = (coeff[1][1]*coeff[2][0]-coeff[1][0]*coeff[2][1]) / nom; - icoeff[2][1] = -(coeff[0][1]*coeff[2][0]-coeff[0][0]*coeff[2][1]) / nom; - icoeff[2][2] = (coeff[0][1]*coeff[1][0]-coeff[0][0]*coeff[1][1]) / nom; -} - -int RawImageSource::load (Glib::ustring fname) { - - fileName = fname; - - if (plistener) { - plistener->setProgressStr ("Decoding..."); - plistener->setProgress (0.0); - } - - ri = new RawImage; - int res = loadRaw (fname.c_str(), ri); - if (res) - return res; - - if(red) { - delete red; - red = 0; - } - if(green) { - delete green; - green = 0; - } - if(blue) { - delete blue; - blue = 0; - } - - W = ri->width; - H = ri->height; - - d1x = !strcmp(ri->model, "D1X"); - fuji = ri->fuji_width; - if (d1x) - border = 8; - - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - coeff[i][j] = ri->coeff[i][j]; - - // compute inverse of the color transformation matrix - inverse33 (coeff, icoeff); - - double cam_r = coeff[0][0]*ri->camwb_red + coeff[0][1]*ri->camwb_green + coeff[0][2]*ri->camwb_blue; - double cam_g = coeff[1][0]*ri->camwb_red + coeff[1][1]*ri->camwb_green + coeff[1][2]*ri->camwb_blue; - double cam_b = coeff[2][0]*ri->camwb_red + coeff[2][1]*ri->camwb_green + coeff[2][2]*ri->camwb_blue; - - wb = ColorTemp (cam_r, cam_g, cam_b); - - double tr = icoeff[0][0] * cam_r + icoeff[0][1] * cam_g + icoeff[0][2] * cam_b; - double tg = icoeff[1][0] * cam_r + icoeff[1][1] * cam_g + icoeff[1][2] * cam_b; - double tb = icoeff[2][0] * cam_r + icoeff[2][1] * cam_g + icoeff[2][2] * cam_b; - - // create profile - memset (cam, 0, sizeof(cam)); - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - for (int k=0; k<3; k++) - cam[i][j] += coeff[k][i] * sRGB_d50[k][j]; - camProfile = iccStore.createFromMatrix (cam, false, "Camera"); - inverse33 (cam, icam); - - if (ri->profile_data) - embProfile = cmsOpenProfileFromMem (ri->profile_data, ri->profile_len); - - defGain = log(ri->defgain) / log(2.0); - - RawMetaDataLocation rml; - rml.exifBase = ri->exifbase; - rml.ciffBase = ri->ciff_base; - rml.ciffLength = ri->ciff_len; - - idata = new ImageData (fname, &rml); - - // check if it is an olympus E camera, if yes, compute G channel pre-compensation factors - if (settings->greenthresh || (((idata->getMake().size()>=7 && idata->getMake().substr(0,7)=="OLYMPUS" && idata->getModel()[0]=='E') || (idata->getMake().size()>=9 && idata->getMake().substr(0,7)=="Panasonic")) && settings->demosaicMethod!="vng4" && ri->filters) ) { - // global correction - int ng1=0, ng2=0; - double avgg1=0, avgg2=0; - for (int i=border; idata[i][j]; - ng1++; - } - else { - avgg2 += ri->data[i][j]; - ng2++; - } - } - double corrg1 = ((double)avgg1/ng1 + (double)avgg2/ng2) / 2.0 / ((double)avgg1/ng1); - double corrg2 = ((double)avgg1/ng1 + (double)avgg2/ng2) / 2.0 / ((double)avgg2/ng2); - for (int i=border; idata[i][j] = CLIP(ri->data[i][j] * (i%2 ? corrg2 : corrg1)); - } - - - // local correction in a 9x9 box - if (settings->greenthresh) { -//Emil's green equilbration - if (plistener) { - plistener->setProgressStr ("Green equilibrate..."); - plistener->setProgress (0.0); - } - green_equilibrate(0.01*(settings->greenthresh)); - -/* unsigned short* corr_alloc = new unsigned short[W*H]; - unsigned short** corr_data = new unsigned short* [H]; - for (int i=0; iallocation, W*H*sizeof(unsigned short)); - for (int i=border; idata[i-4][j-4] + ri->data[i-4][j-2] + ri->data[i-4][j] + ri->data[i-4][j+2] + - ri->data[i-2][j-4] + ri->data[i-2][j-2] + ri->data[i-2][j] + ri->data[i-2][j+2] + - ri->data[i][j-4] + ri->data[i][j-2] + ri->data[i][j] + ri->data[i][j+2] + - ri->data[i+2][j-4] + ri->data[i+2][j-2] + ri->data[i+2][j] + ri->data[i+2][j+2]; - unsigned int ag2 = ri->data[i-3][j-3] + ri->data[i-3][j-1] + ri->data[i-3][j+1] + ri->data[i-3][j+1] + - ri->data[i-1][j-3] + ri->data[i-1][j-1] + ri->data[i-1][j+1] + ri->data[i-1][j+1] + - ri->data[i+1][j-3] + ri->data[i+1][j-1] + ri->data[i+1][j+1] + ri->data[i+1][j+1] + - ri->data[i+3][j-3] + ri->data[i+3][j-1] + ri->data[i+3][j+1] + ri->data[i+3][j+1]; - unsigned int val = (ri->data[i][j] + ri->data[i][j] * ag2 / ag1) / 2; - corr_data[i][j] = CLIP (val); - } - memcpy (ri->allocation, corr_alloc, W*H*sizeof(unsigned short)); - delete corr_alloc; - delete corr_data; */ - - } - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//Emil's CFA hot/dead pixel filter - - if (settings->hotdeadpix_filt) { - if (plistener) { - plistener->setProgressStr ("Hot/Dead Pixel Filter..."); - plistener->setProgress (0.0); - } - - cfa_clean(0.1); - } - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//Emil's line noise filter - - if (settings->linenoise) { - if (plistener) { - plistener->setProgressStr ("Line Denoise..."); - plistener->setProgress (0.0); - } - - cfa_linedn(0.00002*(settings->linenoise)); - } - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//Emil's CA auto correction - - if (settings->ca_autocorrect) { - if (plistener) { - plistener->setProgressStr ("CA Auto Correction..."); - plistener->setProgress (0.0); - } - - CA_correct_RT(); - } -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - if (ri->filters) { - //MyTime t1,t2; - //t1.set(); - // demosaic - if (settings->demosaicMethod=="hphd") - hphd_demosaic (); - else if (settings->demosaicMethod=="vng4") - vng4_demosaic (); - else if (settings->demosaicMethod=="ahd") - ahd_demosaic (); - else if (settings->demosaicMethod=="bilinear") - bilinear_demosaic(); - //else if (settings->demosaicMethod=="ppg") - // ppg_demosaic (); - else if (settings->demosaicMethod=="amaze") - amaze_demosaic_RT ();//Emil's code for AMaZE - else if (settings->demosaicMethod=="dcb") - dcb_demosaic(settings->dcb_iterations, settings->dcb_enhance? 1:0); - else - eahd_demosaic (); - //t2.set(); - //printf("Demosaicing:%d usec\n",t2.etime(t1)); - } - - - if (plistener) { - plistener->setProgressStr ("Ready."); - plistener->setProgress (1.0); - } - - return 0; -} - -int RawImageSource::defTransform (int tran) { - - int deg = ri->rotate_deg; - if ((tran & TR_ROT) == TR_R180) - deg += 180; - else if ((tran & TR_ROT) == TR_R90) - deg += 90; - else if ((tran & TR_ROT) == TR_R270) - deg += 270; - deg %= 360; - - int ret = 0; - if (deg==90) - ret |= TR_R90; - else if (deg==180) - ret |= TR_R180; - else if (deg==270) - ret |= TR_R270; - if (tran & TR_HFLIP) - ret |= TR_HFLIP; - if (tran & TR_VFLIP) - ret |= TR_VFLIP; - return ret; -} - -void RawImageSource::correction_YIQ_LQ_ (Image16* im, int row_from, int row_to) { - - int W = im->width; - - int** rbconv_Y = new int*[3]; - int** rbconv_I = new int*[3]; - int** rbconv_Q = new int*[3]; - int** rbout_I = new int*[3]; - int** rbout_Q = new int*[3]; - for (int i=0; i<3; i++) { - rbconv_Y[i] = new int[W]; - rbconv_I[i] = new int[W]; - rbconv_Q[i] = new int[W]; - rbout_I[i] = new int[W]; - rbout_Q[i] = new int[W]; - } - - int* row_I = new int[W]; - int* row_Q = new int[W]; - - int* pre1_I = new int[3]; - int* pre2_I = new int[3]; - int* post1_I = new int[3]; - int* post2_I = new int[3]; - int middle_I[6]; - int* pre1_Q = new int[3]; - int* pre2_Q = new int[3]; - int* post1_Q = new int[3]; - int* post2_Q = new int[3]; - int middle_Q[6]; - int* tmp; - - int ppx=0, px=(row_from-1)%3, cx=row_from%3, nx=0; - - convert_row_to_YIQ (im->r[row_from-1], im->g[row_from-1], im->b[row_from-1], rbconv_Y[px], rbconv_I[px], rbconv_Q[px], W); - convert_row_to_YIQ (im->r[row_from], im->g[row_from], im->b[row_from], rbconv_Y[cx], rbconv_I[cx], rbconv_Q[cx], W); - - for (int j=0; jr[i+1], im->g[i+1], im->b[i+1], rbconv_Y[nx], rbconv_I[nx], rbconv_Q[nx], W); - - SORT3(rbconv_I[px][0],rbconv_I[cx][0],rbconv_I[nx][0],pre1_I[0],pre1_I[1],pre1_I[2]); - SORT3(rbconv_I[px][1],rbconv_I[cx][1],rbconv_I[nx][1],pre2_I[0],pre2_I[1],pre2_I[2]); - SORT3(rbconv_Q[px][0],rbconv_Q[cx][0],rbconv_Q[nx][0],pre1_Q[0],pre1_Q[1],pre1_Q[2]); - SORT3(rbconv_Q[px][1],rbconv_Q[cx][1],rbconv_Q[nx][1],pre2_Q[0],pre2_Q[1],pre2_Q[2]); - - // median I channel - for (int j=1; jrow_from) { - for (int j=1; jr[i-1], im->g[i-1], im->b[i-1], rbconv_Y[px], row_I, row_Q, W); - } - } - // blur last 3 row and finalize H-1th row - for (int j=1; jr[row_to-1], im->g[row_to-1], im->b[row_to-1], rbconv_Y[cx], row_I, row_Q, W); - - freeArray(rbconv_Y, 3); - freeArray(rbconv_I, 3); - freeArray(rbconv_Q, 3); - freeArray(rbout_I, 3); - freeArray(rbout_Q, 3); - delete [] row_I; - delete [] row_Q; - delete [] pre1_I; - delete [] pre2_I; - delete [] post1_I; - delete [] post2_I; - delete [] pre1_Q; - delete [] pre2_Q; - delete [] post1_Q; - delete [] post2_Q; -} - - -void RawImageSource::correction_YIQ_LQ (Image16* im, int times) { - - if (im->height<4) - return; - - for (int t=0; theight-2)/nthreads; - - if (tidheight - 1); - } -#else - correction_YIQ_LQ_ (im, 1 , im->height - 1); -#endif - } -} - - -void RawImageSource::colorSpaceConversion (Image16* im, ColorManagementParams cmp, cmsHPROFILE embedded, cmsHPROFILE camprofile, double camMatrix[3][3], double& defgain) { - - if (cmp.input == "(none)") - return; - - MyTime t1, t2, t3; - - t1.set (); - - cmsHPROFILE in; - cmsHPROFILE out; - - Glib::ustring inProfile = cmp.input; - - if (inProfile=="(embedded)") { - if (embedded) - in = embedded; - else - in = camprofile; - } - else if (inProfile=="(camera)" || inProfile=="") - in = camprofile; - else { - in = iccStore.getProfile (inProfile); - if (in==NULL) - inProfile = "(camera)"; - } - - - if (inProfile=="(camera)" || inProfile=="" || (inProfile=="(embedded)" && !embedded)) { - // in this case we avoid using the slllllooooooowwww lcms - -// out = iccStore.workingSpace (wProfile); -// hTransform = cmsCreateTransform (in, TYPE_RGB_16_PLANAR, out, TYPE_RGB_16_PLANAR, settings->colorimetricIntent, cmsFLAGS_MATRIXINPUT | cmsFLAGS_MATRIXOUTPUT);//cmsFLAGS_MATRIXINPUT | cmsFLAGS_MATRIXOUTPUT); -// cmsDoTransform (hTransform, im->data, im->data, im->planestride/2); -// cmsDeleteTransform(hTransform); - TMatrix work = iccStore.workingSpaceInverseMatrix (cmp.working); - double mat[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}; - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - for (int k=0; k<3; k++) - mat[i][j] += camMatrix[i][k] * work[k][j]; - - for (int i=0; iheight; i++) - for (int j=0; jwidth; j++) { - - int newr = mat[0][0]*im->r[i][j] + mat[1][0]*im->g[i][j] + mat[2][0]*im->b[i][j]; - int newg = mat[0][1]*im->r[i][j] + mat[1][1]*im->g[i][j] + mat[2][1]*im->b[i][j]; - int newb = mat[0][2]*im->r[i][j] + mat[1][2]*im->g[i][j] + mat[2][2]*im->b[i][j]; - - im->r[i][j] = CLIP(newr); - im->g[i][j] = CLIP(newg); - im->b[i][j] = CLIP(newb); - } - } - else { - out = iccStore.workingSpace (cmp.working); -// out = iccStore.workingSpaceGamma (wProfile); - lcmsMutex->lock (); - cmsHTRANSFORM hTransform = cmsCreateTransform (in, TYPE_RGB_16_PLANAR, out, TYPE_RGB_16_PLANAR, settings->colorimetricIntent, 0); - lcmsMutex->unlock (); - if (hTransform) { - if (cmp.gammaOnInput) { - double gd = pow (2.0, defgain); - defgain = 0.0; - for (int i=0; iheight; i++) - for (int j=0; jwidth; j++) { - im->r[i][j] = CurveFactory::gamma (CLIP(defgain*im->r[i][j])); - im->g[i][j] = CurveFactory::gamma (CLIP(defgain*im->g[i][j])); - im->b[i][j] = CurveFactory::gamma (CLIP(defgain*im->b[i][j])); - } - } - cmsDoTransform (hTransform, im->data, im->data, im->planestride/2); - } - else { - lcmsMutex->lock (); - hTransform = cmsCreateTransform (camprofile, TYPE_RGB_16_PLANAR, out, TYPE_RGB_16_PLANAR, settings->colorimetricIntent, 0); - lcmsMutex->unlock (); - cmsDoTransform (hTransform, im->data, im->data, im->planestride/2); - } - cmsDeleteTransform(hTransform); - } - t3.set (); -// printf ("ICM TIME: %d\n", t3.etime(t1)); -} - -void RawImageSource::eahd_demosaic () { - - if (plistener) { - plistener->setProgressStr ("Demosaicing..."); - plistener->setProgress (0.0); - } - - // prepare chache and constants for cielab conversion - lc00 = (0.412453 * coeff[0][0] + 0.357580 * coeff[1][0] + 0.180423 * coeff[2][0]) / 0.950456; - lc01 = (0.412453 * coeff[0][1] + 0.357580 * coeff[1][1] + 0.180423 * coeff[2][1]) / 0.950456; - lc02 = (0.412453 * coeff[0][2] + 0.357580 * coeff[1][2] + 0.180423 * coeff[2][2]) / 0.950456; - - lc10 = 0.212671 * coeff[0][0] + 0.715160 * coeff[1][0] + 0.072169 * coeff[2][0]; - lc11 = 0.212671 * coeff[0][1] + 0.715160 * coeff[1][1] + 0.072169 * coeff[2][1]; - lc12 = 0.212671 * coeff[0][2] + 0.715160 * coeff[1][2] + 0.072169 * coeff[2][2]; - - lc20 = (0.019334 * coeff[0][0] + 0.119193 * coeff[1][0] + 0.950227 * coeff[2][0]) / 1.088754; - lc21 = (0.019334 * coeff[0][1] + 0.119193 * coeff[1][1] + 0.950227 * coeff[2][1]) / 1.088754; - lc22 = (0.019334 * coeff[0][2] + 0.119193 * coeff[1][2] + 0.950227 * coeff[2][2]) / 1.088754; - - int maxindex = 2*65536; - cache = new double[maxindex]; - threshold = (int)(0.008856*CMAXVAL); - for (int i=0; isv) { // fixate horizontal pixel - dLmaph[dmi] = DIST(lLh[ix][j], lLh[idx][j+y]); - dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j+y]); - dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j+y]); - dLmapv[dmi] = DIST(lLv[ix][j], lLh[idx][j+y]); - dCamapv[dmi] = DIST(lav[ix][j], lah[idx][j+y]); - dCbmapv[dmi] = DIST(lbv[ix][j], lbh[idx][j+y]); - } - else if (shdata[i-1][j]; - else { - hc = homh[imx][j]; - vc = homv[imx][j]; - if (hc > vc) - green[i-1][j] = gh[(i-1)%4][j]; - else if (hc < vc) - green[i-1][j] = gv[(i-1)%4][j]; - else - green[i-1][j] = (gh[(i-1)%4][j] + gv[(i-1)%4][j]) / 2; - } - } - - if (!(i%20) && plistener) - plistener->setProgress ((double)i / (H-2)); - } - // finish H-2th and H-1th row, homogenity value is still valailable - int hc, vc; - for (int i=H-1; i vc) - green[i-1][j] = gh[(i-1)%4][j]; - else if (hc < vc) - green[i-1][j] = gv[(i-1)%4][j]; - else - green[i-1][j] = (gh[(i-1)%4][j] + gv[(i-1)%4][j]) / 2; - } - - freeArray2(rh, 3); - freeArray2(gh, 4); - freeArray2(bh, 3); - freeArray2(rv, 3); - freeArray2(gv, 4); - freeArray2(bv, 3); - freeArray2(lLh, 3); - freeArray2(lah, 3); - freeArray2(lbh, 3); - freeArray2(homh, 3); - freeArray2(lLv, 3); - freeArray2(lav, 3); - freeArray2(lbv, 3); - freeArray2(homv, 3); -} - -void RawImageSource::hphd_vertical (float** hpmap, int col_from, int col_to) { - - float* temp = new float[MAX(W,H)]; - float* avg = new float[MAX(W,H)]; - float* dev = new float[MAX(W,H)]; - - memset (temp, 0, MAX(W,H)*sizeof(float)); - memset (avg, 0, MAX(W,H)*sizeof(float)); - memset (dev, 0, MAX(W,H)*sizeof(float)); - - for (int k=col_from; kdata[i-5][k] - 8*ri->data[i-4][k] + 27*ri->data[i-3][k] - 48*ri->data[i-2][k] + 42*ri->data[i-1][k] - - (ri->data[i+5][k] - 8*ri->data[i+4][k] + 27*ri->data[i+3][k] - 48*ri->data[i+2][k] + 42*ri->data[i+1][k])) / 100.0; - temp[i] = ABS(temp[i]); - } - for (int j=4; jdata[i][j-5] - 8*ri->data[i][j-4] + 27*ri->data[i][j-3] - 48*ri->data[i][j-2] + 42*ri->data[i][j-1] - - (ri->data[i][j+5] - 8*ri->data[i][j+4] + 27*ri->data[i][j+3] - 48*ri->data[i][j+2] + 42*ri->data[i][j+1])) / 100; - temp[j] = ABS(temp[j]); - } - for (int j=4; jhpmap[i][j] = 2; - else if (hpv < 0.8*hpmap[i][j]) - this->hpmap[i][j] = 1; - else - this->hpmap[i][j] = 0; - } - } - delete [] temp; - delete [] avg; - delete [] dev; -} - -void RawImageSource::hphd_green () { - - #pragma omp parallel for - for (int i=3; idata[i][j]; - else { - if (this->hpmap[i][j]==1) { - int g2 = ri->data[i][j+1] + ((ri->data[i][j] - ri->data[i][j+2]) >> 1); - int g4 = ri->data[i][j-1] + ((ri->data[i][j] - ri->data[i][j-2]) >> 1); - - int dx = ri->data[i][j+1] - ri->data[i][j-1]; - int d1 = ri->data[i][j+3] - ri->data[i][j+1]; - int d2 = ri->data[i][j+2] - ri->data[i][j]; - int d3 = (ri->data[i-1][j+2] - ri->data[i-1][j]) >> 1; - int d4 = (ri->data[i+1][j+2] - ri->data[i+1][j]) >> 1; - - double e2 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - - d1 = ri->data[i][j-3] - ri->data[i][j-1]; - d2 = ri->data[i][j-2] - ri->data[i][j]; - d3 = (ri->data[i-1][j-2] - ri->data[i-1][j]) >> 1; - d4 = (ri->data[i+1][j-2] - ri->data[i+1][j]) >> 1; - - double e4 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - - green[i][j] = CLIP((e2 * g2 + e4 * g4) / (e2 + e4)); - } - else if (this->hpmap[i][j]==2) { - int g1 = ri->data[i-1][j] + ((ri->data[i][j] - ri->data[i-2][j]) >> 1); - int g3 = ri->data[i+1][j] + ((ri->data[i][j] - ri->data[i+2][j]) >> 1); - - int dy = ri->data[i+1][j] - ri->data[i-1][j]; - int d1 = ri->data[i-1][j] - ri->data[i-3][j]; - int d2 = ri->data[i][j] - ri->data[i-2][j]; - int d3 = (ri->data[i][j-1] - ri->data[i-2][j-1]) >> 1; - int d4 = (ri->data[i][j+1] - ri->data[i-2][j+1]) >> 1; - - double e1 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - - d1 = ri->data[i+1][j] - ri->data[i+3][j]; - d2 = ri->data[i][j] - ri->data[i+2][j]; - d3 = (ri->data[i][j-1] - ri->data[i+2][j-1]) >> 1; - d4 = (ri->data[i][j+1] - ri->data[i+2][j+1]) >> 1; - - double e3 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - - green[i][j] = CLIP((e1 * g1 + e3 * g3) / (e1 + e3)); - } - else { - int g1 = ri->data[i-1][j] + ((ri->data[i][j] - ri->data[i-2][j]) >> 1); - int g2 = ri->data[i][j+1] + ((ri->data[i][j] - ri->data[i][j+2]) >> 1); - int g3 = ri->data[i+1][j] + ((ri->data[i][j] - ri->data[i+2][j]) >> 1); - int g4 = ri->data[i][j-1] + ((ri->data[i][j] - ri->data[i][j-2]) >> 1); - - int dx = ri->data[i][j+1] - ri->data[i][j-1]; - int dy = ri->data[i+1][j] - ri->data[i-1][j]; - - int d1 = ri->data[i-1][j] - ri->data[i-3][j]; - int d2 = ri->data[i][j] - ri->data[i-2][j]; - int d3 = (ri->data[i][j-1] - ri->data[i-2][j-1]) >> 1; - int d4 = (ri->data[i][j+1] - ri->data[i-2][j+1]) >> 1; - - double e1 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - - d1 = ri->data[i][j+3] - ri->data[i][j+1]; - d2 = ri->data[i][j+2] - ri->data[i][j]; - d3 = (ri->data[i-1][j+2] - ri->data[i-1][j]) >> 1; - d4 = (ri->data[i+1][j+2] - ri->data[i+1][j]) >> 1; - - double e2 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - - d1 = ri->data[i+1][j] - ri->data[i+3][j]; - d2 = ri->data[i][j] - ri->data[i+2][j]; - d3 = (ri->data[i][j-1] - ri->data[i+2][j-1]) >> 1; - d4 = (ri->data[i][j+1] - ri->data[i+2][j+1]) >> 1; - - double e3 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - - d1 = ri->data[i][j-3] - ri->data[i][j-1]; - d2 = ri->data[i][j-2] - ri->data[i][j]; - d3 = (ri->data[i-1][j-2] - ri->data[i-1][j]) >> 1; - d4 = (ri->data[i+1][j-2] - ri->data[i+1][j]) >> 1; - - double e4 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - - green[i][j] = CLIP((e1*g1 + e2*g2 + e3*g3 + e4*g4) / (e1 + e2 + e3 + e4)); - } - } - } - } -} - -void RawImageSource::hphd_demosaic () { - - if (plistener) { - plistener->setProgressStr ("Demosaicing..."); - plistener->setProgress (0.0); - } - - float** hpmap = new float*[H]; - for (int i=0; isetProgress (0.33); - - this->hpmap = allocArray(W, H); - for (int i=0; ihpmap[i], 0, W*sizeof(char)); - -#ifdef _OPENMP - #pragma omp parallel - { - int tid = omp_get_thread_num(); - int nthreads = omp_get_num_threads(); - int blk = H/nthreads; - - if (tid(hpmap, H); - - if (plistener) - plistener->setProgress (0.66); - - green = new unsigned short*[H]; - for (int i=0; isetProgress (1.0); -} - -void RawImageSource::HLRecovery_Luminance (unsigned short* rin, unsigned short* gin, unsigned short* bin, unsigned short* rout, unsigned short* gout, unsigned short* bout, int width, int maxval) { - - for (int i=0; imaxval || g>maxval || b>maxval) { - int ro = MIN (r, maxval); - int go = MIN (g, maxval); - int bo = MIN (b, maxval); - double L = r + g + b; - double C = 1.732050808 * (r - g); - double H = 2 * b - r - g; - double Co = 1.732050808 * (ro - go); - double Ho = 2 * bo - ro - go; - if (r!=g && g!=b) { - double ratio = sqrt ((Co*Co+Ho*Ho) / (C*C+H*H)); - C *= ratio; - H *= ratio; - } - int rr = L / 3.0 - H / 6.0 + C / 3.464101615; - int gr = L / 3.0 - H / 6.0 - C / 3.464101615; - int br = L / 3.0 + H / 3.0; - rout[i] = CLIP(rr); - gout[i] = CLIP(gr); - bout[i] = CLIP(br); - } - else { - rout[i] = rin[i]; - gout[i] = gin[i]; - bout[i] = bin[i]; - } - } -} - -void RawImageSource::HLRecovery_CIELab (unsigned short* rin, unsigned short* gin, unsigned short* bin, unsigned short* rout, unsigned short* gout, unsigned short* bout, int width, int maxval, double cam[3][3], double icam[3][3]) { - - static bool crTableReady = false; - static double fv[0x10000]; - if (!crTableReady) { - for (int ix=0; ix < 0x10000; ix++) { - double rx = ix / 65535.0; - fv[ix] = rx > 0.008856 ? exp(1.0/3 * log(rx)) : 7.787*rx + 16/116.0; - } - crTableReady = true; - } - - for (int i=0; imaxval || g>maxval || b>maxval) { - int ro = MIN (r, maxval); - int go = MIN (g, maxval); - int bo = MIN (b, maxval); - double yy = cam[0][1]*r + cam[1][1]*g + cam[2][1]*b; - double fy = fv[CLIP((int)yy)]; - // compute LCH decompostion of the clipped pixel (only color information, thus C and H will be used) - double x = cam[0][0]*ro + cam[1][0]*go + cam[2][0]*bo; - double y = cam[0][1]*ro + cam[1][1]*go + cam[2][1]*bo; - double z = cam[0][2]*ro + cam[1][2]*go + cam[2][2]*bo; - x = fv[CLIP((int)x)]; - y = fv[CLIP((int)y)]; - z = fv[CLIP((int)z)]; - // convert back to rgb - double fz = fy - y + z; - double fx = fy + x - y; - double zr = (fz<=0.206893) ? ((116.0*fz-16.0)/903.3) : (fz * fz * fz); - double xr = (fx<=0.206893) ? ((116.0*fx-16.0)/903.3) : (fx * fx * fx); - x = xr*65535.0 - 0.5; - y = yy; - z = zr*65535.0 - 0.5; - int rr = icam[0][0]*x + icam[1][0]*y + icam[2][0]*z; - int gr = icam[0][1]*x + icam[1][1]*y + icam[2][1]*z; - int br = icam[0][2]*x + icam[1][2]*y + icam[2][2]*z; - rout[i] = CLIP(rr); - gout[i] = CLIP(gr); - bout[i] = CLIP(br); - } - else { - rout[i] = rin[i]; - gout[i] = gin[i]; - bout[i] = bin[i]; - } - } -} - -void RawImageSource::hlRecovery (std::string method, unsigned short* red, unsigned short* green, unsigned short* blue, int i, int sx1, int width, int skip) { - - if (method=="Luminance") - HLRecovery_Luminance (red, green, blue, red, green, blue, width, 65535 / ri->defgain); - else if (method=="CIELab blending") - HLRecovery_CIELab (red, green, blue, red, green, blue, width, 65535 / ri->defgain, cam, icam); - else if (method=="Color") - HLRecovery_ColorPropagation (red, green, blue, i, sx1, width, skip); -} - -int RawImageSource::getAEHistogram (unsigned int* histogram, int& histcompr) { - - histcompr = 3; - - memset (histogram, 0, (65536>>histcompr)*sizeof(int)); - - for (int i=border; iheight-border; i++) { - int start, end; - if (fuji) { - int fw = ri->fuji_width; - start = ABS(fw-i) + border; - end = MIN( ri->height+ ri->width-fw-i, fw+i) - border; - } - else { - start = border; - end = ri->width-border; - } - if (ri->filters) - for (int j=start; jdata[i][j]>>histcompr]+=2; - else - histogram[ri->data[i][j]>>histcompr]+=4; - else - for (int j=start; j<3*end; j++) { - histogram[ri->data[i][j+0]>>histcompr]++; - histogram[ri->data[i][j+1]>>histcompr]++; - histogram[ri->data[i][j+2]>>histcompr]++; - } - } - return 1; -} - -ColorTemp RawImageSource::getAutoWB () { - - double avg_r = 0; - double avg_g = 0; - double avg_b = 0; - int rn = 0, gn = 0, bn = 0; - - if (fuji) { - for (int i=32; iheight-32; i++) { - int fw = ri->fuji_width; - int start = ABS(fw-i) + 32; - int end = MIN(ri->height+ri->width-fw-i, fw+i) - 32; - for (int j=start; jfilters) { - double d = CLIP(ri->defgain*ri->data[i][3*j]); - if (d>64000) - continue; - avg_r += d*d*d*d*d*d; rn++; - d = CLIP(ri->defgain*ri->data[i][3*j+1]); - if (d>64000) - continue; - avg_g += d*d*d*d*d*d; gn++; - d = CLIP(ri->defgain*ri->data[i][3*j+2]); - if (d>64000) - continue; - avg_b += d*d*d*d*d*d; bn++; - } - else { - double d = CLIP(ri->defgain*ri->data[i][j]); - if (d>64000) - continue; - double dp = d*d*d*d*d*d; - if (ISRED(ri,i,j)) { - avg_r += dp; - rn++; - } - else if (ISGREEN(ri,i,j)) { - avg_g += dp; - gn++; - } - else if (ISBLUE(ri,i,j)) { - avg_b += dp; - bn++; - } - } - } - } - } - else { - for (int i=32; iheight-32; i++) - for (int j=32; jwidth-32; j++) { - if (!ri->filters) { - double d = CLIP(ri->defgain*ri->data[i][3*j]); - if (d>64000) - continue; - avg_r += d*d*d*d*d*d; rn++; - d = CLIP(ri->defgain*ri->data[i][3*j+1]); - if (d>64000) - continue; - avg_g += d*d*d*d*d*d; gn++; - d = CLIP(ri->defgain*ri->data[i][3*j+2]); - if (d>64000) - continue; - avg_b += d*d*d*d*d*d; bn++; - } - else { - double d = CLIP(ri->defgain*ri->data[i][j]); - if (d>64000) - continue; - double dp = d*d*d*d*d*d; - if (ISRED(ri,i,j)) { - avg_r += dp; - rn++; - } - else if (ISGREEN(ri,i,j)) { - avg_g += dp; - gn++; - } - else if (ISBLUE(ri,i,j)) { - avg_b += dp; - bn++; - } - } - } - } - - printf ("AVG: %g %g %g\n", avg_r/rn, avg_g/gn, avg_b/bn); - -// double img_r, img_g, img_b; -// wb.getMultipliers (img_r, img_g, img_b); - -// return ColorTemp (pow(avg_r/rn, 1.0/6.0)*img_r, pow(avg_g/gn, 1.0/6.0)*img_g, pow(avg_b/bn, 1.0/6.0)*img_b); - - double reds = pow (avg_r/rn, 1.0/6.0) * ri->camwb_red; - double greens = pow (avg_g/gn, 1.0/6.0) * ri->camwb_green; - double blues = pow (avg_b/bn, 1.0/6.0) * ri->camwb_blue; - - double rm = coeff[0][0]*reds + coeff[0][1]*greens + coeff[0][2]*blues; - double gm = coeff[1][0]*reds + coeff[1][1]*greens + coeff[1][2]*blues; - double bm = coeff[2][0]*reds + coeff[2][1]*greens + coeff[2][2]*blues; - - return ColorTemp (rm, gm, bm); -} - -void RawImageSource::transformPosition (int x, int y, int tran, int& ttx, int& tty) { - - tran = defTransform (tran); - - x += border; - y += border; - - if (d1x) { - if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) - x /= 2; - else - y /= 2; - } - - int w = W, h = H; - if (fuji) { - w = ri->fuji_width * 2 + 1; - h = (H - ri->fuji_width)*2 + 1; - } - int sw = w, sh = h; - if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) { - sw = h; - sh = w; - } - - int ppx = x, ppy = y; - if (tran & TR_HFLIP) - ppx = sw - 1 - x ; - if (tran & TR_VFLIP) - ppy = sh - 1 - y; - - int tx = ppx; - int ty = ppy; - - if ((tran & TR_ROT) == TR_R180) { - tx = w - 1 - ppx; - ty = h - 1 - ppy; - } - else if ((tran & TR_ROT) == TR_R90) { - tx = ppy; - ty = h - 1 - ppx; - } - else if ((tran & TR_ROT) == TR_R270) { - tx = w - 1 - ppy; - ty = ppx; - } - - if (fuji) { - ttx = (tx+ty) / 2; - tty = (ty-tx) / 2 + ri->fuji_width; - } - else { - ttx = tx; - tty = ty; - } -} - -ColorTemp RawImageSource::getSpotWB (std::vector red, std::vector green, std::vector& blue, int tran) { - - int x; int y; - int d[9][2] = {{0,0}, {-1,-1}, {-1,0}, {-1,1}, {0,-1}, {0,1}, {1,-1}, {1,0}, {1,1}}; - double reds = 0, greens = 0, blues = 0; - int rn = 0, gn = 0, bn = 0; - - if (!ri->filters) { - for (int i=0; i=0 && y>=0 && xdata[y][3*x]; - rn++; - } - transformPosition (green[i].x, green[i].y, tran, x, y); - if (x>=0 && y>=0 && xdata[y][3*x+1]; - gn++; - } - transformPosition (blue[i].x, blue[i].y, tran, x, y); - if (x>=0 && y>=0 && xdata[y][3*x+2]; - bn++; - } - } - } - else { - for (int i=0; i=0 && yv>=0 && xvdata[yv][xv]; - rn++; - break; - } - } - transformPosition (green[i].x, green[i].y, tran, x, y); - for (int k=0; k<9; k++) { - int xv = x + d[k][0]; - int yv = y + d[k][1]; - if (ISGREEN(ri,yv,xv) && xv>=0 && yv>=0 && xvdata[yv][xv]; - gn++; - break; - } - } - transformPosition (blue[i].x, blue[i].y, tran, x, y); - for (int k=0; k<9; k++) { - int xv = x + d[k][0]; - int yv = y + d[k][1]; - if (ISBLUE(ri,yv,xv) && xv>=0 && yv>=0 && xvdata[yv][xv]; - bn++; - break; - } - } - } - } - - reds = reds/rn * ri->camwb_red; - greens = greens/gn * ri->camwb_green; - blues = blues/bn * ri->camwb_blue; - - double rm = coeff[0][0]*reds + coeff[0][1]*greens + coeff[0][2]*blues; - double gm = coeff[1][0]*reds + coeff[1][1]*greens + coeff[1][2]*blues; - double bm = coeff[2][0]*reds + coeff[2][1]*greens + coeff[2][2]*blues; - - return ColorTemp (rm, gm, bm); -} - -#define FORCC for (c=0; c < colors; c++) -#define fc(row,col) \ - (ri->prefilters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3) -typedef unsigned short ushort; -void RawImageSource::vng4_demosaic () { - - static const signed char *cp, terms[] = { - -2,-2,+0,-1,0,0x01, -2,-2,+0,+0,1,0x01, -2,-1,-1,+0,0,0x01, - -2,-1,+0,-1,0,0x02, -2,-1,+0,+0,0,0x03, -2,-1,+0,+1,1,0x01, - -2,+0,+0,-1,0,0x06, -2,+0,+0,+0,1,0x02, -2,+0,+0,+1,0,0x03, - -2,+1,-1,+0,0,0x04, -2,+1,+0,-1,1,0x04, -2,+1,+0,+0,0,0x06, - -2,+1,+0,+1,0,0x02, -2,+2,+0,+0,1,0x04, -2,+2,+0,+1,0,0x04, - -1,-2,-1,+0,0,0x80, -1,-2,+0,-1,0,0x01, -1,-2,+1,-1,0,0x01, - -1,-2,+1,+0,1,0x01, -1,-1,-1,+1,0,0x88, -1,-1,+1,-2,0,0x40, - -1,-1,+1,-1,0,0x22, -1,-1,+1,+0,0,0x33, -1,-1,+1,+1,1,0x11, - -1,+0,-1,+2,0,0x08, -1,+0,+0,-1,0,0x44, -1,+0,+0,+1,0,0x11, - -1,+0,+1,-2,1,0x40, -1,+0,+1,-1,0,0x66, -1,+0,+1,+0,1,0x22, - -1,+0,+1,+1,0,0x33, -1,+0,+1,+2,1,0x10, -1,+1,+1,-1,1,0x44, - -1,+1,+1,+0,0,0x66, -1,+1,+1,+1,0,0x22, -1,+1,+1,+2,0,0x10, - -1,+2,+0,+1,0,0x04, -1,+2,+1,+0,1,0x04, -1,+2,+1,+1,0,0x04, - +0,-2,+0,+0,1,0x80, +0,-1,+0,+1,1,0x88, +0,-1,+1,-2,0,0x40, - +0,-1,+1,+0,0,0x11, +0,-1,+2,-2,0,0x40, +0,-1,+2,-1,0,0x20, - +0,-1,+2,+0,0,0x30, +0,-1,+2,+1,1,0x10, +0,+0,+0,+2,1,0x08, - +0,+0,+2,-2,1,0x40, +0,+0,+2,-1,0,0x60, +0,+0,+2,+0,1,0x20, - +0,+0,+2,+1,0,0x30, +0,+0,+2,+2,1,0x10, +0,+1,+1,+0,0,0x44, - +0,+1,+1,+2,0,0x10, +0,+1,+2,-1,1,0x40, +0,+1,+2,+0,0,0x60, - +0,+1,+2,+1,0,0x20, +0,+1,+2,+2,0,0x10, +1,-2,+1,+0,0,0x80, - +1,-1,+1,+1,0,0x88, +1,+0,+1,+2,0,0x08, +1,+0,+2,-1,0,0x40, - +1,+0,+2,+1,0,0x10 - }, chood[] = { -1,-1, -1,0, -1,+1, 0,+1, +1,+1, +1,0, +1,-1, 0,-1 }; - - if (plistener) { - plistener->setProgressStr ("Demosaicing..."); - plistener->setProgress (0.0); - } - - ushort (*brow[5])[4], *pix; - int prow=7, pcol=1, *ip, *code[16][16], gval[8], gmin, gmax, sum[4]; - int row, col, x, y, x1, x2, y1, y2, t, weight, grads, color, diag; - int g, diff, thold, num, c, width=W, height=H, colors=4; - ushort (*image)[4]; - int lcode[16][16][32], shift, i, j; - - image = (ushort (*)[4]) calloc (H*W, sizeof *image); - for (int ii=0; iidata[ii][jj]; - -// first linear interpolation - for (row=0; row < 16; row++) - for (col=0; col < 16; col++) { - ip = lcode[row][col]; - memset (sum, 0, sizeof sum); - for (y=-1; y <= 1; y++) - for (x=-1; x <= 1; x++) { - shift = (y==0) + (x==0); - if (shift == 2) continue; - color = fc(row+y,col+x); - *ip++ = (width*y + x)*4 + color; - *ip++ = shift; - *ip++ = color; - sum[color] += 1 << shift; - } - FORCC - if (c != fc(row,col)) { - *ip++ = c; - *ip++ = 256 / sum[c]; - } - } - - for (row=1; row < height-1; row++) - for (col=1; col < width-1; col++) { - pix = image[row*width+col]; - ip = lcode[row & 15][col & 15]; - memset (sum, 0, sizeof sum); - for (i=8; i--; ip+=3) - sum[ip[2]] += pix[ip[0]] << ip[1]; - for (i=colors; --i; ip+=2) - pix[ip[0]] = sum[ip[0]] * ip[1] >> 8; - } - -// lin_interpolate(); - - - ip = (int *) calloc ((prow+1)*(pcol+1), 1280); - for (row=0; row <= prow; row++) /* Precalculate for VNG */ - for (col=0; col <= pcol; col++) { - code[row][col] = ip; - for (cp=terms, t=0; t < 64; t++) { - y1 = *cp++; x1 = *cp++; - y2 = *cp++; x2 = *cp++; - weight = *cp++; - grads = *cp++; - color = fc(row+y1,col+x1); - if (fc(row+y2,col+x2) != color) continue; - diag = (fc(row,col+1) == color && fc(row+1,col) == color) ? 2:1; - if (abs(y1-y2) == diag && abs(x1-x2) == diag) continue; - *ip++ = (y1*width + x1)*4 + color; - *ip++ = (y2*width + x2)*4 + color; - *ip++ = weight; - for (g=0; g < 8; g++) - if (grads & 1< gval[g]) gmin = gval[g]; - if (gmax < gval[g]) gmax = gval[g]; - } - if (gmax == 0) { - memcpy (brow[2][col], pix, sizeof *image); - continue; - } - thold = gmin + (gmax >> 1); - memset (sum, 0, sizeof sum); - for (num=g=0; g < 8; g++,ip+=2) { /* Average the neighbors */ - if (gval[g] <= thold) { - FORCC - if (c == color && ip[1]) - sum[c] += (pix[c] + pix[ip[1]]) >> 1; - else - sum[c] += pix[ip[0] + c]; - num++; - } - } - FORCC { /* Save to buffer */ - t = pix[color]; - if (c != color) - t += (sum[c] - sum[color]) / num; - brow[2][col][c] = CLIP(t); - } - } - if (row > 3) /* Write buffer to image */ - memcpy (image[(row-2)*width+2], brow[0]+2, (width-4)*sizeof *image); - for (g=0; g < 4; g++) - brow[(g-1) & 3] = brow[g]; - if (!(row%20) && plistener) - plistener->setProgress ((double)row / (H-2)); - } - memcpy (image[(row-2)*width+2], brow[0]+2, (width-4)*sizeof *image); - memcpy (image[(row-1)*width+2], brow[1]+2, (width-4)*sizeof *image); - free (brow[4]); - free (code[0][0]); - - green = new unsigned short*[H]; - for (int i=0; i> 1; - } - free (image); -} - -//#define ABS(x) (((int)(x) ^ ((int)(x) >> 31)) - ((int)(x) >> 31)) -//#define MIN(a,b) ((a) < (b) ? (a) : (b)) -//#define MAX(a,b) ((a) > (b) ? (a) : (b)) -//#define LIM(x,min,max) MAX(min,MIN(x,max)) -//#define CLIP(x) LIM(x,0,65535) -#undef fc -#define fc(row,col) \ - (ri->filters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3) -#define FC(x,y) fc(x,y) -#define LIM(x,min,max) MAX(min,MIN(x,max)) -#define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y)) - -/* - Patterned Pixel Grouping Interpolation by Alain Desbiolles -*/ -void RawImageSource::ppg_demosaic() -{ - int width=W, height=H; - int dir[5] = { 1, width, -1, -width, 1 }; - int row, col, diff[2], guess[2], c, d, i; - ushort (*pix)[4]; - - ushort (*image)[4]; - int colors = 3; - - if (plistener) { - plistener->setProgressStr ("Demosaicing..."); - plistener->setProgress (0.0); - } - - image = (ushort (*)[4]) calloc (H*W, sizeof *image); - for (int ii=0; iidata[ii][jj]; - - border_interpolate(3, image); - -/* Fill in the green layer with gradients and pattern recognition: */ - for (row=3; row < height-3; row++) { - for (col=3+(FC(row,3) & 1), c=FC(row,col); col < width-3; col+=2) { - pix = image + row*width+col; - for (i=0; (d=dir[i]) > 0; i++) { - guess[i] = (pix[-d][1] + pix[0][c] + pix[d][1]) * 2 - - pix[-2*d][c] - pix[2*d][c]; - diff[i] = ( ABS(pix[-2*d][c] - pix[ 0][c]) + - ABS(pix[ 2*d][c] - pix[ 0][c]) + - ABS(pix[ -d][1] - pix[ d][1]) ) * 3 + - ( ABS(pix[ 3*d][1] - pix[ d][1]) + - ABS(pix[-3*d][1] - pix[-d][1]) ) * 2; - } - d = dir[i = diff[0] > diff[1]]; - pix[0][1] = ULIM(guess[i] >> 2, pix[d][1], pix[-d][1]); - } - if(plistener) plistener->setProgress(0.33*row/(height-3)); - } -/* Calculate red and blue for each green pixel: */ - for (row=1; row < height-1; row++) { - for (col=1+(FC(row,2) & 1), c=FC(row,col+1); col < width-1; col+=2) { - pix = image + row*width+col; - for (i=0; (d=dir[i]) > 0; c=2-c, i++) - pix[0][c] = CLIP((pix[-d][c] + pix[d][c] + 2*pix[0][1] - - pix[-d][1] - pix[d][1]) >> 1); - } - if(plistener) plistener->setProgress(0.33 + 0.33*row/(height-1)); - } -/* Calculate blue for red pixels and vice versa: */ - for (row=1; row < height-1; row++) { - for (col=1+(FC(row,1) & 1), c=2-FC(row,col); col < width-1; col+=2) { - pix = image + row*width+col; - for (i=0; (d=dir[i]+dir[i+1]) > 0; i++) { - diff[i] = ABS(pix[-d][c] - pix[d][c]) + - ABS(pix[-d][1] - pix[0][1]) + - ABS(pix[ d][1] - pix[0][1]); - guess[i] = pix[-d][c] + pix[d][c] + 2*pix[0][1] - - pix[-d][1] - pix[d][1]; - } - if (diff[0] != diff[1]) - pix[0][c] = CLIP(guess[diff[0] > diff[1]] >> 1); - else - pix[0][c] = CLIP((guess[0]+guess[1]) >> 2); - } - if(plistener) plistener->setProgress(0.67 + 0.33*row/(height-1)); - } - - red = new unsigned short*[H]; - for (int i=0; i= border && row < height-border) - col = width-border; - memset (sum, 0, sizeof sum); - for (y=row-1; y != row+2; y++) - for (x=col-1; x != col+2; x++) - if (y < height && x < width) { - f = fc(y,x); - sum[f] += image[y*width+x][f]; - sum[f+4]++; - } - f = fc(row,col); - FORCC if (c != f && sum[c+4]) - image[row*width+col][c] = sum[c] / sum[c+4]; - } -} - -void RawImageSource::bilinear_interpolate_block(ushort (*image)[4], int start, int end) -{ - ushort (*pix); - int i, *ip, sum[4]; - int width=W; - int colors = 3; - - for (int row = start; row < end; row++) - for (int col=1; col < width-1; col++) { - pix = image[row*width+col]; - ip = blcode[row & 15][col & 15]; - memset (sum, 0, sizeof sum); - for (i=8; i--; ip+=3) - sum[ip[2]] += pix[ip[0]] << ip[1]; - for (i=colors; --i; ip+=2) - pix[ip[0]] = sum[ip[0]] * ip[1] >> 8; - } - - -} - -void RawImageSource::bilinear_demosaic() -{ - int width=W, height=H; - int *ip, sum[4]; - int c, x, y, row, col, shift, color; - int colors = 3; - - ushort (*image)[4], *pix; - image = (ushort (*)[4]) calloc (H*W, sizeof *image); - - for (int ii=0; iidata[ii][jj]; - - //if (verbose) fprintf (stderr,_("Bilinear interpolation...\n")); - if (plistener) { - plistener->setProgressStr ("Demosaicing..."); - plistener->setProgress (0.0); - } - - memset(blcode,0,16*16*32); - for (row=0; row < 16; row++) - for (col=0; col < 16; col++) { - ip = blcode[row][col]; - memset (sum, 0, sizeof sum); - for (y=-1; y <= 1; y++) - for (x=-1; x <= 1; x++) { - shift = (y==0) + (x==0); - if (shift == 2) continue; - color = fc(row+y,col+x); - *ip++ = (width*y + x)*4 + color; - *ip++ = shift; - *ip++ = color; - sum[color] += 1 << shift; - } - FORCC - if (c != fc(row,col)) { - *ip++ = c; - *ip++ = 256 / sum[c]; - } - } - -#ifdef _OPENMP - #pragma omp parallel - { - int tid = omp_get_thread_num(); - int nthreads = omp_get_num_threads(); - int blk = H/nthreads; - - int start = 0; - if (tid == 0) start = 1; - if (tidsetProgress (1.0); - free (image); -} - -/* - Adaptive Homogeneity-Directed interpolation is based on - the work of Keigo Hirakawa, Thomas Parks, and Paul Lee. - */ -#define TS 256 /* Tile Size */ -#define FORC(cnt) for (c=0; c < cnt; c++) -#define FORC3 FORC(3) -#define SQR(x) ((x)*(x)) - -void RawImageSource::ahd_demosaic() -{ - int i, j, k, top, left, row, col, tr, tc, c, d, val, hm[2]; - ushort (*pix)[4], (*rix)[3]; - static const int dir[4] = { -1, 1, -TS, TS }; - unsigned ldiff[2][4], abdiff[2][4], leps, abeps; - float r, cbrt[0x10000], xyz[3], xyz_cam[3][4]; - ushort (*rgb)[TS][TS][3]; - short (*lab)[TS][TS][3], (*lix)[3]; - char (*homo)[TS][TS], *buffer; - - int width=W, height=H; - ushort (*image)[4]; - int colors = 3; - - const double xyz_rgb[3][3] = { /* XYZ from RGB */ - { 0.412453, 0.357580, 0.180423 }, - { 0.212671, 0.715160, 0.072169 }, - { 0.019334, 0.119193, 0.950227 } - }; - - const float d65_white[3] = { 0.950456, 1, 1.088754 }; - - if (plistener) { - plistener->setProgressStr ("Demosaicing..."); - plistener->setProgress (0.0); - } - - image = (ushort (*)[4]) calloc (H*W, sizeof *image); - for (int ii=0; iidata[ii][jj]; - - for (i=0; i < 0x10000; i++) { - r = i / 65535.0; - cbrt[i] = r > 0.008856 ? pow(r,1/3.0) : 7.787*r + 16/116.0; - } - - for (i=0; i < 3; i++) - for (j=0; j < colors; j++) - for (xyz_cam[i][j] = k=0; k < 3; k++) - xyz_cam[i][j] += xyz_rgb[i][k] * coeff[k][j] / d65_white[i]; - - border_interpolate(5, image); - buffer = (char *) malloc (26*TS*TS); /* 1664 kB */ - //merror (buffer, "ahd_interpolate()"); - rgb = (ushort(*)[TS][TS][3]) buffer; - lab = (short (*)[TS][TS][3])(buffer + 12*TS*TS); - homo = (char (*)[TS][TS]) (buffer + 24*TS*TS); - - // helper variables for progress indication - int n_tiles = ((height-7 + (TS-7))/(TS-6)) * ((width-7 + (TS-7))/(TS-6)); - int tile = 0; - - for (top=2; top < height-5; top += TS-6) - for (left=2; left < width-5; left += TS-6) { - - /* Interpolate green horizontally and vertically: */ - for (row = top; row < top+TS && row < height-2; row++) { - col = left + (FC(row,left) & 1); - for (c = FC(row,col); col < left+TS && col < width-2; col+=2) { - pix = image + row*width+col; - val = ((pix[-1][1] + pix[0][c] + pix[1][1]) * 2 - - pix[-2][c] - pix[2][c]) >> 2; - rgb[0][row-top][col-left][1] = ULIM(val,pix[-1][1],pix[1][1]); - val = ((pix[-width][1] + pix[0][c] + pix[width][1]) * 2 - - pix[-2*width][c] - pix[2*width][c]) >> 2; - rgb[1][row-top][col-left][1] = ULIM(val,pix[-width][1],pix[width][1]); - } - } - - /* Interpolate red and blue, and convert to CIELab: */ - for (d=0; d < 2; d++) - for (row=top+1; row < top+TS-1 && row < height-3; row++) - for (col=left+1; col < left+TS-1 && col < width-3; col++) { - pix = image + row*width+col; - rix = &rgb[d][row-top][col-left]; - lix = &lab[d][row-top][col-left]; - if ((c = 2 - FC(row,col)) == 1) { - c = FC(row+1,col); - val = pix[0][1] + (( pix[-1][2-c] + pix[1][2-c] - - rix[-1][1] - rix[1][1] ) >> 1); - rix[0][2-c] = CLIP(val); - val = pix[0][1] + (( pix[-width][c] + pix[width][c] - - rix[-TS][1] - rix[TS][1] ) >> 1); - } else - val = rix[0][1] + (( pix[-width-1][c] + pix[-width+1][c] - + pix[+width-1][c] + pix[+width+1][c] - - rix[-TS-1][1] - rix[-TS+1][1] - - rix[+TS-1][1] - rix[+TS+1][1] + 1) >> 2); - rix[0][c] = CLIP(val); - c = FC(row,col); - rix[0][c] = pix[0][c]; - xyz[0] = xyz[1] = xyz[2] = 0.5; - FORCC { - xyz[0] += xyz_cam[0][c] * rix[0][c]; - xyz[1] += xyz_cam[1][c] * rix[0][c]; - xyz[2] += xyz_cam[2][c] * rix[0][c]; - } - xyz[0] = cbrt[CLIP((int) xyz[0])]; - xyz[1] = cbrt[CLIP((int) xyz[1])]; - xyz[2] = cbrt[CLIP((int) xyz[2])]; - lix[0][0] = 64 * (116 * xyz[1] - 16); - lix[0][1] = 64 * 500 * (xyz[0] - xyz[1]); - lix[0][2] = 64 * 200 * (xyz[1] - xyz[2]); - } - - /* Build homogeneity maps from the CIELab images: */ - memset (homo, 0, 2*TS*TS); - for (row=top+2; row < top+TS-2 && row < height-4; row++) { - tr = row-top; - for (col=left+2; col < left+TS-2 && col < width-4; col++) { - tc = col-left; - for (d=0; d < 2; d++) { - lix = &lab[d][tr][tc]; - for (i=0; i < 4; i++) { - ldiff[d][i] = ABS(lix[0][0]-lix[dir[i]][0]); - abdiff[d][i] = SQR(lix[0][1]-lix[dir[i]][1]) - + SQR(lix[0][2]-lix[dir[i]][2]); - } - } - leps = MIN(MAX(ldiff[0][0],ldiff[0][1]), - MAX(ldiff[1][2],ldiff[1][3])); - abeps = MIN(MAX(abdiff[0][0],abdiff[0][1]), - MAX(abdiff[1][2],abdiff[1][3])); - for (d=0; d < 2; d++) - for (i=0; i < 4; i++) - if (ldiff[d][i] <= leps && abdiff[d][i] <= abeps) - homo[d][tr][tc]++; - } - } - - /* Combine the most homogenous pixels for the final result: */ - for (row=top+3; row < top+TS-3 && row < height-5; row++) { - tr = row-top; - for (col=left+3; col < left+TS-3 && col < width-5; col++) { - tc = col-left; - for (d=0; d < 2; d++) - for (hm[d]=0, i=tr-1; i <= tr+1; i++) - for (j=tc-1; j <= tc+1; j++) - hm[d] += homo[d][i][j]; - if (hm[0] != hm[1]) - FORC3 image[row*width+col][c] = rgb[hm[1] > hm[0]][tr][tc][c]; - else - FORC3 image[row*width+col][c] = - (rgb[0][tr][tc][c] + rgb[1][tr][tc][c]) >> 1; - } - } - - tile++; - if(plistener) { - plistener->setProgress((double)tile / n_tiles); - } - } - - if(plistener) plistener->setProgress (1.0); - free (buffer); - red = new unsigned short*[H]; - for (int i=0; i= H-border) rowMax = TILEBORDER+H-border-y0; - if( x0+TILESIZE+TILEBORDER >= W-border) colMax = TILEBORDER+W-border-x0; -} - -void RawImageSource::fill_raw( ushort (*cache )[4], int x0, int y0, ushort** rawData) -{ - int rowMin,colMin,rowMax,colMax; - dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,0); - - for (int row=rowMin,y=y0-TILEBORDER+rowMin; row= border && col < W - border && row >= border && row < H - border){ - col = W - border; - if(col >= x0+TILESIZE+TILEBORDER ) - break; - } - memset(sum, 0, sizeof sum); - for (y = row - 1; y != row + 2; y++) - for (x = col - 1; x != col + 2; x++) - if (y < H && y< y0+TILESIZE+TILEBORDER && x < W && x ( pix[-4] + pix[+4] + pix[-u] + pix[+u])/4 ) - image[indx][3] = ((MIN( pix[-4], pix[+4]) + pix[-4] + pix[+4] ) < (MIN( pix[-u], pix[+u]) + pix[-u] + pix[+u])); - else - image[indx][3] = ((MAX( pix[-4], pix[+4]) + pix[-4] + pix[+4] ) > (MAX( pix[-u], pix[+u]) + pix[-u] + pix[+u])); - } - } -} - - -// interpolated green pixels are corrected using the map -void RawImageSource::dcb_correction(ushort (*image)[4], int x0, int y0) -{ - const int u=CACHESIZE, v=2*CACHESIZE; - int rowMin,colMin,rowMax,colMax; - dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,2); - - for (int row=rowMin; row < rowMax; row++) { - for (int col = colMin+(FC(y0-TILEBORDER+row,x0-TILEBORDER+colMin)&1),indx=row*CACHESIZE+col,c=FC(y0-TILEBORDER+row,x0-TILEBORDER+col); col < colMax; col+=2, indx+=2) { - int current = 4*image[indx][3] + - 2*(image[indx+u][3] + image[indx-u][3] + image[indx+1][3] + image[indx-1][3]) + - image[indx+v][3] + image[indx-v][3] + image[indx+2][3] + image[indx-2][3]; - image[indx][1] = ((16-current)*(image[indx-1][1] + image[indx+1][1])/2 + current*(image[indx-u][1] + image[indx+u][1])/2)/16; - } - } -} - -// R and B smoothing using green contrast, all pixels except 2 pixel wide border -void RawImageSource::dcb_pp(ushort (*image)[4], int x0, int y0) -{ - const int u=CACHESIZE; - int rowMin,colMin,rowMax,colMax; - dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,2); - - for (int row=rowMin; row < rowMax; row++) - for (int col=colMin, indx=row*CACHESIZE+col; col < colMax; col++, indx++) { - //int 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; - //int 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; - //int 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; - ushort (*pix)[4] = image+(indx-u-1); - int r1 = (*pix)[0]; - int g1 = (*pix)[1]; - int b1 = (*pix)[2]; - pix++; - r1 += (*pix)[0]; - g1 += (*pix)[1]; - b1 += (*pix)[2]; - pix++; - r1 += (*pix)[0]; - g1 += (*pix)[1]; - b1 += (*pix)[2]; - pix+=CACHESIZE-2; - r1 += (*pix)[0]; - g1 += (*pix)[1]; - b1 += (*pix)[2]; - pix+=2; - r1 += (*pix)[0]; - g1 += (*pix)[1]; - b1 += (*pix)[2]; - pix+=CACHESIZE-2; - r1 += (*pix)[0]; - g1 += (*pix)[1]; - b1 += (*pix)[2]; - pix++; - r1 += (*pix)[0]; - g1 += (*pix)[1]; - b1 += (*pix)[2]; - pix++; - r1 += (*pix)[0]; - g1 += (*pix)[1]; - b1 += (*pix)[2]; - r1 /=8; - g1 /=8; - b1 /=8; - r1 = r1 + ( image[indx][1] - g1 ); - b1 = b1 + ( image[indx][1] - g1 ); - image[indx][0] = CLIP(r1); - image[indx][2] = CLIP(b1); - } -} - -// interpolated green pixels are corrected using the map -// with correction -void RawImageSource::dcb_correction2(ushort (*image)[4], int x0, int y0) -{ - const int u=CACHESIZE, v=2*CACHESIZE; - int rowMin,colMin,rowMax,colMax; - dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,4); - - for (int row=rowMin; row < rowMax; row++) { - for (int col = colMin+(FC(y0-TILEBORDER+row,x0-TILEBORDER+colMin)&1),indx=row*CACHESIZE+col,c=FC(y0-TILEBORDER+row,x0-TILEBORDER+col); col < colMax; col+=2, indx+=2) { - int current = 4*image[indx][3] + - 2*(image[indx+u][3] + image[indx-u][3] + image[indx+1][3] + image[indx-1][3]) + - image[indx+v][3] + image[indx-v][3] + image[indx+2][3] + image[indx-2][3]; - current = ((16-current)*((image[indx-1][1] + image[indx+1][1])/2 + image[indx][c] - (image[indx+2][c] + image[indx-2][c])/2) + current*((image[indx-u][1] + image[indx+u][1])/2 + image[indx][c] - (image[indx+v][c] + image[indx-v][c])/2))/16; - image[indx][1] = CLIP(current); - } - } -} - -// image refinement -void RawImageSource::dcb_refinement(ushort (*image)[4], int x0, int y0) -{ - const int u=CACHESIZE, v=2*CACHESIZE, w=3*CACHESIZE; - int rowMin,colMin,rowMax,colMax; - dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,4); - - float f[5],g1,g2; - - for (int row=rowMin; row < rowMax; row++) - for (int col=colMin+(FC(y0-TILEBORDER+row,x0-TILEBORDER+colMin)&1),indx=row*CACHESIZE+col,c=FC(y0-TILEBORDER+row,x0-TILEBORDER+col); col < colMax; col+=2,indx+=2){ - int current = 4*image[indx][3] + - 2*(image[indx+u][3] + image[indx-u][3] + image[indx+1][3] + image[indx-1][3]) - +image[indx+v][3] + image[indx-v][3] + image[indx-2][3] + image[indx+2][3]; - - f[0] = (float)(image[indx-u][1] + image[indx+u][1])/(2 + 2*image[indx][c]); - f[1] = 2*(float)image[indx-u][1]/(2 + image[indx-v][c] + image[indx][c]); - f[2] = (float)(image[indx-u][1] + image[indx-w][1])/(2 + 2*image[indx-v][c]); - f[3] = 2*(float)image[indx+u][1]/(2 + image[indx+v][c] + image[indx][c]); - f[4] = (float)(image[indx+u][1] + image[indx+w][1])/(2 + 2*image[indx+v][c]); - - g1 = (f[0] + f[1] + f[2] + f[3] + f[4] - MAX(f[1], MAX(f[2], MAX(f[3], f[4]))) - MIN(f[1], MIN(f[2], MIN(f[3], f[4]))))/3.0; - - f[0] = (float)(image[indx-1][1] + image[indx+1][1])/(2 + 2*image[indx][c]); - f[1] = 2*(float)image[indx-1][1]/(2 + image[indx-2][c] + image[indx][c]); - f[2] = (float)(image[indx-1][1] + image[indx-3][1])/(2 + 2*image[indx-2][c]); - f[3] = 2*(float)image[indx+1][1]/(2 + image[indx+2][c] + image[indx][c]); - f[4] = (float)(image[indx+1][1] + image[indx+3][1])/(2 + 2*image[indx+2][c]); - - g2 = (f[0] + f[1] + f[2] + f[3] + f[4] - MAX(f[1], MAX(f[2], MAX(f[3], f[4]))) - MIN(f[1], MIN(f[2], MIN(f[3], f[4]))))/3.0; - - image[indx][1] = CLIP((2+image[indx][c])*(current*g1 + (16-current)*g2)/16.0); - -// get rid of the overshooted pixels - int 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]))))))); - int 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]))))))); - - image[indx][1] = LIM(image[indx][1], min, max); - - - } -} - -// missing colors are interpolated using high quality algorithm by Luis Sanz Rodríguez -void RawImageSource::dcb_color_full(ushort (*image)[4], int x0, int y0, float (*chroma)[2]) -{ - const int u=CACHESIZE, v=2*CACHESIZE, w=3*CACHESIZE; - int rowMin,colMin,rowMax,colMax; - dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,3); - - int i,j; - float f[4],g[4]; - - for (int row=1; row < CACHESIZE-1; row++) - for (int col=1+(FC(y0-TILEBORDER+row,x0-TILEBORDER+1)&1),indx=row*CACHESIZE+col,c=FC(y0-TILEBORDER+row,x0-TILEBORDER+col),d=c/2; col < CACHESIZE-1; col+=2,indx+=2) - chroma[indx][d]=image[indx][c]-image[indx][1]; - - for (int row=rowMin; rowsetProgressStr ("DCB Demosaicing..."); - plistener->setProgress (currentProgress); - } - - int wTiles = W/TILESIZE + (W%TILESIZE?1:0); - int hTiles = H/TILESIZE + (H%TILESIZE?1:0); - int numTiles = wTiles * hTiles; - int tilesDone=0; -#ifdef _OPENMP - int nthreads = omp_get_max_threads(); - ushort (**image)[4] = (ushort(**)[4]) calloc( nthreads,sizeof( void*) ); - ushort (**image2)[3] = (ushort(**)[3]) calloc( nthreads,sizeof( void*) ); - ushort (**image3)[3] = (ushort(**)[3]) calloc( nthreads,sizeof( void*) ); - float (**chroma)[2] = (float (**)[2]) calloc( nthreads,sizeof( void*) ); - for(int i=0; idata ); - if( !xTile || !yTile || xTile==wTiles-1 || yTile==hTiles-1) - fill_border(tile,6, x0, y0); - dcb_hid(tile,buffer,buffer2,x0,y0); - copy_to_buffer(buffer, tile); - for (int i=iterations; i>0;i--) { - dcb_hid2(tile,x0,y0); - dcb_hid2(tile,x0,y0); - dcb_hid2(tile,x0,y0); - dcb_map(tile,x0,y0); - dcb_correction(tile,x0,y0); - } - dcb_color(tile,x0,y0); - dcb_pp(tile,x0,y0); - dcb_map(tile,x0,y0); - dcb_correction2(tile,x0,y0); - dcb_map(tile,x0,y0); - dcb_correction(tile,x0,y0); - dcb_color(tile,x0,y0); - dcb_map(tile,x0,y0); - dcb_correction(tile,x0,y0); - dcb_map(tile,x0,y0); - dcb_correction(tile,x0,y0); - dcb_map(tile,x0,y0); - restore_from_buffer(tile, buffer); - dcb_color(tile,x0,y0); - if (dcb_enhance) { - dcb_refinement(tile,x0,y0); - dcb_color_full(tile,x0,y0,chrm); - } - - for(int y=0;y currentProgress){ - currentProgress+=0.1; // Show progress each 10% - plistener->setProgress (currentProgress); - } - } -#pragma omp atomic - tilesDone++; - } - -#ifdef _OPENMP - for(int i=0; isetProgress (1.0); -} -#undef TILEBORDER -#undef TILESIZE -#undef CACHESIZE - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//Emil's code for AMaZE -#include "amaze_interpolate_RT.cc"//AMaZE demosaic -#include "CA_correct_RT.cc"//Emil's CA auto correction -#include "cfa_linedn_RT.cc"//Emil's CA auto correction -#include "green_equil_RT.cc"//Emil's green channel equilibration -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -} /* namespace */ - +/* + * This file is part of RawTherapee. + * + * Copyright (c) 2004-2010 Gabor Horvath + * + * RawTherapee 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. + * + * RawTherapee 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 RawTherapee. If not, see . + */ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef _OPENMP +#include +#endif + +namespace rtengine { + +int loadRaw (const char* fname, struct RawImage* ri); + +extern const Settings* settings; + +#undef ABS +#undef MAX +#undef MIN +#undef DIST + +#define ABS(a) ((a)<0?-(a):(a)) +#define MAX(a,b) ((a)<(b)?(b):(a)) +#define MIN(a,b) ((a)>(b)?(b):(a)) +#define DIST(a,b) (ABS(a-b)) + +RawImageSource::RawImageSource () : ImageSource(), plistener(NULL), border(4), cache(NULL), green(NULL), red(NULL), blue(NULL) { + + hrmap[0] = NULL; + hrmap[1] = NULL; + hrmap[2] = NULL; + needhr = NULL; + hpmap = NULL; + oldmethod = "None"; + camProfile = NULL; + embProfile = NULL; +} + +RawImageSource::~RawImageSource () { + + delete idata; + if (ri) { + if (ri->allocation) + free (ri->allocation); + if (ri->data) + free (ri->data); + if (ri->profile_data) + free (ri->profile_data); + delete ri; + } + if (green) + freeArray(green, H); + if (red) + freeArray(red, H); + if (blue) + freeArray(blue, H); + + delete [] cache; + if (hrmap[0]!=NULL) { + int dh = H/HR_SCALE; + freeArray(hrmap[0], dh); + freeArray(hrmap[1], dh); + freeArray(hrmap[2], dh); + } + if (needhr) + freeArray(needhr, H); + if (hpmap) + freeArray(hpmap, H); + if (camProfile) + cmsCloseProfile (camProfile); + if (embProfile) + cmsCloseProfile (embProfile); +} + +void RawImageSource::transformRect (PreviewProps pp, int tran, int &ssx1, int &ssy1, int &width, int &height, int &fw) { + + pp.x += border; + pp.y += border; + + if (d1x) { + if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) { + pp.x /= 2; + pp.w = pp.w/2+1; + } + else { + pp.y /= 2; + pp.h = pp.h/2+1; + } + } + + int w = W, h = H; + if (fuji) { + w = ri->fuji_width * 2 + 1; + h = (H - ri->fuji_width)*2 + 1; + } + + int sw = w, sh = h; + if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) { + sw = h; + sh = w; + } + if( pp.w > sw-2*border) pp.w = sw-2*border; + if( pp.h > sh-2*border) pp.h = sh-2*border; + + int ppx = pp.x, ppy = pp.y; + if (tran & TR_HFLIP) + ppx = sw - pp.x - pp.w; + if (tran & TR_VFLIP) + ppy = sh - pp.y - pp.h; + + int sx1 = ppx; + int sy1 = ppy; + int sx2 = ppx + pp.w; + int sy2 = ppy + pp.h; + + if ((tran & TR_ROT) == TR_R180) { + sx1 = w - ppx - pp.w; + sy1 = h - ppy - pp.h; + sx2 = sx1 + pp.w; + sy2 = sy1 + pp.h; + } + else if ((tran & TR_ROT) == TR_R90) { + sx1 = ppy; + sy1 = h - ppx - pp.w; + sx2 = sx1 + pp.h; + sy2 = sy1 + pp.w; + } + else if ((tran & TR_ROT) == TR_R270) { + sx1 = w - ppy - pp.h; + sy1 = ppx; + sx2 = sx1 + pp.h; + sy2 = sy1 + pp.w; + } + + if (fuji) { + // atszamoljuk a koordinatakat fuji-ra: + ssx1 = (sx1+sy1) / 2; + ssy1 = (sy1 - sx2 ) / 2 + ri->fuji_width; + int ssx2 = (sx2+sy2) / 2 + 1; + int ssy2 = (sy2 - sx1) / 2 + ri->fuji_width; + fw = (sx2 - sx1) / 2 / pp.skip; + width = (ssx2 - ssx1) / pp.skip + ((ssx2 - ssx1) % pp.skip > 0); + height = (ssy2 - ssy1) / pp.skip + ((ssy2 - ssy1) % pp.skip > 0); + } + else { + ssx1 = sx1; + ssy1 = sy1; + width = (sx2 - sx1) / pp.skip + ((sx2 - sx1) % pp.skip > 0); + height = (sy2 - sy1) / pp.skip + ((sy2 - sy1) % pp.skip > 0); + } +} + +void RawImageSource::interpolate_image(Image16* image, HRecParams hrp, double rm, double gm, double bm, int skip, int tran, int fw, int imwidth, int imheight, int sx1, int sy1, int start, int end) +{ + + unsigned short* red = new unsigned short[imwidth]; + unsigned short* grn = new unsigned short[imwidth]; + unsigned short* blue = new unsigned short[imwidth]; + + for (int i=sy1 + start*skip ,ix=start; ixfilters && this->red && this->blue) { + for (int j=0,jx=sx1; jred[i][jx]); + grn[j] = CLIP(gm*green[i][jx]); + blue[j] = CLIP(bm*this->blue[i][jx]); + } + } else if(ri->filters) { + if (i==0) + interpolate_row_rb_mul_pp (red, blue, NULL, green[i], green[i+1], i, rm, gm, bm, sx1, imwidth, skip); + else if (i==H-1) + interpolate_row_rb_mul_pp (red, blue, green[i-1], green[i], NULL, i, rm, gm, bm, sx1, imwidth, skip); + else + interpolate_row_rb_mul_pp (red, blue, green[i-1], green[i], green[i+1], i, rm, gm, bm, sx1, imwidth, skip); + + for (int j=0,jx=sx1; jdata[i][jx*3+0]); + grn[j] = CLIP(gm*ri->data[i][jx*3+1]); + blue[j] = CLIP(bm*ri->data[i][jx*3+2]); + } + } + + if (hrp.enabled) + hlRecovery (hrp.method, red, grn, blue, i, sx1, imwidth, skip); + + transLine (red, grn, blue, ix, image, tran, imwidth, imheight, fw); + + } + delete [] red; + delete [] grn; + delete [] blue; +} + +void RawImageSource::getImage (ColorTemp ctemp, int tran, Image16* image, PreviewProps pp, HRecParams hrp, ColorManagementParams cmp) { + + isrcMutex.lock (); + + tran = defTransform (tran); + + // compute channel multipliers + double r, g, b, rm, gm, bm; + ctemp.getMultipliers (r, g, b); + rm = icoeff[0][0]*r + icoeff[0][1]*g + icoeff[0][2]*b; + gm = icoeff[1][0]*r + icoeff[1][1]*g + icoeff[1][2]*b; + bm = icoeff[2][0]*r + icoeff[2][1]*g + icoeff[2][2]*b; + rm = ri->camwb_red / rm; + gm = ri->camwb_green / gm; + bm = ri->camwb_blue / bm; + double mul_lum = 0.299*rm + 0.587*gm + 0.114*bm; + rm /= mul_lum; + gm /= mul_lum; + bm /= mul_lum; + + if (hrp.enabled) + defGain = log(ri->defgain) / log(2.0); + else { + defGain = 0.0; + rm *= ri->defgain; + gm *= ri->defgain; + bm *= ri->defgain; + } + + if (hrp.enabled==true && hrp.method=="Color" && hrmap[0]==NULL) + updateHLRecoveryMap_ColorPropagation (); + + // compute image area to render in order to provide the requested part of the image + int sx1, sy1, imwidth, imheight, fw; + transformRect (pp, tran, sx1, sy1, imwidth, imheight, fw); + + // check possible overflows + int maximwidth, maximheight; + if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) { + maximwidth = image->height; + maximheight = image->width; + } + else { + maximwidth = image->width; + maximheight = image->height; + } + if (d1x) + maximheight /= 2; + + // correct if overflow (very rare), but not fuji because it is corrected in transline + if (!fuji && imwidth>maximwidth) + imwidth = maximwidth; + if (!fuji && imheight>maximheight) + imheight = maximheight; + +#ifdef _OPENMP +#pragma omp parallel +{ + int tid = omp_get_thread_num(); + int nthreads = omp_get_num_threads(); + int blk = imheight/nthreads; + + + if (tidwidth%2==0) || ((tran & TR_ROT) == TR_R180 && image->height%2+image->width%2==1) || ((tran & TR_ROT) == TR_R270 && image->height%2==0); + // first row + for (int j=1+a; jwidth-1; j+=2) { + image->r[0][j] = (image->r[1][j] + image->r[0][j+1] + image->r[0][j-1]) / 3; + image->g[0][j] = (image->g[1][j] + image->g[0][j+1] + image->g[0][j-1]) / 3; + image->b[0][j] = (image->b[1][j] + image->b[0][j+1] + image->b[0][j-1]) / 3; + } + // other rows + for (int i=1; iheight-1; i++) { + for (int j=2-(a+i+1)%2; jwidth-1; j+=2) { + // edge-adaptive interpolation + double dh = (ABS(image->r[i][j+1] - image->r[i][j-1]) + ABS(image->g[i][j+1] - image->g[i][j-1]) + ABS(image->b[i][j+1] - image->b[i][j-1])) / 1.0; + double dv = (ABS(image->r[i+1][j] - image->r[i-1][j]) + ABS(image->g[i+1][j] - image->g[i-1][j]) + ABS(image->b[i+1][j] - image->b[i-1][j])) / 1.0; + double eh = 1.0 / (1.0 + dh); + double ev = 1.0 / (1.0 + dv); + image->r[i][j] = (eh * (image->r[i][j+1] + image->r[i][j-1]) + ev * (image->r[i+1][j] + image->r[i-1][j])) / (2.0 * (eh + ev)); + image->g[i][j] = (eh * (image->g[i][j+1] + image->g[i][j-1]) + ev * (image->g[i+1][j] + image->g[i-1][j])) / (2.0 * (eh + ev)); + image->b[i][j] = (eh * (image->b[i][j+1] + image->b[i][j-1]) + ev * (image->b[i+1][j] + image->b[i-1][j])) / (2.0 * (eh + ev)); + } + // first pixel + if (2-(a+i+1)%2==2) { + image->r[i][0] = (image->r[i+1][0] + image->r[i-1][0] + image->r[i][1]) / 3; + image->g[i][0] = (image->g[i+1][0] + image->g[i-1][0] + image->g[i][1]) / 3; + image->b[i][0] = (image->b[i+1][0] + image->b[i-1][0] + image->b[i][1]) / 3; + } + // last pixel + if (2-(a+i+image->width)%2==2) { + image->r[i][image->width-1] = (image->r[i+1][image->width-1] + image->r[i-1][image->width-1] + image->r[i][image->width-2]) / 3; + image->g[i][image->width-1] = (image->g[i+1][image->width-1] + image->g[i-1][image->width-1] + image->g[i][image->width-2]) / 3; + image->b[i][image->width-1] = (image->b[i+1][image->width-1] + image->b[i-1][image->width-1] + image->b[i][image->width-2]) / 3; + } + } + // last row + int b = (a==1 && image->height%2) || (a==0 && image->height%2==0); + for (int j=1+b; jwidth-1; j+=2) { + image->r[image->height-1][j] = (image->r[image->height-2][j] + image->r[image->height-1][j+1] + image->r[image->height-1][j-1]) / 3; + image->g[image->height-1][j] = (image->g[image->height-2][j] + image->g[image->height-1][j+1] + image->g[image->height-1][j-1]) / 3; + image->b[image->height-1][j] = (image->b[image->height-2][j] + image->b[image->height-1][j+1] + image->b[image->height-1][j-1]) / 3; + } + } + + // Flip if needed + if (tran & TR_HFLIP) + hflip (image); + if (tran & TR_VFLIP) + vflip (image); + + // Color correction + if (ri->filters && pp.skip==1) + correction_YIQ_LQ (image, settings->colorCorrectionSteps); + + // Applying postmul + colorSpaceConversion (image, cmp, embProfile, camProfile, cam, defGain); + + isrcMutex.unlock (); +} + + +void RawImageSource::cfa_clean(float thresh) //Emil's hot/dead pixel removal -- only filters egregiously impulsive pixels + { + // local variables + int rr, cc; + int gin, g[8]; + + float eps=1.0;//tolerance to avoid dividing by zero + float p[8]; + float pixave, pixratio; + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + //The cleaning algorithm starts here + + for (rr=4; rr < H-4; rr++) + for (cc=4; cc < W-4; cc++) { + + //pixel neighbor average + gin=ri->data[rr][cc]; + g[0]=ri->data[rr-2][cc-2]; + g[1]=ri->data[rr-2][cc]; + g[2]=ri->data[rr-2][cc+2]; + g[3]=ri->data[rr][cc-2]; + g[4]=ri->data[rr][cc+2]; + g[5]=ri->data[rr+2][cc-2]; + g[6]=ri->data[rr+2][cc]; + g[7]=ri->data[rr+2][cc+2]; + + pixave=(float)(g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6]+g[7])/8; + pixratio=MIN(gin,pixave)/(eps+MAX(gin,pixave)); + + if (pixratio > thresh) continue; + + p[0]=1/(eps+fabs(g[0]-gin)+fabs(g[0]-ri->data[rr-4][cc-4])+fabs(ri->data[rr-1][cc-1]-ri->data[rr-3][cc-3])); + p[1]=1/(eps+fabs(g[1]-gin)+fabs(g[1]-ri->data[rr-4][cc])+fabs(ri->data[rr-1][cc]-ri->data[rr-3][cc])); + p[2]=1/(eps+fabs(g[2]-gin)+fabs(g[2]-ri->data[rr-4][cc+4])+fabs(ri->data[rr-1][cc+1]-ri->data[rr-3][cc+3])); + p[3]=1/(eps+fabs(g[3]-gin)+fabs(g[3]-ri->data[rr][cc-4])+fabs(ri->data[rr][cc-1]-ri->data[rr][cc-3])); + p[4]=1/(eps+fabs(g[4]-gin)+fabs(g[4]-ri->data[rr][cc+4])+fabs(ri->data[rr][cc+1]-ri->data[rr][cc+3])); + p[5]=1/(eps+fabs(g[5]-gin)+fabs(g[5]-ri->data[rr+4][cc-4])+fabs(ri->data[rr+1][cc-1]-ri->data[rr+3][cc-3])); + p[6]=1/(eps+fabs(g[6]-gin)+fabs(g[6]-ri->data[rr+4][cc])+fabs(ri->data[rr+1][cc]-ri->data[rr+3][cc])); + p[7]=1/(eps+fabs(g[7]-gin)+fabs(g[7]-ri->data[rr+4][cc+4])+fabs(ri->data[rr+1][cc+1]-ri->data[rr+3][cc+3])); + + ri->data[rr][cc] = (int)((g[0]*p[0]+g[1]*p[1]+g[2]*p[2]+g[3]*p[3]+g[4]*p[4]+ \ + g[5]*p[5]+g[6]*p[6]+g[7]*p[7])/(p[0]+p[1]+p[2]+p[3]+p[4]+p[5]+p[6]+p[7])); + + }//now impulsive values have been corrected + + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + } + + +void RawImageSource::rotateLine (unsigned short* line, unsigned short** channel, int tran, int i, int w, int h) { + + if ((tran & TR_ROT) == TR_R180) + for (int j=0; j=0 && yheight && y>=0 && xwidth) { + image->r[image->height-1-y][image->width-1-x] = red[j]; + image->g[image->height-1-y][image->width-1-x] = green[j]; + image->b[image->height-1-y][image->width-1-x] = blue[j]; + } + } + } + else if ((tran & TR_ROT) == TR_R270) { + int end = MIN(h+fw-i, w-fw+i); + for (int j=start; j=0 && xheight && y>=0 && ywidth) { + image->r[image->height-1-x][y] = red[j]; + image->g[image->height-1-x][y] = green[j]; + image->b[image->height-1-x][y] = blue[j]; + } + } + } + else if ((tran & TR_ROT) == TR_R90) { + int end = MIN(h+fw-i, w-fw+i); + for (int j=start; j=0 && ywidth && y>=0 && xheight) { + image->r[x][image->width-1-y] = red[j]; + image->g[x][image->width-1-y] = green[j]; + image->b[x][image->width-1-y] = blue[j]; + } + } + } + else { + int end = MIN(h+fw-i, w-fw+i); + for (int j=start; j=0 && yheight && y>=0 && xwidth) { + image->r[y][x] = red[j]; + image->g[y][x] = green[j]; + image->b[y][x] = blue[j]; + } + } + } + } + // Nikon D1X vertical interpolation + coarse rotation + else if (d1x) { + // copy new pixels + if ((tran & TR_ROT) == TR_R180) { + for (int j=0; jr[2*imheight-2-2*i][imwidth-1-j] = red[j]; + image->g[2*imheight-2-2*i][imwidth-1-j] = green[j]; + image->b[2*imheight-2-2*i][imwidth-1-j] = blue[j]; + } + + if (i==1 || i==2) { // linear interpolation + int row = 2*imheight-1-2*i; + for (int j=0; jr[row][col] = (red[j] + image->r[row+1][col]) >> 1; + image->g[row][col] = (green[j] + image->g[row+1][col]) >> 1; + image->b[row][col] = (blue[j] + image->b[row+1][col]) >> 1; + } + } + else if (i==imheight-1) { + int row = 2*imheight-1-2*i; + for (int j=0; jr[row][col] = (red[j] + image->r[row+1][col]) >> 1; + image->g[row][col] = (green[j] + image->g[row+1][col]) >> 1; + image->b[row][col] = (blue[j] + image->b[row+1][col]) >> 1; + } + row = 2*imheight-1-2*i+2; + for (int j=0; jr[row][col] = (red[j] + image->r[row+1][col]) >> 1; + image->g[row][col] = (green[j] + image->g[row+1][col]) >> 1; + image->b[row][col] = (blue[j] + image->b[row+1][col]) >> 1; + } + } + else if (i>2 && ir[row][col] = CLIP((int)(-0.0625*red[j] + 0.5625*image->r[row-1][col] + 0.5625*image->r[row+1][col] - 0.0625*image->r[row+3][col])); + image->g[row][col] = CLIP((int)(-0.0625*green[j] + 0.5625*image->g[row-1][col] + 0.5625*image->g[row+1][col] - 0.0625*image->g[row+3][col])); + image->b[row][col] = CLIP((int)(-0.0625*blue[j] + 0.5625*image->b[row-1][col] + 0.5625*image->b[row+1][col] - 0.0625*image->b[row+3][col])); + } + } + } + else if ((tran & TR_ROT) == TR_R90) { + for (int j=0; jr[j][2*imheight-2-2*i] = red[j]; + image->g[j][2*imheight-2-2*i] = green[j]; + image->b[j][2*imheight-2-2*i] = blue[j]; + } + if (i==1 || i==2) { // linear interpolation + int col = 2*imheight-1-2*i; + for (int j=0; jr[j][col] = (red[j] + image->r[j][col+1]) >> 1; + image->g[j][col] = (green[j] + image->g[j][col+1]) >> 1; + image->b[j][col] = (blue[j] + image->b[j][col+1]) >> 1; + } + } + else if (i==imheight-1) { + int col = 2*imheight-1-2*i; + for (int j=0; jr[j][col] = (red[j] + image->r[j][col+1]) >> 1; + image->g[j][col] = (green[j] + image->g[j][col+1]) >> 1; + image->b[j][col] = (blue[j] + image->b[j][col+1]) >> 1; + } + col = 2*imheight-1-2*i+2; + for (int j=0; jr[j][col] = (red[j] + image->r[j][col+1]) >> 1; + image->g[j][col] = (green[j] + image->g[j][col+1]) >> 1; + image->b[j][col] = (blue[j] + image->b[j][col+1]) >> 1; + } + } + else if (i>2 && ir[j][col] = CLIP((int)(-0.0625*red[j] + 0.5625*image->r[j][col-1] + 0.5625*image->r[j][col+1] - 0.0625*image->r[j][col+3])); + image->g[j][col] = CLIP((int)(-0.0625*green[j] + 0.5625*image->g[j][col-1] + 0.5625*image->g[j][col+1] - 0.0625*image->g[j][col+3])); + image->b[j][col] = CLIP((int)(-0.0625*blue[j] + 0.5625*image->b[j][col-1] + 0.5625*image->b[j][col+1] - 0.0625*image->b[j][col+3])); + } + } + } + else if ((tran & TR_ROT) == TR_R270) { + for (int j=0; jr[imwidth-1-j][2*i] = red[j]; + image->g[imwidth-1-j][2*i] = green[j]; + image->b[imwidth-1-j][2*i] = blue[j]; + } + if (i==1 || i==2) { // linear interpolation + for (int j=0; jr[row][2*i-1] = (red[j] + image->r[row][2*i-2]) >> 1; + image->g[row][2*i-1] = (green[j] + image->g[row][2*i-2]) >> 1; + image->b[row][2*i-1] = (blue[j] + image->b[row][2*i-2]) >> 1; + } + } + else if (i==imheight-1) { + for (int j=0; jr[row][2*i-1] = (red[j] + image->r[row][2*i-2]) >> 1; + image->g[row][2*i-1] = (green[j] + image->g[row][2*i-2]) >> 1; + image->b[row][2*i-1] = (blue[j] + image->b[row][2*i-2]) >> 1; + image->r[row][2*i-3] = (image->r[row][2*i-2] + image->r[row][2*i-4]) >> 1; + image->g[row][2*i-3] = (image->g[row][2*i-2] + image->g[row][2*i-4]) >> 1; + image->b[row][2*i-3] = (image->b[row][2*i-2] + image->b[row][2*i-4]) >> 1; + } + } + else if (i>0 && ir[row][2*i-3] = CLIP((int)(-0.0625*red[j] + 0.5625*image->r[row][2*i-2] + 0.5625*image->r[row][2*i-4] - 0.0625*image->r[row][2*i-6])); + image->g[row][2*i-3] = CLIP((int)(-0.0625*green[j] + 0.5625*image->g[row][2*i-2] + 0.5625*image->g[row][2*i-4] - 0.0625*image->g[row][2*i-6])); + image->b[row][2*i-3] = CLIP((int)(-0.0625*blue[j] + 0.5625*image->b[row][2*i-2] + 0.5625*image->b[row][2*i-4] - 0.0625*image->b[row][2*i-6])); + } + } + } + else { + rotateLine (red, image->r, tran, 2*i, imwidth, imheight); + rotateLine (green, image->g, tran, 2*i, imwidth, imheight); + rotateLine (blue, image->b, tran, 2*i, imwidth, imheight); + + if (i==1 || i==2) { // linear interpolation + for (int j=0; jr[2*i-1][j] = (red[j] + image->r[2*i-2][j]) >> 1; + image->g[2*i-1][j] = (green[j] + image->g[2*i-2][j]) >> 1; + image->b[2*i-1][j] = (blue[j] + image->b[2*i-2][j]) >> 1; + } + } + else if (i==imheight-1) { + for (int j=0; jr[2*i-3][j] = (image->r[2*i-4][j] + image->r[2*i-2][j]) >> 1; + image->g[2*i-3][j] = (image->g[2*i-4][j] + image->g[2*i-2][j]) >> 1; + image->b[2*i-3][j] = (image->b[2*i-4][j] + image->b[2*i-2][j]) >> 1; + image->r[2*i-1][j] = (red[j] + image->r[2*i-2][j]) >> 1; + image->g[2*i-1][j] = (green[j] + image->g[2*i-2][j]) >> 1; + image->b[2*i-1][j] = (blue[j] + image->b[2*i-2][j]) >> 1; + } + } + else if (i>2 && ir[2*i-3][j] = CLIP((int)(-0.0625*red[j] + 0.5625*image->r[2*i-2][j] + 0.5625*image->r[2*i-4][j] - 0.0625*image->r[2*i-6][j])); + image->g[2*i-3][j] = CLIP((int)(-0.0625*green[j] + 0.5625*image->g[2*i-2][j] + 0.5625*image->g[2*i-4][j] - 0.0625*image->g[2*i-6][j])); + image->b[2*i-3][j] = CLIP((int)(-0.0625*blue[j] + 0.5625*image->b[2*i-2][j] + 0.5625*image->b[2*i-4][j] - 0.0625*image->b[2*i-6][j])); + } + } + } + } + // other (conventional) CCD coarse rotation + else { + rotateLine (red, image->r, tran, i, imwidth, imheight); + rotateLine (green, image->g, tran, i, imwidth, imheight); + rotateLine (blue, image->b, tran, i, imwidth, imheight); + } +} + +void RawImageSource::getFullSize (int& w, int& h, int tr) { + + tr = defTransform (tr); + + if (fuji) { + w = ri->fuji_width * 2 + 1; + h = (H - ri->fuji_width)*2 + 1; + } + else if (d1x) { + w = W; + h = 2*H-1; + } + else { + w = W; + h = H; + } + + if ((tr & TR_ROT) == TR_R90 || (tr & TR_ROT) == TR_R270) { + int tmp = w; + w = h; + h = tmp; + } + w -= 2 * border; + h -= 2 * border; +} + +void RawImageSource::getSize (int tran, PreviewProps pp, int& w, int& h) { + + tran = defTransform (tran); + +// if (fuji) { +// return; +// } +// else if (d1x) { +// return; +// } +// else { + w = pp.w / pp.skip + (pp.w % pp.skip > 0); + h = pp.h / pp.skip + (pp.h % pp.skip > 0); +// } +} + +void RawImageSource::hflip (Image16* image) { + int width = image->width; + int height = image->height; + + unsigned short* rowr = new unsigned short[width]; + unsigned short* rowg = new unsigned short[width]; + unsigned short* rowb = new unsigned short[width]; + for (int i=0; ir[i][width-1-j]; + rowg[j] = image->g[i][width-1-j]; + rowb[j] = image->b[i][width-1-j]; + } + memcpy (image->r[i], rowr, width*sizeof(unsigned short)); + memcpy (image->g[i], rowg, width*sizeof(unsigned short)); + memcpy (image->b[i], rowb, width*sizeof(unsigned short)); + } + delete [] rowr; + delete [] rowg; + delete [] rowb; +} + +void RawImageSource::vflip (Image16* image) { + int width = image->width; + int height = image->height; + + register unsigned short tmp; + for (int i=0; ir[i][j]; + image->r[i][j] = image->r[height-1-i][j]; + image->r[height-1-i][j] = tmp; + tmp = image->g[i][j]; + image->g[i][j] = image->g[height-1-i][j]; + image->g[height-1-i][j] = tmp; + tmp = image->b[i][j]; + image->b[i][j] = image->b[height-1-i][j]; + image->b[height-1-i][j] = tmp; + } +} + +void RawImageSource::inverse33 (double (*coeff)[3], double (*icoeff)[3]) { + double nom = coeff[0][2]*coeff[1][1]*coeff[2][0] - coeff[0][1]*coeff[1][2]*coeff[2][0] - coeff[0][2]*coeff[1][0]*coeff[2][1] + coeff[0][0]*coeff[1][2]*coeff[2][1] + coeff[0][1]*coeff[1][0]*coeff[2][2] - coeff[0][0]*coeff[1][1]*coeff[2][2]; + icoeff[0][0] = (coeff[1][2]*coeff[2][1]-coeff[1][1]*coeff[2][2]) / nom; + icoeff[0][1] = -(coeff[0][2]*coeff[2][1]-coeff[0][1]*coeff[2][2]) / nom; + icoeff[0][2] = (coeff[0][2]*coeff[1][1]-coeff[0][1]*coeff[1][2]) / nom; + icoeff[1][0] = -(coeff[1][2]*coeff[2][0]-coeff[1][0]*coeff[2][2]) / nom; + icoeff[1][1] = (coeff[0][2]*coeff[2][0]-coeff[0][0]*coeff[2][2]) / nom; + icoeff[1][2] = -(coeff[0][2]*coeff[1][0]-coeff[0][0]*coeff[1][2]) / nom; + icoeff[2][0] = (coeff[1][1]*coeff[2][0]-coeff[1][0]*coeff[2][1]) / nom; + icoeff[2][1] = -(coeff[0][1]*coeff[2][0]-coeff[0][0]*coeff[2][1]) / nom; + icoeff[2][2] = (coeff[0][1]*coeff[1][0]-coeff[0][0]*coeff[1][1]) / nom; +} + +int RawImageSource::load (Glib::ustring fname) { + + fileName = fname; + + if (plistener) { + plistener->setProgressStr ("Decoding..."); + plistener->setProgress (0.0); + } + + ri = new RawImage; + int res = loadRaw (fname.c_str(), ri); + if (res) + return res; + + if(red) { + delete red; + red = 0; + } + if(green) { + delete green; + green = 0; + } + if(blue) { + delete blue; + blue = 0; + } + + W = ri->width; + H = ri->height; + + d1x = !strcmp(ri->model, "D1X"); + fuji = ri->fuji_width; + if (d1x) + border = 8; + + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + coeff[i][j] = ri->coeff[i][j]; + + // compute inverse of the color transformation matrix + inverse33 (coeff, icoeff); + + double cam_r = coeff[0][0]*ri->camwb_red + coeff[0][1]*ri->camwb_green + coeff[0][2]*ri->camwb_blue; + double cam_g = coeff[1][0]*ri->camwb_red + coeff[1][1]*ri->camwb_green + coeff[1][2]*ri->camwb_blue; + double cam_b = coeff[2][0]*ri->camwb_red + coeff[2][1]*ri->camwb_green + coeff[2][2]*ri->camwb_blue; + + wb = ColorTemp (cam_r, cam_g, cam_b); + + double tr = icoeff[0][0] * cam_r + icoeff[0][1] * cam_g + icoeff[0][2] * cam_b; + double tg = icoeff[1][0] * cam_r + icoeff[1][1] * cam_g + icoeff[1][2] * cam_b; + double tb = icoeff[2][0] * cam_r + icoeff[2][1] * cam_g + icoeff[2][2] * cam_b; + + // create profile + memset (cam, 0, sizeof(cam)); + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + for (int k=0; k<3; k++) + cam[i][j] += coeff[k][i] * sRGB_d50[k][j]; + camProfile = iccStore.createFromMatrix (cam, false, "Camera"); + inverse33 (cam, icam); + + if (ri->profile_data) + embProfile = cmsOpenProfileFromMem (ri->profile_data, ri->profile_len); + + defGain = log(ri->defgain) / log(2.0); + + RawMetaDataLocation rml; + rml.exifBase = ri->exifbase; + rml.ciffBase = ri->ciff_base; + rml.ciffLength = ri->ciff_len; + + idata = new ImageData (fname, &rml); + + // check if it is an olympus E camera, if yes, compute G channel pre-compensation factors + if (settings->greenthresh || (((idata->getMake().size()>=7 && idata->getMake().substr(0,7)=="OLYMPUS" && idata->getModel()[0]=='E') || (idata->getMake().size()>=9 && idata->getMake().substr(0,7)=="Panasonic")) && settings->demosaicMethod!="vng4" && ri->filters) ) { + // global correction + int ng1=0, ng2=0; + double avgg1=0, avgg2=0; + for (int i=border; idata[i][j]; + ng1++; + } + else { + avgg2 += ri->data[i][j]; + ng2++; + } + } + double corrg1 = ((double)avgg1/ng1 + (double)avgg2/ng2) / 2.0 / ((double)avgg1/ng1); + double corrg2 = ((double)avgg1/ng1 + (double)avgg2/ng2) / 2.0 / ((double)avgg2/ng2); + for (int i=border; idata[i][j] = CLIP(ri->data[i][j] * (i%2 ? corrg2 : corrg1)); + } + + + // local correction in a 9x9 box + if (settings->greenthresh) { +//Emil's green equilbration + if (plistener) { + plistener->setProgressStr ("Green equilibrate..."); + plistener->setProgress (0.0); + } + green_equilibrate(0.01*(settings->greenthresh)); + +/* unsigned short* corr_alloc = new unsigned short[W*H]; + unsigned short** corr_data = new unsigned short* [H]; + for (int i=0; iallocation, W*H*sizeof(unsigned short)); + for (int i=border; idata[i-4][j-4] + ri->data[i-4][j-2] + ri->data[i-4][j] + ri->data[i-4][j+2] + + ri->data[i-2][j-4] + ri->data[i-2][j-2] + ri->data[i-2][j] + ri->data[i-2][j+2] + + ri->data[i][j-4] + ri->data[i][j-2] + ri->data[i][j] + ri->data[i][j+2] + + ri->data[i+2][j-4] + ri->data[i+2][j-2] + ri->data[i+2][j] + ri->data[i+2][j+2]; + unsigned int ag2 = ri->data[i-3][j-3] + ri->data[i-3][j-1] + ri->data[i-3][j+1] + ri->data[i-3][j+1] + + ri->data[i-1][j-3] + ri->data[i-1][j-1] + ri->data[i-1][j+1] + ri->data[i-1][j+1] + + ri->data[i+1][j-3] + ri->data[i+1][j-1] + ri->data[i+1][j+1] + ri->data[i+1][j+1] + + ri->data[i+3][j-3] + ri->data[i+3][j-1] + ri->data[i+3][j+1] + ri->data[i+3][j+1]; + unsigned int val = (ri->data[i][j] + ri->data[i][j] * ag2 / ag1) / 2; + corr_data[i][j] = CLIP (val); + } + memcpy (ri->allocation, corr_alloc, W*H*sizeof(unsigned short)); + delete corr_alloc; + delete corr_data; */ + + } + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +//Emil's CFA hot/dead pixel filter + + if (settings->hotdeadpix_filt) { + if (plistener) { + plistener->setProgressStr ("Hot/Dead Pixel Filter..."); + plistener->setProgress (0.0); + } + + cfa_clean(0.1); + } + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +//Emil's line noise filter + + if (settings->linenoise) { + if (plistener) { + plistener->setProgressStr ("Line Denoise..."); + plistener->setProgress (0.0); + } + + cfa_linedn(0.00002*(settings->linenoise)); + } + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +//Emil's CA auto correction + + if (settings->ca_autocorrect) { + if (plistener) { + plistener->setProgressStr ("CA Auto Correction..."); + plistener->setProgress (0.0); + } + + CA_correct_RT(); + } +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + if (ri->filters) { + //MyTime t1,t2; + //t1.set(); + // demosaic + if (settings->demosaicMethod=="hphd") + hphd_demosaic (); + else if (settings->demosaicMethod=="vng4") + vng4_demosaic (); + else if (settings->demosaicMethod=="ahd") + ahd_demosaic (); + else if (settings->demosaicMethod=="bilinear") + bilinear_demosaic(); + //else if (settings->demosaicMethod=="ppg") + // ppg_demosaic (); + else if (settings->demosaicMethod=="amaze") + amaze_demosaic_RT ();//Emil's code for AMaZE + else if (settings->demosaicMethod=="dcb") + dcb_demosaic(settings->dcb_iterations, settings->dcb_enhance? 1:0); + else + eahd_demosaic (); + //t2.set(); + //printf("Demosaicing:%d usec\n",t2.etime(t1)); + } + + + if (plistener) { + plistener->setProgressStr ("Ready."); + plistener->setProgress (1.0); + } + + return 0; +} + +int RawImageSource::defTransform (int tran) { + + int deg = ri->rotate_deg; + if ((tran & TR_ROT) == TR_R180) + deg += 180; + else if ((tran & TR_ROT) == TR_R90) + deg += 90; + else if ((tran & TR_ROT) == TR_R270) + deg += 270; + deg %= 360; + + int ret = 0; + if (deg==90) + ret |= TR_R90; + else if (deg==180) + ret |= TR_R180; + else if (deg==270) + ret |= TR_R270; + if (tran & TR_HFLIP) + ret |= TR_HFLIP; + if (tran & TR_VFLIP) + ret |= TR_VFLIP; + return ret; +} + +void RawImageSource::correction_YIQ_LQ_ (Image16* im, int row_from, int row_to) { + + int W = im->width; + + int** rbconv_Y = new int*[3]; + int** rbconv_I = new int*[3]; + int** rbconv_Q = new int*[3]; + int** rbout_I = new int*[3]; + int** rbout_Q = new int*[3]; + for (int i=0; i<3; i++) { + rbconv_Y[i] = new int[W]; + rbconv_I[i] = new int[W]; + rbconv_Q[i] = new int[W]; + rbout_I[i] = new int[W]; + rbout_Q[i] = new int[W]; + } + + int* row_I = new int[W]; + int* row_Q = new int[W]; + + int* pre1_I = new int[3]; + int* pre2_I = new int[3]; + int* post1_I = new int[3]; + int* post2_I = new int[3]; + int middle_I[6]; + int* pre1_Q = new int[3]; + int* pre2_Q = new int[3]; + int* post1_Q = new int[3]; + int* post2_Q = new int[3]; + int middle_Q[6]; + int* tmp; + + int ppx=0, px=(row_from-1)%3, cx=row_from%3, nx=0; + + convert_row_to_YIQ (im->r[row_from-1], im->g[row_from-1], im->b[row_from-1], rbconv_Y[px], rbconv_I[px], rbconv_Q[px], W); + convert_row_to_YIQ (im->r[row_from], im->g[row_from], im->b[row_from], rbconv_Y[cx], rbconv_I[cx], rbconv_Q[cx], W); + + for (int j=0; jr[i+1], im->g[i+1], im->b[i+1], rbconv_Y[nx], rbconv_I[nx], rbconv_Q[nx], W); + + SORT3(rbconv_I[px][0],rbconv_I[cx][0],rbconv_I[nx][0],pre1_I[0],pre1_I[1],pre1_I[2]); + SORT3(rbconv_I[px][1],rbconv_I[cx][1],rbconv_I[nx][1],pre2_I[0],pre2_I[1],pre2_I[2]); + SORT3(rbconv_Q[px][0],rbconv_Q[cx][0],rbconv_Q[nx][0],pre1_Q[0],pre1_Q[1],pre1_Q[2]); + SORT3(rbconv_Q[px][1],rbconv_Q[cx][1],rbconv_Q[nx][1],pre2_Q[0],pre2_Q[1],pre2_Q[2]); + + // median I channel + for (int j=1; jrow_from) { + for (int j=1; jr[i-1], im->g[i-1], im->b[i-1], rbconv_Y[px], row_I, row_Q, W); + } + } + // blur last 3 row and finalize H-1th row + for (int j=1; jr[row_to-1], im->g[row_to-1], im->b[row_to-1], rbconv_Y[cx], row_I, row_Q, W); + + freeArray(rbconv_Y, 3); + freeArray(rbconv_I, 3); + freeArray(rbconv_Q, 3); + freeArray(rbout_I, 3); + freeArray(rbout_Q, 3); + delete [] row_I; + delete [] row_Q; + delete [] pre1_I; + delete [] pre2_I; + delete [] post1_I; + delete [] post2_I; + delete [] pre1_Q; + delete [] pre2_Q; + delete [] post1_Q; + delete [] post2_Q; +} + + +void RawImageSource::correction_YIQ_LQ (Image16* im, int times) { + + if (im->height<4) + return; + + for (int t=0; theight-2)/nthreads; + + if (tidheight - 1); + } +#else + correction_YIQ_LQ_ (im, 1 , im->height - 1); +#endif + } +} + + +void RawImageSource::colorSpaceConversion (Image16* im, ColorManagementParams cmp, cmsHPROFILE embedded, cmsHPROFILE camprofile, double camMatrix[3][3], double& defgain) { + + if (cmp.input == "(none)") + return; + + MyTime t1, t2, t3; + + t1.set (); + + cmsHPROFILE in; + cmsHPROFILE out; + + Glib::ustring inProfile = cmp.input; + + if (inProfile=="(embedded)") { + if (embedded) + in = embedded; + else + in = camprofile; + } + else if (inProfile=="(camera)" || inProfile=="") + in = camprofile; + else { + in = iccStore.getProfile (inProfile); + if (in==NULL) + inProfile = "(camera)"; + } + + + if (inProfile=="(camera)" || inProfile=="" || (inProfile=="(embedded)" && !embedded)) { + // in this case we avoid using the slllllooooooowwww lcms + +// out = iccStore.workingSpace (wProfile); +// hTransform = cmsCreateTransform (in, TYPE_RGB_16_PLANAR, out, TYPE_RGB_16_PLANAR, settings->colorimetricIntent, cmsFLAGS_MATRIXINPUT | cmsFLAGS_MATRIXOUTPUT);//cmsFLAGS_MATRIXINPUT | cmsFLAGS_MATRIXOUTPUT); +// cmsDoTransform (hTransform, im->data, im->data, im->planestride/2); +// cmsDeleteTransform(hTransform); + TMatrix work = iccStore.workingSpaceInverseMatrix (cmp.working); + double mat[3][3] = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}}; + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + for (int k=0; k<3; k++) + mat[i][j] += camMatrix[i][k] * work[k][j]; + + for (int i=0; iheight; i++) + for (int j=0; jwidth; j++) { + + int newr = mat[0][0]*im->r[i][j] + mat[1][0]*im->g[i][j] + mat[2][0]*im->b[i][j]; + int newg = mat[0][1]*im->r[i][j] + mat[1][1]*im->g[i][j] + mat[2][1]*im->b[i][j]; + int newb = mat[0][2]*im->r[i][j] + mat[1][2]*im->g[i][j] + mat[2][2]*im->b[i][j]; + + im->r[i][j] = CLIP(newr); + im->g[i][j] = CLIP(newg); + im->b[i][j] = CLIP(newb); + } + } + else { + out = iccStore.workingSpace (cmp.working); +// out = iccStore.workingSpaceGamma (wProfile); + lcmsMutex->lock (); + cmsHTRANSFORM hTransform = cmsCreateTransform (in, TYPE_RGB_16_PLANAR, out, TYPE_RGB_16_PLANAR, settings->colorimetricIntent, 0); + lcmsMutex->unlock (); + if (hTransform) { + if (cmp.gammaOnInput) { + double gd = pow (2.0, defgain); + defgain = 0.0; + for (int i=0; iheight; i++) + for (int j=0; jwidth; j++) { + im->r[i][j] = CurveFactory::gamma (CLIP(defgain*im->r[i][j])); + im->g[i][j] = CurveFactory::gamma (CLIP(defgain*im->g[i][j])); + im->b[i][j] = CurveFactory::gamma (CLIP(defgain*im->b[i][j])); + } + } + cmsDoTransform (hTransform, im->data, im->data, im->planestride/2); + } + else { + lcmsMutex->lock (); + hTransform = cmsCreateTransform (camprofile, TYPE_RGB_16_PLANAR, out, TYPE_RGB_16_PLANAR, settings->colorimetricIntent, 0); + lcmsMutex->unlock (); + cmsDoTransform (hTransform, im->data, im->data, im->planestride/2); + } + cmsDeleteTransform(hTransform); + } + t3.set (); +// printf ("ICM TIME: %d\n", t3.etime(t1)); +} + +void RawImageSource::eahd_demosaic () { + + if (plistener) { + plistener->setProgressStr ("Demosaicing..."); + plistener->setProgress (0.0); + } + + // prepare chache and constants for cielab conversion + lc00 = (0.412453 * coeff[0][0] + 0.357580 * coeff[1][0] + 0.180423 * coeff[2][0]) / 0.950456; + lc01 = (0.412453 * coeff[0][1] + 0.357580 * coeff[1][1] + 0.180423 * coeff[2][1]) / 0.950456; + lc02 = (0.412453 * coeff[0][2] + 0.357580 * coeff[1][2] + 0.180423 * coeff[2][2]) / 0.950456; + + lc10 = 0.212671 * coeff[0][0] + 0.715160 * coeff[1][0] + 0.072169 * coeff[2][0]; + lc11 = 0.212671 * coeff[0][1] + 0.715160 * coeff[1][1] + 0.072169 * coeff[2][1]; + lc12 = 0.212671 * coeff[0][2] + 0.715160 * coeff[1][2] + 0.072169 * coeff[2][2]; + + lc20 = (0.019334 * coeff[0][0] + 0.119193 * coeff[1][0] + 0.950227 * coeff[2][0]) / 1.088754; + lc21 = (0.019334 * coeff[0][1] + 0.119193 * coeff[1][1] + 0.950227 * coeff[2][1]) / 1.088754; + lc22 = (0.019334 * coeff[0][2] + 0.119193 * coeff[1][2] + 0.950227 * coeff[2][2]) / 1.088754; + + int maxindex = 2*65536; + cache = new double[maxindex]; + threshold = (int)(0.008856*CMAXVAL); + for (int i=0; isv) { // fixate horizontal pixel + dLmaph[dmi] = DIST(lLh[ix][j], lLh[idx][j+y]); + dCamaph[dmi] = DIST(lah[ix][j], lah[idx][j+y]); + dCbmaph[dmi] = DIST(lbh[ix][j], lbh[idx][j+y]); + dLmapv[dmi] = DIST(lLv[ix][j], lLh[idx][j+y]); + dCamapv[dmi] = DIST(lav[ix][j], lah[idx][j+y]); + dCbmapv[dmi] = DIST(lbv[ix][j], lbh[idx][j+y]); + } + else if (shdata[i-1][j]; + else { + hc = homh[imx][j]; + vc = homv[imx][j]; + if (hc > vc) + green[i-1][j] = gh[(i-1)%4][j]; + else if (hc < vc) + green[i-1][j] = gv[(i-1)%4][j]; + else + green[i-1][j] = (gh[(i-1)%4][j] + gv[(i-1)%4][j]) / 2; + } + } + + if (!(i%20) && plistener) + plistener->setProgress ((double)i / (H-2)); + } + // finish H-2th and H-1th row, homogenity value is still valailable + int hc, vc; + for (int i=H-1; i vc) + green[i-1][j] = gh[(i-1)%4][j]; + else if (hc < vc) + green[i-1][j] = gv[(i-1)%4][j]; + else + green[i-1][j] = (gh[(i-1)%4][j] + gv[(i-1)%4][j]) / 2; + } + + freeArray2(rh, 3); + freeArray2(gh, 4); + freeArray2(bh, 3); + freeArray2(rv, 3); + freeArray2(gv, 4); + freeArray2(bv, 3); + freeArray2(lLh, 3); + freeArray2(lah, 3); + freeArray2(lbh, 3); + freeArray2(homh, 3); + freeArray2(lLv, 3); + freeArray2(lav, 3); + freeArray2(lbv, 3); + freeArray2(homv, 3); +} + +void RawImageSource::hphd_vertical (float** hpmap, int col_from, int col_to) { + + float* temp = new float[MAX(W,H)]; + float* avg = new float[MAX(W,H)]; + float* dev = new float[MAX(W,H)]; + + memset (temp, 0, MAX(W,H)*sizeof(float)); + memset (avg, 0, MAX(W,H)*sizeof(float)); + memset (dev, 0, MAX(W,H)*sizeof(float)); + + for (int k=col_from; kdata[i-5][k] - 8*ri->data[i-4][k] + 27*ri->data[i-3][k] - 48*ri->data[i-2][k] + 42*ri->data[i-1][k] - + (ri->data[i+5][k] - 8*ri->data[i+4][k] + 27*ri->data[i+3][k] - 48*ri->data[i+2][k] + 42*ri->data[i+1][k])) / 100.0; + temp[i] = ABS(temp[i]); + } + for (int j=4; jdata[i][j-5] - 8*ri->data[i][j-4] + 27*ri->data[i][j-3] - 48*ri->data[i][j-2] + 42*ri->data[i][j-1] - + (ri->data[i][j+5] - 8*ri->data[i][j+4] + 27*ri->data[i][j+3] - 48*ri->data[i][j+2] + 42*ri->data[i][j+1])) / 100; + temp[j] = ABS(temp[j]); + } + for (int j=4; jhpmap[i][j] = 2; + else if (hpv < 0.8*hpmap[i][j]) + this->hpmap[i][j] = 1; + else + this->hpmap[i][j] = 0; + } + } + delete [] temp; + delete [] avg; + delete [] dev; +} + +void RawImageSource::hphd_green () { + + #pragma omp parallel for + for (int i=3; idata[i][j]; + else { + if (this->hpmap[i][j]==1) { + int g2 = ri->data[i][j+1] + ((ri->data[i][j] - ri->data[i][j+2]) >> 1); + int g4 = ri->data[i][j-1] + ((ri->data[i][j] - ri->data[i][j-2]) >> 1); + + int dx = ri->data[i][j+1] - ri->data[i][j-1]; + int d1 = ri->data[i][j+3] - ri->data[i][j+1]; + int d2 = ri->data[i][j+2] - ri->data[i][j]; + int d3 = (ri->data[i-1][j+2] - ri->data[i-1][j]) >> 1; + int d4 = (ri->data[i+1][j+2] - ri->data[i+1][j]) >> 1; + + double e2 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = ri->data[i][j-3] - ri->data[i][j-1]; + d2 = ri->data[i][j-2] - ri->data[i][j]; + d3 = (ri->data[i-1][j-2] - ri->data[i-1][j]) >> 1; + d4 = (ri->data[i+1][j-2] - ri->data[i+1][j]) >> 1; + + double e4 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + green[i][j] = CLIP((e2 * g2 + e4 * g4) / (e2 + e4)); + } + else if (this->hpmap[i][j]==2) { + int g1 = ri->data[i-1][j] + ((ri->data[i][j] - ri->data[i-2][j]) >> 1); + int g3 = ri->data[i+1][j] + ((ri->data[i][j] - ri->data[i+2][j]) >> 1); + + int dy = ri->data[i+1][j] - ri->data[i-1][j]; + int d1 = ri->data[i-1][j] - ri->data[i-3][j]; + int d2 = ri->data[i][j] - ri->data[i-2][j]; + int d3 = (ri->data[i][j-1] - ri->data[i-2][j-1]) >> 1; + int d4 = (ri->data[i][j+1] - ri->data[i-2][j+1]) >> 1; + + double e1 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = ri->data[i+1][j] - ri->data[i+3][j]; + d2 = ri->data[i][j] - ri->data[i+2][j]; + d3 = (ri->data[i][j-1] - ri->data[i+2][j-1]) >> 1; + d4 = (ri->data[i][j+1] - ri->data[i+2][j+1]) >> 1; + + double e3 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + green[i][j] = CLIP((e1 * g1 + e3 * g3) / (e1 + e3)); + } + else { + int g1 = ri->data[i-1][j] + ((ri->data[i][j] - ri->data[i-2][j]) >> 1); + int g2 = ri->data[i][j+1] + ((ri->data[i][j] - ri->data[i][j+2]) >> 1); + int g3 = ri->data[i+1][j] + ((ri->data[i][j] - ri->data[i+2][j]) >> 1); + int g4 = ri->data[i][j-1] + ((ri->data[i][j] - ri->data[i][j-2]) >> 1); + + int dx = ri->data[i][j+1] - ri->data[i][j-1]; + int dy = ri->data[i+1][j] - ri->data[i-1][j]; + + int d1 = ri->data[i-1][j] - ri->data[i-3][j]; + int d2 = ri->data[i][j] - ri->data[i-2][j]; + int d3 = (ri->data[i][j-1] - ri->data[i-2][j-1]) >> 1; + int d4 = (ri->data[i][j+1] - ri->data[i-2][j+1]) >> 1; + + double e1 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = ri->data[i][j+3] - ri->data[i][j+1]; + d2 = ri->data[i][j+2] - ri->data[i][j]; + d3 = (ri->data[i-1][j+2] - ri->data[i-1][j]) >> 1; + d4 = (ri->data[i+1][j+2] - ri->data[i+1][j]) >> 1; + + double e2 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = ri->data[i+1][j] - ri->data[i+3][j]; + d2 = ri->data[i][j] - ri->data[i+2][j]; + d3 = (ri->data[i][j-1] - ri->data[i+2][j-1]) >> 1; + d4 = (ri->data[i][j+1] - ri->data[i+2][j+1]) >> 1; + + double e3 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + d1 = ri->data[i][j-3] - ri->data[i][j-1]; + d2 = ri->data[i][j-2] - ri->data[i][j]; + d3 = (ri->data[i-1][j-2] - ri->data[i-1][j]) >> 1; + d4 = (ri->data[i+1][j-2] - ri->data[i+1][j]) >> 1; + + double e4 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + green[i][j] = CLIP((e1*g1 + e2*g2 + e3*g3 + e4*g4) / (e1 + e2 + e3 + e4)); + } + } + } + } +} + +void RawImageSource::hphd_demosaic () { + + if (plistener) { + plistener->setProgressStr ("Demosaicing..."); + plistener->setProgress (0.0); + } + + float** hpmap = new float*[H]; + for (int i=0; isetProgress (0.33); + + this->hpmap = allocArray(W, H); + for (int i=0; ihpmap[i], 0, W*sizeof(char)); + +#ifdef _OPENMP + #pragma omp parallel + { + int tid = omp_get_thread_num(); + int nthreads = omp_get_num_threads(); + int blk = H/nthreads; + + if (tid(hpmap, H); + + if (plistener) + plistener->setProgress (0.66); + + green = new unsigned short*[H]; + for (int i=0; isetProgress (1.0); +} + +void RawImageSource::HLRecovery_Luminance (unsigned short* rin, unsigned short* gin, unsigned short* bin, unsigned short* rout, unsigned short* gout, unsigned short* bout, int width, int maxval) { + + for (int i=0; imaxval || g>maxval || b>maxval) { + int ro = MIN (r, maxval); + int go = MIN (g, maxval); + int bo = MIN (b, maxval); + double L = r + g + b; + double C = 1.732050808 * (r - g); + double H = 2 * b - r - g; + double Co = 1.732050808 * (ro - go); + double Ho = 2 * bo - ro - go; + if (r!=g && g!=b) { + double ratio = sqrt ((Co*Co+Ho*Ho) / (C*C+H*H)); + C *= ratio; + H *= ratio; + } + int rr = L / 3.0 - H / 6.0 + C / 3.464101615; + int gr = L / 3.0 - H / 6.0 - C / 3.464101615; + int br = L / 3.0 + H / 3.0; + rout[i] = CLIP(rr); + gout[i] = CLIP(gr); + bout[i] = CLIP(br); + } + else { + rout[i] = rin[i]; + gout[i] = gin[i]; + bout[i] = bin[i]; + } + } +} + +void RawImageSource::HLRecovery_CIELab (unsigned short* rin, unsigned short* gin, unsigned short* bin, unsigned short* rout, unsigned short* gout, unsigned short* bout, int width, int maxval, double cam[3][3], double icam[3][3]) { + + static bool crTableReady = false; + static double fv[0x10000]; + if (!crTableReady) { + for (int ix=0; ix < 0x10000; ix++) { + double rx = ix / 65535.0; + fv[ix] = rx > 0.008856 ? exp(1.0/3 * log(rx)) : 7.787*rx + 16/116.0; + } + crTableReady = true; + } + + for (int i=0; imaxval || g>maxval || b>maxval) { + int ro = MIN (r, maxval); + int go = MIN (g, maxval); + int bo = MIN (b, maxval); + double yy = cam[0][1]*r + cam[1][1]*g + cam[2][1]*b; + double fy = fv[CLIP((int)yy)]; + // compute LCH decompostion of the clipped pixel (only color information, thus C and H will be used) + double x = cam[0][0]*ro + cam[1][0]*go + cam[2][0]*bo; + double y = cam[0][1]*ro + cam[1][1]*go + cam[2][1]*bo; + double z = cam[0][2]*ro + cam[1][2]*go + cam[2][2]*bo; + x = fv[CLIP((int)x)]; + y = fv[CLIP((int)y)]; + z = fv[CLIP((int)z)]; + // convert back to rgb + double fz = fy - y + z; + double fx = fy + x - y; + double zr = (fz<=0.206893) ? ((116.0*fz-16.0)/903.3) : (fz * fz * fz); + double xr = (fx<=0.206893) ? ((116.0*fx-16.0)/903.3) : (fx * fx * fx); + x = xr*65535.0 - 0.5; + y = yy; + z = zr*65535.0 - 0.5; + int rr = icam[0][0]*x + icam[1][0]*y + icam[2][0]*z; + int gr = icam[0][1]*x + icam[1][1]*y + icam[2][1]*z; + int br = icam[0][2]*x + icam[1][2]*y + icam[2][2]*z; + rout[i] = CLIP(rr); + gout[i] = CLIP(gr); + bout[i] = CLIP(br); + } + else { + rout[i] = rin[i]; + gout[i] = gin[i]; + bout[i] = bin[i]; + } + } +} + +void RawImageSource::hlRecovery (std::string method, unsigned short* red, unsigned short* green, unsigned short* blue, int i, int sx1, int width, int skip) { + + if (method=="Luminance") + HLRecovery_Luminance (red, green, blue, red, green, blue, width, 65535 / ri->defgain); + else if (method=="CIELab blending") + HLRecovery_CIELab (red, green, blue, red, green, blue, width, 65535 / ri->defgain, cam, icam); + else if (method=="Color") + HLRecovery_ColorPropagation (red, green, blue, i, sx1, width, skip); +} + +int RawImageSource::getAEHistogram (unsigned int* histogram, int& histcompr) { + + histcompr = 3; + + memset (histogram, 0, (65536>>histcompr)*sizeof(int)); + + for (int i=border; iheight-border; i++) { + int start, end; + if (fuji) { + int fw = ri->fuji_width; + start = ABS(fw-i) + border; + end = MIN( ri->height+ ri->width-fw-i, fw+i) - border; + } + else { + start = border; + end = ri->width-border; + } + if (ri->filters) + for (int j=start; jdata[i][j]>>histcompr]+=2; + else + histogram[ri->data[i][j]>>histcompr]+=4; + else + for (int j=start; j<3*end; j++) { + histogram[ri->data[i][j+0]>>histcompr]++; + histogram[ri->data[i][j+1]>>histcompr]++; + histogram[ri->data[i][j+2]>>histcompr]++; + } + } + return 1; +} + +ColorTemp RawImageSource::getAutoWB () { + + double avg_r = 0; + double avg_g = 0; + double avg_b = 0; + int rn = 0, gn = 0, bn = 0; + + if (fuji) { + for (int i=32; iheight-32; i++) { + int fw = ri->fuji_width; + int start = ABS(fw-i) + 32; + int end = MIN(ri->height+ri->width-fw-i, fw+i) - 32; + for (int j=start; jfilters) { + double d = CLIP(ri->defgain*ri->data[i][3*j]); + if (d>64000) + continue; + avg_r += d*d*d*d*d*d; rn++; + d = CLIP(ri->defgain*ri->data[i][3*j+1]); + if (d>64000) + continue; + avg_g += d*d*d*d*d*d; gn++; + d = CLIP(ri->defgain*ri->data[i][3*j+2]); + if (d>64000) + continue; + avg_b += d*d*d*d*d*d; bn++; + } + else { + double d = CLIP(ri->defgain*ri->data[i][j]); + if (d>64000) + continue; + double dp = d*d*d*d*d*d; + if (ISRED(ri,i,j)) { + avg_r += dp; + rn++; + } + else if (ISGREEN(ri,i,j)) { + avg_g += dp; + gn++; + } + else if (ISBLUE(ri,i,j)) { + avg_b += dp; + bn++; + } + } + } + } + } + else { + for (int i=32; iheight-32; i++) + for (int j=32; jwidth-32; j++) { + if (!ri->filters) { + double d = CLIP(ri->defgain*ri->data[i][3*j]); + if (d>64000) + continue; + avg_r += d*d*d*d*d*d; rn++; + d = CLIP(ri->defgain*ri->data[i][3*j+1]); + if (d>64000) + continue; + avg_g += d*d*d*d*d*d; gn++; + d = CLIP(ri->defgain*ri->data[i][3*j+2]); + if (d>64000) + continue; + avg_b += d*d*d*d*d*d; bn++; + } + else { + double d = CLIP(ri->defgain*ri->data[i][j]); + if (d>64000) + continue; + double dp = d*d*d*d*d*d; + if (ISRED(ri,i,j)) { + avg_r += dp; + rn++; + } + else if (ISGREEN(ri,i,j)) { + avg_g += dp; + gn++; + } + else if (ISBLUE(ri,i,j)) { + avg_b += dp; + bn++; + } + } + } + } + + printf ("AVG: %g %g %g\n", avg_r/rn, avg_g/gn, avg_b/bn); + +// double img_r, img_g, img_b; +// wb.getMultipliers (img_r, img_g, img_b); + +// return ColorTemp (pow(avg_r/rn, 1.0/6.0)*img_r, pow(avg_g/gn, 1.0/6.0)*img_g, pow(avg_b/bn, 1.0/6.0)*img_b); + + double reds = pow (avg_r/rn, 1.0/6.0) * ri->camwb_red; + double greens = pow (avg_g/gn, 1.0/6.0) * ri->camwb_green; + double blues = pow (avg_b/bn, 1.0/6.0) * ri->camwb_blue; + + double rm = coeff[0][0]*reds + coeff[0][1]*greens + coeff[0][2]*blues; + double gm = coeff[1][0]*reds + coeff[1][1]*greens + coeff[1][2]*blues; + double bm = coeff[2][0]*reds + coeff[2][1]*greens + coeff[2][2]*blues; + + return ColorTemp (rm, gm, bm); +} + +void RawImageSource::transformPosition (int x, int y, int tran, int& ttx, int& tty) { + + tran = defTransform (tran); + + x += border; + y += border; + + if (d1x) { + if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) + x /= 2; + else + y /= 2; + } + + int w = W, h = H; + if (fuji) { + w = ri->fuji_width * 2 + 1; + h = (H - ri->fuji_width)*2 + 1; + } + int sw = w, sh = h; + if ((tran & TR_ROT) == TR_R90 || (tran & TR_ROT) == TR_R270) { + sw = h; + sh = w; + } + + int ppx = x, ppy = y; + if (tran & TR_HFLIP) + ppx = sw - 1 - x ; + if (tran & TR_VFLIP) + ppy = sh - 1 - y; + + int tx = ppx; + int ty = ppy; + + if ((tran & TR_ROT) == TR_R180) { + tx = w - 1 - ppx; + ty = h - 1 - ppy; + } + else if ((tran & TR_ROT) == TR_R90) { + tx = ppy; + ty = h - 1 - ppx; + } + else if ((tran & TR_ROT) == TR_R270) { + tx = w - 1 - ppy; + ty = ppx; + } + + if (fuji) { + ttx = (tx+ty) / 2; + tty = (ty-tx) / 2 + ri->fuji_width; + } + else { + ttx = tx; + tty = ty; + } +} + +ColorTemp RawImageSource::getSpotWB (std::vector red, std::vector green, std::vector& blue, int tran) { + + int x; int y; + int d[9][2] = {{0,0}, {-1,-1}, {-1,0}, {-1,1}, {0,-1}, {0,1}, {1,-1}, {1,0}, {1,1}}; + double reds = 0, greens = 0, blues = 0; + int rn = 0, gn = 0, bn = 0; + + if (!ri->filters) { + for (int i=0; i=0 && y>=0 && xdata[y][3*x]; + rn++; + } + transformPosition (green[i].x, green[i].y, tran, x, y); + if (x>=0 && y>=0 && xdata[y][3*x+1]; + gn++; + } + transformPosition (blue[i].x, blue[i].y, tran, x, y); + if (x>=0 && y>=0 && xdata[y][3*x+2]; + bn++; + } + } + } + else { + for (int i=0; i=0 && yv>=0 && xvdata[yv][xv]; + rn++; + break; + } + } + transformPosition (green[i].x, green[i].y, tran, x, y); + for (int k=0; k<9; k++) { + int xv = x + d[k][0]; + int yv = y + d[k][1]; + if (ISGREEN(ri,yv,xv) && xv>=0 && yv>=0 && xvdata[yv][xv]; + gn++; + break; + } + } + transformPosition (blue[i].x, blue[i].y, tran, x, y); + for (int k=0; k<9; k++) { + int xv = x + d[k][0]; + int yv = y + d[k][1]; + if (ISBLUE(ri,yv,xv) && xv>=0 && yv>=0 && xvdata[yv][xv]; + bn++; + break; + } + } + } + } + + reds = reds/rn * ri->camwb_red; + greens = greens/gn * ri->camwb_green; + blues = blues/bn * ri->camwb_blue; + + double rm = coeff[0][0]*reds + coeff[0][1]*greens + coeff[0][2]*blues; + double gm = coeff[1][0]*reds + coeff[1][1]*greens + coeff[1][2]*blues; + double bm = coeff[2][0]*reds + coeff[2][1]*greens + coeff[2][2]*blues; + + return ColorTemp (rm, gm, bm); +} + +#define FORCC for (c=0; c < colors; c++) +#define fc(row,col) \ + (ri->prefilters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3) +typedef unsigned short ushort; +void RawImageSource::vng4_demosaic () { + + static const signed char *cp, terms[] = { + -2,-2,+0,-1,0,0x01, -2,-2,+0,+0,1,0x01, -2,-1,-1,+0,0,0x01, + -2,-1,+0,-1,0,0x02, -2,-1,+0,+0,0,0x03, -2,-1,+0,+1,1,0x01, + -2,+0,+0,-1,0,0x06, -2,+0,+0,+0,1,0x02, -2,+0,+0,+1,0,0x03, + -2,+1,-1,+0,0,0x04, -2,+1,+0,-1,1,0x04, -2,+1,+0,+0,0,0x06, + -2,+1,+0,+1,0,0x02, -2,+2,+0,+0,1,0x04, -2,+2,+0,+1,0,0x04, + -1,-2,-1,+0,0,0x80, -1,-2,+0,-1,0,0x01, -1,-2,+1,-1,0,0x01, + -1,-2,+1,+0,1,0x01, -1,-1,-1,+1,0,0x88, -1,-1,+1,-2,0,0x40, + -1,-1,+1,-1,0,0x22, -1,-1,+1,+0,0,0x33, -1,-1,+1,+1,1,0x11, + -1,+0,-1,+2,0,0x08, -1,+0,+0,-1,0,0x44, -1,+0,+0,+1,0,0x11, + -1,+0,+1,-2,1,0x40, -1,+0,+1,-1,0,0x66, -1,+0,+1,+0,1,0x22, + -1,+0,+1,+1,0,0x33, -1,+0,+1,+2,1,0x10, -1,+1,+1,-1,1,0x44, + -1,+1,+1,+0,0,0x66, -1,+1,+1,+1,0,0x22, -1,+1,+1,+2,0,0x10, + -1,+2,+0,+1,0,0x04, -1,+2,+1,+0,1,0x04, -1,+2,+1,+1,0,0x04, + +0,-2,+0,+0,1,0x80, +0,-1,+0,+1,1,0x88, +0,-1,+1,-2,0,0x40, + +0,-1,+1,+0,0,0x11, +0,-1,+2,-2,0,0x40, +0,-1,+2,-1,0,0x20, + +0,-1,+2,+0,0,0x30, +0,-1,+2,+1,1,0x10, +0,+0,+0,+2,1,0x08, + +0,+0,+2,-2,1,0x40, +0,+0,+2,-1,0,0x60, +0,+0,+2,+0,1,0x20, + +0,+0,+2,+1,0,0x30, +0,+0,+2,+2,1,0x10, +0,+1,+1,+0,0,0x44, + +0,+1,+1,+2,0,0x10, +0,+1,+2,-1,1,0x40, +0,+1,+2,+0,0,0x60, + +0,+1,+2,+1,0,0x20, +0,+1,+2,+2,0,0x10, +1,-2,+1,+0,0,0x80, + +1,-1,+1,+1,0,0x88, +1,+0,+1,+2,0,0x08, +1,+0,+2,-1,0,0x40, + +1,+0,+2,+1,0,0x10 + }, chood[] = { -1,-1, -1,0, -1,+1, 0,+1, +1,+1, +1,0, +1,-1, 0,-1 }; + + if (plistener) { + plistener->setProgressStr ("Demosaicing..."); + plistener->setProgress (0.0); + } + + ushort (*brow[5])[4], *pix; + int prow=7, pcol=1, *ip, *code[16][16], gval[8], gmin, gmax, sum[4]; + int row, col, x, y, x1, x2, y1, y2, t, weight, grads, color, diag; + int g, diff, thold, num, c, width=W, height=H, colors=4; + ushort (*image)[4]; + int lcode[16][16][32], shift, i, j; + + image = (ushort (*)[4]) calloc (H*W, sizeof *image); + for (int ii=0; iidata[ii][jj]; + +// first linear interpolation + for (row=0; row < 16; row++) + for (col=0; col < 16; col++) { + ip = lcode[row][col]; + memset (sum, 0, sizeof sum); + for (y=-1; y <= 1; y++) + for (x=-1; x <= 1; x++) { + shift = (y==0) + (x==0); + if (shift == 2) continue; + color = fc(row+y,col+x); + *ip++ = (width*y + x)*4 + color; + *ip++ = shift; + *ip++ = color; + sum[color] += 1 << shift; + } + FORCC + if (c != fc(row,col)) { + *ip++ = c; + *ip++ = 256 / sum[c]; + } + } + + for (row=1; row < height-1; row++) + for (col=1; col < width-1; col++) { + pix = image[row*width+col]; + ip = lcode[row & 15][col & 15]; + memset (sum, 0, sizeof sum); + for (i=8; i--; ip+=3) + sum[ip[2]] += pix[ip[0]] << ip[1]; + for (i=colors; --i; ip+=2) + pix[ip[0]] = sum[ip[0]] * ip[1] >> 8; + } + +// lin_interpolate(); + + + ip = (int *) calloc ((prow+1)*(pcol+1), 1280); + for (row=0; row <= prow; row++) /* Precalculate for VNG */ + for (col=0; col <= pcol; col++) { + code[row][col] = ip; + for (cp=terms, t=0; t < 64; t++) { + y1 = *cp++; x1 = *cp++; + y2 = *cp++; x2 = *cp++; + weight = *cp++; + grads = *cp++; + color = fc(row+y1,col+x1); + if (fc(row+y2,col+x2) != color) continue; + diag = (fc(row,col+1) == color && fc(row+1,col) == color) ? 2:1; + if (abs(y1-y2) == diag && abs(x1-x2) == diag) continue; + *ip++ = (y1*width + x1)*4 + color; + *ip++ = (y2*width + x2)*4 + color; + *ip++ = weight; + for (g=0; g < 8; g++) + if (grads & 1< gval[g]) gmin = gval[g]; + if (gmax < gval[g]) gmax = gval[g]; + } + if (gmax == 0) { + memcpy (brow[2][col], pix, sizeof *image); + continue; + } + thold = gmin + (gmax >> 1); + memset (sum, 0, sizeof sum); + for (num=g=0; g < 8; g++,ip+=2) { /* Average the neighbors */ + if (gval[g] <= thold) { + FORCC + if (c == color && ip[1]) + sum[c] += (pix[c] + pix[ip[1]]) >> 1; + else + sum[c] += pix[ip[0] + c]; + num++; + } + } + FORCC { /* Save to buffer */ + t = pix[color]; + if (c != color) + t += (sum[c] - sum[color]) / num; + brow[2][col][c] = CLIP(t); + } + } + if (row > 3) /* Write buffer to image */ + memcpy (image[(row-2)*width+2], brow[0]+2, (width-4)*sizeof *image); + for (g=0; g < 4; g++) + brow[(g-1) & 3] = brow[g]; + if (!(row%20) && plistener) + plistener->setProgress ((double)row / (H-2)); + } + memcpy (image[(row-2)*width+2], brow[0]+2, (width-4)*sizeof *image); + memcpy (image[(row-1)*width+2], brow[1]+2, (width-4)*sizeof *image); + free (brow[4]); + free (code[0][0]); + + green = new unsigned short*[H]; + for (int i=0; i> 1; + } + free (image); +} + +//#define ABS(x) (((int)(x) ^ ((int)(x) >> 31)) - ((int)(x) >> 31)) +//#define MIN(a,b) ((a) < (b) ? (a) : (b)) +//#define MAX(a,b) ((a) > (b) ? (a) : (b)) +//#define LIM(x,min,max) MAX(min,MIN(x,max)) +//#define CLIP(x) LIM(x,0,65535) +#undef fc +#define fc(row,col) \ + (ri->filters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3) +#define FC(x,y) fc(x,y) +#define LIM(x,min,max) MAX(min,MIN(x,max)) +#define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y)) + +/* + Patterned Pixel Grouping Interpolation by Alain Desbiolles +*/ +void RawImageSource::ppg_demosaic() +{ + int width=W, height=H; + int dir[5] = { 1, width, -1, -width, 1 }; + int row, col, diff[2], guess[2], c, d, i; + ushort (*pix)[4]; + + ushort (*image)[4]; + int colors = 3; + + if (plistener) { + plistener->setProgressStr ("Demosaicing..."); + plistener->setProgress (0.0); + } + + image = (ushort (*)[4]) calloc (H*W, sizeof *image); + for (int ii=0; iidata[ii][jj]; + + border_interpolate(3, image); + +/* Fill in the green layer with gradients and pattern recognition: */ + for (row=3; row < height-3; row++) { + for (col=3+(FC(row,3) & 1), c=FC(row,col); col < width-3; col+=2) { + pix = image + row*width+col; + for (i=0; (d=dir[i]) > 0; i++) { + guess[i] = (pix[-d][1] + pix[0][c] + pix[d][1]) * 2 + - pix[-2*d][c] - pix[2*d][c]; + diff[i] = ( ABS(pix[-2*d][c] - pix[ 0][c]) + + ABS(pix[ 2*d][c] - pix[ 0][c]) + + ABS(pix[ -d][1] - pix[ d][1]) ) * 3 + + ( ABS(pix[ 3*d][1] - pix[ d][1]) + + ABS(pix[-3*d][1] - pix[-d][1]) ) * 2; + } + d = dir[i = diff[0] > diff[1]]; + pix[0][1] = ULIM(guess[i] >> 2, pix[d][1], pix[-d][1]); + } + if(plistener) plistener->setProgress(0.33*row/(height-3)); + } +/* Calculate red and blue for each green pixel: */ + for (row=1; row < height-1; row++) { + for (col=1+(FC(row,2) & 1), c=FC(row,col+1); col < width-1; col+=2) { + pix = image + row*width+col; + for (i=0; (d=dir[i]) > 0; c=2-c, i++) + pix[0][c] = CLIP((pix[-d][c] + pix[d][c] + 2*pix[0][1] + - pix[-d][1] - pix[d][1]) >> 1); + } + if(plistener) plistener->setProgress(0.33 + 0.33*row/(height-1)); + } +/* Calculate blue for red pixels and vice versa: */ + for (row=1; row < height-1; row++) { + for (col=1+(FC(row,1) & 1), c=2-FC(row,col); col < width-1; col+=2) { + pix = image + row*width+col; + for (i=0; (d=dir[i]+dir[i+1]) > 0; i++) { + diff[i] = ABS(pix[-d][c] - pix[d][c]) + + ABS(pix[-d][1] - pix[0][1]) + + ABS(pix[ d][1] - pix[0][1]); + guess[i] = pix[-d][c] + pix[d][c] + 2*pix[0][1] + - pix[-d][1] - pix[d][1]; + } + if (diff[0] != diff[1]) + pix[0][c] = CLIP(guess[diff[0] > diff[1]] >> 1); + else + pix[0][c] = CLIP((guess[0]+guess[1]) >> 2); + } + if(plistener) plistener->setProgress(0.67 + 0.33*row/(height-1)); + } + + red = new unsigned short*[H]; + for (int i=0; i= border && row < height-border) + col = width-border; + memset (sum, 0, sizeof sum); + for (y=row-1; y != row+2; y++) + for (x=col-1; x != col+2; x++) + if (y < height && x < width) { + f = fc(y,x); + sum[f] += image[y*width+x][f]; + sum[f+4]++; + } + f = fc(row,col); + FORCC if (c != f && sum[c+4]) + image[row*width+col][c] = sum[c] / sum[c+4]; + } +} + +void RawImageSource::bilinear_interpolate_block(ushort (*image)[4], int start, int end) +{ + ushort (*pix); + int i, *ip, sum[4]; + int width=W; + int colors = 3; + + for (int row = start; row < end; row++) + for (int col=1; col < width-1; col++) { + pix = image[row*width+col]; + ip = blcode[row & 15][col & 15]; + memset (sum, 0, sizeof sum); + for (i=8; i--; ip+=3) + sum[ip[2]] += pix[ip[0]] << ip[1]; + for (i=colors; --i; ip+=2) + pix[ip[0]] = sum[ip[0]] * ip[1] >> 8; + } + + +} + +void RawImageSource::bilinear_demosaic() +{ + int width=W, height=H; + int *ip, sum[4]; + int c, x, y, row, col, shift, color; + int colors = 3; + + ushort (*image)[4], *pix; + image = (ushort (*)[4]) calloc (H*W, sizeof *image); + + for (int ii=0; iidata[ii][jj]; + + //if (verbose) fprintf (stderr,_("Bilinear interpolation...\n")); + if (plistener) { + plistener->setProgressStr ("Demosaicing..."); + plistener->setProgress (0.0); + } + + memset(blcode,0,16*16*32); + for (row=0; row < 16; row++) + for (col=0; col < 16; col++) { + ip = blcode[row][col]; + memset (sum, 0, sizeof sum); + for (y=-1; y <= 1; y++) + for (x=-1; x <= 1; x++) { + shift = (y==0) + (x==0); + if (shift == 2) continue; + color = fc(row+y,col+x); + *ip++ = (width*y + x)*4 + color; + *ip++ = shift; + *ip++ = color; + sum[color] += 1 << shift; + } + FORCC + if (c != fc(row,col)) { + *ip++ = c; + *ip++ = 256 / sum[c]; + } + } + +#ifdef _OPENMP + #pragma omp parallel + { + int tid = omp_get_thread_num(); + int nthreads = omp_get_num_threads(); + int blk = H/nthreads; + + int start = 0; + if (tid == 0) start = 1; + if (tidsetProgress (1.0); + free (image); +} + +/* + Adaptive Homogeneity-Directed interpolation is based on + the work of Keigo Hirakawa, Thomas Parks, and Paul Lee. + */ +#define TS 256 /* Tile Size */ +#define FORC(cnt) for (c=0; c < cnt; c++) +#define FORC3 FORC(3) +#define SQR(x) ((x)*(x)) + +void RawImageSource::ahd_demosaic() +{ + int i, j, k, top, left, row, col, tr, tc, c, d, val, hm[2]; + ushort (*pix)[4], (*rix)[3]; + static const int dir[4] = { -1, 1, -TS, TS }; + unsigned ldiff[2][4], abdiff[2][4], leps, abeps; + float r, cbrt[0x10000], xyz[3], xyz_cam[3][4]; + ushort (*rgb)[TS][TS][3]; + short (*lab)[TS][TS][3], (*lix)[3]; + char (*homo)[TS][TS], *buffer; + + int width=W, height=H; + ushort (*image)[4]; + int colors = 3; + + const double xyz_rgb[3][3] = { /* XYZ from RGB */ + { 0.412453, 0.357580, 0.180423 }, + { 0.212671, 0.715160, 0.072169 }, + { 0.019334, 0.119193, 0.950227 } + }; + + const float d65_white[3] = { 0.950456, 1, 1.088754 }; + + if (plistener) { + plistener->setProgressStr ("Demosaicing..."); + plistener->setProgress (0.0); + } + + image = (ushort (*)[4]) calloc (H*W, sizeof *image); + for (int ii=0; iidata[ii][jj]; + + for (i=0; i < 0x10000; i++) { + r = i / 65535.0; + cbrt[i] = r > 0.008856 ? pow(r,1/3.0) : 7.787*r + 16/116.0; + } + + for (i=0; i < 3; i++) + for (j=0; j < colors; j++) + for (xyz_cam[i][j] = k=0; k < 3; k++) + xyz_cam[i][j] += xyz_rgb[i][k] * coeff[k][j] / d65_white[i]; + + border_interpolate(5, image); + buffer = (char *) malloc (26*TS*TS); /* 1664 kB */ + //merror (buffer, "ahd_interpolate()"); + rgb = (ushort(*)[TS][TS][3]) buffer; + lab = (short (*)[TS][TS][3])(buffer + 12*TS*TS); + homo = (char (*)[TS][TS]) (buffer + 24*TS*TS); + + // helper variables for progress indication + int n_tiles = ((height-7 + (TS-7))/(TS-6)) * ((width-7 + (TS-7))/(TS-6)); + int tile = 0; + + for (top=2; top < height-5; top += TS-6) + for (left=2; left < width-5; left += TS-6) { + + /* Interpolate green horizontally and vertically: */ + for (row = top; row < top+TS && row < height-2; row++) { + col = left + (FC(row,left) & 1); + for (c = FC(row,col); col < left+TS && col < width-2; col+=2) { + pix = image + row*width+col; + val = ((pix[-1][1] + pix[0][c] + pix[1][1]) * 2 + - pix[-2][c] - pix[2][c]) >> 2; + rgb[0][row-top][col-left][1] = ULIM(val,pix[-1][1],pix[1][1]); + val = ((pix[-width][1] + pix[0][c] + pix[width][1]) * 2 + - pix[-2*width][c] - pix[2*width][c]) >> 2; + rgb[1][row-top][col-left][1] = ULIM(val,pix[-width][1],pix[width][1]); + } + } + + /* Interpolate red and blue, and convert to CIELab: */ + for (d=0; d < 2; d++) + for (row=top+1; row < top+TS-1 && row < height-3; row++) + for (col=left+1; col < left+TS-1 && col < width-3; col++) { + pix = image + row*width+col; + rix = &rgb[d][row-top][col-left]; + lix = &lab[d][row-top][col-left]; + if ((c = 2 - FC(row,col)) == 1) { + c = FC(row+1,col); + val = pix[0][1] + (( pix[-1][2-c] + pix[1][2-c] + - rix[-1][1] - rix[1][1] ) >> 1); + rix[0][2-c] = CLIP(val); + val = pix[0][1] + (( pix[-width][c] + pix[width][c] + - rix[-TS][1] - rix[TS][1] ) >> 1); + } else + val = rix[0][1] + (( pix[-width-1][c] + pix[-width+1][c] + + pix[+width-1][c] + pix[+width+1][c] + - rix[-TS-1][1] - rix[-TS+1][1] + - rix[+TS-1][1] - rix[+TS+1][1] + 1) >> 2); + rix[0][c] = CLIP(val); + c = FC(row,col); + rix[0][c] = pix[0][c]; + xyz[0] = xyz[1] = xyz[2] = 0.5; + FORCC { + xyz[0] += xyz_cam[0][c] * rix[0][c]; + xyz[1] += xyz_cam[1][c] * rix[0][c]; + xyz[2] += xyz_cam[2][c] * rix[0][c]; + } + xyz[0] = cbrt[CLIP((int) xyz[0])]; + xyz[1] = cbrt[CLIP((int) xyz[1])]; + xyz[2] = cbrt[CLIP((int) xyz[2])]; + lix[0][0] = 64 * (116 * xyz[1] - 16); + lix[0][1] = 64 * 500 * (xyz[0] - xyz[1]); + lix[0][2] = 64 * 200 * (xyz[1] - xyz[2]); + } + + /* Build homogeneity maps from the CIELab images: */ + memset (homo, 0, 2*TS*TS); + for (row=top+2; row < top+TS-2 && row < height-4; row++) { + tr = row-top; + for (col=left+2; col < left+TS-2 && col < width-4; col++) { + tc = col-left; + for (d=0; d < 2; d++) { + lix = &lab[d][tr][tc]; + for (i=0; i < 4; i++) { + ldiff[d][i] = ABS(lix[0][0]-lix[dir[i]][0]); + abdiff[d][i] = SQR(lix[0][1]-lix[dir[i]][1]) + + SQR(lix[0][2]-lix[dir[i]][2]); + } + } + leps = MIN(MAX(ldiff[0][0],ldiff[0][1]), + MAX(ldiff[1][2],ldiff[1][3])); + abeps = MIN(MAX(abdiff[0][0],abdiff[0][1]), + MAX(abdiff[1][2],abdiff[1][3])); + for (d=0; d < 2; d++) + for (i=0; i < 4; i++) + if (ldiff[d][i] <= leps && abdiff[d][i] <= abeps) + homo[d][tr][tc]++; + } + } + + /* Combine the most homogenous pixels for the final result: */ + for (row=top+3; row < top+TS-3 && row < height-5; row++) { + tr = row-top; + for (col=left+3; col < left+TS-3 && col < width-5; col++) { + tc = col-left; + for (d=0; d < 2; d++) + for (hm[d]=0, i=tr-1; i <= tr+1; i++) + for (j=tc-1; j <= tc+1; j++) + hm[d] += homo[d][i][j]; + if (hm[0] != hm[1]) + FORC3 image[row*width+col][c] = rgb[hm[1] > hm[0]][tr][tc][c]; + else + FORC3 image[row*width+col][c] = + (rgb[0][tr][tc][c] + rgb[1][tr][tc][c]) >> 1; + } + } + + tile++; + if(plistener) { + plistener->setProgress((double)tile / n_tiles); + } + } + + if(plistener) plistener->setProgress (1.0); + free (buffer); + red = new unsigned short*[H]; + for (int i=0; i= H-border) rowMax = TILEBORDER+H-border-y0; + if( x0+TILESIZE+TILEBORDER >= W-border) colMax = TILEBORDER+W-border-x0; +} + +void RawImageSource::fill_raw( ushort (*cache )[4], int x0, int y0, ushort** rawData) +{ + int rowMin,colMin,rowMax,colMax; + dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,0); + + for (int row=rowMin,y=y0-TILEBORDER+rowMin; row= border && col < W - border && row >= border && row < H - border){ + col = W - border; + if(col >= x0+TILESIZE+TILEBORDER ) + break; + } + memset(sum, 0, sizeof sum); + for (y = row - 1; y != row + 2; y++) + for (x = col - 1; x != col + 2; x++) + if (y < H && y< y0+TILESIZE+TILEBORDER && x < W && x ( pix[-4] + pix[+4] + pix[-u] + pix[+u])/4 ) + image[indx][3] = ((MIN( pix[-4], pix[+4]) + pix[-4] + pix[+4] ) < (MIN( pix[-u], pix[+u]) + pix[-u] + pix[+u])); + else + image[indx][3] = ((MAX( pix[-4], pix[+4]) + pix[-4] + pix[+4] ) > (MAX( pix[-u], pix[+u]) + pix[-u] + pix[+u])); + } + } +} + + +// interpolated green pixels are corrected using the map +void RawImageSource::dcb_correction(ushort (*image)[4], int x0, int y0) +{ + const int u=CACHESIZE, v=2*CACHESIZE; + int rowMin,colMin,rowMax,colMax; + dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,2); + + for (int row=rowMin; row < rowMax; row++) { + for (int col = colMin+(FC(y0-TILEBORDER+row,x0-TILEBORDER+colMin)&1),indx=row*CACHESIZE+col,c=FC(y0-TILEBORDER+row,x0-TILEBORDER+col); col < colMax; col+=2, indx+=2) { + int current = 4*image[indx][3] + + 2*(image[indx+u][3] + image[indx-u][3] + image[indx+1][3] + image[indx-1][3]) + + image[indx+v][3] + image[indx-v][3] + image[indx+2][3] + image[indx-2][3]; + image[indx][1] = ((16-current)*(image[indx-1][1] + image[indx+1][1])/2 + current*(image[indx-u][1] + image[indx+u][1])/2)/16; + } + } +} + +// R and B smoothing using green contrast, all pixels except 2 pixel wide border +void RawImageSource::dcb_pp(ushort (*image)[4], int x0, int y0) +{ + const int u=CACHESIZE; + int rowMin,colMin,rowMax,colMax; + dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,2); + + for (int row=rowMin; row < rowMax; row++) + for (int col=colMin, indx=row*CACHESIZE+col; col < colMax; col++, indx++) { + //int 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; + //int 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; + //int 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; + ushort (*pix)[4] = image+(indx-u-1); + int r1 = (*pix)[0]; + int g1 = (*pix)[1]; + int b1 = (*pix)[2]; + pix++; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix++; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix+=CACHESIZE-2; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix+=2; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix+=CACHESIZE-2; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix++; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + pix++; + r1 += (*pix)[0]; + g1 += (*pix)[1]; + b1 += (*pix)[2]; + r1 /=8; + g1 /=8; + b1 /=8; + r1 = r1 + ( image[indx][1] - g1 ); + b1 = b1 + ( image[indx][1] - g1 ); + image[indx][0] = CLIP(r1); + image[indx][2] = CLIP(b1); + } +} + +// interpolated green pixels are corrected using the map +// with correction +void RawImageSource::dcb_correction2(ushort (*image)[4], int x0, int y0) +{ + const int u=CACHESIZE, v=2*CACHESIZE; + int rowMin,colMin,rowMax,colMax; + dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,4); + + for (int row=rowMin; row < rowMax; row++) { + for (int col = colMin+(FC(y0-TILEBORDER+row,x0-TILEBORDER+colMin)&1),indx=row*CACHESIZE+col,c=FC(y0-TILEBORDER+row,x0-TILEBORDER+col); col < colMax; col+=2, indx+=2) { + int current = 4*image[indx][3] + + 2*(image[indx+u][3] + image[indx-u][3] + image[indx+1][3] + image[indx-1][3]) + + image[indx+v][3] + image[indx-v][3] + image[indx+2][3] + image[indx-2][3]; + current = ((16-current)*((image[indx-1][1] + image[indx+1][1])/2 + image[indx][c] - (image[indx+2][c] + image[indx-2][c])/2) + current*((image[indx-u][1] + image[indx+u][1])/2 + image[indx][c] - (image[indx+v][c] + image[indx-v][c])/2))/16; + image[indx][1] = CLIP(current); + } + } +} + +// image refinement +void RawImageSource::dcb_refinement(ushort (*image)[4], int x0, int y0) +{ + const int u=CACHESIZE, v=2*CACHESIZE, w=3*CACHESIZE; + int rowMin,colMin,rowMax,colMax; + dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,4); + + float f[5],g1,g2; + + for (int row=rowMin; row < rowMax; row++) + for (int col=colMin+(FC(y0-TILEBORDER+row,x0-TILEBORDER+colMin)&1),indx=row*CACHESIZE+col,c=FC(y0-TILEBORDER+row,x0-TILEBORDER+col); col < colMax; col+=2,indx+=2){ + int current = 4*image[indx][3] + + 2*(image[indx+u][3] + image[indx-u][3] + image[indx+1][3] + image[indx-1][3]) + +image[indx+v][3] + image[indx-v][3] + image[indx-2][3] + image[indx+2][3]; + + f[0] = (float)(image[indx-u][1] + image[indx+u][1])/(2 + 2*image[indx][c]); + f[1] = 2*(float)image[indx-u][1]/(2 + image[indx-v][c] + image[indx][c]); + f[2] = (float)(image[indx-u][1] + image[indx-w][1])/(2 + 2*image[indx-v][c]); + f[3] = 2*(float)image[indx+u][1]/(2 + image[indx+v][c] + image[indx][c]); + f[4] = (float)(image[indx+u][1] + image[indx+w][1])/(2 + 2*image[indx+v][c]); + + g1 = (f[0] + f[1] + f[2] + f[3] + f[4] - MAX(f[1], MAX(f[2], MAX(f[3], f[4]))) - MIN(f[1], MIN(f[2], MIN(f[3], f[4]))))/3.0; + + f[0] = (float)(image[indx-1][1] + image[indx+1][1])/(2 + 2*image[indx][c]); + f[1] = 2*(float)image[indx-1][1]/(2 + image[indx-2][c] + image[indx][c]); + f[2] = (float)(image[indx-1][1] + image[indx-3][1])/(2 + 2*image[indx-2][c]); + f[3] = 2*(float)image[indx+1][1]/(2 + image[indx+2][c] + image[indx][c]); + f[4] = (float)(image[indx+1][1] + image[indx+3][1])/(2 + 2*image[indx+2][c]); + + g2 = (f[0] + f[1] + f[2] + f[3] + f[4] - MAX(f[1], MAX(f[2], MAX(f[3], f[4]))) - MIN(f[1], MIN(f[2], MIN(f[3], f[4]))))/3.0; + + image[indx][1] = CLIP((2+image[indx][c])*(current*g1 + (16-current)*g2)/16.0); + +// get rid of the overshooted pixels + int 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]))))))); + int 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]))))))); + + image[indx][1] = LIM(image[indx][1], min, max); + + + } +} + +// missing colors are interpolated using high quality algorithm by Luis Sanz Rodríguez +void RawImageSource::dcb_color_full(ushort (*image)[4], int x0, int y0, float (*chroma)[2]) +{ + const int u=CACHESIZE, v=2*CACHESIZE, w=3*CACHESIZE; + int rowMin,colMin,rowMax,colMax; + dcb_initTileLimits(colMin,rowMin,colMax,rowMax,x0,y0,3); + + int i,j; + float f[4],g[4]; + + for (int row=1; row < CACHESIZE-1; row++) + for (int col=1+(FC(y0-TILEBORDER+row,x0-TILEBORDER+1)&1),indx=row*CACHESIZE+col,c=FC(y0-TILEBORDER+row,x0-TILEBORDER+col),d=c/2; col < CACHESIZE-1; col+=2,indx+=2) + chroma[indx][d]=image[indx][c]-image[indx][1]; + + for (int row=rowMin; rowsetProgressStr ("DCB Demosaicing..."); + plistener->setProgress (currentProgress); + } + + int wTiles = W/TILESIZE + (W%TILESIZE?1:0); + int hTiles = H/TILESIZE + (H%TILESIZE?1:0); + int numTiles = wTiles * hTiles; + int tilesDone=0; +#ifdef _OPENMP + int nthreads = omp_get_max_threads(); + ushort (**image)[4] = (ushort(**)[4]) calloc( nthreads,sizeof( void*) ); + ushort (**image2)[3] = (ushort(**)[3]) calloc( nthreads,sizeof( void*) ); + ushort (**image3)[3] = (ushort(**)[3]) calloc( nthreads,sizeof( void*) ); + float (**chroma)[2] = (float (**)[2]) calloc( nthreads,sizeof( void*) ); + for(int i=0; idata ); + if( !xTile || !yTile || xTile==wTiles-1 || yTile==hTiles-1) + fill_border(tile,6, x0, y0); + dcb_hid(tile,buffer,buffer2,x0,y0); + copy_to_buffer(buffer, tile); + for (int i=iterations; i>0;i--) { + dcb_hid2(tile,x0,y0); + dcb_hid2(tile,x0,y0); + dcb_hid2(tile,x0,y0); + dcb_map(tile,x0,y0); + dcb_correction(tile,x0,y0); + } + dcb_color(tile,x0,y0); + dcb_pp(tile,x0,y0); + dcb_map(tile,x0,y0); + dcb_correction2(tile,x0,y0); + dcb_map(tile,x0,y0); + dcb_correction(tile,x0,y0); + dcb_color(tile,x0,y0); + dcb_map(tile,x0,y0); + dcb_correction(tile,x0,y0); + dcb_map(tile,x0,y0); + dcb_correction(tile,x0,y0); + dcb_map(tile,x0,y0); + restore_from_buffer(tile, buffer); + dcb_color(tile,x0,y0); + if (dcb_enhance) { + dcb_refinement(tile,x0,y0); + dcb_color_full(tile,x0,y0,chrm); + } + + for(int y=0;y currentProgress){ + currentProgress+=0.1; // Show progress each 10% + plistener->setProgress (currentProgress); + } + } +#pragma omp atomic + tilesDone++; + } + +#ifdef _OPENMP + for(int i=0; isetProgress (1.0); +} +#undef TILEBORDER +#undef TILESIZE +#undef CACHESIZE + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +//Emil's code for AMaZE +#include "amaze_interpolate_RT.cc"//AMaZE demosaic +#include "CA_correct_RT.cc"//Emil's CA auto correction +#include "cfa_linedn_RT.cc"//Emil's CA auto correction +#include "green_equil_RT.cc"//Emil's green channel equilibration +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +} /* namespace */ + diff --git a/rtengine/simpleprocess.cc b/rtengine/simpleprocess.cc index b64238e3b..de7ed1584 100644 --- a/rtengine/simpleprocess.cc +++ b/rtengine/simpleprocess.cc @@ -149,7 +149,8 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p // luminance processing CurveFactory::complexCurve (0.0, 0.0, 0.0, 0.0, params.lumaCurve.brightness, params.lumaCurve.contrast, 0.0, 0.0, false, params.lumaCurve.curve, hist16, curve, NULL); ipf.luminanceCurve (labView, labView, curve, 0, fh); - ipf.lumadenoise (labView, buffer); + ipf.impulsedenoise (labView); + ipf.lumadenoise (labView, buffer); ipf.sharpening (labView, (unsigned short**)buffer); delete [] curve; @@ -158,6 +159,7 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p // color processing ipf.colorCurve (labView, labView); ipf.colordenoise (labView, buffer); + ipf.dirpyrdenoise (labView); // wavelet equalizer ipf.waveletEqualizer (labView, true, true); diff --git a/rtgui/cropwindow.cc b/rtgui/cropwindow.cc index 206758e9b..e42834067 100644 --- a/rtgui/cropwindow.cc +++ b/rtgui/cropwindow.cc @@ -31,13 +31,13 @@ struct ZoomStep { int czoom; }; -ZoomStep zoomSteps[] = {{"10%", 0.1, 10}, +ZoomStep zoomSteps[] = {{" 10%", 0.1, 10}, {"12.5%", 0.125, 8}, {"16.6%", 1.0/6.0, 6}, - {"20%", 0.2, 5}, - {"25%", 0.25, 4}, - {"33%", 1.0/3.0, 3}, - {"50%", 0.5, 2}, + {" 20%", 0.2, 5}, + {" 25%", 0.25, 4}, + {" 33%", 1.0/3.0, 3}, + {" 50%", 0.5, 2}, {"100%", 1.0, 1000}, {"200%", 2.0, 2000}, {"300%", 3.0, 3000}, diff --git a/rtgui/zoompanel.cc b/rtgui/zoompanel.cc index df70ca6c2..29eef8fa0 100644 --- a/rtgui/zoompanel.cc +++ b/rtgui/zoompanel.cc @@ -111,7 +111,11 @@ void ZoomPanel::refreshZoomLabel () { if (iarea->mainCropWindow) { int z = (int)(iarea->mainCropWindow->getZoom () * 100); - zoomLabel->set_text (Glib::ustring::compose("%1%%", z)); + if (z<100) { + zoomLabel->set_text (Glib::ustring::compose(" %1%%", z)); + } else { + zoomLabel->set_text (Glib::ustring::compose("%1%%", z)); + } } }