diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc index eab5f1b85..d4e053a85 100644 --- a/rtengine/demosaic_algos.cc +++ b/rtengine/demosaic_algos.cc @@ -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 @@ -28,29 +28,26 @@ #include #include #include - - +#include #ifdef _OPENMP #include #endif namespace rtengine { - #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)) - + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% void RawImageSource::eahd_demosaic () { - if (plistener) { plistener->setProgressStr ("Demosaicing..."); plistener->setProgress (0.0); @@ -75,7 +72,7 @@ void RawImageSource::eahd_demosaic () { threshold = (int)(0.008856*MAXVAL); for (int i=0; iISGREEN(i-1,j)) green[i-1][j] = rawData[i-1][j]; - else { + else { hc = homh[imx][j]; vc = homv[imx][j]; - if (hc > vc) + if (hc > vc) green[i-1][j] = gh[(i-1)%4][j]; - else if (hc < vc) + else if (hc < vc) green[i-1][j] = gv[(i-1)%4][j]; - else + else green[i-1][j] = (gh[(i-1)%4][j] + gv[(i-1)%4][j]) / 2; } } - if (!(i%20) && plistener) + if (!(i%20) && plistener) plistener->setProgress ((double)i / (H-2)); } // finish H-2th and H-1th row, homogenity value is still valailable @@ -281,10 +275,10 @@ void RawImageSource::eahd_demosaic () { green[i-1][j] = gh[(i-1)%4][j]; else if (hc < vc) green[i-1][j] = gv[(i-1)%4][j]; - else + 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); @@ -310,22 +304,20 @@ void RawImageSource::eahd_demosaic () { interpolate_row_rb_mul_pp (red[i], blue[i], green[i-1], green[i], NULL, i, 1.0, 1.0, 1.0, 0, W, 1); else interpolate_row_rb_mul_pp (red[i], blue[i], green[i-1], green[i], green[i+1], i, 1.0, 1.0, 1.0, 0, W, 1); - } } - + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 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; kISGREEN(i,j)) green[i][j] = rawData[i][j]; else { - if (hpmap[i][j]==1) { + if (hpmap[i][j]==1) { int g2 = rawData[i][j+1] + ((rawData[i][j] - rawData[i][j+2]) /2); int g4 = rawData[i][j-1] + ((rawData[i][j] - rawData[i][j-2]) /2); @@ -415,19 +405,19 @@ void RawImageSource::hphd_green (float** hpmap) { int d2 = rawData[i][j+2] - rawData[i][j]; int d3 = (rawData[i-1][j+2] - rawData[i-1][j]) /2; int d4 = (rawData[i+1][j+2] - rawData[i+1][j]) /2; - + double e2 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - + d1 = rawData[i][j-3] - rawData[i][j-1]; d2 = rawData[i][j-2] - rawData[i][j]; d3 = (rawData[i-1][j-2] - rawData[i-1][j]) /2; d4 = (rawData[i+1][j-2] - rawData[i+1][j]) /2; - + double e4 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); green[i][j] = (e2 * g2 + e4 * g4) / (e2 + e4); } - else if (hpmap[i][j]==2) { + else if (hpmap[i][j]==2) { int g1 = rawData[i-1][j] + ((rawData[i][j] - rawData[i-2][j]) /2); int g3 = rawData[i+1][j] + ((rawData[i][j] - rawData[i+2][j]) /2); @@ -436,16 +426,16 @@ void RawImageSource::hphd_green (float** hpmap) { int d2 = rawData[i][j] - rawData[i-2][j]; int d3 = (rawData[i][j-1] - rawData[i-2][j-1]) /2; int d4 = (rawData[i][j+1] - rawData[i-2][j+1]) /2; - + double e1 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - + d1 = rawData[i+1][j] - rawData[i+3][j]; d2 = rawData[i][j] - rawData[i+2][j]; d3 = (rawData[i][j-1] - rawData[i+2][j-1]) /2; d4 = (rawData[i][j+1] - rawData[i+2][j+1]) /2; - + double e3 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - + green[i][j] = (e1 * g1 + e3 * g3) / (e1 + e3); } else { @@ -453,7 +443,7 @@ void RawImageSource::hphd_green (float** hpmap) { int g2 = rawData[i][j+1] + ((rawData[i][j] - rawData[i][j+2]) /2); int g3 = rawData[i+1][j] + ((rawData[i][j] - rawData[i+2][j]) /2); int g4 = rawData[i][j-1] + ((rawData[i][j] - rawData[i][j-2]) /2); - + int dx = rawData[i][j+1] - rawData[i][j-1]; int dy = rawData[i+1][j] - rawData[i-1][j]; @@ -461,29 +451,29 @@ void RawImageSource::hphd_green (float** hpmap) { int d2 = rawData[i][j] - rawData[i-2][j]; int d3 = (rawData[i][j-1] - rawData[i-2][j-1]) /2; int d4 = (rawData[i][j+1] - rawData[i-2][j+1]) /2; - + double e1 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); d1 = rawData[i][j+3] - rawData[i][j+1]; d2 = rawData[i][j+2] - rawData[i][j]; d3 = (rawData[i-1][j+2] - rawData[i-1][j]) /2; d4 = (rawData[i+1][j+2] - rawData[i+1][j]) /2; - + double e2 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); d1 = rawData[i+1][j] - rawData[i+3][j]; d2 = rawData[i][j] - rawData[i+2][j]; d3 = (rawData[i][j-1] - rawData[i+2][j-1]) /2; d4 = (rawData[i][j+1] - rawData[i+2][j+1]) /2; - + double e3 = 1.0 / (1.0 + ABS(dy) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); - + d1 = rawData[i][j-3] - rawData[i][j-1]; d2 = rawData[i][j-2] - rawData[i][j]; d3 = (rawData[i-1][j-2] - rawData[i-1][j]) /2; d4 = (rawData[i+1][j-2] - rawData[i+1][j]) /2; - - double e4 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); + + double e4 = 1.0 / (1.0 + ABS(dx) + ABS(d1) + ABS(d2) + ABS(d3) + ABS(d4)); green[i][j] = (e1*g1 + e2*g2 + e3*g3 + e4*g4) / (e1 + e2 + e3 + e4); } @@ -491,18 +481,17 @@ void RawImageSource::hphd_green (float** hpmap) { } } } - + //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% void RawImageSource::hphd_demosaic () { - if (plistener) { plistener->setProgressStr ("Demosaicing..."); plistener->setProgress (0.0); } float** hpmap = allocArray< float >(W,H, true); - + #ifdef _OPENMP #pragma omp parallel { @@ -518,7 +507,7 @@ void RawImageSource::hphd_demosaic () { #else hphd_vertical (hpmap, 0, W); #endif - if (plistener) + if (plistener) plistener->setProgress (0.33); for (int i=0; i(hpmap, H); - - if (plistener) + + if (plistener) plistener->setProgress (0.66); for (int i=0; isetProgress (1.0); } - -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #define FORCC for (c=0; c < colors; c++) @@ -569,7 +557,6 @@ void RawImageSource::hphd_demosaic () { (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, @@ -647,7 +634,6 @@ void RawImageSource::vng4_demosaic () { // 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++) { @@ -730,7 +716,7 @@ void RawImageSource::vng4_demosaic () { 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) + if (!(row%20) && plistener) plistener->setProgress ((double)row / (H-2)); } memcpy (image[(row-2)*width+2], brow[0]+2, (width-4)*sizeof *image); @@ -756,9 +742,8 @@ void RawImageSource::vng4_demosaic () { free (image); } -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #undef fc #define fc(row,col) \ @@ -783,7 +768,7 @@ void RawImageSource::ppg_demosaic() plistener->setProgressStr ("Demosaicing..."); plistener->setProgress (0.0); } - + image = (float (*)[4]) calloc (H*W, sizeof *image); for (int ii=0; iisetProgressStr ("AHD Demosaicing..."); plistener->setProgress (0.0); } - + image = (float (*)[4]) calloc (H*W, sizeof *image); for (int ii=0; ii 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++) @@ -947,14 +931,13 @@ void RawImageSource::ahd_demosaic(int winx, int winy, int winw, int winh) rgb = (float(*)[TS][TS][3]) buffer; lab = (float(*)[TS][TS][3])(buffer + 6*TS*TS*sizeof(float)); homo = (char (*)[TS][TS]) (buffer + 12*TS*TS*sizeof(float)); - + // 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); @@ -997,11 +980,11 @@ void RawImageSource::ahd_demosaic(int winx, int winy, int winw, int winh) xyz[1] += xyz_cam[1][c] * rix[0][c]; xyz[2] += xyz_cam[2][c] * rix[0][c]; } - + xyz[0] = CurveFactory::flinterp(cbrt,xyz[0]); xyz[1] = CurveFactory::flinterp(cbrt,xyz[1]); xyz[2] = CurveFactory::flinterp(cbrt,xyz[2]); - + //xyz[0] = xyz[0] > 0.008856 ? pow(xyz[0]/65535,1/3.0) : 7.787*xyz[0] + 16/116.0; //xyz[1] = xyz[1] > 0.008856 ? pow(xyz[1]/65535,1/3.0) : 7.787*xyz[1] + 16/116.0; //xyz[2] = xyz[2] > 0.008856 ? pow(xyz[2]/65535,1/3.0) : 7.787*xyz[2] + 16/116.0; @@ -1052,13 +1035,13 @@ void RawImageSource::ahd_demosaic(int winx, int winy, int winw, int winh) 0.5*(rgb[0][tr][tc][c] + rgb[1][tr][tc][c]) ; } } - + tile++; if(plistener) { plistener->setProgress((double)tile / n_tiles); } } - + if(plistener) plistener->setProgress (1.0); free (buffer); for (int i=0; i0) cache[(row-y0+TILEBORDER) * CACHESIZE +TILEBORDER + col-x0][c] = sum[c] / sum[c + 4]; } } @@ -1213,14 +1197,17 @@ void RawImageSource::dcb_hid(float (*image)[4],float (*bufferH)[3], float (*buff // green pixels 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) { - bufferH[indx][1] = ( image[indx-1][1] + image[indx+1][1]) * 0.5; - bufferV[indx][1] = (image[indx+u][1] + image[indx-u][1])*0.5; + assert(indx>=0 && indx=0 && indx=0 && c<3); + + bufferH[indx][c] = ( 4*bufferH[indx][1] - bufferH[indx+u+1][1] - bufferH[indx+u-1][1] - bufferH[indx-u+1][1] - bufferH[indx-u-1][1] + image[indx+u+1][c] + image[indx+u-1][c] + image[indx-u+1][c] + image[indx-u-1][c] ) * 0.25; bufferV[indx][c] = ( 4*bufferV[indx][1] @@ -1231,6 +1218,7 @@ void RawImageSource::dcb_hid(float (*image)[4],float (*bufferH)[3], float (*buff // red or blue in green pixels for (int row=rowMin; row=0 && indx=0 && c<3 && d>=0 && d<3); bufferH[indx][c] = ( image[indx+1][c] + image[indx-1][c]) * 0.5; bufferH[indx][d] = (2*bufferH[indx][1] - bufferH[indx+u][1] - bufferH[indx-u][1] + image[indx+u][d] + image[indx-u][d]) * 0.5; bufferV[indx][c] = (2*bufferV[indx][1] - bufferV[indx+1][1] - bufferV[indx-1][1] + image[indx+1][c] + image[indx-1][c]) * 0.5; @@ -1240,30 +1228,29 @@ void RawImageSource::dcb_hid(float (*image)[4],float (*bufferH)[3], float (*buff // Decide green pixels 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),d=2-c; col < colMax; col+=2, indx+=2) { - int current = MAX(image[indx+v][c], MAX(image[indx-v][c], MAX(image[indx-2][c], image[indx+2][c]))) - + float current = MAX(image[indx+v][c], MAX(image[indx-v][c], MAX(image[indx-2][c], image[indx+2][c]))) - MIN(image[indx+v][c], MIN(image[indx-v][c], MIN(image[indx-2][c], image[indx+2][c]))) + MAX(image[indx+1+u][d], MAX(image[indx+1-u][d], MAX(image[indx-1+u][d], image[indx-1-u][d]))) - MIN(image[indx+1+u][d], MIN(image[indx+1-u][d], MIN(image[indx-1+u][d], image[indx-1-u][d]))); - int currentH = MAX(bufferH[indx+v][d], MAX(bufferH[indx-v][d], MAX(bufferH[indx-2][d], bufferH[indx+2][d]))) - + float currentH = MAX(bufferH[indx+v][d], MAX(bufferH[indx-v][d], MAX(bufferH[indx-2][d], bufferH[indx+2][d]))) - MIN(bufferH[indx+v][d], MIN(bufferH[indx-v][d], MIN(bufferH[indx-2][d], bufferH[indx+2][d]))) + MAX(bufferH[indx+1+u][c], MAX(bufferH[indx+1-u][c], MAX(bufferH[indx-1+u][c], bufferH[indx-1-u][c]))) - MIN(bufferH[indx+1+u][c], MIN(bufferH[indx+1-u][c], MIN(bufferH[indx-1+u][c], bufferH[indx-1-u][c]))); - int currentV = MAX(bufferV[indx+v][d], MAX(bufferV[indx-v][d], MAX(bufferV[indx-2][d], bufferV[indx+2][d]))) - + float currentV = MAX(bufferV[indx+v][d], MAX(bufferV[indx-v][d], MAX(bufferV[indx-2][d], bufferV[indx+2][d]))) - MIN(bufferV[indx+v][d], MIN(bufferV[indx-v][d], MIN(bufferV[indx-2][d], bufferV[indx+2][d]))) + MAX(bufferV[indx+1+u][c], MAX(bufferV[indx+1-u][c], MAX(bufferV[indx-1+u][c], bufferV[indx-1-u][c]))) - MIN(bufferV[indx+1+u][c], MIN(bufferV[indx+1-u][c], MIN(bufferV[indx-1+u][c], bufferV[indx-1-u][c]))); + assert(indx>=0 && indx=0 && indx=0 && c<4); + image[indx][c] = ( 4.0 * image[indx][1] - image[indx+u+1][1] - image[indx+u-1][1] - image[indx-u+1][1] - image[indx-u-1][1] + image[indx+u+1][c] + image[indx+u-1][c] + image[indx-u+1][c] + image[indx-u-1][c] ) * 0.25; } @@ -1282,9 +1270,10 @@ void RawImageSource::dcb_color(float (*image)[4], int x0, int y0) // red or blue in green pixels for (int row=rowMin; row=0 && indx=0 && c<4); + image[indx][c] = (2.0 * image[indx][1] - image[indx+1][1] - image[indx-1][1] + image[indx+1][c] + image[indx-1][c]) * 0.5; + image[indx][d] = (2.0 * image[indx][1] - image[indx+u][1] - image[indx-u][1] + image[indx+u][d] + image[indx-u][d]) * 0.5; + } } // green correction @@ -1293,22 +1282,23 @@ void RawImageSource::dcb_hid2(float (*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) { - image[indx][1] = (image[indx+v][1] + image[indx-v][1] + image[indx-2][1] + image[indx+2][1])/4 + + assert(indx>=0 && indx=0 && indx ( 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 @@ -1324,20 +1316,21 @@ void RawImageSource::dcb_map(float (*image)[4], int x0, int y0) } } - // interpolated green pixels are corrected using the map void RawImageSource::dcb_correction(float (*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]) + + register float current = 4.0 * image[indx][3] + + 2.0 * (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; + + assert(indx>=0 && indx=0 && indx=0 && indx=0 && indx=0 && indx=0 && c<4); chroma[indx][d]=image[indx][c]-image[indx][1]; + } for (int row=rowMin; row=0 && indx=0 && c<2); chroma[indx][c]=(f[0]*g[0]+f[1]*g[1]+f[2]*g[2]+f[3]*g[3])/(f[0]+f[1]+f[2]+f[3]); } for (int row=rowMin; row=0 && indx=0 && c<2); chroma[indx][c]=(f[0]*g[0]+f[1]*g[1]+f[2]*g[2]+f[3]*g[3])/(f[0]+f[1]+f[2]+f[3]); } for(int row=rowMin; row=0 && indx