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

rcd demosaic: small speedup and reduction of memory usage

Browse Source
This commit is contained in:
Ingo Weyrich
2021-01-04 13:39:29 +01:00
parent cc7f1ecdab
commit 064a79971b
1 changed files with 44 additions and 41 deletions
Download Patch File Download Diff File Expand all files Collapse all files

85
rtengine/rcd_demosaic.cc
View File

@@ -84,17 +84,18 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
//Tolerance to avoid dividing by zero
constexpr float eps = 1e-5f;
constexpr float epssq = 1e-10f;
constexpr float scale = 65536.f;
#ifdef _OPENMP
#pragma omp parallel
#endif
{
int progresscounter = 0;
float *cfa = (float*) calloc(tileSize * tileSize, sizeof *cfa);
float (*rgb)[tileSize * tileSize] = (float (*)[tileSize * tileSize])malloc(3 * sizeof *rgb);
float *VH_Dir = (float*) calloc(tileSize * tileSize, sizeof *VH_Dir);
float *PQ_Dir = (float*) calloc(tileSize * tileSize, sizeof *PQ_Dir);
float *lpf = PQ_Dir; // reuse buffer, they don't overlap in usage
float *const cfa = (float*) calloc(tileSize * tileSize, sizeof *cfa);
float (*const rgb)[tileSize * tileSize] = (float (*)[tileSize * tileSize])malloc(3 * sizeof *rgb);
float *const VH_Dir = (float*) calloc(tileSize * tileSize, sizeof *VH_Dir);
float *const PQ_Dir = (float*) calloc(tileSize * tileSize / 2, sizeof *PQ_Dir);
float *const lpf = PQ_Dir; // reuse buffer, they don't overlap in usage
#ifdef _OPENMP
#pragma omp for schedule(dynamic, chunkSize) collapse(2) nowait
@@ -117,16 +118,15 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
for (int row = rowStart; row < rowEnd; row++) {
int indx = (row - rowStart) * tileSize;
int c0 = fc(cfarray, row, colStart);
int c1 = fc(cfarray, row, colStart + 1);
const int c0 = fc(cfarray, row, colStart);
const int c1 = fc(cfarray, row, colStart + 1);
int col = colStart;
for (; col < colEnd - 1; col+=2, indx+=2) {
cfa[indx] = rgb[c0][indx] = LIM01(rawData[row][col] / 65535.f);
cfa[indx + 1] = rgb[c1][indx + 1] = LIM01(rawData[row][col + 1] / 65535.f);
cfa[indx] = rgb[c0][indx] = LIM01(rawData[row][col] / scale);
cfa[indx + 1] = rgb[c1][indx + 1] = LIM01(rawData[row][col + 1] / scale);
}
if (col < colEnd) {
cfa[indx] = rgb[c0][indx] = LIM01(rawData[row][col] / 65535.f);
cfa[indx] = rgb[c0][indx] = LIM01(rawData[row][col] / scale);
}
}
@@ -138,8 +138,8 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
for (int col = 4, indx = row * tileSize + col; col < tilecols - 4; col++, indx++) {
const float cfai = cfa[indx];
//Calculate h/v local discrimination
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] - 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.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] - 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]);
float V_Stat = std::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] - 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 = std::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] - 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);
}
@@ -151,8 +151,10 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
// Step 2.1: Low pass filter incorporating green, red and blue local samples from the raw data
for (int row = 2; row < tileRows - 2; row++) {
for (int col = 2 + (fc(cfarray, row, 0) & 1), indx = row * tileSize + col; col < tilecols - 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]);
for (int col = 2 + (fc(cfarray, row, 0) & 1), indx = row * tileSize + col, lpindx = indx / 2; col < tilecols - 2; col += 2, indx += 2, ++lpindx) {
lpf[lpindx] = 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]);
}
}
@@ -163,11 +165,12 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
for (int row = 4; row < tileRows - 4; row++) {
int col = 4 + (fc(cfarray, row, 0) & 1);
int indx = row * tileSize + col;
int lpindx = indx / 2;
#ifdef __SSE2__
const vfloat zd5v = F2V(0.5f);
const vfloat zd25v = F2V(0.25f);
const vfloat epsv = F2V(eps);
for (; col < tilecols - 7; col += 8, indx += 8) {
for (; col < tilecols - 7; col += 8, indx += 8, lpindx += 4) {
// Cardinal gradients
const vfloat cfai = LC2VFU(cfa[indx]);
const vfloat N_Grad = epsv + (vabsf(LC2VFU(cfa[indx - w1]) - LC2VFU(cfa[indx + w1])) + vabsf(cfai - LC2VFU(cfa[indx - w2]))) + (vabsf(LC2VFU(cfa[indx - w1]) - LC2VFU(cfa[indx - w3])) + vabsf(LC2VFU(cfa[indx - w2]) - LC2VFU(cfa[indx - w4])));
@@ -176,11 +179,11 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
const vfloat E_Grad = epsv + (vabsf(LC2VFU(cfa[indx - 1]) - LC2VFU(cfa[indx + 1])) + vabsf(cfai - LC2VFU(cfa[indx + 2]))) + (vabsf(LC2VFU(cfa[indx + 1]) - LC2VFU(cfa[indx + 3])) + vabsf(LC2VFU(cfa[indx + 2]) - LC2VFU(cfa[indx + 4])));
// Cardinal pixel estimations
const vfloat lpfi = LVFU(lpf[indx>>1]);
const vfloat N_Est = LC2VFU(cfa[indx - w1]) + (LC2VFU(cfa[indx - w1]) * (lpfi - LVFU(lpf[(indx - w2)>>1])) / (epsv + lpfi + LVFU(lpf[(indx - w2)>>1])));
const vfloat S_Est = LC2VFU(cfa[indx + w1]) + (LC2VFU(cfa[indx + w1]) * (lpfi - LVFU(lpf[(indx + w2)>>1])) / (epsv + lpfi + LVFU(lpf[(indx + w2)>>1])));
const vfloat W_Est = LC2VFU(cfa[indx - 1]) + (LC2VFU(cfa[indx - 1]) * (lpfi - LVFU(lpf[(indx - 2)>>1])) / (epsv + lpfi + LVFU(lpf[(indx - 2)>>1])));
const vfloat E_Est = LC2VFU(cfa[indx + 1]) + (LC2VFU(cfa[indx + 1]) * (lpfi - LVFU(lpf[(indx + 2)>>1])) / (epsv + lpfi + LVFU(lpf[(indx + 2)>>1])));
const vfloat lpfi = LVFU(lpf[lpindx]);
const vfloat N_Est = LC2VFU(cfa[indx - w1]) + (LC2VFU(cfa[indx - w1]) * (lpfi - LVFU(lpf[lpindx - w1])) / (epsv + lpfi + LVFU(lpf[lpindx - w1])));
const vfloat S_Est = LC2VFU(cfa[indx + w1]) + (LC2VFU(cfa[indx + w1]) * (lpfi - LVFU(lpf[lpindx + w1])) / (epsv + lpfi + LVFU(lpf[lpindx + w1])));
const vfloat W_Est = LC2VFU(cfa[indx - 1]) + (LC2VFU(cfa[indx - 1]) * (lpfi - LVFU(lpf[lpindx - 1])) / (epsv + lpfi + LVFU(lpf[lpindx - 1])));
const vfloat E_Est = LC2VFU(cfa[indx + 1]) + (LC2VFU(cfa[indx + 1]) * (lpfi - LVFU(lpf[lpindx + 1])) / (epsv + lpfi + LVFU(lpf[lpindx + 1])));
// Vertical and horizontal estimations
const vfloat V_Est = (S_Grad * N_Est + N_Grad * S_Est) / (N_Grad + S_Grad);
@@ -200,7 +203,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
STC2VFU(rgb[1][indx], result);
}
#endif
for (; col < tilecols - 4; col += 2, indx += 2) {
for (; col < tilecols - 4; col += 2, indx += 2, ++lpindx) {
// Cardinal gradients
const float cfai = cfa[indx];
const float N_Grad = eps + (std::fabs(cfa[indx - w1] - cfa[indx + w1]) + std::fabs(cfai - cfa[indx - w2])) + (std::fabs(cfa[indx - w1] - cfa[indx - w3]) + std::fabs(cfa[indx - w2] - cfa[indx - w4]));
@@ -209,11 +212,11 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
const float E_Grad = eps + (std::fabs(cfa[indx - 1] - cfa[indx + 1]) + std::fabs(cfai - cfa[indx + 2])) + (std::fabs(cfa[indx + 1] - cfa[indx + 3]) + std::fabs(cfa[indx + 2] - cfa[indx + 4]));
// Cardinal pixel estimations
const float lpfi = lpf[indx>>1];
const float N_Est = cfa[indx - w1] * (1.f + (lpfi - lpf[(indx - w2)>>1]) / (eps + lpfi + lpf[(indx - w2)>>1]));
const float S_Est = cfa[indx + w1] * (1.f + (lpfi - lpf[(indx + w2)>>1]) / (eps + lpfi + lpf[(indx + w2)>>1]));
const float W_Est = cfa[indx - 1] * (1.f + (lpfi - lpf[(indx - 2)>>1]) / (eps + lpfi + lpf[(indx - 2)>>1]));
const float E_Est = cfa[indx + 1] * (1.f + (lpfi - lpf[(indx + 2)>>1]) / (eps + lpfi + lpf[(indx + 2)>>1]));
const float lpfi = lpf[lpindx];
const float N_Est = cfa[indx - w1] * (1.f + (lpfi - lpf[lpindx - w1]) / (eps + lpfi + lpf[lpindx - w1]));
const float S_Est = cfa[indx + w1] * (1.f + (lpfi - lpf[lpindx + w1]) / (eps + lpfi + lpf[lpindx + w1]));
const float W_Est = cfa[indx - 1] * (1.f + (lpfi - lpf[lpindx - 1]) / (eps + lpfi + lpf[lpindx - 1]));
const float E_Est = cfa[indx + 1] * (1.f + (lpfi - lpf[lpindx + 1]) / (eps + lpfi + lpf[lpindx + 1]));
// Vertical and horizontal estimations
const float V_Est = (S_Grad * N_Est + N_Grad * S_Est) / (N_Grad + S_Grad);
@@ -225,7 +228,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
const 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]));
const float VH_Disc = std::fabs(0.5f - VH_Central_Value) < std::fabs(0.5f - VH_Neighbourhood_Value) ? VH_Neighbourhood_Value : VH_Central_Value;
rgb[1][indx] = VH_Disc * H_Est + (1.f - VH_Disc) * V_Est;
rgb[1][indx] = intp(VH_Disc, H_Est, V_Est);
}
}
@@ -235,23 +238,23 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
// Step 4.1: Calculate P/Q diagonal local discrimination
for (int row = rcdBorder - 4; row < tileRows - rcdBorder + 4; row++) {
for (int col = rcdBorder - 4 + (fc(cfarray, row, rcdBorder) & 1), indx = row * tileSize + col; col < tilecols - rcdBorder + 4; col += 2, indx += 2) {
for (int col = rcdBorder - 4 + (fc(cfarray, row, rcdBorder) & 1), indx = row * tileSize + col, pqindx = indx / 2; col < tilecols - rcdBorder + 4; col += 2, indx += 2, ++pqindx) {
const float cfai = cfa[indx];
float P_Stat = max(epssq, - 18.f * cfai * (cfa[indx - w1 - 1] + cfa[indx + w1 + 1] + 2.f * (cfa[indx - w2 - 2] + cfa[indx + w2 + 2]) - cfa[indx - w3 - 3] - cfa[indx + w3 + 3]) - 2.f * cfai * (cfa[indx - w4 - 4] + cfa[indx + w4 + 4] - 19.f * cfai) - cfa[indx - w1 - 1] * (70.f * cfa[indx + w1 + 1] - 12.f * cfa[indx - w2 - 2] + 24.f * cfa[indx + w2 + 2] - 38.f * cfa[indx - w3 - 3] + 16.f * cfa[indx + w3 + 3] + 12.f * cfa[indx - w4 - 4] - 6.f * cfa[indx + w4 + 4] + 46.f * cfa[indx - w1 - 1]) + cfa[indx + w1 + 1] * (24.f * cfa[indx - w2 - 2] - 12.f * cfa[indx + w2 + 2] + 16.f * cfa[indx - w3 - 3] - 38.f * cfa[indx + w3 + 3] - 6.f * cfa[indx - w4 - 4] + 12.f * cfa[indx + w4 + 4] + 46.f * cfa[indx + w1 + 1]) + cfa[indx - w2 - 2] * (14.f * cfa[indx + w2 + 2] - 12.f * cfa[indx + w3 + 3] - 2.f * (cfa[indx - w4 - 4] - cfa[indx + w4 + 4]) + 11.f * cfa[indx - w2 - 2]) - cfa[indx + w2 + 2] * (12.f * cfa[indx - w3 - 3] + 2.f * (cfa[indx - w4 - 4] - cfa[indx + w4 + 4]) + 11.f * cfa[indx + w2 + 2]) + cfa[indx - w3 - 3] * (2.f * cfa[indx + w3 + 3] - 6.f * cfa[indx - w4 - 4] + 10.f * cfa[indx - w3 - 3]) - cfa[indx + w3 + 3] * (6.f * cfa[indx + w4 + 4] + 10.f * cfa[indx + w3 + 3]) + cfa[indx - w4 - 4] * cfa[indx - w4 - 4] + cfa[indx + w4 + 4] * cfa[indx + w4 + 4]);
float Q_Stat = max(epssq, - 18.f * cfai * (cfa[indx + w1 - 1] + cfa[indx - w1 + 1] + 2.f * (cfa[indx + w2 - 2] + cfa[indx - w2 + 2]) - cfa[indx + w3 - 3] - cfa[indx - w3 + 3]) - 2.f * cfai * (cfa[indx + w4 - 4] + cfa[indx - w4 + 4] - 19.f * cfai) - cfa[indx + w1 - 1] * (70.f * cfa[indx - w1 + 1] - 12.f * cfa[indx + w2 - 2] + 24.f * cfa[indx - w2 + 2] - 38.f * cfa[indx + w3 - 3] + 16.f * cfa[indx - w3 + 3] + 12.f * cfa[indx + w4 - 4] - 6.f * cfa[indx - w4 + 4] + 46.f * cfa[indx + w1 - 1]) + cfa[indx - w1 + 1] * (24.f * cfa[indx + w2 - 2] - 12.f * cfa[indx - w2 + 2] + 16.f * cfa[indx + w3 - 3] - 38.f * cfa[indx - w3 + 3] - 6.f * cfa[indx + w4 - 4] + 12.f * cfa[indx - w4 + 4] + 46.f * cfa[indx - w1 + 1]) + cfa[indx + w2 - 2] * (14.f * cfa[indx - w2 + 2] - 12.f * cfa[indx - w3 + 3] - 2.f * (cfa[indx + w4 - 4] - cfa[indx - w4 + 4]) + 11.f * cfa[indx + w2 - 2]) - cfa[indx - w2 + 2] * (12.f * cfa[indx + w3 - 3] + 2.f * (cfa[indx + w4 - 4] - cfa[indx - w4 + 4]) + 11.f * cfa[indx - w2 + 2]) + cfa[indx + w3 - 3] * (2.f * cfa[indx - w3 + 3] - 6.f * cfa[indx + w4 - 4] + 10.f * cfa[indx + w3 - 3]) - cfa[indx - w3 + 3] * (6.f * cfa[indx - w4 + 4] + 10.f * cfa[indx - w3 + 3]) + cfa[indx + w4 - 4] * cfa[indx + w4 - 4] + cfa[indx - w4 + 4] * cfa[indx - w4 + 4]);
float P_Stat = std::max(epssq, - 18.f * cfai * (cfa[indx - w1 - 1] + cfa[indx + w1 + 1] + 2.f * (cfa[indx - w2 - 2] + cfa[indx + w2 + 2]) - cfa[indx - w3 - 3] - cfa[indx + w3 + 3]) - 2.f * cfai * (cfa[indx - w4 - 4] + cfa[indx + w4 + 4] - 19.f * cfai) - cfa[indx - w1 - 1] * (70.f * cfa[indx + w1 + 1] - 12.f * cfa[indx - w2 - 2] + 24.f * cfa[indx + w2 + 2] - 38.f * cfa[indx - w3 - 3] + 16.f * cfa[indx + w3 + 3] + 12.f * cfa[indx - w4 - 4] - 6.f * cfa[indx + w4 + 4] + 46.f * cfa[indx - w1 - 1]) + cfa[indx + w1 + 1] * (24.f * cfa[indx - w2 - 2] - 12.f * cfa[indx + w2 + 2] + 16.f * cfa[indx - w3 - 3] - 38.f * cfa[indx + w3 + 3] - 6.f * cfa[indx - w4 - 4] + 12.f * cfa[indx + w4 + 4] + 46.f * cfa[indx + w1 + 1]) + cfa[indx - w2 - 2] * (14.f * cfa[indx + w2 + 2] - 12.f * cfa[indx + w3 + 3] - 2.f * (cfa[indx - w4 - 4] - cfa[indx + w4 + 4]) + 11.f * cfa[indx - w2 - 2]) - cfa[indx + w2 + 2] * (12.f * cfa[indx - w3 - 3] + 2.f * (cfa[indx - w4 - 4] - cfa[indx + w4 + 4]) + 11.f * cfa[indx + w2 + 2]) + cfa[indx - w3 - 3] * (2.f * cfa[indx + w3 + 3] - 6.f * cfa[indx - w4 - 4] + 10.f * cfa[indx - w3 - 3]) - cfa[indx + w3 + 3] * (6.f * cfa[indx + w4 + 4] + 10.f * cfa[indx + w3 + 3]) + cfa[indx - w4 - 4] * cfa[indx - w4 - 4] + cfa[indx + w4 + 4] * cfa[indx + w4 + 4]);
float Q_Stat = std::max(epssq, - 18.f * cfai * (cfa[indx + w1 - 1] + cfa[indx - w1 + 1] + 2.f * (cfa[indx + w2 - 2] + cfa[indx - w2 + 2]) - cfa[indx + w3 - 3] - cfa[indx - w3 + 3]) - 2.f * cfai * (cfa[indx + w4 - 4] + cfa[indx - w4 + 4] - 19.f * cfai) - cfa[indx + w1 - 1] * (70.f * cfa[indx - w1 + 1] - 12.f * cfa[indx + w2 - 2] + 24.f * cfa[indx - w2 + 2] - 38.f * cfa[indx + w3 - 3] + 16.f * cfa[indx - w3 + 3] + 12.f * cfa[indx + w4 - 4] - 6.f * cfa[indx - w4 + 4] + 46.f * cfa[indx + w1 - 1]) + cfa[indx - w1 + 1] * (24.f * cfa[indx + w2 - 2] - 12.f * cfa[indx - w2 + 2] + 16.f * cfa[indx + w3 - 3] - 38.f * cfa[indx - w3 + 3] - 6.f * cfa[indx + w4 - 4] + 12.f * cfa[indx - w4 + 4] + 46.f * cfa[indx - w1 + 1]) + cfa[indx + w2 - 2] * (14.f * cfa[indx - w2 + 2] - 12.f * cfa[indx - w3 + 3] - 2.f * (cfa[indx + w4 - 4] - cfa[indx - w4 + 4]) + 11.f * cfa[indx + w2 - 2]) - cfa[indx - w2 + 2] * (12.f * cfa[indx + w3 - 3] + 2.f * (cfa[indx + w4 - 4] - cfa[indx - w4 + 4]) + 11.f * cfa[indx - w2 + 2]) + cfa[indx + w3 - 3] * (2.f * cfa[indx - w3 + 3] - 6.f * cfa[indx + w4 - 4] + 10.f * cfa[indx + w3 - 3]) - cfa[indx - w3 + 3] * (6.f * cfa[indx - w4 + 4] + 10.f * cfa[indx - w3 + 3]) + cfa[indx + w4 - 4] * cfa[indx + w4 - 4] + cfa[indx - w4 + 4] * cfa[indx - w4 + 4]);
PQ_Dir[indx] = P_Stat / (P_Stat + Q_Stat);
PQ_Dir[pqindx] = P_Stat / (P_Stat + Q_Stat);
}
}
// Step 4.2: Populate the red and blue channels at blue and red CFA positions
for (int row = rcdBorder - 3; row < tileRows - rcdBorder + 3; row++) {
for (int col = rcdBorder - 3 + (fc(cfarray, row, rcdBorder - 1) & 1), indx = row * tileSize + col, c = 2 - fc(cfarray, row, col); col < tilecols - rcdBorder + 3; col += 2, indx += 2) {
for (int col = rcdBorder - 3 + (fc(cfarray, row, rcdBorder - 1) & 1), indx = row * tileSize + col, c = 2 - fc(cfarray, row, col), pqindx = indx / 2, pqindx2 = (indx - w1 - 1) / 2, pqindx3 = (indx + w1 - 1) / 2; col < tilecols - rcdBorder + 3; col += 2, indx += 2, ++pqindx, ++pqindx2, ++pqindx3) {
// Refined P/Q diagonal local discrimination
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_Central_Value = PQ_Dir[pqindx];
float PQ_Neighbourhood_Value = 0.25f * (PQ_Dir[pqindx2] + PQ_Dir[pqindx2 + 1] + PQ_Dir[pqindx3] + PQ_Dir[pqindx3 + 1]);
float PQ_Disc = (std::fabs(0.5f - PQ_Central_Value) < std::fabs(0.5f - PQ_Neighbourhood_Value)) ? PQ_Neighbourhood_Value : PQ_Central_Value;
@@ -272,7 +275,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
float Q_Est = (NE_Grad * SW_Est + SW_Grad * NE_Est) / (NE_Grad + SW_Grad);
// 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] + intp(PQ_Disc, Q_Est, P_Est);
}
}
@@ -281,7 +284,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
for (int col = rcdBorder + (fc(cfarray, row, rcdBorder - 1) & 1), indx = row * tileSize + col; col < tilecols - rcdBorder; col += 2, indx += 2) {
// 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_Disc = (std::fabs(0.5f - VH_Central_Value) < std::fabs(0.5f - VH_Neighbourhood_Value)) ? VH_Neighbourhood_Value : VH_Central_Value;
@@ -315,7 +318,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
float H_Est = (E_Grad * W_Est + W_Grad * E_Est) / (E_Grad + W_Grad);
// R@G and B@G interpolation
rgb[c][indx] = rgb1 + (1.f - VH_Disc) * V_Est + VH_Disc * H_Est;
rgb[c][indx] = rgb1 + intp(VH_Disc, H_Est, V_Est);
}
}
}
@@ -323,9 +326,9 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
for (int row = rowStart + rcdBorder; row < rowEnd - rcdBorder; ++row) {
for (int col = colStart + rcdBorder; col < colEnd - rcdBorder; ++col) {
int idx = (row - rowStart) * tileSize + col - colStart ;
red[row][col] = std::max(0.f, rgb[0][idx] * 65535.f);
green[row][col] = std::max(0.f, rgb[1][idx] * 65535.f);
blue[row][col] = std::max(0.f, rgb[2][idx] * 65535.f);
red[row][col] = std::max(0.f, rgb[0][idx] * scale);
green[row][col] = std::max(0.f, rgb[1][idx] * scale);
blue[row][col] = std::max(0.f, rgb[2][idx] * scale);
}
}