liz/rawTherapee
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

rcd_demosaic: further speedup

Browse Source
This commit is contained in:
heckflosse
2018-02-28 14:53:02 +01:00
parent bad28bb0ad
commit 2fb4403823
1 changed files with 57 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files

114
rtengine/rcd_demosaic.cc
View File

@@ -51,7 +51,7 @@ void RawImageSource::rcd_demosaic()
const int width = W, height = H; const int width = W, height = H;
constexpr int tileBorder = 8; constexpr int tileBorder = 8;
constexpr int tileSize = 228; constexpr int tileSize = 220;
constexpr int tileSizeN = tileSize - 2 * tileBorder; constexpr int tileSizeN = tileSize - 2 * tileBorder;
const int numTh = H / (tileSizeN) + ((H % (tileSizeN)) ? 1 : 0); const int numTh = H / (tileSizeN) + ((H % (tileSizeN)) ? 1 : 0);
const int numTw = W / (tileSizeN) + ((W % (tileSizeN)) ? 1 : 0); const int numTw = W / (tileSizeN) + ((W % (tileSizeN)) ? 1 : 0);
@@ -64,10 +64,10 @@ void RawImageSource::rcd_demosaic()
#pragma omp parallel #pragma omp parallel
#endif #endif
{ {
float *cfa = (float*) calloc( tileSize * tileSize, sizeof *cfa ); float *cfa = (float*) calloc(tileSize * tileSize, sizeof *cfa);
float (*rgb)[tileSize * tileSize] = (float (*)[tileSize * tileSize])malloc(3 * sizeof *rgb); float (*rgb)[tileSize * tileSize] = (float (*)[tileSize * tileSize])malloc(3 * sizeof *rgb);
float *VH_Dir = (float*) calloc( tileSize * tileSize, sizeof *VH_Dir ); float *VH_Dir = (float*) calloc(tileSize * tileSize, sizeof *VH_Dir);
float *PQ_Dir = (float*) calloc( tileSize * tileSize, sizeof *PQ_Dir ); float *PQ_Dir = (float*) calloc(tileSize * tileSize, sizeof *PQ_Dir);
float *lpf = PQ_Dir; // reuse buffer, they don't overlap in usage float *lpf = PQ_Dir; // reuse buffer, they don't overlap in usage
#ifdef _OPENMP #ifdef _OPENMP
@@ -99,12 +99,13 @@ void RawImageSource::rcd_demosaic()
* STEP 1: Find cardinal and diagonal interpolation directions * STEP 1: Find cardinal and diagonal interpolation directions
*/ */
for (int row = 4; row < tileSize - 4; row++ ) { for (int row = 4; row < tileSize - 4; row++) {
for (int col = 4, indx = row * tileSize + col; col < tileSize - 4; col++, indx++ ) { for (int col = 4, indx = row * tileSize + col; col < tileSize - 4; col++, indx++) {
const float cfai = cfa[indx];
//Calculate h/v local discrimination //Calculate h/v local discrimination
float V_Stat = max(epssq, - 18.0f * cfa[indx] * cfa[indx - w1] - 18.0f * cfa[indx] * cfa[indx + w1] - 36.0f * cfa[indx] * cfa[indx - w2] - 36.0f * cfa[indx] * cfa[indx + w2] + 18.0f * cfa[indx] * cfa[indx - w3] + 18.0f * cfa[indx] * cfa[indx + w3] - 2.0f * cfa[indx] * cfa[indx - w4] - 2.0f * cfa[indx] * cfa[indx + w4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx - w1] * cfa[indx + w1] - 12.0f * cfa[indx - w1] * cfa[indx - w2] + 24.0f * cfa[indx - w1] * cfa[indx + w2] - 38.0f * cfa[indx - w1] * cfa[indx - w3] + 16.0f * cfa[indx - w1] * cfa[indx + w3] + 12.0f * cfa[indx - w1] * cfa[indx - w4] - 6.0f * cfa[indx - w1] * cfa[indx + w4] + 46.0f * cfa[indx - w1] * cfa[indx - w1] + 24.0f * cfa[indx + w1] * cfa[indx - w2] - 12.0f * cfa[indx + w1] * cfa[indx + w2] + 16.0f * cfa[indx + w1] * cfa[indx - w3] - 38.0f * cfa[indx + w1] * cfa[indx + w3] - 6.0f * cfa[indx + w1] * cfa[indx - w4] + 12.0f * cfa[indx + w1] * cfa[indx + w4] + 46.0f * cfa[indx + w1] * cfa[indx + w1] + 14.0f * cfa[indx - w2] * cfa[indx + w2] - 12.0f * cfa[indx - w2] * cfa[indx + w3] - 2.0f * cfa[indx - w2] * cfa[indx - w4] + 2.0f * cfa[indx - w2] * cfa[indx + w4] + 11.0f * cfa[indx - w2] * cfa[indx - w2] - 12.0f * cfa[indx + w2] * cfa[indx - w3] + 2.0f * cfa[indx + w2] * cfa[indx - w4] - 2.0f * cfa[indx + w2] * cfa[indx + w4] + 11.0f * cfa[indx + w2] * cfa[indx + w2] + 2.0f * cfa[indx - w3] * cfa[indx + w3] - 6.0f * cfa[indx - w3] * cfa[indx - w4] + 10.0f * cfa[indx - w3] * cfa[indx - w3] - 6.0f * cfa[indx + w3] * cfa[indx + w4] + 10.0f * cfa[indx + w3] * cfa[indx + w3] + 1.0f * cfa[indx - w4] * cfa[indx - w4] + 1.0f * cfa[indx + w4] * cfa[indx + w4]); float V_Stat = max(epssq, - 18.f * cfai * (cfa[indx - w1] + cfa[indx + w1] + 2.f * (cfa[indx - w2] + cfa[indx + w2]) - cfa[indx - w3] - cfa[indx + w3]) - 2.f * cfai * (cfa[indx - w4] + cfa[indx + w4] - 19.f * cfai) - cfa[indx - w1] * (70.f * cfa[indx + w1] + 12.f * cfa[indx - w2] - 24.f * cfa[indx + w2] + 38.f * cfa[indx - w3] - 16.f * cfa[indx + w3] - 12.f * cfa[indx - w4] + 6.f * cfa[indx + w4] - 46.f * cfa[indx - w1]) + cfa[indx + w1] * (24.f * cfa[indx - w2] - 12.f * cfa[indx + w2] + 16.f * cfa[indx - w3] - 38.f * cfa[indx + w3] - 6.f * cfa[indx - w4] + 12.f * cfa[indx + w4] + 46.f * cfa[indx + w1]) + cfa[indx - w2] * (14.f * cfa[indx + w2] - 12.f * cfa[indx + w3] - 2.f * cfa[indx - w4] + 2.f * cfa[indx + w4] + 11.f * cfa[indx - w2]) + cfa[indx + w2] * (-12.f * cfa[indx - w3] + 2.f * cfa[indx - w4] - 2.f * cfa[indx + w4] + 11.f * cfa[indx + w2]) + cfa[indx - w3] * (2.f * cfa[indx + w3] - 6.f * cfa[indx - w4] + 10.f * cfa[indx - w3]) + cfa[indx + w3] * (-6.f * cfa[indx + w4] + 10.f * cfa[indx + w3]) + cfa[indx - w4] * cfa[indx - w4] + cfa[indx + w4] * cfa[indx + w4]);
float H_Stat = max(epssq, - 18.0f * cfa[indx] * cfa[indx - 1] - 18.0f * cfa[indx] * cfa[indx + 1] - 36.0f * cfa[indx] * cfa[indx - 2] - 36.0f * cfa[indx] * cfa[indx + 2] + 18.0f * cfa[indx] * cfa[indx - 3] + 18.0f * cfa[indx] * cfa[indx + 3] - 2.0f * cfa[indx] * cfa[indx - 4] - 2.0f * cfa[indx] * cfa[indx + 4] + 38.0f * cfa[indx] * cfa[indx] - 70.0f * cfa[indx - 1] * cfa[indx + 1] - 12.0f * cfa[indx - 1] * cfa[indx - 2] + 24.0f * cfa[indx - 1] * cfa[indx + 2] - 38.0f * cfa[indx - 1] * cfa[indx - 3] + 16.0f * cfa[indx - 1] * cfa[indx + 3] + 12.0f * cfa[indx - 1] * cfa[indx - 4] - 6.0f * cfa[indx - 1] * cfa[indx + 4] + 46.0f * cfa[indx - 1] * cfa[indx - 1] + 24.0f * cfa[indx + 1] * cfa[indx - 2] - 12.0f * cfa[indx + 1] * cfa[indx + 2] + 16.0f * cfa[indx + 1] * cfa[indx - 3] - 38.0f * cfa[indx + 1] * cfa[indx + 3] - 6.0f * cfa[indx + 1] * cfa[indx - 4] + 12.0f * cfa[indx + 1] * cfa[indx + 4] + 46.0f * cfa[indx + 1] * cfa[indx + 1] + 14.0f * cfa[indx - 2] * cfa[indx + 2] - 12.0f * cfa[indx - 2] * cfa[indx + 3] - 2.0f * cfa[indx - 2] * cfa[indx - 4] + 2.0f * cfa[indx - 2] * cfa[indx + 4] + 11.0f * cfa[indx - 2] * cfa[indx - 2] - 12.0f * cfa[indx + 2] * cfa[indx - 3] + 2.0f * cfa[indx + 2] * cfa[indx - 4] - 2.0f * cfa[indx + 2] * cfa[indx + 4] + 11.0f * cfa[indx + 2] * cfa[indx + 2] + 2.0f * cfa[indx - 3] * cfa[indx + 3] - 6.0f * cfa[indx - 3] * cfa[indx - 4] + 10.0f * cfa[indx - 3] * cfa[indx - 3] - 6.0f * cfa[indx + 3] * cfa[indx + 4] + 10.0f * cfa[indx + 3] * cfa[indx + 3] + 1.0f * cfa[indx - 4] * cfa[indx - 4] + 1.0f * cfa[indx + 4] * cfa[indx + 4]); float H_Stat = max(epssq, - 18.f * cfai * (cfa[indx - 1] + cfa[indx + 1] + 2.f * (cfa[indx - 2] + cfa[indx + 2]) - cfa[indx - 3] - cfa[indx + 3]) - 2.f * cfai * (cfa[indx - 4] + cfa[indx + 4] - 19.f * cfai) - cfa[indx - 1] * (70.f * cfa[indx + 1] + 12.f * cfa[indx - 2] - 24.f * cfa[indx + 2] + 38.f * cfa[indx - 3] - 16.f * cfa[indx + 3] - 12.f * cfa[indx - 4] + 6.f * cfa[indx + 4] - 46.f * cfa[indx - 1]) + cfa[indx + 1] * (24.f * cfa[indx - 2] - 12.f * cfa[indx + 2] + 16.f * cfa[indx - 3] - 38.f * cfa[indx + 3] - 6.f * cfa[indx - 4] + 12.f * cfa[indx + 4] + 46.f * cfa[indx + 1]) + cfa[indx - 2] * (14.f * cfa[indx + 2] - 12.f * cfa[indx + 3] - 2.f * cfa[indx - 4] + 2.f * cfa[indx + 4] + 11.f * cfa[indx - 2]) + cfa[indx + 2] * (-12.f * cfa[indx - 3] + 2.f * cfa[indx - 4] - 2.f * cfa[indx + 4] + 11.f * cfa[indx + 2]) + cfa[indx - 3] * (2.f * cfa[indx + 3] - 6.f * cfa[indx - 4] + 10.f * cfa[indx - 3]) + cfa[indx + 3] * (-6.f * cfa[indx + 4] + 10.f * cfa[indx + 3]) + cfa[indx - 4] * cfa[indx - 4] + cfa[indx + 4] * cfa[indx + 4]);
VH_Dir[indx] = V_Stat / (V_Stat + H_Stat); VH_Dir[indx] = V_Stat / (V_Stat + H_Stat);
} }
@@ -115,9 +116,9 @@ void RawImageSource::rcd_demosaic()
*/ */
// Step 2.1: Low pass filter incorporating green, red and blue local samples from the raw data // Step 2.1: Low pass filter incorporating green, red and blue local samples from the raw data
for ( int row = 2; row < tileSize - 2; row++ ) { for (int row = 2; row < tileSize - 2; row++) {
for ( int col = 2 + (FC( row, 0 ) & 1), indx = row * tileSize + col; col < tileSize - 2; col += 2, indx += 2 ) { for (int col = 2 + (FC(row, 0) & 1), indx = row * tileSize + col; col < tileSize - 2; col += 2, indx += 2) {
lpf[indx>>1] = 0.25f * cfa[indx] + 0.125f * ( cfa[indx - w1] + cfa[indx + w1] + cfa[indx - 1] + cfa[indx + 1] ) + 0.0625f * ( cfa[indx - w1 - 1] + cfa[indx - w1 + 1] + cfa[indx + w1 - 1] + cfa[indx + w1 + 1] ); lpf[indx>>1] = 0.25f * cfa[indx] + 0.125f * (cfa[indx - w1] + cfa[indx + w1] + cfa[indx - 1] + cfa[indx + 1]) + 0.0625f * (cfa[indx - w1 - 1] + cfa[indx - w1 + 1] + cfa[indx + w1 - 1] + cfa[indx + w1 + 1]);
} }
} }
@@ -125,8 +126,8 @@ void RawImageSource::rcd_demosaic()
* STEP 3: Populate the green channel * STEP 3: Populate the green channel
*/ */
// Step 3.1: Populate the green channel at blue and red CFA positions // Step 3.1: Populate the green channel at blue and red CFA positions
for ( int row = 4; row < tileSize - 4; row++ ) { for (int row = 4; row < tileSize - 4; row++) {
for ( int col = 4 + (FC( row, 0 ) & 1), indx = row * tileSize + col; col < tileSize - 4; col += 2, indx += 2 ) { for (int col = 4 + (FC(row, 0) & 1), indx = row * tileSize + col; col < tileSize - 4; col += 2, indx += 2) {
// Refined vertical and horizontal local discrimination // Refined vertical and horizontal local discrimination
float VH_Central_Value = VH_Dir[indx]; float VH_Central_Value = VH_Dir[indx];
@@ -135,23 +136,23 @@ void RawImageSource::rcd_demosaic()
float VH_Disc = std::fabs(0.5f - VH_Central_Value) < std::fabs(0.5f - VH_Neighbourhood_Value) ? VH_Neighbourhood_Value : VH_Central_Value; float VH_Disc = std::fabs(0.5f - VH_Central_Value) < std::fabs(0.5f - VH_Neighbourhood_Value) ? VH_Neighbourhood_Value : VH_Central_Value;
// Cardinal gradients // Cardinal gradients
float N_Grad = eps + std::fabs( cfa[indx - w1] - cfa[indx + w1] ) + std::fabs( cfa[indx] - cfa[indx - w2] ) + std::fabs( cfa[indx - w1] - cfa[indx - w3] ) + std::fabs( cfa[indx - w2] - cfa[indx - w4] ); float N_Grad = eps + std::fabs(cfa[indx - w1] - cfa[indx + w1]) + std::fabs(cfa[indx] - cfa[indx - w2]) + std::fabs(cfa[indx - w1] - cfa[indx - w3]) + std::fabs(cfa[indx - w2] - cfa[indx - w4]);
float S_Grad = eps + std::fabs( cfa[indx - w1] - cfa[indx + w1] ) + std::fabs( cfa[indx] - cfa[indx + w2] ) + std::fabs( cfa[indx + w1] - cfa[indx + w3] ) + std::fabs( cfa[indx + w2] - cfa[indx + w4] ); float S_Grad = eps + std::fabs(cfa[indx - w1] - cfa[indx + w1]) + std::fabs(cfa[indx] - cfa[indx + w2]) + std::fabs(cfa[indx + w1] - cfa[indx + w3]) + std::fabs(cfa[indx + w2] - cfa[indx + w4]);
float W_Grad = eps + std::fabs( cfa[indx - 1] - cfa[indx + 1] ) + std::fabs( cfa[indx] - cfa[indx - 2] ) + std::fabs( cfa[indx - 1] - cfa[indx - 3] ) + std::fabs( cfa[indx - 2] - cfa[indx - 4] ); float W_Grad = eps + std::fabs(cfa[indx - 1] - cfa[indx + 1]) + std::fabs(cfa[indx] - cfa[indx - 2]) + std::fabs(cfa[indx - 1] - cfa[indx - 3]) + std::fabs(cfa[indx - 2] - cfa[indx - 4]);
float E_Grad = eps + std::fabs( cfa[indx - 1] - cfa[indx + 1] ) + std::fabs( cfa[indx] - cfa[indx + 2] ) + std::fabs( cfa[indx + 1] - cfa[indx + 3] ) + std::fabs( cfa[indx + 2] - cfa[indx + 4] ); float E_Grad = eps + std::fabs(cfa[indx - 1] - cfa[indx + 1]) + std::fabs(cfa[indx] - cfa[indx + 2]) + std::fabs(cfa[indx + 1] - cfa[indx + 3]) + std::fabs(cfa[indx + 2] - cfa[indx + 4]);
// Cardinal pixel estimations // Cardinal pixel estimations
float N_Est = cfa[indx - w1] * ( 1.f + ( lpf[indx>>1] - lpf[(indx - w2)>>1] ) / ( eps + lpf[indx>>1] + lpf[(indx - w2)>>1] ) ); float N_Est = cfa[indx - w1] * (1.f + (lpf[indx>>1] - lpf[(indx - w2)>>1]) / (eps + lpf[indx>>1] + lpf[(indx - w2)>>1]));
float S_Est = cfa[indx + w1] * ( 1.f + ( lpf[indx>>1] - lpf[(indx + w2)>>1] ) / ( eps + lpf[indx>>1] + lpf[(indx + w2)>>1] ) ); float S_Est = cfa[indx + w1] * (1.f + (lpf[indx>>1] - lpf[(indx + w2)>>1]) / (eps + lpf[indx>>1] + lpf[(indx + w2)>>1]));
float W_Est = cfa[indx - 1] * ( 1.f + ( lpf[indx>>1] - lpf[(indx - 2)>>1] ) / ( eps + lpf[indx>>1] + lpf[(indx - 2)>>1] ) ); float W_Est = cfa[indx - 1] * (1.f + (lpf[indx>>1] - lpf[(indx - 2)>>1]) / (eps + lpf[indx>>1] + lpf[(indx - 2)>>1]));
float E_Est = cfa[indx + 1] * ( 1.f + ( lpf[indx>>1] - lpf[(indx + 2)>>1] ) / ( eps + lpf[indx>>1] + lpf[(indx + 2)>>1] ) ); float E_Est = cfa[indx + 1] * (1.f + (lpf[indx>>1] - lpf[(indx + 2)>>1]) / (eps + lpf[indx>>1] + lpf[(indx + 2)>>1]));
// Vertical and horizontal estimations // Vertical and horizontal estimations
float V_Est = ( S_Grad * N_Est + N_Grad * S_Est ) / (N_Grad + S_Grad ); float V_Est = (S_Grad * N_Est + N_Grad * S_Est) / (N_Grad + S_Grad);
float H_Est = ( W_Grad * E_Est + E_Grad * W_Est ) / (E_Grad + W_Grad ); float H_Est = (W_Grad * E_Est + E_Grad * W_Est) / (E_Grad + W_Grad);
// G@B and G@R interpolation // G@B and G@R interpolation
rgb[1][indx] = VH_Disc * H_Est + ( 1.f - VH_Disc ) * V_Est; rgb[1][indx] = VH_Disc * H_Est + (1.f - VH_Disc) * V_Est;
} }
} }
@@ -160,32 +161,32 @@ void RawImageSource::rcd_demosaic()
*/ */
// Step 4.1: Calculate P/Q diagonal local discrimination // Step 4.1: Calculate P/Q diagonal local discrimination
for ( int row = 4; row < tileSize - 4; row++ ) { for (int row = 4; row < tileSize - 4; row++) {
for ( int col = 4 + (FC( row, 0 ) & 1), indx = row * tileSize + col; col < tileSize - 4; col += 2, indx += 2 ) { for (int col = 4 + (FC(row, 0) & 1), indx = row * tileSize + col; col < tileSize - 4; col += 2, indx += 2) {
float P_Stat = max( - 18.f * cfa[indx] * cfa[indx - w1 - 1] - 18.f * cfa[indx] * cfa[indx + w1 + 1] - 36.f * cfa[indx] * cfa[indx - w2 - 2] - 36.f * cfa[indx] * cfa[indx + w2 + 2] + 18.f * cfa[indx] * cfa[indx - w3 - 3] + 18.f * cfa[indx] * cfa[indx + w3 + 3] - 2.f * cfa[indx] * cfa[indx - w4 - 4] - 2.f * cfa[indx] * cfa[indx + w4 + 4] + 38.f * cfa[indx] * cfa[indx] - 70.f * cfa[indx - w1 - 1] * cfa[indx + w1 + 1] - 12.f * cfa[indx - w1 - 1] * cfa[indx - w2 - 2] + 24.f * cfa[indx - w1 - 1] * cfa[indx + w2 + 2] - 38.f * cfa[indx - w1 - 1] * cfa[indx - w3 - 3] + 16.f * cfa[indx - w1 - 1] * cfa[indx + w3 + 3] + 12.f * cfa[indx - w1 - 1] * cfa[indx - w4 - 4] - 6.f * cfa[indx - w1 - 1] * cfa[indx + w4 + 4] + 46.f * cfa[indx - w1 - 1] * cfa[indx - w1 - 1] + 24.f * cfa[indx + w1 + 1] * cfa[indx - w2 - 2] - 12.f * cfa[indx + w1 + 1] * cfa[indx + w2 + 2] + 16.f * cfa[indx + w1 + 1] * cfa[indx - w3 - 3] - 38.f * cfa[indx + w1 + 1] * cfa[indx + w3 + 3] - 6.f * cfa[indx + w1 + 1] * cfa[indx - w4 - 4] + 12.f * cfa[indx + w1 + 1] * cfa[indx + w4 + 4] + 46.f * cfa[indx + w1 + 1] * cfa[indx + w1 + 1] + 14.f * cfa[indx - w2 - 2] * cfa[indx + w2 + 2] - 12.f * cfa[indx - w2 - 2] * cfa[indx + w3 + 3] - 2.f * cfa[indx - w2 - 2] * cfa[indx - w4 - 4] + 2.f * cfa[indx - w2 - 2] * cfa[indx + w4 + 4] + 11.f * cfa[indx - w2 - 2] * cfa[indx - w2 - 2] - 12.f * cfa[indx + w2 + 2] * cfa[indx - w3 - 3] + 2 * cfa[indx + w2 + 2] * cfa[indx - w4 - 4] - 2.f * cfa[indx + w2 + 2] * cfa[indx + w4 + 4] + 11.f * cfa[indx + w2 + 2] * cfa[indx + w2 + 2] + 2.f * cfa[indx - w3 - 3] * cfa[indx + w3 + 3] - 6.f * cfa[indx - w3 - 3] * cfa[indx - w4 - 4] + 10.f * cfa[indx - w3 - 3] * cfa[indx - w3 - 3] - 6.f * cfa[indx + w3 + 3] * cfa[indx + w4 + 4] + 10.f * cfa[indx + w3 + 3] * cfa[indx + w3 + 3] + 1.f * cfa[indx - w4 - 4] * cfa[indx - w4 - 4] + 1.f * cfa[indx + w4 + 4] * cfa[indx + w4 + 4], epssq ); float P_Stat = max(- 18.f * cfa[indx] * cfa[indx - w1 - 1] - 18.f * cfa[indx] * cfa[indx + w1 + 1] - 36.f * cfa[indx] * cfa[indx - w2 - 2] - 36.f * cfa[indx] * cfa[indx + w2 + 2] + 18.f * cfa[indx] * cfa[indx - w3 - 3] + 18.f * cfa[indx] * cfa[indx + w3 + 3] - 2.f * cfa[indx] * cfa[indx - w4 - 4] - 2.f * cfa[indx] * cfa[indx + w4 + 4] + 38.f * cfa[indx] * cfa[indx] - 70.f * cfa[indx - w1 - 1] * cfa[indx + w1 + 1] - 12.f * cfa[indx - w1 - 1] * cfa[indx - w2 - 2] + 24.f * cfa[indx - w1 - 1] * cfa[indx + w2 + 2] - 38.f * cfa[indx - w1 - 1] * cfa[indx - w3 - 3] + 16.f * cfa[indx - w1 - 1] * cfa[indx + w3 + 3] + 12.f * cfa[indx - w1 - 1] * cfa[indx - w4 - 4] - 6.f * cfa[indx - w1 - 1] * cfa[indx + w4 + 4] + 46.f * cfa[indx - w1 - 1] * cfa[indx - w1 - 1] + 24.f * cfa[indx + w1 + 1] * cfa[indx - w2 - 2] - 12.f * cfa[indx + w1 + 1] * cfa[indx + w2 + 2] + 16.f * cfa[indx + w1 + 1] * cfa[indx - w3 - 3] - 38.f * cfa[indx + w1 + 1] * cfa[indx + w3 + 3] - 6.f * cfa[indx + w1 + 1] * cfa[indx - w4 - 4] + 12.f * cfa[indx + w1 + 1] * cfa[indx + w4 + 4] + 46.f * cfa[indx + w1 + 1] * cfa[indx + w1 + 1] + 14.f * cfa[indx - w2 - 2] * cfa[indx + w2 + 2] - 12.f * cfa[indx - w2 - 2] * cfa[indx + w3 + 3] - 2.f * cfa[indx - w2 - 2] * cfa[indx - w4 - 4] + 2.f * cfa[indx - w2 - 2] * cfa[indx + w4 + 4] + 11.f * cfa[indx - w2 - 2] * cfa[indx - w2 - 2] - 12.f * cfa[indx + w2 + 2] * cfa[indx - w3 - 3] + 2 * cfa[indx + w2 + 2] * cfa[indx - w4 - 4] - 2.f * cfa[indx + w2 + 2] * cfa[indx + w4 + 4] + 11.f * cfa[indx + w2 + 2] * cfa[indx + w2 + 2] + 2.f * cfa[indx - w3 - 3] * cfa[indx + w3 + 3] - 6.f * cfa[indx - w3 - 3] * cfa[indx - w4 - 4] + 10.f * cfa[indx - w3 - 3] * cfa[indx - w3 - 3] - 6.f * cfa[indx + w3 + 3] * cfa[indx + w4 + 4] + 10.f * cfa[indx + w3 + 3] * cfa[indx + w3 + 3] + 1.f * cfa[indx - w4 - 4] * cfa[indx - w4 - 4] + 1.f * cfa[indx + w4 + 4] * cfa[indx + w4 + 4], epssq);
float Q_Stat = max( - 18.f * cfa[indx] * cfa[indx + w1 - 1] - 18.f * cfa[indx] * cfa[indx - w1 + 1] - 36.f * cfa[indx] * cfa[indx + w2 - 2] - 36.f * cfa[indx] * cfa[indx - w2 + 2] + 18.f * cfa[indx] * cfa[indx + w3 - 3] + 18.f * cfa[indx] * cfa[indx - w3 + 3] - 2.f * cfa[indx] * cfa[indx + w4 - 4] - 2.f * cfa[indx] * cfa[indx - w4 + 4] + 38.f * cfa[indx] * cfa[indx] - 70.f * cfa[indx + w1 - 1] * cfa[indx - w1 + 1] - 12.f * cfa[indx + w1 - 1] * cfa[indx + w2 - 2] + 24.f * cfa[indx + w1 - 1] * cfa[indx - w2 + 2] - 38.f * cfa[indx + w1 - 1] * cfa[indx + w3 - 3] + 16.f * cfa[indx + w1 - 1] * cfa[indx - w3 + 3] + 12.f * cfa[indx + w1 - 1] * cfa[indx + w4 - 4] - 6.f * cfa[indx + w1 - 1] * cfa[indx - w4 + 4] + 46.f * cfa[indx + w1 - 1] * cfa[indx + w1 - 1] + 24.f * cfa[indx - w1 + 1] * cfa[indx + w2 - 2] - 12.f * cfa[indx - w1 + 1] * cfa[indx - w2 + 2] + 16.f * cfa[indx - w1 + 1] * cfa[indx + w3 - 3] - 38.f * cfa[indx - w1 + 1] * cfa[indx - w3 + 3] - 6.f * cfa[indx - w1 + 1] * cfa[indx + w4 - 4] + 12.f * cfa[indx - w1 + 1] * cfa[indx - w4 + 4] + 46.f * cfa[indx - w1 + 1] * cfa[indx - w1 + 1] + 14.f * cfa[indx + w2 - 2] * cfa[indx - w2 + 2] - 12.f * cfa[indx + w2 - 2] * cfa[indx - w3 + 3] - 2.f * cfa[indx + w2 - 2] * cfa[indx + w4 - 4] + 2.f * cfa[indx + w2 - 2] * cfa[indx - w4 + 4] + 11.f * cfa[indx + w2 - 2] * cfa[indx + w2 - 2] - 12.f * cfa[indx - w2 + 2] * cfa[indx + w3 - 3] + 2 * cfa[indx - w2 + 2] * cfa[indx + w4 - 4] - 2.f * cfa[indx - w2 + 2] * cfa[indx - w4 + 4] + 11.f * cfa[indx - w2 + 2] * cfa[indx - w2 + 2] + 2.f * cfa[indx + w3 - 3] * cfa[indx - w3 + 3] - 6.f * cfa[indx + w3 - 3] * cfa[indx + w4 - 4] + 10.f * cfa[indx + w3 - 3] * cfa[indx + w3 - 3] - 6.f * cfa[indx - w3 + 3] * cfa[indx - w4 + 4] + 10.f * cfa[indx - w3 + 3] * cfa[indx - w3 + 3] + 1.f * cfa[indx + w4 - 4] * cfa[indx + w4 - 4] + 1.f * cfa[indx - w4 + 4] * cfa[indx - w4 + 4], epssq ); float Q_Stat = max(- 18.f * cfa[indx] * cfa[indx + w1 - 1] - 18.f * cfa[indx] * cfa[indx - w1 + 1] - 36.f * cfa[indx] * cfa[indx + w2 - 2] - 36.f * cfa[indx] * cfa[indx - w2 + 2] + 18.f * cfa[indx] * cfa[indx + w3 - 3] + 18.f * cfa[indx] * cfa[indx - w3 + 3] - 2.f * cfa[indx] * cfa[indx + w4 - 4] - 2.f * cfa[indx] * cfa[indx - w4 + 4] + 38.f * cfa[indx] * cfa[indx] - 70.f * cfa[indx + w1 - 1] * cfa[indx - w1 + 1] - 12.f * cfa[indx + w1 - 1] * cfa[indx + w2 - 2] + 24.f * cfa[indx + w1 - 1] * cfa[indx - w2 + 2] - 38.f * cfa[indx + w1 - 1] * cfa[indx + w3 - 3] + 16.f * cfa[indx + w1 - 1] * cfa[indx - w3 + 3] + 12.f * cfa[indx + w1 - 1] * cfa[indx + w4 - 4] - 6.f * cfa[indx + w1 - 1] * cfa[indx - w4 + 4] + 46.f * cfa[indx + w1 - 1] * cfa[indx + w1 - 1] + 24.f * cfa[indx - w1 + 1] * cfa[indx + w2 - 2] - 12.f * cfa[indx - w1 + 1] * cfa[indx - w2 + 2] + 16.f * cfa[indx - w1 + 1] * cfa[indx + w3 - 3] - 38.f * cfa[indx - w1 + 1] * cfa[indx - w3 + 3] - 6.f * cfa[indx - w1 + 1] * cfa[indx + w4 - 4] + 12.f * cfa[indx - w1 + 1] * cfa[indx - w4 + 4] + 46.f * cfa[indx - w1 + 1] * cfa[indx - w1 + 1] + 14.f * cfa[indx + w2 - 2] * cfa[indx - w2 + 2] - 12.f * cfa[indx + w2 - 2] * cfa[indx - w3 + 3] - 2.f * cfa[indx + w2 - 2] * cfa[indx + w4 - 4] + 2.f * cfa[indx + w2 - 2] * cfa[indx - w4 + 4] + 11.f * cfa[indx + w2 - 2] * cfa[indx + w2 - 2] - 12.f * cfa[indx - w2 + 2] * cfa[indx + w3 - 3] + 2 * cfa[indx - w2 + 2] * cfa[indx + w4 - 4] - 2.f * cfa[indx - w2 + 2] * cfa[indx - w4 + 4] + 11.f * cfa[indx - w2 + 2] * cfa[indx - w2 + 2] + 2.f * cfa[indx + w3 - 3] * cfa[indx - w3 + 3] - 6.f * cfa[indx + w3 - 3] * cfa[indx + w4 - 4] + 10.f * cfa[indx + w3 - 3] * cfa[indx + w3 - 3] - 6.f * cfa[indx - w3 + 3] * cfa[indx - w4 + 4] + 10.f * cfa[indx - w3 + 3] * cfa[indx - w3 + 3] + 1.f * cfa[indx + w4 - 4] * cfa[indx + w4 - 4] + 1.f * cfa[indx - w4 + 4] * cfa[indx - w4 + 4], epssq);
PQ_Dir[indx] = P_Stat / ( P_Stat + Q_Stat ); PQ_Dir[indx] = P_Stat / (P_Stat + Q_Stat);
} }
} }
// Step 4.2: Populate the red and blue channels at blue and red CFA positions // Step 4.2: Populate the red and blue channels at blue and red CFA positions
for ( int row = 4; row < tileSize - 4; row++ ) { for (int row = 4; row < tileSize - 4; row++) {
for ( int col = 4 + (FC( row, 0 ) & 1), indx = row * tileSize + col, c = 2 - FC( row, col ); col < tileSize - 4; col += 2, indx += 2 ) { for (int col = 4 + (FC(row, 0) & 1), indx = row * tileSize + col, c = 2 - FC(row, col); col < tileSize - 4; col += 2, indx += 2) {
// Refined P/Q diagonal local discrimination // Refined P/Q diagonal local discrimination
float PQ_Central_Value = PQ_Dir[indx]; float PQ_Central_Value = PQ_Dir[indx];
float PQ_Neighbourhood_Value = 0.25f * ( PQ_Dir[indx - w1 - 1] + PQ_Dir[indx - w1 + 1] + PQ_Dir[indx + w1 - 1] + PQ_Dir[indx + w1 + 1] ); float PQ_Neighbourhood_Value = 0.25f * (PQ_Dir[indx - w1 - 1] + PQ_Dir[indx - w1 + 1] + PQ_Dir[indx + w1 - 1] + PQ_Dir[indx + w1 + 1]);
float PQ_Disc = ( std::fabs( 0.5f - PQ_Central_Value ) < std::fabs( 0.5f - PQ_Neighbourhood_Value ) ) ? PQ_Neighbourhood_Value : PQ_Central_Value; float PQ_Disc = (std::fabs(0.5f - PQ_Central_Value) < std::fabs(0.5f - PQ_Neighbourhood_Value)) ? PQ_Neighbourhood_Value : PQ_Central_Value;
// Diagonal gradients // Diagonal gradients
float NW_Grad = eps + std::fabs( rgb[c][indx - w1 - 1] - rgb[c][indx + w1 + 1]) + std::fabs( rgb[c][indx - w1 - 1] - rgb[c][indx - w3 - 3] ) + std::fabs( rgb[1][indx] - rgb[1][indx - w2 - 2] ); float NW_Grad = eps + std::fabs(rgb[c][indx - w1 - 1] - rgb[c][indx + w1 + 1]) + std::fabs(rgb[c][indx - w1 - 1] - rgb[c][indx - w3 - 3]) + std::fabs(rgb[1][indx] - rgb[1][indx - w2 - 2]);
float NE_Grad = eps + std::fabs( rgb[c][indx - w1 + 1] - rgb[c][indx + w1 - 1]) + std::fabs( rgb[c][indx - w1 + 1] - rgb[c][indx - w3 + 3] ) + std::fabs( rgb[1][indx] - rgb[1][indx - w2 + 2] ); float NE_Grad = eps + std::fabs(rgb[c][indx - w1 + 1] - rgb[c][indx + w1 - 1]) + std::fabs(rgb[c][indx - w1 + 1] - rgb[c][indx - w3 + 3]) + std::fabs(rgb[1][indx] - rgb[1][indx - w2 + 2]);
float SW_Grad = eps + std::fabs( rgb[c][indx - w1 + 1] - rgb[c][indx + w1 - 1]) + std::fabs( rgb[c][indx + w1 - 1] - rgb[c][indx + w3 - 3] ) + std::fabs( rgb[1][indx] - rgb[1][indx + w2 - 2] ); float SW_Grad = eps + std::fabs(rgb[c][indx - w1 + 1] - rgb[c][indx + w1 - 1]) + std::fabs(rgb[c][indx + w1 - 1] - rgb[c][indx + w3 - 3]) + std::fabs(rgb[1][indx] - rgb[1][indx + w2 - 2]);
float SE_Grad = eps + std::fabs( rgb[c][indx - w1 - 1] - rgb[c][indx + w1 + 1]) + std::fabs( rgb[c][indx + w1 + 1] - rgb[c][indx + w3 + 3] ) + std::fabs( rgb[1][indx] - rgb[1][indx + w2 + 2] ); float SE_Grad = eps + std::fabs(rgb[c][indx - w1 - 1] - rgb[c][indx + w1 + 1]) + std::fabs(rgb[c][indx + w1 + 1] - rgb[c][indx + w3 + 3]) + std::fabs(rgb[1][indx] - rgb[1][indx + w2 + 2]);
// Diagonal colour differences // Diagonal colour differences
float NW_Est = rgb[c][indx - w1 - 1] - rgb[1][indx - w1 - 1]; float NW_Est = rgb[c][indx - w1 - 1] - rgb[1][indx - w1 - 1];
@@ -194,36 +195,35 @@ void RawImageSource::rcd_demosaic()
float SE_Est = rgb[c][indx + w1 + 1] - rgb[1][indx + w1 + 1]; float SE_Est = rgb[c][indx + w1 + 1] - rgb[1][indx + w1 + 1];
// P/Q estimations // P/Q estimations
float P_Est = ( NW_Grad * SE_Est + SE_Grad * NW_Est ) / (NW_Grad + SE_Grad ); float P_Est = (NW_Grad * SE_Est + SE_Grad * NW_Est) / (NW_Grad + SE_Grad);
float Q_Est = ( NE_Grad * SW_Est + SW_Grad * NE_Est ) / (NE_Grad + SW_Grad ); float Q_Est = (NE_Grad * SW_Est + SW_Grad * NE_Est) / (NE_Grad + SW_Grad);
// R@B and B@R interpolation // R@B and B@R interpolation
rgb[c][indx] = rgb[1][indx] + ( 1.f - PQ_Disc ) * P_Est + PQ_Disc * Q_Est; rgb[c][indx] = rgb[1][indx] + (1.f - PQ_Disc) * P_Est + PQ_Disc * Q_Est;
} }
} }
// Step 4.3: Populate the red and blue channels at green CFA positions // Step 4.3: Populate the red and blue channels at green CFA positions
for ( int row = 4; row < tileSize - 4; row++ ) { for (int row = 4; row < tileSize - 4; row++) {
for ( int col = 4 + (FC( row, 1 ) & 1), indx = row * tileSize + col; col < tileSize - 4; col += 2, indx += 2 ) { for (int col = 4 + (FC(row, 1) & 1), indx = row * tileSize + col; col < tileSize - 4; col += 2, indx += 2) {
// Refined vertical and horizontal local discrimination // Refined vertical and horizontal local discrimination
float VH_Central_Value = VH_Dir[indx]; float VH_Central_Value = VH_Dir[indx];
float VH_Neighbourhood_Value = 0.25f * ( (VH_Dir[indx - w1 - 1] + VH_Dir[indx - w1 + 1]) + (VH_Dir[indx + w1 - 1] + VH_Dir[indx + w1 + 1]) ); float VH_Neighbourhood_Value = 0.25f * ((VH_Dir[indx - w1 - 1] + VH_Dir[indx - w1 + 1]) + (VH_Dir[indx + w1 - 1] + VH_Dir[indx + w1 + 1]));
float VH_Disc = ( std::fabs( 0.5f - VH_Central_Value ) < std::fabs( 0.5f - VH_Neighbourhood_Value ) ) ? VH_Neighbourhood_Value : VH_Central_Value;
float N1 = eps + std::fabs( rgb[1][indx] - rgb[1][indx - w2] );
float S1 = eps + std::fabs( rgb[1][indx] - rgb[1][indx + w2] );
float W1 = eps + std::fabs( rgb[1][indx] - rgb[1][indx - 2] );
float E1 = eps + std::fabs( rgb[1][indx] - rgb[1][indx + 2] );
for ( int c = 0; c <= 2; c += 2 ) {
float VH_Disc = (std::fabs(0.5f - VH_Central_Value) < std::fabs(0.5f - VH_Neighbourhood_Value)) ? VH_Neighbourhood_Value : VH_Central_Value;
float N1 = eps + std::fabs(rgb[1][indx] - rgb[1][indx - w2]);
float S1 = eps + std::fabs(rgb[1][indx] - rgb[1][indx + w2]);
float W1 = eps + std::fabs(rgb[1][indx] - rgb[1][indx - 2]);
float E1 = eps + std::fabs(rgb[1][indx] - rgb[1][indx + 2]);
for (int c = 0; c <= 2; c += 2) {
// Cardinal gradients // Cardinal gradients
float N_Grad = N1 + std::fabs( rgb[c][indx - w1] - rgb[c][indx + w1] ) + std::fabs( rgb[c][indx - w1] - rgb[c][indx - w3] ); float N_Grad = N1 + std::fabs(rgb[c][indx - w1] - rgb[c][indx + w1]) + std::fabs(rgb[c][indx - w1] - rgb[c][indx - w3]);
float S_Grad = S1 + std::fabs( rgb[c][indx - w1] - rgb[c][indx + w1] ) + std::fabs( rgb[c][indx + w1] - rgb[c][indx + w3] ); float S_Grad = S1 + std::fabs(rgb[c][indx - w1] - rgb[c][indx + w1]) + std::fabs(rgb[c][indx + w1] - rgb[c][indx + w3]);
float W_Grad = W1 + std::fabs( rgb[c][indx - 1] - rgb[c][indx + 1] ) + std::fabs( rgb[c][indx - 1] - rgb[c][indx - 3] ); float W_Grad = W1 + std::fabs(rgb[c][indx - 1] - rgb[c][indx + 1]) + std::fabs(rgb[c][indx - 1] - rgb[c][indx - 3]);
float E_Grad = E1 + std::fabs( rgb[c][indx - 1] - rgb[c][indx + 1] ) + std::fabs( rgb[c][indx + 1] - rgb[c][indx + 3] ); float E_Grad = E1 + std::fabs(rgb[c][indx - 1] - rgb[c][indx + 1]) + std::fabs(rgb[c][indx + 1] - rgb[c][indx + 3]);
// Cardinal colour differences // Cardinal colour differences
float N_Est = rgb[c][indx - w1] - rgb[1][indx - w1]; float N_Est = rgb[c][indx - w1] - rgb[1][indx - w1];
@@ -232,11 +232,11 @@ void RawImageSource::rcd_demosaic()
float E_Est = rgb[c][indx + 1] - rgb[1][indx + 1]; float E_Est = rgb[c][indx + 1] - rgb[1][indx + 1];
// Vertical and horizontal estimations // Vertical and horizontal estimations
float V_Est = ( N_Grad * S_Est + S_Grad * N_Est ) / (N_Grad + S_Grad ); float V_Est = (N_Grad * S_Est + S_Grad * N_Est) / (N_Grad + S_Grad);
float H_Est = ( E_Grad * W_Est + W_Grad * E_Est ) / (E_Grad + W_Grad ); float H_Est = (E_Grad * W_Est + W_Grad * E_Est) / (E_Grad + W_Grad);
// R@G and B@G interpolation // R@G and B@G interpolation
rgb[c][indx] = rgb[1][indx] + ( 1.f - VH_Disc ) * V_Est + VH_Disc * H_Est; rgb[c][indx] = rgb[1][indx] + (1.f - VH_Disc) * V_Est + VH_Disc * H_Est;
} }
} }