From 33708560ce834f18d67f427dbd7e84274d6fb7ff Mon Sep 17 00:00:00 2001 From: Ingo Date: Tue, 17 Dec 2013 11:51:28 +0100 Subject: [PATCH] Improvement to the raw Auto CA correction, Issue 2128 --- rtengine/CA_correct_RT.cc | 378 ++++++++++++++++++++++---------------- rtengine/rawimagesource.h | 6 +- 2 files changed, 219 insertions(+), 165 deletions(-) diff --git a/rtengine/CA_correct_RT.cc b/rtengine/CA_correct_RT.cc index c79db2ea7..d43fc41bb 100644 --- a/rtengine/CA_correct_RT.cc +++ b/rtengine/CA_correct_RT.cc @@ -30,7 +30,7 @@ using namespace std; using namespace rtengine; -int RawImageSource::LinEqSolve(int nDim, float* pfMatr, float* pfVect, float* pfSolution) +int RawImageSource::LinEqSolve(int nDim, double* pfMatr, double* pfVect, double* pfSolution) { //============================================================================== // return 1 if system not solving, 0 if system solved @@ -45,11 +45,11 @@ int RawImageSource::LinEqSolve(int nDim, float* pfMatr, float* pfVect, float* pf // //============================================================================== - float fMaxElem; - float fAcc; - + double fMaxElem; + double fAcc; + int i, j, k, m; - + for(k=0; k<(nDim-1); k++) {// base row of matrix // search of line with max element fMaxElem = fabsf( pfMatr[k*nDim + k] ); @@ -60,7 +60,7 @@ int RawImageSource::LinEqSolve(int nDim, float* pfMatr, float* pfVect, float* pf m = i; } } - + // permutation of base line (index k) and max element line(index m) if(m != k) { for(i=k; i=0; k--) { pfSolution[k] = pfVect[k]; for(i=(k+1); i(b)) {temp=(a);(a)=(b);(b)=temp;} } volatile double progress = 0.0; if(plistener) plistener->setProgress (progress); - + bool autoCA = (cared==0 && cablue==0); // local variables int width=W, height=H; //temporary array to store simple interpolation of G @@ -130,7 +130,8 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { //block CA shift values and weight assigned to block char *buffer1; // vblsz*hblsz*(3*2+1) float (*blockwt); // vblsz*hblsz - float (*blockshifts)[3][2]; // vblsz*hblsz*3*2 + float (*blockshifts)[3][2]; // vblsz*hblsz*3*2 + const int border=8; const int border2=16; @@ -149,16 +150,21 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { // block CA shifts blockwt = (float (*)) (buffer1); blockshifts = (float (*)[3][2]) (buffer1+(vblsz*hblsz*sizeof(float))); - - float polymat[3][2][256], shiftmat[3][2][16], fitparams[3][2][16]; -#pragma omp parallel shared(Gtmp,width,height,blockave,blocksqave,blockdenom,blockvar,blockwt,blockshifts,polymat,shiftmat,fitparams) -{ + double polymat[3][2][256], shiftmat[3][2][16], fitparams[3][2][16]; + for (int i=0; i<256; i++) {polymat[0][0][i] = polymat[0][1][i] = polymat[2][0][i] = polymat[2][1][i] = 0;} + for (int i=0; i<16; i++) {shiftmat[0][0][i] = shiftmat[0][1][i] = shiftmat[2][0][i] = shiftmat[2][1][i] = 0;} + //order of 2d polynomial fit (polyord), and numpar=polyord^2 int polyord=4, numpar=16; - //number of blocks used in the fit int numblox[3]={0,0,0}; - + +#pragma omp parallel shared(Gtmp,width,height,blockave,blocksqave,blockdenom,blockvar,blockwt,blockshifts,polymat,shiftmat,fitparams,polyord,numpar) +{ + int progresscounter = 0; + //number of blocks used in the fit + int numbloxthr[3]={0,0,0}; + int rrmin, rrmax, ccmin, ccmax; int top, left, row, col; int rr, cc, c, indx, indx1, i, j, k, m, n, dir; @@ -176,9 +182,9 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { int res; //shifts to location of vertical and diagonal neighbors const int v1=TS, v2=2*TS, v3=3*TS, v4=4*TS;//, p1=-TS+1, p2=-2*TS+2, p3=-3*TS+3, m1=TS+1, m2=2*TS+2, m3=3*TS+3; - + float eps=1e-5f, eps2=1e-10f; //tolerance to avoid dividing by zero - + //adaptive weights for green interpolation float wtu, wtd, wtl, wtr; //local quadratic fit to shift data within a tile @@ -201,15 +207,15 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { //low and high pass 1D filters of G in vertical/horizontal directions float glpfh, glpfv; - + //max allowed CA shift const float bslim = 3.99; //gaussians for low pass filtering of G and R/B //static const float gaussg[5] = {0.171582, 0.15839, 0.124594, 0.083518, 0.0477063};//sig=2.5 //static const float gaussrb[3] = {0.332406, 0.241376, 0.0924212};//sig=1.25 - + //block CA shift values and weight assigned to block - + char *buffer; // TS*TS*16 //rgb data in a tile float* rgb[3]; @@ -230,38 +236,38 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { //low pass filter for color differences in vertical direction float (*grblpfv); // TS*TS*4 - + /* assign working space; this would not be necessary if the algorithm is part of the larger pre-interpolation processing */ buffer = (char *) malloc(3*sizeof(float)*TS*TS + 8*sizeof(float)*TS*TSH + 10*64 + 64); //merror(buffer,"CA_correct()"); memset(buffer,0,3*sizeof(float)*TS*TS + 8*sizeof(float)*TS*TSH + 10*64 + 64); - + char *data; data = buffer; // buffers aligned to size of cacheline // data = (char*)( ( uintptr_t(buffer) + uintptr_t(63)) / 64 * 64); - - + + // shift the beginning of all arrays but the first by 64 bytes to avoid cache miss conflicts on CPUs which have <=4-way associative L1-Cache rgb[0] = (float (*)) data; rgb[1] = (float (*)) (data + 1*sizeof(float)*TS*TS + 1*64); rgb[2] = (float (*)) (data + 2*sizeof(float)*TS*TS + 2*64); - grbdiff = (float (*)) (data + 3*sizeof(float)*TS*TS + 3*64); - gshift = (float (*)) (data + 3*sizeof(float)*TS*TS + sizeof(float)*TS*TSH + 4*64); - rbhpfh = (float (*)) (data + 4*sizeof(float)*TS*TS + 5*64); - rbhpfv = (float (*)) (data + 4*sizeof(float)*TS*TS + sizeof(float)*TS*TSH + 6*64); - rblpfh = (float (*)) (data + 5*sizeof(float)*TS*TS + 7*64); - rblpfv = (float (*)) (data + 5*sizeof(float)*TS*TS + sizeof(float)*TS*TSH + 8*64); - grblpfh = (float (*)) (data + 6*sizeof(float)*TS*TS + 9*64); - grblpfv = (float (*)) (data + 6*sizeof(float)*TS*TS + sizeof(float)*TS*TSH + 10*64); - - - if (cared==0 && cablue==0) { + grbdiff = (float (*)) (data + 3*sizeof(float)*TS*TS + 3*64); + gshift = (float (*)) (data + 3*sizeof(float)*TS*TS + sizeof(float)*TS*TSH + 4*64); + rbhpfh = (float (*)) (data + 4*sizeof(float)*TS*TS + 5*64); + rbhpfv = (float (*)) (data + 4*sizeof(float)*TS*TS + sizeof(float)*TS*TSH + 6*64); + rblpfh = (float (*)) (data + 5*sizeof(float)*TS*TS + 7*64); + rblpfv = (float (*)) (data + 5*sizeof(float)*TS*TS + sizeof(float)*TS*TSH + 8*64); + grblpfh = (float (*)) (data + 6*sizeof(float)*TS*TS + 9*64); + grblpfv = (float (*)) (data + 6*sizeof(float)*TS*TS + sizeof(float)*TS*TSH + 10*64); + + + if (autoCA) { // Main algorithm: Tile loop #pragma omp for collapse(2) schedule(dynamic) nowait - for (top=-border ; top < height; top += TS-border2) + for (top=-border ; top < height; top += TS-border2) for (left=-border; left < width; left += TS-border2) { vblock = ((top+border)/(TS-border2))+1; hblock = ((left+border)/(TS-border2))+1; @@ -276,9 +282,9 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { if (right>width) {ccmax=width-left;} else {ccmax=cc1;} // rgb from input CFA data - // rgb values should be floating point number between 0 and 1 - // after white balance multipliers are applied - + // rgb values should be floating point number between 0 and 1 + // after white balance multipliers are applied + for (rr=rrmin; rr < rrmax; rr++) for (row=rr+top, cc=ccmin; cc < ccmax; cc++) { col = cc+left; @@ -288,18 +294,18 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { rgb[c][indx1] = (rawData[row][col])/65535.0f; //rgb[indx1][c] = image[indx][c]/65535.0f;//for dcraw implementation } - + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //fill borders if (rrmin>0) { - for (rr=0; rr0) { - for (rr=rrmin; rr0 && ccmin>0) { - for (rr=0; rr0 && ccmax0) { - for (rr=0; rr-1 && row-1 && col>1] = fabsf(fabsf((rgb[1][indx]-rgb[c][indx])-(rgb[1][indx+v4]-rgb[c][indx+v4])) + fabsf((rgb[1][indx-v4]-rgb[c][indx-v4])-(rgb[1][indx]-rgb[c][indx])) - fabsf((rgb[1][indx-v4]-rgb[c][indx-v4])-(rgb[1][indx+v4]-rgb[c][indx+v4]))); rbhpfh[indx>>1] = fabsf(fabsf((rgb[1][indx]-rgb[c][indx])-(rgb[1][indx+4]-rgb[c][indx+4])) + fabsf((rgb[1][indx-4]-rgb[c][indx-4])-(rgb[1][indx]-rgb[c][indx])) - fabsf((rgb[1][indx-4]-rgb[c][indx-4])-(rgb[1][indx+4]-rgb[c][indx+4]))); - + /*ghpfv = fabsf(fabsf(rgb[indx][1]-rgb[indx+v4][1])+fabsf(rgb[indx][1]-rgb[indx-v4][1]) - fabsf(rgb[indx+v4][1]-rgb[indx-v4][1])); ghpfh = fabsf(fabsf(rgb[indx][1]-rgb[indx+4][1])+fabsf(rgb[indx][1]-rgb[indx-4][1]) - @@ -408,7 +414,7 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { fabsf(rgb[indx+v4][c]-rgb[indx-v4][c]))); rbhpfh[indx] = fabsf(ghpfh - fabsf(fabsf(rgb[indx][c]-rgb[indx+4][c])+fabsf(rgb[indx][c]-rgb[indx-4][c]) - fabsf(rgb[indx+4][c]-rgb[indx-4][c])));*/ - + glpfv = 0.25*(2.0*rgb[1][indx]+rgb[1][indx+v2]+rgb[1][indx-v2]); glpfh = 0.25*(2.0*rgb[1][indx]+rgb[1][indx+2]+rgb[1][indx-2]); rblpfv[indx>>1] = eps+fabsf(glpfv - 0.25*(2.0*rgb[c][indx]+rgb[c][indx+v2]+rgb[c][indx-v2])); @@ -416,27 +422,29 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { grblpfv[indx>>1] = glpfv + 0.25*(2.0*rgb[c][indx]+rgb[c][indx+v2]+rgb[c][indx-v2]); grblpfh[indx>>1] = glpfh + 0.25*(2.0*rgb[c][indx]+rgb[c][indx+2]+rgb[c][indx-2]); } - + areawt[0][0]=areawt[1][0]=1; + areawt[0][2]=areawt[1][2]=1; + // along line segments, find the point along each segment that minimizes the color variance - // averaged over the tile; evaluate for up/down and left/right away from R/B grid point + // averaged over the tile; evaluate for up/down and left/right away from R/B grid point for (rr=8; rr < rr1-8; rr++) for (cc=8+(FC(rr,2)&1), indx=rr*TS+cc, c = FC(rr,cc); cc < cc1-8; cc+=2, indx+=2) { - - areawt[0][c]=areawt[1][c]=0; - + +// areawt[0][c]=areawt[1][c]=0; + //in linear interpolation, color differences are a quadratic function of interpolation position; //solve for the interpolation position that minimizes color difference variance over the tile - + //vertical gdiff=0.3125*(rgb[1][indx+TS]-rgb[1][indx-TS])+0.09375*(rgb[1][indx+TS+1]-rgb[1][indx-TS+1]+rgb[1][indx+TS-1]-rgb[1][indx-TS-1]); deltgrb=(rgb[c][indx]-rgb[1][indx]); - + gradwt=fabsf(0.25*rbhpfv[indx>>1]+0.125*(rbhpfv[(indx>>1)+1]+rbhpfv[(indx>>1)-1]) )*(grblpfv[(indx>>1)-v1]+grblpfv[(indx>>1)+v1])/(eps+0.1*grblpfv[(indx>>1)-v1]+rblpfv[(indx>>1)-v1]+0.1*grblpfv[(indx>>1)+v1]+rblpfv[(indx>>1)+v1]); - + coeff[0][0][c] += gradwt*deltgrb*deltgrb; coeff[0][1][c] += gradwt*gdiff*deltgrb; coeff[0][2][c] += gradwt*gdiff*gdiff; - areawt[0][c]+=1; +// areawt[0][c]+=1; //horizontal gdiff=0.3125*(rgb[1][indx+1]-rgb[1][indx-1])+0.09375*(rgb[1][indx+1+TS]-rgb[1][indx-1+TS]+rgb[1][indx+1-TS]-rgb[1][indx-1-TS]); @@ -447,19 +455,19 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { coeff[1][0][c] += gradwt*deltgrb*deltgrb; coeff[1][1][c] += gradwt*gdiff*deltgrb; coeff[1][2][c] += gradwt*gdiff*gdiff; - areawt[1][c]+=1; +// areawt[1][c]+=1; // In Mathematica, // f[x_]=Expand[Total[Flatten[ // ((1-x) RotateLeft[Gint,shift1]+x RotateLeft[Gint,shift2]-cfapad)^2[[dv;;-1;;2,dh;;-1;;2]]]]]; - // extremum = -.5Coefficient[f[x],x]/Coefficient[f[x],x^2] - } - + // extremum = -.5Coefficient[f[x],x]/Coefficient[f[x],x^2] + } + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% /* for (rr=4; rr < rr1-4; rr++) for (cc=4+(FC(rr,2)&1), indx=rr*TS+cc, c = FC(rr,cc); cc < cc1-4; cc+=2, indx+=2) { - + rbhpfv[indx] = SQR(fabs((rgb[indx][1]-rgb[indx][c])-(rgb[indx+v4][1]-rgb[indx+v4][c])) + fabs((rgb[indx-v4][1]-rgb[indx-v4][c])-(rgb[indx][1]-rgb[indx][c])) - @@ -467,8 +475,8 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { rbhpfh[indx] = SQR(fabs((rgb[indx][1]-rgb[indx][c])-(rgb[indx+4][1]-rgb[indx+4][c])) + fabs((rgb[indx-4][1]-rgb[indx-4][c])-(rgb[indx][1]-rgb[indx][c])) - fabs((rgb[indx-4][1]-rgb[indx-4][c])-(rgb[indx+4][1]-rgb[indx+4][c]))); - - + + glpfv = 0.25*(2*rgb[indx][1]+rgb[indx+v2][1]+rgb[indx-v2][1]); glpfh = 0.25*(2*rgb[indx][1]+rgb[indx+2][1]+rgb[indx-2][1]); rblpfv[indx] = eps+fabs(glpfv - 0.25*(2*rgb[indx][c]+rgb[indx+v2][c]+rgb[indx-v2][c])); @@ -476,23 +484,23 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { grblpfv[indx] = glpfv + 0.25*(2*rgb[indx][c]+rgb[indx+v2][c]+rgb[indx-v2][c]); grblpfh[indx] = glpfh + 0.25*(2*rgb[indx][c]+rgb[indx+2][c]+rgb[indx-2][c]); } - + for (c=0;c<3;c++) {areawt[0][c]=areawt[1][c]=0;} - + // along line segments, find the point along each segment that minimizes the color variance - // averaged over the tile; evaluate for up/down and left/right away from R/B grid point + // averaged over the tile; evaluate for up/down and left/right away from R/B grid point for (rr=rrmin+8; rr < rrmax-8; rr++) for (cc=ccmin+8+(FC(rr,2)&1), indx=rr*TS+cc, c = FC(rr,cc); cc < ccmax-8; cc+=2, indx+=2) { - + if (rgb[indx][c]>0.8*clip_pt || Gtmp[indx]>0.8*clip_pt) continue; //in linear interpolation, color differences are a quadratic function of interpolation position; //solve for the interpolation position that minimizes color difference variance over the tile - + //vertical gdiff=0.3125*(rgb[indx+TS][1]-rgb[indx-TS][1])+0.09375*(rgb[indx+TS+1][1]-rgb[indx-TS+1][1]+rgb[indx+TS-1][1]-rgb[indx-TS-1][1]); deltgrb=(rgb[indx][c]-rgb[indx][1])-0.5*((rgb[indx-v4][c]-rgb[indx-v4][1])+(rgb[indx+v4][c]-rgb[indx+v4][1])); - + gradwt=fabs(0.25*rbhpfv[indx]+0.125*(rbhpfv[indx+2]+rbhpfv[indx-2]) );// *(grblpfv[indx-v2]+grblpfv[indx+v2])/(eps+0.1*grblpfv[indx-v2]+rblpfv[indx-v2]+0.1*grblpfv[indx+v2]+rblpfv[indx+v2]); if (gradwt>eps) { coeff[0][0][c] += gradwt*deltgrb*deltgrb; @@ -500,11 +508,11 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { coeff[0][2][c] += gradwt*gdiff*gdiff; areawt[0][c]++; } - + //horizontal gdiff=0.3125*(rgb[indx+1][1]-rgb[indx-1][1])+0.09375*(rgb[indx+1+TS][1]-rgb[indx-1+TS][1]+rgb[indx+1-TS][1]-rgb[indx-1-TS][1]); deltgrb=(rgb[indx][c]-rgb[indx][1])-0.5*((rgb[indx-4][c]-rgb[indx-4][1])+(rgb[indx+4][c]-rgb[indx+4][1])); - + gradwt=fabs(0.25*rbhpfh[indx]+0.125*(rbhpfh[indx+v2]+rbhpfh[indx-v2]) );// *(grblpfh[indx-2]+grblpfh[indx+2])/(eps+0.1*grblpfh[indx-2]+rblpfh[indx-2]+0.1*grblpfh[indx+2]+rblpfh[indx+2]); if (gradwt>eps) { coeff[1][0][c] += gradwt*deltgrb*deltgrb; @@ -512,11 +520,11 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { coeff[1][2][c] += gradwt*gdiff*gdiff; areawt[1][c]++; } - + // In Mathematica, // f[x_]=Expand[Total[Flatten[ // ((1-x) RotateLeft[Gint,shift1]+x RotateLeft[Gint,shift2]-cfapad)^2[[dv;;-1;;2,dh;;-1;;2]]]]]; - // extremum = -.5Coefficient[f[x],x]/Coefficient[f[x],x^2] + // extremum = -.5Coefficient[f[x],x]/Coefficient[f[x],x^2] }*/ for (c=0; c<3; c+=2){ for (j=0; j<2; j++) {// vert/hor @@ -538,7 +546,7 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { //j=0=vert, 1=hor //offset gives NW corner of square containing the min; j=0=vert, 1=hor - + if (fabsf(CAshift[j][c])<2.0f) { blockavethr[j][c] += CAshift[j][c]; blocksqavethr[j][c] += SQR(CAshift[j][c]); @@ -546,9 +554,9 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { } }//vert/hor }//color - - /* CAshift[j][c] are the locations - that minimize color difference variances; + + /* CAshift[j][c] are the locations + that minimize color difference variances; This is the approximate _optical_ location of the R/B pixels */ for (c=0; c<3; c+=2) { @@ -558,19 +566,24 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { //data structure: blockshifts[blocknum][R/B][v/h] //if (c==0) printf("vblock= %d hblock= %d blockshiftsmedian= %f \n",vblock,hblock,blockshifts[(vblock)*hblsz+hblock][c][0]); } - if(plistener) { - progress+=(double)((TS-border2)*(TS-border2))/(2*height*width); - if (progress>1.0) - { - progress=1.0; + progresscounter++; + if(progresscounter % 8 == 0) +#pragma omp critical + { + progress+=(double)(8.0*(TS-border2)*(TS-border2))/(2*height*width); + if (progress>1.0) + { + progress=1.0; + } + plistener->setProgress(progress); } - plistener->setProgress(progress); } + } //end of diagnostic pass -#pragma omp critical +#pragma omp critical { for (j=0; j<2; j++) for (c=0; c<3; c+=2) { @@ -593,15 +606,18 @@ void RawImageSource::CA_correct_RT(double cared, double cablue) { break; } } - +} //printf ("tile variances %f %f %f %f \n",blockvar[0][0],blockvar[1][0],blockvar[0][2],blockvar[1][2] ); - - + + // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - + //now prepare for CA correction pass //first, fill border blocks of blockshift array if(processpasstwo) { +#pragma omp sections +{ +#pragma omp section for (vblock=1; vblock4.0*blockvar[0][c] || SQR(blockshifts[(vblock)*hblsz+hblock][c][1])>4.0*blockvar[1][c]) continue; - numblox[c] += 1; + if (SQR(bstemp[c][0])>4.0*blockvar[0][c] || SQR(bstemp[c][1])>4.0*blockvar[1][c]) + continue; + numbloxthr[c]++; for (dir=0; dir<2; dir++) { for (i=0; iwidth) {ccmax=width-left;} else {ccmax=cc1;} // rgb from input CFA data - // rgb values should be floating point number between 0 and 1 - // after white balance multipliers are applied + // rgb values should be floating point number between 0 and 1 + // after white balance multipliers are applied for (rr=rrmin; rr < rrmax; rr++) for (row=rr+top, cc=ccmin; cc < ccmax; cc++) { col = cc+left; c = FC(rr,cc); indx=row*width+col; - indx1=rr*TS+cc; + indx1=rr*TS+cc; //rgb[indx1][c] = image[indx][c]/65535.0f; rgb[c][indx1] = (rawData[row][col])/65535.0f; //rgb[indx1][c] = image[indx][c]/65535.0f;//for dcraw implementation @@ -732,7 +778,7 @@ if(processpasstwo) { // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //fill borders if (rrmin>0) { - for (rr=0; rr0) { - for (rr=rrmin; rr0 && ccmin>0) { - for (rr=0; rr0 && ccmax0) { - for (rr=0; rr0) { GRBdir[0][c] = 1; } else { @@ -881,28 +927,28 @@ if(processpasstwo) { } else { GRBdir[1][c] = -1; } - + } - - + + for (rr=4; rr < rr1-4; rr++) for (cc=4+(FC(rr,2)&1), c = FC(rr,cc); cc < cc1-4; cc+=2) { //perform CA correction using color ratios or color differences - + Ginthfloor=(1-shifthfrac[c])*rgb[1][(rr+shiftvfloor[c])*TS+cc+shifthfloor[c]]+(shifthfrac[c])*rgb[1][(rr+shiftvfloor[c])*TS+cc+shifthceil[c]]; Ginthceil=(1-shifthfrac[c])*rgb[1][(rr+shiftvceil[c])*TS+cc+shifthfloor[c]]+(shifthfrac[c])*rgb[1][(rr+shiftvceil[c])*TS+cc+shifthceil[c]]; //Gint is blinear interpolation of G at CA shift point Gint=(1-shiftvfrac[c])*Ginthfloor+(shiftvfrac[c])*Ginthceil; - + //determine R/B at grid points using color differences at shift point plus interpolated G value at grid point //but first we need to interpolate G-R/G-B to grid points... grbdiff[((rr)*TS+cc)>>1]=Gint-rgb[c][(rr)*TS+cc]; gshift[((rr)*TS+cc)>>1]=Gint; } - + for (rr=8; rr < rr1-8; rr++) for (cc=8+(FC(rr,2)&1), c = FC(rr,cc), indx=rr*TS+cc; cc < cc1-8; cc+=2, indx+=2) { - + //if (rgb[indx][c]>clip_pt || Gtmp[indx]>clip_pt) continue; grbdiffold = rgb[1][indx]-rgb[c][indx]; @@ -912,37 +958,37 @@ if(processpasstwo) { grbdiffinthceil=(1.0f-shifthfrac[c]/2.0f)*grbdiff[((rr-2*GRBdir[0][c])*TS+cc)>>1]+(shifthfrac[c]/2.0f)*grbdiff[((rr-2*GRBdir[0][c])*TS+cc-2*GRBdir[1][c])>>1]; //grbdiffint is bilinear interpolation of G-R/G-B at grid point grbdiffint=(1.0f-shiftvfrac[c]/2.0f)*grbdiffinthfloor+(shiftvfrac[c]/2.0f)*grbdiffinthceil; - + //now determine R/B at grid points using interpolated color differences and interpolated G value at grid point RBint=rgb[1][indx]-grbdiffint; - + if (fabsf(RBint-rgb[c][indx])<0.25f*(RBint+rgb[c][indx])) { if (fabsf(grbdiffold)>fabsf(grbdiffint) ) { rgb[c][indx]=RBint; } } else { - + //gradient weights using difference from G at CA shift points and G at grid points p[0]=1.0f/(eps+fabsf(rgb[1][indx]-gshift[indx>>1])); p[1]=1.0f/(eps+fabsf(rgb[1][indx]-gshift[(indx-2*GRBdir[1][c])>>1])); p[2]=1.0f/(eps+fabsf(rgb[1][indx]-gshift[((rr-2*GRBdir[0][c])*TS+cc)>>1])); p[3]=1.0f/(eps+fabsf(rgb[1][indx]-gshift[((rr-2*GRBdir[0][c])*TS+cc-2*GRBdir[1][c])>>1])); - + grbdiffint = (p[0]*grbdiff[indx>>1]+p[1]*grbdiff[(indx-2*GRBdir[1][c])>>1]+ p[2]*grbdiff[((rr-2*GRBdir[0][c])*TS+cc)>>1]+p[3]*grbdiff[((rr-2*GRBdir[0][c])*TS+cc-2*GRBdir[1][c])>>1])/(p[0]+p[1]+p[2]+p[3]); - + //now determine R/B at grid points using interpolated color differences and interpolated G value at grid point if (fabsf(grbdiffold)>fabsf(grbdiffint) ) { rgb[c][indx]=rgb[1][indx]-grbdiffint; } } - + //if color difference interpolation overshot the correction, just desaturate if (grbdiffold*grbdiffint<0) { rgb[c][indx]=rgb[1][indx]-0.5f*(grbdiffold+grbdiffint); } } - + // copy CA corrected results to temporary image matrix for (rr=border; rr < rr1-border; rr++){ c = FC(rr+top, left + border+FC(rr+top,2)&1); @@ -954,13 +1000,19 @@ if(processpasstwo) { } if(plistener) { - progress+=(double)((TS-border2)*(TS-border2))/(2*height*width); - if (progress>1.0) - { - progress=1.0; + progresscounter++; + if(progresscounter % 8 == 0) +#pragma omp critical + { + progress+=(double)(8.0*(TS-border2)*(TS-border2))/(2*height*width); + if (progress>1.0) + { + progress=1.0; + } + plistener->setProgress(progress); } - plistener->setProgress(progress); } + } #pragma omp barrier @@ -973,14 +1025,16 @@ if(processpasstwo) { } // clean up free(buffer); - + } - + free(Gtmp); free(buffer1); free(RawDataTmp); - + if(plistener) + plistener->setProgress(1.0); + #undef TS #undef TSH #undef PIX_SORT diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h index 34994e720..8e21cc145 100644 --- a/rtengine/rawimagesource.h +++ b/rtengine/rawimagesource.h @@ -7,7 +7,7 @@ * 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 @@ -109,7 +109,7 @@ class RawImageSource : public ImageSource { array2D rawData; // holds preprocessed pixel values, rowData[i][j] corresponds to the ith row and jth column // the interpolated green plane: - array2D green; + array2D green; // the interpolated red plane: array2D red; // the interpolated blue plane: @@ -207,7 +207,7 @@ class RawImageSource : public ImageSource { inline void interpolate_row_rb (float* ar, float* ab, float* pg, float* cg, float* ng, int i); inline void interpolate_row_rb_mul_pp (float* ar, float* ab, float* pg, float* cg, float* ng, int i, double r_mul, double g_mul, double b_mul, int x1, int width, int skip); - int LinEqSolve( int nDim, float* pfMatr, float* pfVect, float* pfSolution);//Emil's CA auto correction + int LinEqSolve( int nDim, double* pfMatr, double* pfVect, double* pfSolution);//Emil's CA auto correction void CA_correct_RT (double cared, double cablue); void ddct8x8s(int isgn, float a[8][8]); void processRawWhitepoint (float expos, float preser); // exposure before interpolation