diff --git a/rtengine/array2D.h b/rtengine/array2D.h index 93734a421..f9d92ac36 100644 --- a/rtengine/array2D.h +++ b/rtengine/array2D.h @@ -23,7 +23,7 @@ * Usage: * * array2D name (X-size,Y-size); - * array2D name (X-size,Y-size type ** data); + * array2D name (X-size,Y-size,type ** data); * * creates an array which is valid within the normal C/C++ scope "{ ... }" * @@ -75,7 +75,7 @@ private: T ** ptr; T * data; bool lock; // useful lock to ensure data is not changed anymore. - void ar_realloc(int w, int h) + void ar_realloc(int w, int h, int offset = 0) { if ((ptr) && ((h > y) || (4 * h < y))) { delete[] ptr; @@ -92,14 +92,14 @@ private: } if (data == NULL) { - data = new T[h * w]; + data = new T[h * w + offset]; } x = w; y = h; for (int i = 0; i < h; i++) { - ptr[i] = data + w * i; + ptr[i] = data + offset + w * i; } owner = 1; @@ -184,6 +184,19 @@ public: } } + void free() + { + if ((owner) && (data)) { + delete[] data; + data = NULL; + } + + if (ptr) { + delete [] ptr; + ptr = NULL; + } + } + // use with indices T * operator[](int index) { @@ -207,7 +220,7 @@ public: // useful within init of parent object // or use as resize of 2D array - void operator()(int w, int h, unsigned int flgs = 0) + void operator()(int w, int h, unsigned int flgs = 0, int offset = 0) { flags = flgs; @@ -223,10 +236,10 @@ public: lock = flags & ARRAY2D_LOCK_DATA; - ar_realloc(w, h); + ar_realloc(w, h, offset); if (flags & ARRAY2D_CLEAR_DATA) { - memset(data, 0, w * h * sizeof(T)); + memset(data + offset, 0, w * h * sizeof(T)); } } @@ -298,10 +311,10 @@ private: array2D list[num]; public: - multi_array2D(int x, int y, int flags = 0) + multi_array2D(int x, int y, int flags = 0, int offset = 0) { for (size_t i = 0; i < num; i++) { - list[i](x, y, flags); + list[i](x, y, flags, (i + 1) * offset); } } @@ -312,11 +325,7 @@ public: array2D & operator[](int index) { - if (static_cast(index) >= num) { - printf("index %0u is out of range[0..%0lu]", index, num - 1); - raise( SIGSEGV); - } - + assert(static_cast(index) < num); return list[index]; } }; diff --git a/rtengine/hilite_recon.cc b/rtengine/hilite_recon.cc index bc5e973e2..46c666c80 100644 --- a/rtengine/hilite_recon.cc +++ b/rtengine/hilite_recon.cc @@ -29,25 +29,15 @@ #include "rt_math.h" #include "opthelper.h" -#ifdef _OPENMP -#include -#endif - - #define FOREACHCOLOR for (int c=0; c < ColorCount; c++) -//#include "RGBdefringe.cc" - namespace rtengine { - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W, int box ) +SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, float** temp, int H, int W, int box ) { - array2D temp(W, H); - //box blur image channel; box size = 2*box+1 //horizontal blur #ifdef _OPENMP @@ -84,10 +74,10 @@ SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W #endif { float len = box + 1; - __m128 lenv = _mm_set1_ps( len ); - __m128 lenp1v = _mm_set1_ps( len + 1.0f ); - __m128 onev = _mm_set1_ps( 1.0f ); - __m128 tempv, temp2v; + vfloat lenv = F2V( len ); + vfloat lenp1v = F2V( len + 1.0f ); + vfloat onev = F2V( 1.0f ); + vfloat tempv, temp2v; #ifdef _OPENMP #pragma omp for nowait #endif @@ -130,7 +120,9 @@ SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W } } +#ifdef _OPENMP #pragma omp single +#endif { for (int col = W - (W % 8); col < W - 3; col += 4) { tempv = LVFU(temp[0][col]) / lenv; @@ -160,29 +152,28 @@ SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W _mm_storeu_ps( &dst[row][col], tempv ); } } - } - } + for (int col = W - (W % 4); col < W; col++) { + int len = box + 1; + dst[0][col] = temp[0][col] / len; - for (int col = W - (W % 4); col < W; col++) { - int len = box + 1; - dst[0][col] = temp[0][col] / len; + for (int i = 1; i <= box; i++) { + dst[0][col] += temp[i][col] / len; + } - for (int i = 1; i <= box; i++) { - dst[0][col] += temp[i][col] / len; - } + for (int row = 1; row <= box; row++) { + dst[row][col] = (dst[(row - 1)][col] * len + temp[(row + box)][col]) / (len + 1); + len ++; + } - for (int row = 1; row <= box; row++) { - dst[row][col] = (dst[(row - 1)][col] * len + temp[(row + box)][col]) / (len + 1); - len ++; - } + for (int row = box + 1; row < H - box; row++) { + dst[row][col] = dst[(row - 1)][col] + (temp[(row + box)][col] - temp[(row - box - 1)][col]) / len; + } - for (int row = box + 1; row < H - box; row++) { - dst[row][col] = dst[(row - 1)][col] + (temp[(row + box)][col] - temp[(row - box - 1)][col]) / len; - } - - for (int row = H - box; row < H; row++) { - dst[row][col] = (dst[(row - 1)][col] * len - temp[(row - box - 1)][col]) / (len - 1); - len --; + for (int row = H - box; row < H; row++) { + dst[row][col] = (dst[(row - 1)][col] * len - temp[(row - box - 1)][col]) / (len - 1); + len --; + } + } } } @@ -219,11 +210,9 @@ SSEFUNCTION void RawImageSource::boxblur2(float** src, float** dst, int H, int W } -void RawImageSource::boxblur_resamp(float **src, float **dst, int H, int W, int box, int samp ) +void RawImageSource::boxblur_resamp(float **src, float **dst, float ** temp, int H, int W, int box, int samp ) { - array2D temp((W / samp) + ((W % samp) == 0 ? 0 : 1), H); - #ifdef _OPENMP #pragma omp parallel #endif @@ -235,7 +224,8 @@ void RawImageSource::boxblur_resamp(float **src, float **dst, int H, int W, int //box blur image channel; box size = 2*box+1 //horizontal blur - for (int row = 0; row < H; row++) { + for (int row = 0; row < H; row++) + { int len = box + 1; tempval = src[row][0] / len; @@ -255,8 +245,10 @@ void RawImageSource::boxblur_resamp(float **src, float **dst, int H, int W, int len ++; } + float oneByLen = 1.f / (float)len; + for (int col = box + 1; col < W - box; col++) { - tempval = tempval + (src[row][col + box] - src[row][col - box - 1]) / len; + tempval = tempval + (src[row][col + box] - src[row][col - box - 1]) * oneByLen; if(col % samp == 0) { temp[row][col / samp] = tempval; @@ -275,63 +267,125 @@ void RawImageSource::boxblur_resamp(float **src, float **dst, int H, int W, int } } + static const int numCols = 8; // process numCols columns at once for better L1 CPU cache usage #ifdef _OPENMP #pragma omp parallel #endif { - float tempval; + float tempvalN[numCols] ALIGNED16; #ifdef _OPENMP - #pragma omp for + #pragma omp for nowait #endif //vertical blur - for (int col = 0; col < W / samp; col++) { + for (int col = 0; col < (W / samp) - (numCols - 1); col += numCols) { int len = box + 1; - tempval = temp[0][col] / len; - for (int i = 1; i <= box; i++) { - tempval += temp[i][col] / len; + for(int n = 0; n < numCols; n++) { + tempvalN[n] = temp[0][col + n] / len; } - dst[0][col] = tempval; + for (int i = 1; i <= box; i++) { + for(int n = 0; n < numCols; n++) { + tempvalN[n] += temp[i][col + n] / len; + } + } + + for(int n = 0; n < numCols; n++) { + dst[0][col + n] = tempvalN[n]; + } for (int row = 1; row <= box; row++) { - tempval = (tempval * len + temp[(row + box)][col]) / (len + 1); + for(int n = 0; n < numCols; n++) { + tempvalN[n] = (tempvalN[n] * len + temp[(row + box)][col + n]) / (len + 1); + } if(row % samp == 0) { - dst[row / samp][col] = tempval; + for(int n = 0; n < numCols; n++) { + dst[row / samp][col + n] = tempvalN[n]; + } } len ++; } for (int row = box + 1; row < H - box; row++) { - tempval = tempval + (temp[(row + box)][col] - temp[(row - box - 1)][col]) / len; + for(int n = 0; n < numCols; n++) { + tempvalN[n] = tempvalN[n] + (temp[(row + box)][col + n] - temp[(row - box - 1)][col + n]) / len; + } if(row % samp == 0) { - dst[row / samp][col] = tempval; + for(int n = 0; n < numCols; n++) { + dst[row / samp][col + n] = tempvalN[n]; + } } } for (int row = H - box; row < H; row++) { - tempval = (tempval * len - temp[(row - box - 1)][col]) / (len - 1); + for(int n = 0; n < numCols; n++) { + tempvalN[n] = (tempvalN[n] * len - temp[(row - box - 1)][col + n]) / (len - 1); + } if(row % samp == 0) { - dst[row / samp][col] = tempval; + for(int n = 0; n < numCols; n++) { + dst[row / samp][col + n] = tempvalN[n]; + } } len --; } } + + // process remaining columns + #pragma omp single + { + float tempval; + + //vertical blur + for (int col = (W / samp) - ((W / samp) % numCols); col < W / samp; col++) { + int len = box + 1; + tempval = temp[0][col] / len; + + for (int i = 1; i <= box; i++) { + tempval += temp[i][col] / len; + } + + dst[0][col] = tempval; + + for (int row = 1; row <= box; row++) { + tempval = (tempval * len + temp[(row + box)][col]) / (len + 1); + + if(row % samp == 0) { + dst[row / samp][col] = tempval; + } + + len ++; + } + + for (int row = box + 1; row < H - box; row++) { + tempval = tempval + (temp[(row + box)][col] - temp[(row - box - 1)][col]) / len; + + if(row % samp == 0) { + dst[row / samp][col] = tempval; + } + } + + for (int row = H - box; row < H; row++) { + tempval = (tempval * len - temp[(row - box - 1)][col]) / (len - 1); + + if(row % samp == 0) { + dst[row / samp][col] = tempval; + } + + len --; + } + } + } } } -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** blue) { double progress = 0.0; @@ -344,12 +398,8 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b int height = H; int width = W; - const int range = 2; - const int pitch = 4; - - - int hfh = (height - (height % pitch)) / pitch; - int hfw = (width - (width % pitch)) / pitch; + static const int range = 2; + static const int pitch = 4; static const int numdirs = 4; @@ -360,7 +410,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b //%%%%%%%%%%%%%%%%%%%% //for blend algorithm: static const float blendthresh = 1.0; - const int ColorCount = 3; + static const int ColorCount = 3; // Transform matrixes rgb>lab and back static const float trans[2][ColorCount][ColorCount] = { { { 1, 1, 1 }, { 1.7320508, -1.7320508, 0 }, { -1, -1, 2 } }, @@ -370,49 +420,12 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b { { 1, 0.8660254, -0.5 }, { 1, -0.8660254, -0.5 }, { 1, 0, 1 } }, { { 1, 1, 1 }, { 1, -1, 1 }, { 1, 1, -1 } } }; - //%%%%%%%%%%%%%%%%%%%% - - float max_f[3], thresh[3], fixthresh[3], norm[3]; - - //float red1, green1, blue1;//diagnostic -// float chmaxalt[4]={0,0,0,0};//diagnostic - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //halfsize demosaic - - /* - multi_array2D hfsize (hfw+1,hfh+1,ARRAY2D_CLEAR_DATA); - - boxblur_resamp(red,hfsize[0],chmaxalt[0],height,width,range,pitch); - if(plistener){ - progress += 0.05; - plistener->setProgress(progress); - } - boxblur_resamp(green,hfsize[1],chmaxalt[1],height,width,range,pitch); - if(plistener){ - progress += 0.05; - plistener->setProgress(progress); - } - boxblur_resamp(blue,hfsize[2],chmaxalt[2],height,width,range,pitch); - if(plistener){ - progress += 0.05; - plistener->setProgress(progress); - } - - //blur image - //for (int m=0; m<3; m++) - // boxblur2(hfsize[m],hfsizeblur[m],hfh,hfw,3); - - */ - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + float max_f[3], thresh[3]; for (int c = 0; c < 3; c++) { thresh[c] = chmax[c] * threshpct; - fixthresh[c] = chmax[c] * fixthreshpct; max_f[c] = chmax[c] * maxpct; //min(chmax[0],chmax[1],chmax[2])*maxpct; - norm[c] = 1.0 / (max_f[c] - fixthresh[c]); } float whitept = max(max_f[0], max_f[1], max_f[2]); @@ -420,43 +433,57 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b float medpt = max_f[0] + max_f[1] + max_f[2] - whitept - clippt; float blendpt = blendthresh * clippt; - float camwb[4]; - - for (int c = 0; c < 4; c++) { - camwb[c] = ri->get_cam_mul(c); - } - - multi_array2D channelblur(width, height, ARRAY2D_CLEAR_DATA); - multi_array2D hilite_full(width, height, ARRAY2D_CLEAR_DATA); - - if(plistener) { - progress += 0.05; - plistener->setProgress(progress); - } + multi_array2D channelblur(width, height, 0, 48); + array2D temp(width, height); // allocate temporary buffer // blur RGB channels - boxblur2(red , channelblur[0], height, width, 4); + + boxblur2(red , channelblur[0], temp, height, width, 4); + if(plistener) { + progress += 0.05; + plistener->setProgress(progress); + } + + boxblur2(green, channelblur[1], temp, height, width, 4); + if(plistener) { + progress += 0.05; + plistener->setProgress(progress); + } + + boxblur2(blue , channelblur[2], temp, height, width, 4); + if(plistener) { + progress += 0.05; + plistener->setProgress(progress); + } + + // reduce channel blur to one array +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for(int i = 0; i < height; i++) + for(int j = 0; j < width; j++) { + channelblur[0][i][j] = fabsf(channelblur[0][i][j] - red[i][j]) + fabsf(channelblur[1][i][j] - green[i][j]) + fabsf(channelblur[2][i][j] - blue[i][j]); + } + + for (int c = 1; c < 3; c++) { + channelblur[c].free(); //free up some memory + } if(plistener) { progress += 0.05; plistener->setProgress(progress); } - boxblur2(green, channelblur[1], height, width, 4); - + multi_array2D hilite_full(width, height, ARRAY2D_CLEAR_DATA, 32); if(plistener) { - progress += 0.05; + progress += 0.10; plistener->setProgress(progress); } - boxblur2(blue , channelblur[2], height, width, 4); - if(plistener) { - progress += 0.05; - plistener->setProgress(progress); - } - - float hipass_sum = 0, hipass_norm = 0.00; + double hipass_sum = 0.f; + int hipass_norm = 0; // set up which pixels are clipped or near clipping #ifdef _OPENMP @@ -465,40 +492,34 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) { - //if one or more channels is highlight but none are blown, add to highlight accumulator - if ((red[i][j] > thresh[0] || green[i][j] > thresh[1] || blue[i][j] > thresh[2]) && (red[i][j] < max_f[0] && green[i][j] < max_f[1] && blue[i][j] < max_f[2])) { - hipass_sum += fabs(channelblur[0][i][j] - red[i][j]) + fabs(channelblur[1][i][j] - green[i][j]) + fabs(channelblur[2][i][j] - blue[i][j]); - hipass_norm++; + hipass_sum += channelblur[0][i][j]; + hipass_norm ++; hilite_full[0][i][j] = red[i][j]; hilite_full[1][i][j] = green[i][j]; hilite_full[2][i][j] = blue[i][j]; - hilite_full[3][i][j] = 1; - hilite_full[4][i][j] = 1; + hilite_full[3][i][j] = 1.f; } - - //if (i%100==0 && j%100==0) - // printf("row=%d col=%d r=%f g=%f b=%f hilite=%f \n",i,j,hilite_full[0][i][j],hilite_full[1][i][j],hilite_full[2][i][j],hilite_full[3][i][j]); } }//end of filling highlight array - hipass_norm += 0.01; - - float hipass_ave = (hipass_sum / hipass_norm); - + float hipass_ave = 2.f * hipass_sum / (hipass_norm + 0.01); if(plistener) { progress += 0.05; plistener->setProgress(progress); } + array2D hilite_full4(width, height); //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //blur highlight data - boxblur2(hilite_full[4], hilite_full[4], height, width, 1); + boxblur2(hilite_full[3], hilite_full4, temp, height, width, 1); + + temp.free(); // free temporary buffer if(plistener) { progress += 0.05; @@ -511,32 +532,34 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int i = 0; i < height; i++) { for (int j = 0; j < width; j++) { - - float hipass = fabs(channelblur[0][i][j] - red[i][j]) + fabs(channelblur[1][i][j] - green[i][j]) + fabs(channelblur[2][i][j] - blue[i][j]); - - if (hipass > 2 * hipass_ave) { + if (channelblur[0][i][j] > hipass_ave) { //too much variation - hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0; + hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0.f; continue; } - if (hilite_full[4][i][j] > 0.00001 && hilite_full[4][i][j] < 0.95) { + if (hilite_full4[i][j] > 0.00001 && hilite_full4[i][j] < 0.95) { //too near an edge, could risk using CA affected pixels, therefore omit - hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0; + hilite_full[0][i][j] = hilite_full[1][i][j] = hilite_full[2][i][j] = hilite_full[3][i][j] = 0.f; } } } - for (int c = 0; c < 3; c++) { - channelblur[c](1, 1); //free up some memory - } + channelblur[0].free(); //free up some memory + hilite_full4.free(); //free up some memory - multi_array2D hilite(hfw + 1, hfh + 1, ARRAY2D_CLEAR_DATA); + int hfh = (height - (height % pitch)) / pitch; + int hfw = (width - (width % pitch)) / pitch; + + multi_array2D hilite(hfw + 1, hfh + 1, ARRAY2D_CLEAR_DATA, 48); //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // blur and resample highlight data; range=size of blur, pitch=sample spacing + + array2D temp2((width / pitch) + ((width % pitch) == 0 ? 0 : 1), height); + for (int m = 0; m < 4; m++) { - boxblur_resamp(hilite_full[m], hilite[m], height, width, range, pitch); + boxblur_resamp(hilite_full[m], hilite[m], temp2, height, width, range, pitch); if(plistener) { progress += 0.05; @@ -544,101 +567,235 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } } - for (int c = 0; c < 5; c++) { - hilite_full[c](1, 1); //free up some memory + temp2.free(); + + for (int c = 0; c < 4; c++) { + hilite_full[c].free(); //free up some memory } - multi_array2D hilite_dir(hfw, hfh, ARRAY2D_CLEAR_DATA); - + multi_array2D hilite_dir(hfw, hfh, ARRAY2D_CLEAR_DATA, 64); + // for faster processing we create two buffers using (height,width) instead of (width,height) + multi_array2D hilite_dir0(hfh, hfw, ARRAY2D_CLEAR_DATA, 64); + multi_array2D hilite_dir4(hfh, hfw, ARRAY2D_CLEAR_DATA, 64); //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //blur highlights - //for (int m=0; m<4; m++) - // boxblur2(hilite[m],hilite[m],hfh,hfw,4); - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if(plistener) { progress += 0.05; plistener->setProgress(progress); } - LUTf invfn(0x10000); - - for (int i = 0; i < 0x10000; i++) { - invfn[i] = 1.0 / (1.0 + i); - } - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //fill gaps in highlight map by directional extension //raster scan from four corners - for (int j = 1; j < hfw - 1; j++) + for (int j = 1; j < hfw - 1; j++) { for (int i = 2; i < hfh - 2; i++) { //from left if (hilite[3][i][j] > 0.01) { - for (int c = 0; c < 4; c++) { - hilite_dir[c][i][j] = hilite[c][i][j] / hilite[3][i][j]; - } + hilite_dir0[3][j][i] = 1.f; } else { - for (int c = 0; c < 4; c++) { - hilite_dir[c][i][j] = 0.1 * ((hilite_dir[0 + c][i - 2][j - 1] + hilite_dir[0 + c][i - 1][j - 1] + hilite_dir[0 + c][i][j - 1] + hilite_dir[0 + c][i + 1][j - 1] + hilite_dir[0 + c][i + 2][j - 1]) / - (hilite_dir[0 + 3][i - 2][j - 1] + hilite_dir[0 + 3][i - 1][j - 1] + hilite_dir[0 + 3][i][j - 1] + hilite_dir[0 + 3][i + 1][j - 1] + hilite_dir[0 + 3][i + 2][j - 1] + 0.00001)); - hilite_dir[4 + c][i][j + 1] += hilite_dir[c][i][j]; - hilite_dir[8 + c][i - 2][j] += hilite_dir[c][i][j]; - hilite_dir[12 + c][i + 2][j] += hilite_dir[c][i][j]; - } + hilite_dir0[3][j][i] = (hilite_dir0[0 + 3][j - 1][i - 2] + hilite_dir0[0 + 3][j - 1][i - 1] + hilite_dir0[0 + 3][j - 1][i] + hilite_dir0[0 + 3][j - 1][i + 1] + hilite_dir0[0 + 3][j - 1][i + 2]) == 0.f ? 0.f : 0.1f; } } + if(hilite[3][2][j] <= 0.01) { + hilite_dir[0 + 3][0][j] = hilite_dir0[3][j][2]; + } + + if(hilite[3][3][j] <= 0.01) { + hilite_dir[0 + 3][1][j] = hilite_dir0[3][j][3]; + } + + if(hilite[3][hfh - 3][j] <= 0.01) { + hilite_dir[4 + 3][hfh - 1][j] = hilite_dir0[3][j][hfh - 3]; + } + + if(hilite[3][hfh - 4][j] <= 0.01) { + hilite_dir[4 + 3][hfh - 2][j] = hilite_dir0[3][j][hfh - 4]; + } + } + + for (int i = 2; i < hfh - 2; i++) { + if(hilite[3][i][hfw - 2] <= 0.01) { + hilite_dir4[3][hfw - 1][i] = hilite_dir0[3][hfw - 2][i]; + } + } + + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int c = 0; c < 3; c++) { + for (int j = 1; j < hfw - 1; j++) { + for (int i = 2; i < hfh - 2; i++) { + //from left + if (hilite[3][i][j] > 0.01f) { + hilite_dir0[c][j][i] = hilite[c][i][j] / hilite[3][i][j]; + } else { + hilite_dir0[c][j][i] = 0.1f * ((hilite_dir0[0 + c][j - 1][i - 2] + hilite_dir0[0 + c][j - 1][i - 1] + hilite_dir0[0 + c][j - 1][i] + hilite_dir0[0 + c][j - 1][i + 1] + hilite_dir0[0 + c][j - 1][i + 2]) / + (hilite_dir0[0 + 3][j - 1][i - 2] + hilite_dir0[0 + 3][j - 1][i - 1] + hilite_dir0[0 + 3][j - 1][i] + hilite_dir0[0 + 3][j - 1][i + 1] + hilite_dir0[0 + 3][j - 1][i + 2] + 0.00001f)); + } + } + + if(hilite[3][2][j] <= 0.01f) { + hilite_dir[0 + c][0][j] = hilite_dir0[c][j][2]; + } + + if(hilite[3][3][j] <= 0.01f) { + hilite_dir[0 + c][1][j] = hilite_dir0[c][j][3]; + } + + if(hilite[3][hfh - 3][j] <= 0.01f) { + hilite_dir[4 + c][hfh - 1][j] = hilite_dir0[c][j][hfh - 3]; + } + + if(hilite[3][hfh - 4][j] <= 0.01f) { + hilite_dir[4 + c][hfh - 2][j] = hilite_dir0[c][j][hfh - 4]; + } + } + + for (int i = 2; i < hfh - 2; i++) { + if(hilite[3][i][hfw - 2] <= 0.01f) { + hilite_dir4[c][hfw - 1][i] = hilite_dir0[c][hfw - 2][i]; + } + } + } + if(plistener) { progress += 0.05; plistener->setProgress(progress); } - for (int j = hfw - 2; j > 0; j--) + for (int j = hfw - 2; j > 0; j--) { for (int i = 2; i < hfh - 2; i++) { //from right if (hilite[3][i][j] > 0.01) { - for (int c = 0; c < 4; c++) { - hilite_dir[4 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; - } + hilite_dir4[3][j][i] = 1.f; } else { - for (int c = 0; c < 4; c++) { - hilite_dir[4 + c][i][j] = 0.1 * ((hilite_dir[4 + c][(i - 2)][(j + 1)] + hilite_dir[4 + c][(i - 1)][(j + 1)] + hilite_dir[4 + c][(i)][(j + 1)] + hilite_dir[4 + c][(i + 1)][(j + 1)] + hilite_dir[4 + c][(i + 2)][(j + 1)]) / - (hilite_dir[4 + 3][(i - 2)][(j + 1)] + hilite_dir[4 + 3][(i - 1)][(j + 1)] + hilite_dir[4 + 3][(i)][(j + 1)] + hilite_dir[4 + 3][(i + 1)][(j + 1)] + hilite_dir[4 + 3][(i + 2)][(j + 1)] + 0.00001)); - hilite_dir[8 + c][i - 2][j] += hilite_dir[4 + c][i][j]; - hilite_dir[12 + c][i + 2][j] += hilite_dir[4 + c][i][j]; - } + hilite_dir4[3][j][i] = (hilite_dir4[3][(j + 1)][(i - 2)] + hilite_dir4[3][(j + 1)][(i - 1)] + hilite_dir4[3][(j + 1)][(i)] + hilite_dir4[3][(j + 1)][(i + 1)] + hilite_dir4[3][(j + 1)][(i + 2)]) == 0.f ? 0.f : 0.1f; } } + if(hilite[3][2][j] <= 0.01) { + hilite_dir[0 + 3][0][j] += hilite_dir4[3][j][2]; + } + + if(hilite[3][hfh - 3][j] <= 0.01) { + hilite_dir[4 + 3][hfh - 1][j] += hilite_dir4[3][j][hfh - 3]; + } + } + + for (int i = 2; i < hfh - 2; i++) { + if(hilite[3][i][0] <= 0.01) { + hilite_dir[0 + 3][i - 2][0] += hilite_dir4[3][0][i]; + hilite_dir[4 + 3][i + 2][0] += hilite_dir4[3][0][i]; + } + + if(hilite[3][i][1] <= 0.01) { + hilite_dir[0 + 3][i - 2][1] += hilite_dir4[3][1][i]; + hilite_dir[4 + 3][i + 2][1] += hilite_dir4[3][1][i]; + } + + if(hilite[3][i][hfw - 2] <= 0.01) { + hilite_dir[0 + 3][i - 2][hfw - 2] += hilite_dir4[3][hfw - 2][i]; + hilite_dir[4 + 3][i + 2][hfw - 2] += hilite_dir4[3][hfw - 2][i]; + } + } + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int c = 0; c < 3; c++) { + for (int j = hfw - 2; j > 0; j--) { + for (int i = 2; i < hfh - 2; i++) { + //from right + if (hilite[3][i][j] > 0.01) { + hilite_dir4[c][j][i] = hilite[c][i][j] / hilite[3][i][j]; + } else { + hilite_dir4[c][j][i] = 0.1 * ((hilite_dir4[c][(j + 1)][(i - 2)] + hilite_dir4[c][(j + 1)][(i - 1)] + hilite_dir4[c][(j + 1)][(i)] + hilite_dir4[c][(j + 1)][(i + 1)] + hilite_dir4[c][(j + 1)][(i + 2)]) / + (hilite_dir4[3][(j + 1)][(i - 2)] + hilite_dir4[3][(j + 1)][(i - 1)] + hilite_dir4[3][(j + 1)][(i)] + hilite_dir4[3][(j + 1)][(i + 1)] + hilite_dir4[3][(j + 1)][(i + 2)] + 0.00001)); + } + } + + if(hilite[3][2][j] <= 0.01) { + hilite_dir[0 + c][0][j] += hilite_dir4[c][j][2]; + } + + if(hilite[3][hfh - 3][j] <= 0.01) { + hilite_dir[4 + c][hfh - 1][j] += hilite_dir4[c][j][hfh - 3]; + } + } + + for (int i = 2; i < hfh - 2; i++) { + if(hilite[3][i][0] <= 0.01) { + hilite_dir[0 + c][i - 2][0] += hilite_dir4[c][0][i]; + hilite_dir[4 + c][i + 2][0] += hilite_dir4[c][0][i]; + } + + if(hilite[3][i][1] <= 0.01) { + hilite_dir[0 + c][i - 2][1] += hilite_dir4[c][1][i]; + hilite_dir[4 + c][i + 2][1] += hilite_dir4[c][1][i]; + } + + if(hilite[3][i][hfw - 2] <= 0.01) { + hilite_dir[0 + c][i - 2][hfw - 2] += hilite_dir4[c][hfw - 2][i]; + hilite_dir[4 + c][i + 2][hfw - 2] += hilite_dir4[c][hfw - 2][i]; + } + } + } + if(plistener) { progress += 0.05; plistener->setProgress(progress); } + for (int i = 1; i < hfh - 1; i++) for (int j = 2; j < hfw - 2; j++) { - //if (i%100==0 && j%100==0) - // printf("row=%d col=%d r=%f g=%f b=%f hilite=%f \n",i,j,hilite[0][i][j],hilite[1][i][j],hilite[2][i][j],hilite[3][i][j]); - //from top if (hilite[3][i][j] > 0.01) { - for (int c = 0; c < 4; c++) { - hilite_dir[8 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; - } + hilite_dir[0 + 3][i][j] = 1.f; } else { - for (int c = 0; c < 4; c++) { - hilite_dir[8 + c][i][j] = 0.1 * ((hilite_dir[8 + c][i - 1][j - 2] + hilite_dir[8 + c][i - 1][j - 1] + hilite_dir[8 + c][i - 1][j] + hilite_dir[8 + c][i - 1][j + 1] + hilite_dir[8 + c][i - 1][j + 2]) / - (hilite_dir[8 + 3][i - 1][j - 2] + hilite_dir[8 + 3][i - 1][j - 1] + hilite_dir[8 + 3][i - 1][j] + hilite_dir[8 + 3][i - 1][j + 1] + hilite_dir[8 + 3][i - 1][j + 2] + 0.00001)); - hilite_dir[12 + c][i + 1][j] += hilite_dir[8 + c][i][j]; + hilite_dir[0 + 3][i][j] = (hilite_dir[0 + 3][i - 1][j - 2] + hilite_dir[0 + 3][i - 1][j - 1] + hilite_dir[0 + 3][i - 1][j] + hilite_dir[0 + 3][i - 1][j + 1] + hilite_dir[0 + 3][i - 1][j + 2]) == 0.f ? 0.f : 0.1f; + } + } + + for (int j = 2; j < hfw - 2; j++) { + if(hilite[3][hfh - 2][j] <= 0.01) { + hilite_dir[4 + 3][hfh - 1][j] += hilite_dir[0 + 3][hfh - 2][j]; + } + } + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int c = 0; c < 3; c++) { + for (int i = 1; i < hfh - 1; i++) { + for (int j = 2; j < hfw - 2; j++) { + //from top + if (hilite[3][i][j] > 0.01) { + hilite_dir[0 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; + } else { + hilite_dir[0 + c][i][j] = 0.1 * ((hilite_dir[0 + c][i - 1][j - 2] + hilite_dir[0 + c][i - 1][j - 1] + hilite_dir[0 + c][i - 1][j] + hilite_dir[0 + c][i - 1][j + 1] + hilite_dir[0 + c][i - 1][j + 2]) / + (hilite_dir[0 + 3][i - 1][j - 2] + hilite_dir[0 + 3][i - 1][j - 1] + hilite_dir[0 + 3][i - 1][j] + hilite_dir[0 + 3][i - 1][j + 1] + hilite_dir[0 + 3][i - 1][j + 2] + 0.00001)); } } } + for (int j = 2; j < hfw - 2; j++) { + if(hilite[3][hfh - 2][j] <= 0.01) { + hilite_dir[4 + c][hfh - 1][j] += hilite_dir[0 + c][hfh - 2][j]; + } + } + + } + if(plistener) { progress += 0.05; plistener->setProgress(progress); @@ -648,17 +805,29 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b for (int j = 2; j < hfw - 2; j++) { //from bottom if (hilite[3][i][j] > 0.01) { - for (int c = 0; c < 4; c++) { - hilite_dir[12 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; - } + hilite_dir[4 + 3][i][j] = 1.f; } else { - for (int c = 0; c < 4; c++) { - hilite_dir[12 + c][i][j] = 0.1 * ((hilite_dir[12 + c][(i + 1)][(j - 2)] + hilite_dir[12 + c][(i + 1)][(j - 1)] + hilite_dir[12 + c][(i + 1)][(j)] + hilite_dir[12 + c][(i + 1)][(j + 1)] + hilite_dir[12 + c][(i + 1)][(j + 2)]) / - (hilite_dir[12 + 3][(i + 1)][(j - 2)] + hilite_dir[12 + 3][(i + 1)][(j - 1)] + hilite_dir[12 + 3][(i + 1)][(j)] + hilite_dir[12 + 3][(i + 1)][(j + 1)] + hilite_dir[12 + 3][(i + 1)][(j + 2)] + 0.00001)); + hilite_dir[4 + 3][i][j] = (hilite_dir[4 + 3][(i + 1)][(j - 2)] + hilite_dir[4 + 3][(i + 1)][(j - 1)] + hilite_dir[4 + 3][(i + 1)][(j)] + hilite_dir[4 + 3][(i + 1)][(j + 1)] + hilite_dir[4 + 3][(i + 1)][(j + 2)]) == 0.f ? 0.f : 0.1f; + } + } + +#ifdef _OPENMP + #pragma omp parallel for +#endif + + for (int c = 0; c < 4; c++) { + for (int i = hfh - 2; i > 0; i--) { + for (int j = 2; j < hfw - 2; j++) { + //from bottom + if (hilite[3][i][j] > 0.01) { + hilite_dir[4 + c][i][j] = hilite[c][i][j] / hilite[3][i][j]; + } else { + hilite_dir[4 + c][i][j] = 0.1 * ((hilite_dir[4 + c][(i + 1)][(j - 2)] + hilite_dir[4 + c][(i + 1)][(j - 1)] + hilite_dir[4 + c][(i + 1)][(j)] + hilite_dir[4 + c][(i + 1)][(j + 1)] + hilite_dir[4 + c][(i + 1)][(j + 2)]) / + (hilite_dir[4 + 3][(i + 1)][(j - 2)] + hilite_dir[4 + 3][(i + 1)][(j - 1)] + hilite_dir[4 + 3][(i + 1)][(j)] + hilite_dir[4 + 3][(i + 1)][(j + 1)] + hilite_dir[4 + 3][(i + 1)][(j + 2)] + 0.00001)); } } - } + } if(plistener) { progress += 0.05; @@ -666,7 +835,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } //fill in edges - for (int dir = 0; dir < numdirs; dir++) { + for (int dir = 0; dir < 2; dir++) { for (int i = 1; i < hfh - 1; i++) for (int c = 0; c < 4; c++) { hilite_dir[dir * 4 + c][i][0] = hilite_dir[dir * 4 + c][i][1]; @@ -687,27 +856,63 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } } + for (int i = 1; i < hfh - 1; i++) + for (int c = 0; c < 4; c++) { + hilite_dir0[c][0][i] = hilite_dir0[c][1][i]; + hilite_dir0[c][hfw - 1][i] = hilite_dir0[c][hfw - 2][i]; + } + + for (int j = 1; j < hfw - 1; j++) + for (int c = 0; c < 4; c++) { + hilite_dir0[c][j][0] = hilite_dir0[c][j][1]; + hilite_dir0[c][j][hfh - 1] = hilite_dir0[c][j][hfh - 2]; + } + + for (int c = 0; c < 4; c++) { + hilite_dir0[c][0][0] = hilite_dir0[c][0][1] = hilite_dir0[c][1][0] = hilite_dir0[c][1][1] = hilite_dir0[c][2][2]; + hilite_dir0[c][hfw - 1][0] = hilite_dir0[c][hfw - 1][1] = hilite_dir0[c][hfw - 2][0] = hilite_dir0[c][hfw - 2][1] = hilite_dir0[c][hfw - 3][2]; + hilite_dir0[c][0][hfh - 1] = hilite_dir0[c][0][hfh - 2] = hilite_dir0[c][1][hfh - 1] = hilite_dir0[c][1][hfh - 2] = hilite_dir0[c][2][hfh - 3]; + hilite_dir0[c][hfw - 1][hfh - 1] = hilite_dir0[c][hfw - 1][hfh - 2] = hilite_dir0[c][hfw - 2][hfh - 1] = hilite_dir0[c][hfw - 2][hfh - 2] = hilite_dir0[c][hfw - 3][hfh - 3]; + } + + for (int i = 1; i < hfh - 1; i++) + for (int c = 0; c < 4; c++) { + hilite_dir4[c][0][i] = hilite_dir4[c][1][i]; + hilite_dir4[c][hfw - 1][i] = hilite_dir4[c][hfw - 2][i]; + } + + for (int j = 1; j < hfw - 1; j++) + for (int c = 0; c < 4; c++) { + hilite_dir4[c][j][0] = hilite_dir4[c][j][1]; + hilite_dir4[c][j][hfh - 1] = hilite_dir4[c][j][hfh - 2]; + } + + for (int c = 0; c < 4; c++) { + hilite_dir4[c][0][0] = hilite_dir4[c][0][1] = hilite_dir4[c][1][0] = hilite_dir4[c][1][1] = hilite_dir4[c][2][2]; + hilite_dir4[c][hfw - 1][0] = hilite_dir4[c][hfw - 1][1] = hilite_dir4[c][hfw - 2][0] = hilite_dir4[c][hfw - 2][1] = hilite_dir4[c][hfw - 3][2]; + hilite_dir4[c][0][hfh - 1] = hilite_dir4[c][0][hfh - 2] = hilite_dir4[c][1][hfh - 1] = hilite_dir4[c][1][hfh - 2] = hilite_dir4[c][2][hfh - 3]; + hilite_dir4[c][hfw - 1][hfh - 1] = hilite_dir4[c][hfw - 1][hfh - 2] = hilite_dir4[c][hfw - 2][hfh - 1] = hilite_dir4[c][hfw - 2][hfh - 2] = hilite_dir4[c][hfw - 3][hfh - 3]; + } + + + if(plistener) { progress += 0.05; plistener->setProgress(progress); } - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - /*for (int m=0; m<4*numdirs; m++) { - boxblur2(hilite_dir[m],hilite_dir[m],hfh,hfw,4); - }*/ - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - //now we have highlight data extended in various directions - //next step is to build back local data by averaging - - //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% // now reconstruct clipped channels using color ratios - //const float Yclip = 0.3333*(max[0] + max[1] + max[2]); - //float sumwt=0, counts=0; + for(int c = 0; c < 4; c++) { + hilite[c].free(); + } + + LUTf invfn(0x10000); + + for (int i = 0; i < 0x10000; i++) { + invfn[i] = 1.0 / (1.0 + i); + } #ifdef _OPENMP #pragma omp parallel for schedule(dynamic,16) @@ -725,8 +930,6 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b continue; //pixel not clipped } - //if (pixel[0] 0.5) { + dirwt = invfn[65535 * (SQR(rgb_blend[0] / Y - hilite_dir0[0][j1][i1] / Yhi) + + SQR(rgb_blend[1] / Y - hilite_dir0[1][j1][i1] / Yhi) + + SQR(rgb_blend[2] / Y - hilite_dir0[2][j1][i1] / Yhi))]; + totwt += dirwt; + clipfix[0] += dirwt * hilite_dir0[0][j1][i1] / (hilite_dir0[3][j1][i1] + 0.00001); + clipfix[1] += dirwt * hilite_dir0[1][j1][i1] / (hilite_dir0[3][j1][i1] + 0.00001); + clipfix[2] += dirwt * hilite_dir0[2][j1][i1] / (hilite_dir0[3][j1][i1] + 0.00001); + } + + for (int dir = 0; dir < 2; dir++) { float Yhi = 0.001 + (hilite_dir[dir * 4 + 0][i1][j1] + hilite_dir[dir * 4 + 1][i1][j1] + hilite_dir[dir * 4 + 2][i1][j1]); float Y = 0.001 + (rgb_blend[0] + rgb_blend[1] + rgb_blend[2]); @@ -827,6 +1043,19 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } } + Yhi = 0.001 + (hilite_dir4[0][j1][i1] + hilite_dir4[1][j1][i1] + hilite_dir4[2][j1][i1]); + Y = 0.001 + (rgb_blend[0] + rgb_blend[1] + rgb_blend[2]); + + if (hilite_dir4[0][j1][i1] + hilite_dir4[1][j1][i1] + hilite_dir4[2][j1][i1] > 0.5) { + dirwt = invfn[65535 * (SQR(rgb_blend[0] / Y - hilite_dir4[0][j1][i1] / Yhi) + + SQR(rgb_blend[1] / Y - hilite_dir4[1][j1][i1] / Yhi) + + SQR(rgb_blend[2] / Y - hilite_dir4[2][j1][i1] / Yhi))]; + totwt += dirwt; + clipfix[0] += dirwt * hilite_dir4[0][j1][i1] / (hilite_dir4[3][j1][i1] + 0.00001); + clipfix[1] += dirwt * hilite_dir4[1][j1][i1] / (hilite_dir4[3][j1][i1] + 0.00001); + clipfix[2] += dirwt * hilite_dir4[2][j1][i1] / (hilite_dir4[3][j1][i1] + 0.00001); + } + if(totwt == 0.f) { continue; } @@ -841,7 +1070,7 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b if (pixel[0] > max_f[0] && pixel[1] > max_f[1] && pixel[2] > max_f[2]) { //all channels clipped float Y = (0.299 * clipfix[0] + 0.587 * clipfix[1] + 0.114 * clipfix[2]); - //float Y = (clipfix[0] + clipfix[1] + clipfix[2]); + factor = whitept / Y; red[i][j] = clipfix[0] * factor; green[i][j] = clipfix[1] * factor; @@ -866,28 +1095,11 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b } } - /*if (hilite[3][i1][j1]>0.01) { - red[i][j] = (red[i][j] + hilite[0][i1][j1])/(1+hilite[3][i1][j1]); - green[i][j] = (green[i][j]+ hilite[1][i1][j1])/(1+hilite[3][i1][j1]); - blue[i][j] = (blue[i][j] + hilite[2][i1][j1])/(1+hilite[3][i1][j1]); - }*/ - Y = (0.299 * red[i][j] + 0.587 * green[i][j] + 0.114 * blue[i][j]); if (Y > whitept) { factor = whitept / Y; - /*I = (0.596 * red[i][j] - 0.275 * green[i][j] - 0.321 * blue[i][j]); - Q = (0.212 * red[i][j] - 0.523 * green[i][j] + 0.311 * blue[i][j]); - - Y *= factor; - I *= factor;//max(0,min(1,(whitept-Y)/(whitept-clippt))); - Q *= factor;//max(0,min(1,(whitept-Y)/(whitept-clippt))); - - red[i][j] = Y + 0.956*I + 0.621*Q; - green[i][j] = Y - 0.272*I - 0.647*Q; - blue[i][j] = Y - 1.105*I + 1.702*Q;*/ - red[i][j] *= factor; green[i][j] *= factor; blue[i][j] *= factor; @@ -899,76 +1111,8 @@ void RawImageSource :: HLRecovery_inpaint (float** red, float** green, float** b plistener->setProgress(1.00); } - //printf("ave wt=%f\n",sumwt/counts); - - - // diagnostic output - /*for (int i=0; i & hfsize, multi_array2D & hilite, int box );