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

Some cleanups, also fix some fallouts from reduce_include_dependencies branch

Browse Source
This commit is contained in:
Ingo Weyrich 2019-11-05 20:31:17 +01:00
parent d3e7f6a591
commit 03c0d6c86c
16 changed files with 175 additions and 410 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

102
rtengine/CA_correct_RT.cc
View File

@ -29,6 +29,13 @@
#include "median.h" #include "median.h"
#include "StopWatch.h" #include "StopWatch.h"
namespace
{
unsigned fc(const unsigned int cfa[2][2], int r, int c) {
return cfa[r & 1][c & 1];
}
}
namespace { namespace {
bool LinEqSolve(int nDim, double* pfMatr, double* pfVect, double* pfSolution) bool LinEqSolve(int nDim, double* pfMatr, double* pfVect, double* pfSolution)
@ -135,6 +142,7 @@ float* RawImageSource::CA_correct_RT(
} }
// multithreaded and vectorized by Ingo Weyrich // multithreaded and vectorized by Ingo Weyrich
const unsigned int cfa[2][2] = {{FC(0,0), FC(0,1)}, {FC(1,0), FC(1,1)}};
constexpr int ts = 128; constexpr int ts = 128;
constexpr int tsh = ts / 2; constexpr int tsh = ts / 2;
constexpr int cb = 2; // 2 pixels border will be excluded from correction constexpr int cb = 2; // 2 pixels border will be excluded from correction
@ -145,7 +153,7 @@ float* RawImageSource::CA_correct_RT(
// Test for RGB cfa // Test for RGB cfa
for (int i = 0; i < 2; i++) { for (int i = 0; i < 2; i++) {
for (int j = 0; j < 2; j++) { for (int j = 0; j < 2; j++) {
if (FC(i, j) == 3) { if (fc(cfa, i, j) == 3) {
std::cout << "CA correction supports only RGB Colour filter arrays" << std::endl; std::cout << "CA correction supports only RGB Colour filter arrays" << std::endl;
return buffer; return buffer;
} }
@ -164,7 +172,7 @@ float* RawImageSource::CA_correct_RT(
#pragma omp parallel for #pragma omp parallel for
#endif #endif
for (int i = cb; i < H - cb; ++i) { for (int i = cb; i < H - cb; ++i) {
for (int j = cb + (FC(i, 0) & 1); j < W - cb; j += 2) { for (int j = cb + (fc(cfa, i, 0) & 1); j < W - cb; j += 2) {
(*oldraw)[i - cb][(j - cb) / 2] = rawData[i][j]; (*oldraw)[i - cb][(j - cb) / 2] = rawData[i][j];
} }
} }
@ -317,12 +325,12 @@ float* RawImageSource::CA_correct_RT(
int cc = ccmin; int cc = ccmin;
int col = cc + left; int col = cc + left;
#ifdef __SSE2__ #ifdef __SSE2__
int c0 = FC(rr, cc); int c0 = fc(cfa, rr, cc);
if (c0 == 1) { if (c0 == 1) {
rgb[c0][rr * ts + cc] = rawData[row][col] / 65535.f; rgb[c0][rr * ts + cc] = rawData[row][col] / 65535.f;
cc++; cc++;
col++; col++;
c0 = FC(rr, cc); c0 = fc(cfa, rr, cc);
} }
int indx1 = rr * ts + cc; int indx1 = rr * ts + cc;
for (; cc < ccmax - 7; cc+=8, col+=8, indx1 += 8) { for (; cc < ccmax - 7; cc+=8, col+=8, indx1 += 8) {
@ -335,7 +343,7 @@ float* RawImageSource::CA_correct_RT(
} }
#endif #endif
for (; cc < ccmax; cc++, col++) { for (; cc < ccmax; cc++, col++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
int indx1 = rr * ts + cc; int indx1 = rr * ts + cc;
rgb[c][indx1 >> ((c & 1) ^ 1)] = rawData[row][col] / 65535.f; rgb[c][indx1 >> ((c & 1) ^ 1)] = rawData[row][col] / 65535.f;
} }
@ -345,7 +353,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmin > 0) { if (rrmin > 0) {
for (int rr = 0; rr < border; rr++) { for (int rr = 0; rr < border; rr++) {
for (int cc = ccmin; cc < ccmax; cc++) { for (int cc = ccmin; cc < ccmax; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rgb[c][((border2 - rr) * ts + cc) >> ((c & 1) ^ 1)]; rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rgb[c][((border2 - rr) * ts + cc) >> ((c & 1) ^ 1)];
} }
} }
@ -354,7 +362,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmax < rr1) { if (rrmax < rr1) {
for (int rr = 0; rr < border; rr++) { for (int rr = 0; rr < border; rr++) {
for (int cc = ccmin; cc < ccmax; cc++) { for (int cc = ccmin; cc < ccmax; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][((rrmax + rr)*ts + cc) >> ((c & 1) ^ 1)] = rawData[(height - rr - 2)][left + cc] / 65535.f; rgb[c][((rrmax + rr)*ts + cc) >> ((c & 1) ^ 1)] = rawData[(height - rr - 2)][left + cc] / 65535.f;
} }
} }
@ -363,7 +371,7 @@ float* RawImageSource::CA_correct_RT(
if (ccmin > 0) { if (ccmin > 0) {
for (int rr = rrmin; rr < rrmax; rr++) { for (int rr = rrmin; rr < rrmax; rr++) {
for (int cc = 0; cc < border; cc++) { for (int cc = 0; cc < border; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rgb[c][(rr * ts + border2 - cc) >> ((c & 1) ^ 1)]; rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rgb[c][(rr * ts + border2 - cc) >> ((c & 1) ^ 1)];
} }
} }
@ -372,7 +380,7 @@ float* RawImageSource::CA_correct_RT(
if (ccmax < cc1) { if (ccmax < cc1) {
for (int rr = rrmin; rr < rrmax; rr++) { for (int rr = rrmin; rr < rrmax; rr++) {
for (int cc = 0; cc < border; cc++) { for (int cc = 0; cc < border; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + ccmax + cc) >> ((c & 1) ^ 1)] = rawData[(top + rr)][(width - cc - 2)] / 65535.f; rgb[c][(rr * ts + ccmax + cc) >> ((c & 1) ^ 1)] = rawData[(top + rr)][(width - cc - 2)] / 65535.f;
} }
} }
@ -382,7 +390,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmin > 0 && ccmin > 0) { if (rrmin > 0 && ccmin > 0) {
for (int rr = 0; rr < border; rr++) { for (int rr = 0; rr < border; rr++) {
for (int cc = 0; cc < border; cc++) { for (int cc = 0; cc < border; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rawData[border2 - rr][border2 - cc] / 65535.f; rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rawData[border2 - rr][border2 - cc] / 65535.f;
} }
} }
@ -391,7 +399,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmax < rr1 && ccmax < cc1) { if (rrmax < rr1 && ccmax < cc1) {
for (int rr = 0; rr < border; rr++) { for (int rr = 0; rr < border; rr++) {
for (int cc = 0; cc < border; cc++) { for (int cc = 0; cc < border; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][((rrmax + rr)*ts + ccmax + cc) >> ((c & 1) ^ 1)] = rawData[(height - rr - 2)][(width - cc - 2)] / 65535.f; rgb[c][((rrmax + rr)*ts + ccmax + cc) >> ((c & 1) ^ 1)] = rawData[(height - rr - 2)][(width - cc - 2)] / 65535.f;
} }
} }
@ -400,7 +408,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmin > 0 && ccmax < cc1) { if (rrmin > 0 && ccmax < cc1) {
for (int rr = 0; rr < border; rr++) { for (int rr = 0; rr < border; rr++) {
for (int cc = 0; cc < border; cc++) { for (int cc = 0; cc < border; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + ccmax + cc) >> ((c & 1) ^ 1)] = rawData[(border2 - rr)][(width - cc - 2)] / 65535.f; rgb[c][(rr * ts + ccmax + cc) >> ((c & 1) ^ 1)] = rawData[(border2 - rr)][(width - cc - 2)] / 65535.f;
} }
} }
@ -409,7 +417,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmax < rr1 && ccmin > 0) { if (rrmax < rr1 && ccmin > 0) {
for (int rr = 0; rr < border; rr++) { for (int rr = 0; rr < border; rr++) {
for (int cc = 0; cc < border; cc++) { for (int cc = 0; cc < border; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][((rrmax + rr)*ts + cc) >> ((c & 1) ^ 1)] = rawData[(height - rr - 2)][(border2 - cc)] / 65535.f; rgb[c][((rrmax + rr)*ts + cc) >> ((c & 1) ^ 1)] = rawData[(height - rr - 2)][(border2 - cc)] / 65535.f;
} }
} }
@ -424,9 +432,9 @@ float* RawImageSource::CA_correct_RT(
#endif #endif
for (int rr = 3; rr < rr1 - 3; rr++) { for (int rr = 3; rr < rr1 - 3; rr++) {
int row = rr + top; int row = rr + top;
int cc = 3 + (FC(rr,3) & 1); int cc = 3 + (fc(cfa, rr,3) & 1);
int indx = rr * ts + cc; int indx = rr * ts + cc;
int c = FC(rr,cc); int c = fc(cfa, rr,cc);
#ifdef __SSE2__ #ifdef __SSE2__
for (; cc < cc1 - 9; cc+=8, indx+=8) { for (; cc < cc1 - 9; cc+=8, indx+=8) {
//compute directional weights using image gradients //compute directional weights using image gradients
@ -460,7 +468,7 @@ float* RawImageSource::CA_correct_RT(
} }
if (row > -1 && row < height) { if (row > -1 && row < height) {
int offset = (FC(row,max(left + 3, 0)) & 1); int offset = (fc(cfa, row,max(left + 3, 0)) & 1);
int col = max(left + 3, 0) + offset; int col = max(left + 3, 0) + offset;
int indx = rr * ts + 3 - (left < 0 ? (left+3) : 0) + offset; int indx = rr * ts + 3 - (left < 0 ? (left+3) : 0) + offset;
#ifdef __SSE2__ #ifdef __SSE2__
@ -478,9 +486,9 @@ float* RawImageSource::CA_correct_RT(
vfloat zd25v = F2V(0.25f); vfloat zd25v = F2V(0.25f);
#endif #endif
for (int rr = 4; rr < rr1 - 4; rr++) { for (int rr = 4; rr < rr1 - 4; rr++) {
int cc = 4 + (FC(rr, 2) & 1); int cc = 4 + (fc(cfa, rr, 2) & 1);
int indx = rr * ts + cc; int indx = rr * ts + cc;
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
#ifdef __SSE2__ #ifdef __SSE2__
for (; cc < cc1 - 10; cc += 8, indx += 8) { for (; cc < cc1 - 10; cc += 8, indx += 8) {
vfloat rgb1v = LC2VFU(rgb[1][indx]); vfloat rgb1v = LC2VFU(rgb[1][indx]);
@ -544,9 +552,9 @@ float* RawImageSource::CA_correct_RT(
// along line segments, find the point along each segment that minimizes the colour variance // along line segments, find the point along each segment that minimizes the colour variance
// averaged over the tile; evaluate for up/down and left/right away from R/B grid point // averaged over the tile; evaluate for up/down and left/right away from R/B grid point
for (int rr = 8; rr < rr1 - 8; rr++) { for (int rr = 8; rr < rr1 - 8; rr++) {
int cc = 8 + (FC(rr, 2) & 1); int cc = 8 + (fc(cfa, rr, 2) & 1);
int indx = rr * ts + cc; int indx = rr * ts + cc;
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
#ifdef __SSE2__ #ifdef __SSE2__
vfloat coeff00v = ZEROV; vfloat coeff00v = ZEROV;
vfloat coeff01v = ZEROV; vfloat coeff01v = ZEROV;
@ -868,14 +876,14 @@ float* RawImageSource::CA_correct_RT(
int indx = row * width + col; int indx = row * width + col;
int indx1 = rr * ts + cc; int indx1 = rr * ts + cc;
#ifdef __SSE2__ #ifdef __SSE2__
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
if (c & 1) { if (c & 1) {
rgb[1][indx1] = rawData[row][col] / 65535.f; rgb[1][indx1] = rawData[row][col] / 65535.f;
indx++; indx++;
indx1++; indx1++;
cc++; cc++;
col++; col++;
c = FC(rr, cc); c = fc(cfa, rr, cc);
} }
for (; cc < ccmax - 7; cc += 8, col += 8, indx += 8, indx1 += 8) { for (; cc < ccmax - 7; cc += 8, col += 8, indx += 8, indx1 += 8) {
vfloat val1v = LVFU(rawData[row][col]) / c65535v; vfloat val1v = LVFU(rawData[row][col]) / c65535v;
@ -887,7 +895,7 @@ float* RawImageSource::CA_correct_RT(
} }
#endif #endif
for (; cc < ccmax; cc++, col++, indx++, indx1++) { for (; cc < ccmax; cc++, col++, indx++, indx1++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][indx1 >> ((c & 1) ^ 1)] = rawData[row][col] / 65535.f; rgb[c][indx1 >> ((c & 1) ^ 1)] = rawData[row][col] / 65535.f;
if ((c & 1) == 0) { if ((c & 1) == 0) {
@ -900,7 +908,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmin > 0) { if (rrmin > 0) {
for (int rr = 0; rr < border; rr++) { for (int rr = 0; rr < border; rr++) {
for (int cc = ccmin; cc < ccmax; cc++) { for (int cc = ccmin; cc < ccmax; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rgb[c][((border2 - rr) * ts + cc) >> ((c & 1) ^ 1)]; rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rgb[c][((border2 - rr) * ts + cc) >> ((c & 1) ^ 1)];
rgb[1][rr * ts + cc] = rgb[1][(border2 - rr) * ts + cc]; rgb[1][rr * ts + cc] = rgb[1][(border2 - rr) * ts + cc];
} }
@ -910,7 +918,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmax < rr1) { if (rrmax < rr1) {
for (int rr = 0; rr < std::min(border, rr1 - rrmax); rr++) { for (int rr = 0; rr < std::min(border, rr1 - rrmax); rr++) {
for (int cc = ccmin; cc < ccmax; cc++) { for (int cc = ccmin; cc < ccmax; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][((rrmax + rr)*ts + cc) >> ((c & 1) ^ 1)] = (rawData[(height - rr - 2)][left + cc]) / 65535.f; rgb[c][((rrmax + rr)*ts + cc) >> ((c & 1) ^ 1)] = (rawData[(height - rr - 2)][left + cc]) / 65535.f;
if ((c & 1) == 0) { if ((c & 1) == 0) {
rgb[1][(rrmax + rr)*ts + cc] = Gtmp[((height - rr - 2) * width + left + cc) >> 1]; rgb[1][(rrmax + rr)*ts + cc] = Gtmp[((height - rr - 2) * width + left + cc) >> 1];
@ -922,7 +930,7 @@ float* RawImageSource::CA_correct_RT(
if (ccmin > 0) { if (ccmin > 0) {
for (int rr = rrmin; rr < rrmax; rr++) { for (int rr = rrmin; rr < rrmax; rr++) {
for (int cc = 0; cc < border; cc++) { for (int cc = 0; cc < border; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rgb[c][(rr * ts + border2 - cc) >> ((c & 1) ^ 1)]; rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = rgb[c][(rr * ts + border2 - cc) >> ((c & 1) ^ 1)];
rgb[1][rr * ts + cc] = rgb[1][rr * ts + border2 - cc]; rgb[1][rr * ts + cc] = rgb[1][rr * ts + border2 - cc];
} }
@ -932,7 +940,7 @@ float* RawImageSource::CA_correct_RT(
if (ccmax < cc1) { if (ccmax < cc1) {
for (int rr = rrmin; rr < rrmax; rr++) { for (int rr = rrmin; rr < rrmax; rr++) {
for (int cc = 0; cc < std::min(border, cc1 - ccmax); cc++) { for (int cc = 0; cc < std::min(border, cc1 - ccmax); cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + ccmax + cc) >> ((c & 1) ^ 1)] = (rawData[(top + rr)][(width - cc - 2)]) / 65535.f; rgb[c][(rr * ts + ccmax + cc) >> ((c & 1) ^ 1)] = (rawData[(top + rr)][(width - cc - 2)]) / 65535.f;
if ((c & 1) == 0) { if ((c & 1) == 0) {
rgb[1][rr * ts + ccmax + cc] = Gtmp[((top + rr) * width + (width - cc - 2)) >> 1]; rgb[1][rr * ts + ccmax + cc] = Gtmp[((top + rr) * width + (width - cc - 2)) >> 1];
@ -945,7 +953,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmin > 0 && ccmin > 0) { if (rrmin > 0 && ccmin > 0) {
for (int rr = 0; rr < border; rr++) { for (int rr = 0; rr < border; rr++) {
for (int cc = 0; cc < border; cc++) { for (int cc = 0; cc < border; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = (rawData[border2 - rr][border2 - cc]) / 65535.f; rgb[c][(rr * ts + cc) >> ((c & 1) ^ 1)] = (rawData[border2 - rr][border2 - cc]) / 65535.f;
if ((c & 1) == 0) { if ((c & 1) == 0) {
rgb[1][rr * ts + cc] = Gtmp[((border2 - rr) * width + border2 - cc) >> 1]; rgb[1][rr * ts + cc] = Gtmp[((border2 - rr) * width + border2 - cc) >> 1];
@ -957,7 +965,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmax < rr1 && ccmax < cc1) { if (rrmax < rr1 && ccmax < cc1) {
for (int rr = 0; rr < std::min(border, rr1 - rrmax); rr++) { for (int rr = 0; rr < std::min(border, rr1 - rrmax); rr++) {
for (int cc = 0; cc < std::min(border, cc1 - ccmax); cc++) { for (int cc = 0; cc < std::min(border, cc1 - ccmax); cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][((rrmax + rr)*ts + ccmax + cc) >> ((c & 1) ^ 1)] = (rawData[(height - rr - 2)][(width - cc - 2)]) / 65535.f; rgb[c][((rrmax + rr)*ts + ccmax + cc) >> ((c & 1) ^ 1)] = (rawData[(height - rr - 2)][(width - cc - 2)]) / 65535.f;
if ((c & 1) == 0) { if ((c & 1) == 0) {
rgb[1][(rrmax + rr)*ts + ccmax + cc] = Gtmp[((height - rr - 2) * width + (width - cc - 2)) >> 1]; rgb[1][(rrmax + rr)*ts + ccmax + cc] = Gtmp[((height - rr - 2) * width + (width - cc - 2)) >> 1];
@ -969,7 +977,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmin > 0 && ccmax < cc1) { if (rrmin > 0 && ccmax < cc1) {
for (int rr = 0; rr < border; rr++) { for (int rr = 0; rr < border; rr++) {
for (int cc = 0; cc < std::min(border, cc1 - ccmax); cc++) { for (int cc = 0; cc < std::min(border, cc1 - ccmax); cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][(rr * ts + ccmax + cc) >> ((c & 1) ^ 1)] = (rawData[(border2 - rr)][(width - cc - 2)]) / 65535.f; rgb[c][(rr * ts + ccmax + cc) >> ((c & 1) ^ 1)] = (rawData[(border2 - rr)][(width - cc - 2)]) / 65535.f;
if ((c & 1) == 0) { if ((c & 1) == 0) {
rgb[1][rr * ts + ccmax + cc] = Gtmp[((border2 - rr) * width + (width - cc - 2)) >> 1]; rgb[1][rr * ts + ccmax + cc] = Gtmp[((border2 - rr) * width + (width - cc - 2)) >> 1];
@ -981,7 +989,7 @@ float* RawImageSource::CA_correct_RT(
if (rrmax < rr1 && ccmin > 0) { if (rrmax < rr1 && ccmin > 0) {
for (int rr = 0; rr < std::min(border, rr1 - rrmax); rr++) { for (int rr = 0; rr < std::min(border, rr1 - rrmax); rr++) {
for (int cc = 0; cc < border; cc++) { for (int cc = 0; cc < border; cc++) {
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
rgb[c][((rrmax + rr)*ts + cc) >> ((c & 1) ^ 1)] = (rawData[(height - rr - 2)][(border2 - cc)]) / 65535.f; rgb[c][((rrmax + rr)*ts + cc) >> ((c & 1) ^ 1)] = (rawData[(height - rr - 2)][(border2 - cc)]) / 65535.f;
if ((c & 1) == 0) { if ((c & 1) == 0) {
rgb[1][(rrmax + rr)*ts + cc] = Gtmp[((height - rr - 2) * width + (border2 - cc)) >> 1]; rgb[1][(rrmax + rr)*ts + cc] = Gtmp[((height - rr - 2) * width + (border2 - cc)) >> 1];
@ -998,7 +1006,7 @@ float* RawImageSource::CA_correct_RT(
#endif #endif
//manual CA correction; use red/blue slider values to set CA shift parameters //manual CA correction; use red/blue slider values to set CA shift parameters
for (int rr = 3; rr < rr1 - 3; rr++) { for (int rr = 3; rr < rr1 - 3; rr++) {
int cc = 3 + FC(rr, 1), c = FC(rr,cc), indx = rr * ts + cc; int cc = 3 + fc(cfa, rr, 1), c = fc(cfa, rr,cc), indx = rr * ts + cc;
#ifdef __SSE2__ #ifdef __SSE2__
for (; cc < cc1 - 10; cc += 8, indx += 8) { for (; cc < cc1 - 10; cc += 8, indx += 8) {
//compute directional weights using image gradients //compute directional weights using image gradients
@ -1079,8 +1087,8 @@ float* RawImageSource::CA_correct_RT(
} }
for (int rr = 4; rr < rr1 - 4; rr++) { for (int rr = 4; rr < rr1 - 4; rr++) {
int cc = 4 + (FC(rr, 2) & 1); int cc = 4 + (fc(cfa, rr, 2) & 1);
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
int indx = (rr * ts + cc) >> 1; int indx = (rr * ts + cc) >> 1;
int indxfc = (rr + shiftvfloor[c]) * ts + cc + shifthceil[c]; int indxfc = (rr + shiftvfloor[c]) * ts + cc + shifthceil[c];
int indxff = (rr + shiftvfloor[c]) * ts + cc + shifthfloor[c]; int indxff = (rr + shiftvfloor[c]) * ts + cc + shifthfloor[c];
@ -1129,8 +1137,8 @@ float* RawImageSource::CA_correct_RT(
vfloat epsv = F2V(eps); vfloat epsv = F2V(eps);
#endif #endif
for (int rr = 8; rr < rr1 - 8; rr++) { for (int rr = 8; rr < rr1 - 8; rr++) {
int cc = 8 + (FC(rr, 2) & 1); int cc = 8 + (fc(cfa, rr, 2) & 1);
int c = FC(rr, cc); int c = fc(cfa, rr, cc);
int GRBdir0 = GRBdir[0][c]; int GRBdir0 = GRBdir[0][c];
int GRBdir1 = GRBdir[1][c]; int GRBdir1 = GRBdir[1][c];
#ifdef __SSE2__ #ifdef __SSE2__
@ -1167,7 +1175,7 @@ float* RawImageSource::CA_correct_RT(
STVFU(rgb[c][indx >> 1], RBint); STVFU(rgb[c][indx >> 1], RBint);
} }
#endif #endif
for (int c = FC(rr, cc), indx = rr * ts + cc; cc < cc1 - 8; cc += 2, indx += 2) { for (int c = fc(cfa, rr, cc), indx = rr * ts + cc; cc < cc1 - 8; cc += 2, indx += 2) {
float grbdiffold = rgb[1][indx] - rgb[c][indx >> 1]; float grbdiffold = rgb[1][indx] - rgb[c][indx >> 1];
//interpolate colour difference from optical R/B locations to grid locations //interpolate colour difference from optical R/B locations to grid locations
@ -1209,9 +1217,9 @@ float* RawImageSource::CA_correct_RT(
// copy CA corrected results to temporary image matrix // copy CA corrected results to temporary image matrix
for (int rr = border; rr < rr1 - border; rr++) { for (int rr = border; rr < rr1 - border; rr++) {
int c = FC(rr + top, left + border + (FC(rr + top, 2) & 1)); int c = fc(cfa, rr + top, left + border + (fc(cfa, rr + top, 2) & 1));
int row = rr + top; int row = rr + top;
int cc = border + (FC(rr, 2) & 1); int cc = border + (fc(cfa, rr, 2) & 1);
int indx = (row * width + cc + left) >> 1; int indx = (row * width + cc + left) >> 1;
int indx1 = (rr * ts + cc) >> 1; int indx1 = (rr * ts + cc) >> 1;
#ifdef __SSE2__ #ifdef __SSE2__
@ -1246,7 +1254,7 @@ float* RawImageSource::CA_correct_RT(
#endif #endif
for (int row = cb; row < height - cb; row++) { for (int row = cb; row < height - cb; row++) {
int col = cb + (FC(row, 0) & 1); int col = cb + (fc(cfa, row, 0) & 1);
int indx = (row * width + col) >> 1; int indx = (row * width + col) >> 1;
#ifdef __SSE2__ #ifdef __SSE2__
for (; col < width - 7 - cb; col += 8, indx += 4) { for (; col < width - 7 - cb; col += 8, indx += 4) {
@ -1281,8 +1289,8 @@ float* RawImageSource::CA_correct_RT(
#pragma omp for #pragma omp for
#endif #endif
for (int i = 0; i < H - 2 * cb; ++i) { for (int i = 0; i < H - 2 * cb; ++i) {
const int firstCol = FC(i, 0) & 1; const int firstCol = fc(cfa, i, 0) & 1;
const int colour = FC(i, firstCol); const int colour = fc(cfa, i, firstCol);
const array2D<float>* nonGreen = colour == 0 ? redFactor : blueFactor; const array2D<float>* nonGreen = colour == 0 ? redFactor : blueFactor;
int j = firstCol; int j = firstCol;
#ifdef __SSE2__ #ifdef __SSE2__
@ -1314,9 +1322,9 @@ float* RawImageSource::CA_correct_RT(
if (W % 2) { if (W % 2) {
// odd width => factors for one channel are not set in last column => use value of preceding column // odd width => factors for one channel are not set in last column => use value of preceding column
const int ngRow = 1 - (FC(0, 0) & 1); const int ngRow = 1 - (fc(cfa, 0, 0) & 1);
const int ngCol = FC(ngRow, 0) & 1; const int ngCol = fc(cfa, ngRow, 0) & 1;
const int colour = FC(ngRow, ngCol); const int colour = fc(cfa, ngRow, ngCol);
const array2D<float>* nonGreen = colour == 0 ? redFactor : blueFactor; const array2D<float>* nonGreen = colour == 0 ? redFactor : blueFactor;
for (int i = 0; i < (H + 1 - 2 * cb) / 2; ++i) { for (int i = 0; i < (H + 1 - 2 * cb) / 2; ++i) {
(*nonGreen)[i][(W - 2 * cb + 1) / 2 - 1] = (*nonGreen)[i][(W - 2* cb + 1) / 2 - 2]; (*nonGreen)[i][(W - 2 * cb + 1) / 2 - 1] = (*nonGreen)[i][(W - 2* cb + 1) / 2 - 2];
@ -1333,8 +1341,8 @@ float* RawImageSource::CA_correct_RT(
#pragma omp for #pragma omp for
#endif #endif
for (int i = 0; i < H - 2 * cb; ++i) { for (int i = 0; i < H - 2 * cb; ++i) {
const int firstCol = FC(i, 0) & 1; const int firstCol = fc(cfa, i, 0) & 1;
const int colour = FC(i, firstCol); const int colour = fc(cfa, i, firstCol);
const array2D<float>* nonGreen = colour == 0 ? redFactor : blueFactor; const array2D<float>* nonGreen = colour == 0 ? redFactor : blueFactor;
for (int j = firstCol; j < W - 2 * cb; j += 2) { for (int j = firstCol; j < W - 2 * cb; j += 2) {
rawData[i + cb][j + cb] *= (*nonGreen)[i / 2][j / 2]; rawData[i + cb][j + cb] *= (*nonGreen)[i / 2][j / 2];

15
rtengine/ahd_demosaic_RT.cc
View File

@ -32,6 +32,12 @@
//#define BENCHMARK //#define BENCHMARK
#include "StopWatch.h" #include "StopWatch.h"
namespace
{
unsigned fc(const unsigned int cfa[2][2], int r, int c) {
return cfa[r & 1][c & 1];
}
}
namespace rtengine namespace rtengine
{ {
#define TS 144 #define TS 144
@ -39,6 +45,7 @@ void RawImageSource::ahd_demosaic()
{ {
BENCHFUN BENCHFUN
const unsigned int cfa[2][2] = {{FC(0,0), FC(0,1)}, {FC(1,0), FC(1,1)}};
constexpr int dirs[4] = { -1, 1, -TS, TS }; constexpr int dirs[4] = { -1, 1, -TS, TS };
float xyz_cam[3][3]; float xyz_cam[3][3];
LUTf cbrt(65536); LUTf cbrt(65536);
@ -73,7 +80,7 @@ void RawImageSource::ahd_demosaic()
} }
} }
} }
border_interpolate2(W, H, 5, rawData, red, green, blue); border_interpolate(W, H, 5, rawData, red, green, blue);
#ifdef _OPENMP #ifdef _OPENMP
@ -93,7 +100,7 @@ void RawImageSource::ahd_demosaic()
for (int left = 2; left < width - 5; left += TS - 6) { for (int left = 2; left < width - 5; left += TS - 6) {
// Interpolate green horizontally and vertically: // Interpolate green horizontally and vertically:
for (int row = top; row < top + TS && row < height - 2; row++) { for (int row = top; row < top + TS && row < height - 2; row++) {
for (int col = left + (FC(row, left) & 1); col < std::min(left + TS, width - 2); col += 2) { for (int col = left + (fc(cfa, row, left) & 1); col < std::min(left + TS, width - 2); col += 2) {
auto pix = &rawData[row][col]; auto pix = &rawData[row][col];
float val0 = 0.25f * ((pix[-1] + pix[0] + pix[1]) * 2 float val0 = 0.25f * ((pix[-1] + pix[0] + pix[1]) * 2
- pix[-2] - pix[2]) ; - pix[-2] - pix[2]) ;
@ -107,12 +114,12 @@ void RawImageSource::ahd_demosaic()
// Interpolate red and blue, and convert to CIELab: // Interpolate red and blue, and convert to CIELab:
for (int d = 0; d < 2; d++) for (int d = 0; d < 2; d++)
for (int row = top + 1; row < top + TS - 1 && row < height - 3; row++) { for (int row = top + 1; row < top + TS - 1 && row < height - 3; row++) {
int cng = FC(row + 1, FC(row + 1, 0) & 1); int cng = fc(cfa, row + 1, fc(cfa, row + 1, 0) & 1);
for (int col = left + 1; col < std::min(left + TS - 1, width - 3); col++) { for (int col = left + 1; col < std::min(left + TS - 1, width - 3); col++) {
auto pix = &rawData[row][col]; auto pix = &rawData[row][col];
auto rix = &rgb[d][row - top][col - left]; auto rix = &rgb[d][row - top][col - left];
auto lix = lab[d][row - top][col - left]; auto lix = lab[d][row - top][col - left];
if (FC(row, col) == 1) { if (fc(cfa, row, col) == 1) {
rix[0][2 - cng] = CLIP(pix[0] + (0.5f * (pix[-1] + pix[1] rix[0][2 - cng] = CLIP(pix[0] + (0.5f * (pix[-1] + pix[1]
- rix[-1][1] - rix[1][1] ) )); - rix[-1][1] - rix[1][1] ) ));
rix[0][cng] = CLIP(pix[0] + (0.5f * (pix[-width] + pix[width] rix[0][cng] = CLIP(pix[0] + (0.5f * (pix[-width] + pix[width]

64
rtengine/amaze_demosaic_RT.cc
View File

@ -35,6 +35,13 @@
#include "median.h" #include "median.h"
#include "StopWatch.h" #include "StopWatch.h"
namespace
{
unsigned fc(const unsigned int cfa[2][2], int r, int c) {
return cfa[r & 1][c & 1];
}
}
namespace rtengine namespace rtengine
{ {
@ -55,6 +62,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
plistener->setProgress(progress); plistener->setProgress(progress);
} }
const unsigned int cfarray[2][2] = {{FC(0,0), FC(0,1)}, {FC(1,0), FC(1,1)}};
const int width = winw, height = winh; const int width = winw, height = winh;
const float clip_pt = 1.0 / initialGain; const float clip_pt = 1.0 / initialGain;
const float clip_pt8 = 0.8 / initialGain; const float clip_pt8 = 0.8 / initialGain;
@ -73,8 +81,8 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
int ex, ey; int ex, ey;
//determine GRBG coset; (ey,ex) is the offset of the R subarray //determine GRBG coset; (ey,ex) is the offset of the R subarray
if (FC(0, 0) == 1) { //first pixel is G if (fc(cfarray, 0, 0) == 1) { //first pixel is G
if (FC(0, 1) == 0) { if (fc(cfarray, 0, 1) == 0) {
ey = 0; ey = 0;
ex = 1; ex = 1;
} else { } else {
@ -82,7 +90,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
ex = 0; ex = 0;
} }
} else {//first pixel is R or B } else {//first pixel is R or B
if (FC(0, 0) == 0) { if (fc(cfarray, 0, 0) == 0) {
ey = 0; ey = 0;
ex = 0; ex = 0;
} else { } else {
@ -372,7 +380,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#ifdef __SSE2__ #ifdef __SSE2__
vfloat sgnv; vfloat sgnv;
if( !(FC(4, 4) & 1) ) { if( !(fc(cfarray, 4, 4) & 1) ) {
sgnv = _mm_set_ps( 1.f, -1.f, 1.f, -1.f ); sgnv = _mm_set_ps( 1.f, -1.f, 1.f, -1.f );
} else { } else {
sgnv = _mm_set_ps( -1.f, 1.f, -1.f, 1.f ); sgnv = _mm_set_ps( -1.f, 1.f, -1.f, 1.f );
@ -439,7 +447,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
for (int rr = 4; rr < rr1 - 4; rr++) { for (int rr = 4; rr < rr1 - 4; rr++) {
bool fcswitch = FC(rr, 4) & 1; bool fcswitch = fc(cfarray, rr, 4) & 1;
for (int cc = 4, indx = rr * ts + cc; cc < cc1 - 4; cc++, indx++) { for (int cc = 4, indx = rr * ts + cc; cc < cc1 - 4; cc++, indx++) {
@ -531,7 +539,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
vfloat clip_ptv = F2V( clip_pt ); vfloat clip_ptv = F2V( clip_pt );
vfloat sgn3v; vfloat sgn3v;
if( !(FC(4, 4) & 1) ) { if( !(fc(cfarray, 4, 4) & 1) ) {
sgnv = _mm_set_ps( 1.f, -1.f, 1.f, -1.f ); sgnv = _mm_set_ps( 1.f, -1.f, 1.f, -1.f );
} else { } else {
sgnv = _mm_set_ps( -1.f, 1.f, -1.f, 1.f ); sgnv = _mm_set_ps( -1.f, 1.f, -1.f, 1.f );
@ -589,7 +597,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
for (int rr = 4; rr < rr1 - 4; rr++) { for (int rr = 4; rr < rr1 - 4; rr++) {
for (int cc = 4, indx = rr * ts + cc, c = FC(rr, cc) & 1; cc < cc1 - 4; cc++, indx++) { for (int cc = 4, indx = rr * ts + cc, c = fc(cfarray, rr, cc) & 1; cc < cc1 - 4; cc++, indx++) {
float hcdvar = 3.f * (SQR(hcd[indx - 2]) + SQR(hcd[indx]) + SQR(hcd[indx + 2])) - SQR(hcd[indx - 2] + hcd[indx] + hcd[indx + 2]); float hcdvar = 3.f * (SQR(hcd[indx - 2]) + SQR(hcd[indx]) + SQR(hcd[indx + 2])) - SQR(hcd[indx - 2] + hcd[indx] + hcd[indx + 2]);
float hcdaltvar = 3.f * (SQR(hcdalt[indx - 2]) + SQR(hcdalt[indx]) + SQR(hcdalt[indx + 2])) - SQR(hcdalt[indx - 2] + hcdalt[indx] + hcdalt[indx + 2]); float hcdaltvar = 3.f * (SQR(hcdalt[indx - 2]) + SQR(hcdalt[indx]) + SQR(hcdalt[indx + 2])) - SQR(hcdalt[indx - 2] + hcdalt[indx] + hcdalt[indx + 2]);
float vcdvar = 3.f * (SQR(vcd[indx - v2]) + SQR(vcd[indx]) + SQR(vcd[indx + v2])) - SQR(vcd[indx - v2] + vcd[indx] + vcd[indx + v2]); float vcdvar = 3.f * (SQR(vcd[indx - v2]) + SQR(vcd[indx]) + SQR(vcd[indx + v2])) - SQR(vcd[indx - v2] + vcd[indx] + vcd[indx + v2]);
@ -685,7 +693,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
vfloat epssqv = F2V( epssq ); vfloat epssqv = F2V( epssq );
for (int rr = 6; rr < rr1 - 6; rr++) { for (int rr = 6; rr < rr1 - 6; rr++) {
for (int indx = rr * ts + 6 + (FC(rr, 2) & 1); indx < rr * ts + cc1 - 6; indx += 8) { for (int indx = rr * ts + 6 + (fc(cfarray, rr, 2) & 1); indx < rr * ts + cc1 - 6; indx += 8) {
//compute colour difference variances in cardinal directions //compute colour difference variances in cardinal directions
vfloat tempv = LC2VFU(vcd[indx]); vfloat tempv = LC2VFU(vcd[indx]);
vfloat uavev = tempv + LC2VFU(vcd[indx - v1]) + LC2VFU(vcd[indx - v2]) + LC2VFU(vcd[indx - v3]); vfloat uavev = tempv + LC2VFU(vcd[indx - v1]) + LC2VFU(vcd[indx - v2]) + LC2VFU(vcd[indx - v3]);
@ -731,7 +739,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
for (int rr = 6; rr < rr1 - 6; rr++) { for (int rr = 6; rr < rr1 - 6; rr++) {
for (int cc = 6 + (FC(rr, 2) & 1), indx = rr * ts + cc; cc < cc1 - 6; cc += 2, indx += 2) { for (int cc = 6 + (fc(cfarray, rr, 2) & 1), indx = rr * ts + cc; cc < cc1 - 6; cc += 2, indx += 2) {
//compute colour difference variances in cardinal directions //compute colour difference variances in cardinal directions
@ -793,7 +801,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
// precompute nyquist // precompute nyquist
for (int rr = 6; rr < rr1 - 6; rr++) { for (int rr = 6; rr < rr1 - 6; rr++) {
int cc = 6 + (FC(rr, 2) & 1); int cc = 6 + (fc(cfarray, rr, 2) & 1);
int indx = rr * ts + cc; int indx = rr * ts + cc;
#ifdef __SSE2__ #ifdef __SSE2__
@ -856,7 +864,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
int nyendcol = 0; int nyendcol = 0;
for (int rr = 6; rr < rr1 - 6; rr++) { for (int rr = 6; rr < rr1 - 6; rr++) {
for (int cc = 6 + (FC(rr, 2) & 1), indx = rr * ts + cc; cc < cc1 - 6; cc += 2, indx += 2) { for (int cc = 6 + (fc(cfarray, rr, 2) & 1), indx = rr * ts + cc; cc < cc1 - 6; cc += 2, indx += 2) {
//nyquist texture test: ask if difference of vcd compared to hcd is larger or smaller than RGGB gradients //nyquist texture test: ask if difference of vcd compared to hcd is larger or smaller than RGGB gradients
if(nyqutest[indx >> 1] > 0.f) { if(nyqutest[indx >> 1] > 0.f) {
@ -907,7 +915,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
for (int indx = rr * ts + nystartcol + (FC(rr, 2) & 1); indx < rr * ts + nyendcol; indx += 2) { for (int indx = rr * ts + nystartcol + (fc(cfarray, rr, 2) & 1); indx < rr * ts + nyendcol; indx += 2) {
unsigned int nyquisttemp = (nyquist[(indx - v2) >> 1] + nyquist[(indx - m1) >> 1] + nyquist[(indx + p1) >> 1] + unsigned int nyquisttemp = (nyquist[(indx - v2) >> 1] + nyquist[(indx - m1) >> 1] + nyquist[(indx + p1) >> 1] +
nyquist[(indx - 2) >> 1] + nyquist[(indx + 2) >> 1] + nyquist[(indx - 2) >> 1] + nyquist[(indx + 2) >> 1] +
nyquist[(indx - p1) >> 1] + nyquist[(indx + m1) >> 1] + nyquist[(indx + v2) >> 1]); nyquist[(indx - p1) >> 1] + nyquist[(indx + m1) >> 1] + nyquist[(indx + v2) >> 1]);
@ -922,7 +930,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
// in areas of Nyquist texture, do area interpolation // in areas of Nyquist texture, do area interpolation
for (int rr = nystartrow; rr < nyendrow; rr++) for (int rr = nystartrow; rr < nyendrow; rr++)
for (int indx = rr * ts + nystartcol + (FC(rr, 2) & 1); indx < rr * ts + nyendcol; indx += 2) { for (int indx = rr * ts + nystartcol + (fc(cfarray, rr, 2) & 1); indx < rr * ts + nyendcol; indx += 2) {
if (nyquist2[indx >> 1]) { if (nyquist2[indx >> 1]) {
// area interpolation // area interpolation
@ -962,7 +970,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
//populate G at R/B sites //populate G at R/B sites
for (int rr = 8; rr < rr1 - 8; rr++) for (int rr = 8; rr < rr1 - 8; rr++)
for (int indx = rr * ts + 8 + (FC(rr, 2) & 1); indx < rr * ts + cc1 - 8; indx += 2) { for (int indx = rr * ts + 8 + (fc(cfarray, rr, 2) & 1); indx < rr * ts + cc1 - 8; indx += 2) {
//first ask if one gets more directional discrimination from nearby B/R sites //first ask if one gets more directional discrimination from nearby B/R sites
float hvwtalt = xdivf(hvwt[(indx - m1) >> 1] + hvwt[(indx + p1) >> 1] + hvwt[(indx - p1) >> 1] + hvwt[(indx + m1) >> 1], 2); float hvwtalt = xdivf(hvwt[(indx - m1) >> 1] + hvwt[(indx + p1) >> 1] + hvwt[(indx - p1) >> 1] + hvwt[(indx + m1) >> 1], 2);
@ -985,7 +993,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
// refine Nyquist areas using G curvatures // refine Nyquist areas using G curvatures
if(doNyquist) { if(doNyquist) {
for (int rr = nystartrow; rr < nyendrow; rr++) for (int rr = nystartrow; rr < nyendrow; rr++)
for (int indx = rr * ts + nystartcol + (FC(rr, 2) & 1); indx < rr * ts + nyendcol; indx += 2) { for (int indx = rr * ts + nystartcol + (fc(cfarray, rr, 2) & 1); indx < rr * ts + nyendcol; indx += 2) {
if (nyquist2[indx >> 1]) { if (nyquist2[indx >> 1]) {
//local averages (over Nyquist pixels only) of G curvature squared //local averages (over Nyquist pixels only) of G curvature squared
@ -1008,7 +1016,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#ifdef __SSE2__ #ifdef __SSE2__
for (int rr = 6; rr < rr1 - 6; rr++) { for (int rr = 6; rr < rr1 - 6; rr++) {
if((FC(rr, 2) & 1) == 0) { if((fc(cfarray, rr, 2) & 1) == 0) {
for (int cc = 6, indx = rr * ts + cc; cc < cc1 - 6; cc += 8, indx += 8) { for (int cc = 6, indx = rr * ts + cc; cc < cc1 - 6; cc += 8, indx += 8) {
vfloat tempv = LC2VFU(cfa[indx + 1]); vfloat tempv = LC2VFU(cfa[indx + 1]);
vfloat Dgrbsq1pv = (SQRV(tempv - LC2VFU(cfa[indx + 1 - p1])) + SQRV(tempv - LC2VFU(cfa[indx + 1 + p1]))); vfloat Dgrbsq1pv = (SQRV(tempv - LC2VFU(cfa[indx + 1 - p1])) + SQRV(tempv - LC2VFU(cfa[indx + 1 + p1])));
@ -1034,7 +1042,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
for (int rr = 6; rr < rr1 - 6; rr++) { for (int rr = 6; rr < rr1 - 6; rr++) {
if((FC(rr, 2) & 1) == 0) { if((fc(cfarray, rr, 2) & 1) == 0) {
for (int cc = 6, indx = rr * ts + cc; cc < cc1 - 6; cc += 2, indx += 2) { for (int cc = 6, indx = rr * ts + cc; cc < cc1 - 6; cc += 2, indx += 2) {
delp[indx >> 1] = fabsf(cfa[indx + p1] - cfa[indx - p1]); delp[indx >> 1] = fabsf(cfa[indx + p1] - cfa[indx - p1]);
delm[indx >> 1] = fabsf(cfa[indx + m1] - cfa[indx - m1]); delm[indx >> 1] = fabsf(cfa[indx + m1] - cfa[indx - m1]);
@ -1063,7 +1071,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
for (int rr = 8; rr < rr1 - 8; rr++) { for (int rr = 8; rr < rr1 - 8; rr++) {
#ifdef __SSE2__ #ifdef __SSE2__
for (int indx = rr * ts + 8 + (FC(rr, 2) & 1), indx1 = indx >> 1; indx < rr * ts + cc1 - 8; indx += 8, indx1 += 4) { for (int indx = rr * ts + 8 + (fc(cfarray, rr, 2) & 1), indx1 = indx >> 1; indx < rr * ts + cc1 - 8; indx += 8, indx1 += 4) {
//diagonal colour ratios //diagonal colour ratios
vfloat cfav = LC2VFU(cfa[indx]); vfloat cfav = LC2VFU(cfa[indx]);
@ -1128,7 +1136,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
for (int cc = 8 + (FC(rr, 2) & 1), indx = rr * ts + cc, indx1 = indx >> 1; cc < cc1 - 8; cc += 2, indx += 2, indx1++) { for (int cc = 8 + (fc(cfarray, rr, 2) & 1), indx = rr * ts + cc, indx1 = indx >> 1; cc < cc1 - 8; cc += 2, indx += 2, indx1++) {
//diagonal colour ratios //diagonal colour ratios
float crse = xmul2f(cfa[indx + m1]) / (eps + cfa[indx] + (cfa[indx + m2])); float crse = xmul2f(cfa[indx + m1]) / (eps + cfa[indx] + (cfa[indx + m2]));
@ -1218,7 +1226,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
for (int rr = 10; rr < rr1 - 10; rr++) for (int rr = 10; rr < rr1 - 10; rr++)
#ifdef __SSE2__ #ifdef __SSE2__
for (int indx = rr * ts + 10 + (FC(rr, 2) & 1), indx1 = indx >> 1; indx < rr * ts + cc1 - 10; indx += 8, indx1 += 4) { for (int indx = rr * ts + 10 + (fc(cfarray, rr, 2) & 1), indx1 = indx >> 1; indx < rr * ts + cc1 - 10; indx += 8, indx1 += 4) {
//first ask if one gets more directional discrimination from nearby B/R sites //first ask if one gets more directional discrimination from nearby B/R sites
vfloat pmwtaltv = zd25v * (LVFU(pmwt[(indx - m1) >> 1]) + LVFU(pmwt[(indx + p1) >> 1]) + LVFU(pmwt[(indx - p1) >> 1]) + LVFU(pmwt[(indx + m1) >> 1])); vfloat pmwtaltv = zd25v * (LVFU(pmwt[(indx - m1) >> 1]) + LVFU(pmwt[(indx + p1) >> 1]) + LVFU(pmwt[(indx - p1) >> 1]) + LVFU(pmwt[(indx + m1) >> 1]));
@ -1230,7 +1238,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
for (int cc = 10 + (FC(rr, 2) & 1), indx = rr * ts + cc, indx1 = indx >> 1; cc < cc1 - 10; cc += 2, indx += 2, indx1++) { for (int cc = 10 + (fc(cfarray, rr, 2) & 1), indx = rr * ts + cc, indx1 = indx >> 1; cc < cc1 - 10; cc += 2, indx += 2, indx1++) {
//first ask if one gets more directional discrimination from nearby B/R sites //first ask if one gets more directional discrimination from nearby B/R sites
float pmwtalt = xdivf(pmwt[(indx - m1) >> 1] + pmwt[(indx + p1) >> 1] + pmwt[(indx - p1) >> 1] + pmwt[(indx + m1) >> 1], 2); float pmwtalt = xdivf(pmwt[(indx - m1) >> 1] + pmwt[(indx + p1) >> 1] + pmwt[(indx - p1) >> 1] + pmwt[(indx + m1) >> 1], 2);
@ -1246,7 +1254,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
for (int rr = 12; rr < rr1 - 12; rr++) for (int rr = 12; rr < rr1 - 12; rr++)
#ifdef __SSE2__ #ifdef __SSE2__
for (int indx = rr * ts + 12 + (FC(rr, 2) & 1), indx1 = indx >> 1; indx < rr * ts + cc1 - 12; indx += 8, indx1 += 4) { for (int indx = rr * ts + 12 + (fc(cfarray, rr, 2) & 1), indx1 = indx >> 1; indx < rr * ts + cc1 - 12; indx += 8, indx1 += 4) {
vmask copymask = vmaskf_ge(vabsf(zd5v - LVFU(pmwt[indx1])), vabsf(zd5v - LVFU(hvwt[indx1]))); vmask copymask = vmaskf_ge(vabsf(zd5v - LVFU(pmwt[indx1])), vabsf(zd5v - LVFU(hvwt[indx1])));
if(_mm_movemask_ps((vfloat)copymask)) { // if for any of the 4 pixels the condition is true, do the maths for all 4 pixels and mask the unused out at the end if(_mm_movemask_ps((vfloat)copymask)) { // if for any of the 4 pixels the condition is true, do the maths for all 4 pixels and mask the unused out at the end
@ -1301,7 +1309,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
for (int cc = 12 + (FC(rr, 2) & 1), indx = rr * ts + cc, indx1 = indx >> 1; cc < cc1 - 12; cc += 2, indx += 2, indx1++) { for (int cc = 12 + (fc(cfarray, rr, 2) & 1), indx = rr * ts + cc, indx1 = indx >> 1; cc < cc1 - 12; cc += 2, indx += 2, indx1++) {
if (fabsf(0.5f - pmwt[indx >> 1]) < fabsf(0.5f - hvwt[indx >> 1]) ) { if (fabsf(0.5f - pmwt[indx >> 1]) < fabsf(0.5f - hvwt[indx >> 1]) ) {
continue; continue;
@ -1399,7 +1407,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
for (int rr = 14; rr < rr1 - 14; rr++) for (int rr = 14; rr < rr1 - 14; rr++)
#ifdef __SSE2__ #ifdef __SSE2__
for (int cc = 14 + (FC(rr, 2) & 1), indx = rr * ts + cc, c = 1 - FC(rr, cc) / 2; cc < cc1 - 14; cc += 8, indx += 8) { for (int cc = 14 + (fc(cfarray, rr, 2) & 1), indx = rr * ts + cc, c = 1 - fc(cfarray, rr, cc) / 2; cc < cc1 - 14; cc += 8, indx += 8) {
vfloat tempv = epsv + vabsf(LVFU(Dgrb[c][(indx - m1) >> 1]) - LVFU(Dgrb[c][(indx + m1) >> 1])); vfloat tempv = epsv + vabsf(LVFU(Dgrb[c][(indx - m1) >> 1]) - LVFU(Dgrb[c][(indx + m1) >> 1]));
vfloat temp2v = epsv + vabsf(LVFU(Dgrb[c][(indx + p1) >> 1]) - LVFU(Dgrb[c][(indx - p1) >> 1])); vfloat temp2v = epsv + vabsf(LVFU(Dgrb[c][(indx + p1) >> 1]) - LVFU(Dgrb[c][(indx - p1) >> 1]));
vfloat wtnwv = onev / (tempv + vabsf(LVFU(Dgrb[c][(indx - m1) >> 1]) - LVFU(Dgrb[c][(indx - m3) >> 1])) + vabsf(LVFU(Dgrb[c][(indx + m1) >> 1]) - LVFU(Dgrb[c][(indx - m3) >> 1]))); vfloat wtnwv = onev / (tempv + vabsf(LVFU(Dgrb[c][(indx - m1) >> 1]) - LVFU(Dgrb[c][(indx - m3) >> 1])) + vabsf(LVFU(Dgrb[c][(indx + m1) >> 1]) - LVFU(Dgrb[c][(indx - m3) >> 1])));
@ -1415,7 +1423,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
for (int cc = 14 + (FC(rr, 2) & 1), indx = rr * ts + cc, c = 1 - FC(rr, cc) / 2; cc < cc1 - 14; cc += 2, indx += 2) { for (int cc = 14 + (fc(cfarray, rr, 2) & 1), indx = rr * ts + cc, c = 1 - fc(cfarray, rr, cc) / 2; cc < cc1 - 14; cc += 2, indx += 2) {
float wtnw = 1.f / (eps + fabsf(Dgrb[c][(indx - m1) >> 1] - Dgrb[c][(indx + m1) >> 1]) + fabsf(Dgrb[c][(indx - m1) >> 1] - Dgrb[c][(indx - m3) >> 1]) + fabsf(Dgrb[c][(indx + m1) >> 1] - Dgrb[c][(indx - m3) >> 1])); float wtnw = 1.f / (eps + fabsf(Dgrb[c][(indx - m1) >> 1] - Dgrb[c][(indx + m1) >> 1]) + fabsf(Dgrb[c][(indx - m1) >> 1] - Dgrb[c][(indx - m3) >> 1]) + fabsf(Dgrb[c][(indx + m1) >> 1] - Dgrb[c][(indx - m3) >> 1]));
float wtne = 1.f / (eps + fabsf(Dgrb[c][(indx + p1) >> 1] - Dgrb[c][(indx - p1) >> 1]) + fabsf(Dgrb[c][(indx + p1) >> 1] - Dgrb[c][(indx + p3) >> 1]) + fabsf(Dgrb[c][(indx - p1) >> 1] - Dgrb[c][(indx + p3) >> 1])); float wtne = 1.f / (eps + fabsf(Dgrb[c][(indx + p1) >> 1] - Dgrb[c][(indx - p1) >> 1]) + fabsf(Dgrb[c][(indx + p1) >> 1] - Dgrb[c][(indx + p3) >> 1]) + fabsf(Dgrb[c][(indx - p1) >> 1] - Dgrb[c][(indx + p3) >> 1]));
float wtsw = 1.f / (eps + fabsf(Dgrb[c][(indx - p1) >> 1] - Dgrb[c][(indx + p1) >> 1]) + fabsf(Dgrb[c][(indx - p1) >> 1] - Dgrb[c][(indx + m3) >> 1]) + fabsf(Dgrb[c][(indx + p1) >> 1] - Dgrb[c][(indx - p3) >> 1])); float wtsw = 1.f / (eps + fabsf(Dgrb[c][(indx - p1) >> 1] - Dgrb[c][(indx + p1) >> 1]) + fabsf(Dgrb[c][(indx - p1) >> 1] - Dgrb[c][(indx + m3) >> 1]) + fabsf(Dgrb[c][(indx + p1) >> 1] - Dgrb[c][(indx - p3) >> 1]));
@ -1434,7 +1442,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
vfloat twov = F2V(2.f); vfloat twov = F2V(2.f);
vmask selmask; vmask selmask;
if((FC(16, 2) & 1) == 1) { if((fc(cfarray, 16, 2) & 1) == 1) {
selmask = _mm_set_epi32(0xffffffff, 0, 0xffffffff, 0); selmask = _mm_set_epi32(0xffffffff, 0, 0xffffffff, 0);
offset = 1; offset = 1;
} else { } else {
@ -1509,7 +1517,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
#else #else
if((FC(rr, 2) & 1) == 1) { if((fc(cfarray, rr, 2) & 1) == 1) {
for (; indx < rr * ts + cc1 - 16 - (cc1 & 1); indx++, col++) { for (; indx < rr * ts + cc1 - 16 - (cc1 & 1); indx++, col++) {
float temp = 1.f / (hvwt[(indx - v1) >> 1] + 2.f - hvwt[(indx + 1) >> 1] - hvwt[(indx - 1) >> 1] + hvwt[(indx + v1) >> 1]); float temp = 1.f / (hvwt[(indx - v1) >> 1] + 2.f - hvwt[(indx + 1) >> 1] - hvwt[(indx - 1) >> 1] + hvwt[(indx + v1) >> 1]);
red[row][col] = 65535.f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[0][(indx - v1) >> 1] + (1.f - hvwt[(indx + 1) >> 1]) * Dgrb[0][(indx + 1) >> 1] + (1.f - hvwt[(indx - 1) >> 1]) * Dgrb[0][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[0][(indx + v1) >> 1]) * red[row][col] = 65535.f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[0][(indx - v1) >> 1] + (1.f - hvwt[(indx + 1) >> 1]) * Dgrb[0][(indx + 1) >> 1] + (1.f - hvwt[(indx - 1) >> 1]) * Dgrb[0][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[0][(indx + v1) >> 1]) *
@ -1591,7 +1599,7 @@ void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, int winh, c
free(buffer); free(buffer);
} }
if(border < 4) { if(border < 4) {
border_interpolate2(W, H, 3, rawData, red, green, blue); border_interpolate(W, H, 3, rawData, red, green, blue);
} }
if(plistener) { if(plistener) {

10
rtengine/badpixels.cc
View File

@ -23,6 +23,13 @@
#include "rawimage.h" #include "rawimage.h"
#include "rawimagesource.h" #include "rawimagesource.h"
namespace
{
unsigned fc(const unsigned int cfa[2][2], int r, int c) {
return cfa[r & 1][c & 1];
}
}
namespace rtengine namespace rtengine
{ {
@ -31,6 +38,7 @@ namespace rtengine
*/ */
int RawImageSource::interpolateBadPixelsBayer(const PixelsMap &bitmapBads, array2D<float> &rawData) int RawImageSource::interpolateBadPixelsBayer(const PixelsMap &bitmapBads, array2D<float> &rawData)
{ {
const unsigned int cfarray[2][2] = {{FC(0,0), FC(0,1)}, {FC(1,0), FC(1,1)}};
constexpr float eps = 1.f; constexpr float eps = 1.f;
int counter = 0; int counter = 0;
@ -54,7 +62,7 @@ int RawImageSource::interpolateBadPixelsBayer(const PixelsMap &bitmapBads, array
float wtdsum = 0.f, norm = 0.f; float wtdsum = 0.f, norm = 0.f;
// diagonal interpolation // diagonal interpolation
if (FC(row, col) == 1) { if (fc(cfarray, row, col) == 1) {
// green channel. We can use closer pixels than for red or blue channel. Distance to center pixel is sqrt(2) => weighting is 0.70710678 // green channel. We can use closer pixels than for red or blue channel. Distance to center pixel is sqrt(2) => weighting is 0.70710678
// For green channel following pixels will be used for interpolation. Pixel to be interpolated is in center. // For green channel following pixels will be used for interpolation. Pixel to be interpolated is in center.
// 1 means that pixel is used in this step, if itself and his counterpart are not marked bad // 1 means that pixel is used in this step, if itself and his counterpart are not marked bad

287
rtengine/demosaic_algos.cc
View File

@ -55,158 +55,7 @@ namespace rtengine
#define FORCC for (unsigned int c=0; c < colors; c++) #define FORCC for (unsigned int c=0; c < colors; c++)
/* void RawImageSource::border_interpolate( int winw, int winh, int lborders, const array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue)
Patterned Pixel Grouping Interpolation by Alain Desbiolles
*/
void RawImageSource::ppg_demosaic()
{
int width = W, height = H;
int dir[5] = { 1, width, -1, -width, 1 };
int row, col, diff[2] = {}, guess[2], c, d, i;
float (*pix)[4];
float (*image)[4];
if (plistener) {
// looks like ppg isn't supported anymore
//plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::ppg)));
plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), M("GENERAL_NA")));
plistener->setProgress (0.0);
}
image = (float (*)[4]) calloc (static_cast<size_t>(H) * W, sizeof * image);
for (int ii = 0; ii < H; ii++)
for (int jj = 0; jj < W; jj++) {
image[ii * W + jj][fc(ii, jj)] = rawData[ii][jj];
}
border_interpolate(3, image);
/* Fill in the green layer with gradients and pattern recognition: */
for (row = 3; row < height - 3; row++) {
for (col = 3 + (FC(row, 3) & 1), c = FC(row, col); col < width - 3; col += 2) {
pix = image + row * width + col;
for (i = 0; (d = dir[i]) > 0; i++) {
guess[i] = (pix[-d][1] + pix[0][c] + pix[d][1]) * 2
- pix[-2 * d][c] - pix[2 * d][c];
diff[i] = ( ABS(pix[-2 * d][c] - pix[ 0][c]) +
ABS(pix[ 2 * d][c] - pix[ 0][c]) +
ABS(pix[ -d][1] - pix[ d][1]) ) * 3 +
( ABS(pix[ 3 * d][1] - pix[ d][1]) +
ABS(pix[-3 * d][1] - pix[-d][1]) ) * 2;
}
d = dir[i = diff[0] > diff[1]];
pix[0][1] = median(static_cast<float>(guess[i] >> 2), pix[d][1], pix[-d][1]);
}
if(plistener) {
plistener->setProgress(0.33 * row / (height - 3));
}
}
/* Calculate red and blue for each green pixel: */
for (row = 1; row < height - 1; row++) {
for (col = 1 + (FC(row, 2) & 1), c = FC(row, col + 1); col < width - 1; col += 2) {
pix = image + row * width + col;
for (i = 0; (d = dir[i]) > 0; c = 2 - c, i++)
pix[0][c] = CLIP(0.5 * (pix[-d][c] + pix[d][c] + 2 * pix[0][1]
- pix[-d][1] - pix[d][1]) );
}
if(plistener) {
plistener->setProgress(0.33 + 0.33 * row / (height - 1));
}
}
/* Calculate blue for red pixels and vice versa: */
for (row = 1; row < height - 1; row++) {
for (col = 1 + (FC(row, 1) & 1), c = 2 - FC(row, col); col < width - 1; col += 2) {
pix = image + row * width + col;
for (i = 0; (d = dir[i] + dir[i + 1]) > 0; i++) {
diff[i] = ABS(pix[-d][c] - pix[d][c]) +
ABS(pix[-d][1] - pix[0][1]) +
ABS(pix[ d][1] - pix[0][1]);
guess[i] = pix[-d][c] + pix[d][c] + 2 * pix[0][1]
- pix[-d][1] - pix[d][1];
}
if (diff[0] != diff[1]) {
pix[0][c] = CLIP(guess[diff[0] > diff[1]] / 2);
} else {
pix[0][c] = CLIP((guess[0] + guess[1]) / 4);
}
}
if(plistener) {
plistener->setProgress(0.67 + 0.33 * row / (height - 1));
}
}
red(W, H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++) {
red[i][j] = image[i * W + j][0];
}
green(W, H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++) {
green[i][j] = image[i * W + j][1];
}
blue(W, H);
for (int i = 0; i < H; i++)
for (int j = 0; j < W; j++) {
blue[i][j] = image[i * W + j][2];
}
free (image);
}
void RawImageSource::border_interpolate(unsigned int border, float (*image)[4], unsigned int start, unsigned int end)
{
unsigned row, col, y, x, f;
float sum[8];
unsigned int width = W, height = H;
unsigned int colors = 3;
if (end == 0 ) {
end = H;
}
for (row = start; row < end; row++)
for (col = 0; col < width; col++) {
if (col == border && row >= border && row < height - border) {
col = width - border;
}
memset (sum, 0, sizeof sum);
for (y = row - 1; y != row + 2; y++)
for (x = col - 1; x != col + 2; x++)
if (y < height && x < width) {
f = fc(y, x);
sum[f] += image[y * width + x][f];
sum[f + 4]++;
}
f = fc(row, col);
FORCC if (c != f && sum[c + 4]) {
image[row * width + col][c] = sum[c] / sum[c + 4];
}
}
}
void RawImageSource::border_interpolate2( int winw, int winh, int lborders, const array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue)
{ {
int bord = lborders; int bord = lborders;
int width = winw; int width = winw;
@ -361,129 +210,6 @@ void RawImageSource::border_interpolate2( int winw, int winh, int lborders, cons
} }
// Joint Demosaicing and Denoising using High Order Interpolation Techniques
// Revision 0.9.1a - 09/02/2010 - Contact info: luis.sanz.rodriguez@gmail.com
// Copyright Luis Sanz Rodriguez 2010
// Adapted to RawTherapee by Jacques Desmis 3/2013
void RawImageSource::jdl_interpolate_omp() // from "Lassus"
{
int width = W, height = H;
int row, col, c, d, i, u = width, v = 2 * u, w = 3 * u, x = 4 * u, y = 5 * u, z = 6 * u, indx, (*dif)[2], (*chr)[2];
float f[4], g[4];
float (*image)[4];
image = (float (*)[4]) calloc (static_cast<size_t>(width) * height, sizeof * image);
dif = (int (*)[2]) calloc(static_cast<size_t>(width) * height, sizeof * dif);
chr = (int (*)[2]) calloc(static_cast<size_t>(width) * height, sizeof * chr);
if (plistener) {
// this function seems to be unused
//plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), RAWParams::BayerSensor::getMethodString(RAWParams::BayerSensor::jdl)));
plistener->setProgressStr (Glib::ustring::compose(M("TP_RAW_DMETHOD_PROGRESSBAR"), M("GENERAL_NA")));
plistener->setProgress (0.0);
}
#ifdef _OPENMP
#pragma omp parallel shared(image,width,height,u,w,v,y,x,z,dif,chr) private(row,col,f,g,indx,c,d,i)
#endif
{
#ifdef _OPENMP
#pragma omp for
#endif
for (int ii = 0; ii < height; ii++)
for (int jj = 0; jj < width; jj++) {
image[ii * width + jj][fc(ii, jj)] = rawData[ii][jj];
}
border_interpolate(6, image);
#ifdef _OPENMP
#pragma omp for
#endif
for (row = 5; row < height - 5; row++)
for (col = 5 + (FC(row, 1) & 1), indx = row * width + col, c = FC(row, col); col < u - 5; col += 2, indx += 2) {
f[0] = 1.f + abs(image[indx - u][1] - image[indx - w][1]) + abs(image[indx - u][1] - image[indx + u][1]) + abs(image[indx][c] - image[indx - v][c]) + abs(image[indx - v][c] - image[indx - x][c]);
f[1] = 1.f + abs(image[indx + 1][1] - image[indx + 3][1]) + abs(image[indx + 1][1] - image[indx - 1][1]) + abs(image[indx][c] - image[indx + 2][c]) + abs(image[indx + 2][c] - image[indx + 4][c]);
f[2] = 1.f + abs(image[indx - 1][1] - image[indx - 3][1]) + abs(image[indx - 1][1] - image[indx + 1][1]) + abs(image[indx][c] - image[indx - 2][c]) + abs(image[indx - 2][c] - image[indx - 4][c]);
f[3] = 1.f + abs(image[indx + u][1] - image[indx + w][1]) + abs(image[indx + u][1] - image[indx - u][1]) + abs(image[indx][c] - image[indx + v][c]) + abs(image[indx + v][c] - image[indx + x][c]);
g[0] = CLIP((22.f * image[indx - u][1] + 22.f * image[indx - w][1] + 2.f * image[indx - y][1] + 2.f * image[indx + u][1] + 40.f * image[indx][c] - 32.f * image[indx - v][c] - 8.f * image[indx - x][c]) / 48.f);
g[1] = CLIP((22.f * image[indx + 1][1] + 22.f * image[indx + 3][1] + 2.f * image[indx + 5][1] + 2.f * image[indx - 1][1] + 40.f * image[indx][c] - 32.f * image[indx + 2][c] - 8.f * image[indx + 4][c]) / 48.f);
g[2] = CLIP((22.f * image[indx - 1][1] + 22.f * image[indx - 3][1] + 2.f * image[indx - 5][1] + 2.f * image[indx + 1][1] + 40.f * image[indx][c] - 32.f * image[indx - 2][c] - 8.f * image[indx - 4][c]) / 48.f);
g[3] = CLIP((22.f * image[indx + u][1] + 22.f * image[indx + w][1] + 2.f * image[indx + y][1] + 2.f * image[indx - u][1] + 40.f * image[indx][c] - 32.f * image[indx + v][c] - 8.f * image[indx + x][c]) / 48.f);
dif[indx][0] = CLIP((f[3] * g[0] + f[0] * g[3]) / (f[0] + f[3])) - image[indx][c];
dif[indx][1] = CLIP((f[2] * g[1] + f[1] * g[2]) / (f[1] + f[2])) - image[indx][c];
}
#ifdef _OPENMP
#pragma omp for
#endif
for (row = 6; row < height - 6; row++)
for (col = 6 + (FC(row, 2) & 1), indx = row * width + col, c = FC(row, col) / 2; col < u - 6; col += 2, indx += 2) {
f[0] = 1.f + 78.f * SQR((float)dif[indx][0]) + 69.f * (SQR((float) dif[indx - v][0]) + SQR((float)dif[indx + v][0])) + 51.f * (SQR((float)dif[indx - x][0]) + SQR((float)dif[indx + x][0])) + 21.f * (SQR((float)dif[indx - z][0]) + SQR((float)dif[indx + z][0])) - 6.f * SQR((float)dif[indx - v][0] + dif[indx][0] + dif[indx + v][0]) - 10.f * (SQR((float)dif[indx - x][0] + dif[indx - v][0] + dif[indx][0]) + SQR((float)dif[indx][0] + dif[indx + v][0] + dif[indx + x][0])) - 7.f * (SQR((float)dif[indx - z][0] + dif[indx - x][0] + dif[indx - v][0]) + SQR((float)dif[indx + v][0] + dif[indx + x][0] + dif[indx + z][0]));
f[1] = 1.f + 78.f * SQR((float)dif[indx][1]) + 69.f * (SQR((float)dif[indx - 2][1]) + SQR((float)dif[indx + 2][1])) + 51.f * (SQR((float)dif[indx - 4][1]) + SQR((float)dif[indx + 4][1])) + 21.f * (SQR((float)dif[indx - 6][1]) + SQR((float)dif[indx + 6][1])) - 6.f * SQR((float)dif[indx - 2][1] + dif[indx][1] + dif[indx + 2][1]) - 10.f * (SQR((float)dif[indx - 4][1] + dif[indx - 2][1] + dif[indx][1]) + SQR((float)dif[indx][1] + dif[indx + 2][1] + dif[indx + 4][1])) - 7.f * (SQR((float)dif[indx - 6][1] + dif[indx - 4][1] + dif[indx - 2][1]) + SQR((float)dif[indx + 2][1] + dif[indx + 4][1] + dif[indx + 6][1]));
g[0] = median(0.725f * dif[indx][0] + 0.1375f * dif[indx - v][0] + 0.1375f * dif[indx + v][0], static_cast<float>(dif[indx - v][0]), static_cast<float>(dif[indx + v][0]));
g[1] = median(0.725f * dif[indx][1] + 0.1375f * dif[indx - 2][1] + 0.1375f * dif[indx + 2][1], static_cast<float>(dif[indx - 2][1]), static_cast<float>(dif[indx + 2][1]));
chr[indx][c] = (f[1] * g[0] + f[0] * g[1]) / (f[0] + f[1]);
}
#ifdef _OPENMP
#pragma omp for
#endif
for (row = 6; row < height - 6; row++)
for (col = 6 + (FC(row, 2) & 1), indx = row * width + col, c = 1 - FC(row, col) / 2, d = 2 * c; col < u - 6; col += 2, indx += 2) {
f[0] = 1.f / (float)(1.f + fabs((float)chr[indx - u - 1][c] - chr[indx + u + 1][c]) + fabs((float)chr[indx - u - 1][c] - chr[indx - w - 3][c]) + fabs((float)chr[indx + u + 1][c] - chr[indx - w - 3][c]));
f[1] = 1.f / (float)(1.f + fabs((float)chr[indx - u + 1][c] - chr[indx + u - 1][c]) + fabs((float)chr[indx - u + 1][c] - chr[indx - w + 3][c]) + fabs((float)chr[indx + u - 1][c] - chr[indx - w + 3][c]));
f[2] = 1.f / (float)(1.f + fabs((float)chr[indx + u - 1][c] - chr[indx - u + 1][c]) + fabs((float)chr[indx + u - 1][c] - chr[indx + w + 3][c]) + fabs((float)chr[indx - u + 1][c] - chr[indx + w - 3][c]));
f[3] = 1.f / (float)(1.f + fabs((float)chr[indx + u + 1][c] - chr[indx - u - 1][c]) + fabs((float)chr[indx + u + 1][c] - chr[indx + w - 3][c]) + fabs((float)chr[indx - u - 1][c] - chr[indx + w + 3][c]));
g[0] = median(chr[indx - u - 1][c], chr[indx - w - 1][c], chr[indx - u - 3][c]);
g[1] = median(chr[indx - u + 1][c], chr[indx - w + 1][c], chr[indx - u + 3][c]);
g[2] = median(chr[indx + u - 1][c], chr[indx + w - 1][c], chr[indx + u - 3][c]);
g[3] = median(chr[indx + u + 1][c], chr[indx + w + 1][c], chr[indx + u + 3][c]);
chr[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]);
image[indx][1] = CLIP(image[indx][2 - d] + chr[indx][1 - c]);
image[indx][d] = CLIP(image[indx][1] - chr[indx][c]);
}
#ifdef _OPENMP
#pragma omp for
#endif
for (row = 6; row < height - 6; row++)
for (col = 6 + (FC(row, 1) & 1), indx = row * width + col, c = FC(row, col + 1) / 2, d = 2 * c; col < u - 6; col += 2, indx += 2)
for(i = 0; i <= 1; c = 1 - c, d = 2 * c, i++) {
f[0] = 1.f / (float)(1.f + fabs((float)chr[indx - u][c] - chr[indx + u][c]) + fabs((float)chr[indx - u][c] - chr[indx - w][c]) + fabs((float)chr[indx + u][c] - chr[indx - w][c]));
f[1] = 1.f / (float)(1.f + fabs((float)chr[indx + 1][c] - chr[indx - 1][c]) + fabs((float)chr[indx + 1][c] - chr[indx + 3][c]) + fabs((float)chr[indx - 1][c] - chr[indx + 3][c]));
f[2] = 1.f / (float)(1.f + fabs((float)chr[indx - 1][c] - chr[indx + 1][c]) + fabs((float)chr[indx - 1][c] - chr[indx - 3][c]) + fabs((float)chr[indx + 1][c] - chr[indx - 3][c]));
f[3] = 1.f / (float)(1.f + fabs((float)chr[indx + u][c] - chr[indx - u][c]) + fabs((float)chr[indx + u][c] - chr[indx + w][c]) + fabs((float)chr[indx - u][c] - chr[indx + w][c]));
g[0] = 0.875f * chr[indx - u][c] + 0.125f * chr[indx - w][c];
g[1] = 0.875f * chr[indx + 1][c] + 0.125f * chr[indx + 3][c];
g[2] = 0.875f * chr[indx - 1][c] + 0.125f * chr[indx - 3][c];
g[3] = 0.875f * chr[indx + u][c] + 0.125f * chr[indx + w][c];
image[indx][d] = CLIP(image[indx][1] - (f[0] * g[0] + f[1] * g[1] + f[2] * g[2] + f[3] * g[3]) / (f[0] + f[1] + f[2] + f[3]));
}
#ifdef _OPENMP
#pragma omp for
#endif
for (int ii = 0; ii < height; ii++) {
for (int jj = 0; jj < width; jj++) {
red[ii][jj] = CLIP(image[ii * width + jj][0]);
green[ii][jj] = CLIP(image[ii * width + jj][1]);
blue[ii][jj] = CLIP(image[ii * width + jj][2]);
}
}
} // End of parallelization
free (image);
free(dif);
free(chr);
//RawImageSource::refinement_lassus();
}
// LSMME demosaicing algorithm // LSMME demosaicing algorithm
// L. Zhang and X. Wu, // L. Zhang and X. Wu,
// Color demozaicing via directional Linear Minimum Mean Square-error Estimation, // Color demozaicing via directional Linear Minimum Mean Square-error Estimation,
@ -1503,7 +1229,7 @@ void RawImageSource::igv_interpolate(int winw, int winh)
} }
} }
}// End of parallelization }// End of parallelization
border_interpolate2(winw, winh, 8, rawData, red, green, blue); border_interpolate(winw, winh, 8, rawData, red, green, blue);
if (plistener) { if (plistener) {
plistener->setProgress (1.0); plistener->setProgress (1.0);
@ -1560,8 +1286,6 @@ void RawImageSource::igv_interpolate(int winw, int winh)
rgb[c][indx] = CLIP(rawData[row][col]); //rawData = RT data rgb[c][indx] = CLIP(rawData[row][col]); //rawData = RT data
} }
// border_interpolate2(7, rgb);
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp single #pragma omp single
#endif #endif
@ -1736,9 +1460,6 @@ void RawImageSource::igv_interpolate(int winw, int winh)
if (plistener) { if (plistener) {
plistener->setProgress (0.91); plistener->setProgress (0.91);
} }
//Interpolate borders
// border_interpolate2(7, rgb);
} }
/* /*
#ifdef _OPENMP #ifdef _OPENMP
@ -1766,7 +1487,7 @@ void RawImageSource::igv_interpolate(int winw, int winh)
blue [row][col] = CLIP(rgb[1][indx] - 65535.f * chr[1][indx]); blue [row][col] = CLIP(rgb[1][indx] - 65535.f * chr[1][indx]);
} }
}// End of parallelization }// End of parallelization
border_interpolate2(winw, winh, 8, rawData, red, green, blue); border_interpolate(winw, winh, 8, rawData, red, green, blue);
if (plistener) { if (plistener) {
@ -2809,7 +2530,7 @@ BENCHFUN
free(buffer0); free(buffer0);
} }
border_interpolate2(W, H, 1, rawData, red, green, blue); border_interpolate(W, H, 1, rawData, red, green, blue);
if(plistener) { if(plistener) {
plistener->setProgress (1.0); plistener->setProgress (1.0);
} }

40
rtengine/fast_demo.cc
View File

@ -30,6 +30,13 @@
using namespace std; using namespace std;
using namespace rtengine; using namespace rtengine;
namespace
{
unsigned fc(const unsigned int cfa[2][2], int r, int c) {
return cfa[r & 1][c & 1];
}
}
#define TS 224 #define TS 224
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -66,6 +73,7 @@ void RawImageSource::fast_demosaic()
} }
const unsigned int cfarray[2][2] = {{FC(0,0), FC(0,1)}, {FC(1,0), FC(1,1)}};
const int bord = 5; const int bord = 5;
float clip_pt = 4 * 65535 * initialGain; float clip_pt = 4 * 65535 * initialGain;
@ -116,12 +124,12 @@ void RawImageSource::fast_demosaic()
for (int i1 = imin; i1 < imax; i1++) for (int i1 = imin; i1 < imax; i1++)
for (int j1 = jmin; j1 < j + 2; j1++) { for (int j1 = jmin; j1 < j + 2; j1++) {
int c = FC(i1, j1); int c = fc(cfarray, i1, j1);
sum[c] += rawData[i1][j1]; sum[c] += rawData[i1][j1];
sum[c + 3]++; sum[c + 3]++;
} }
int c = FC(i, j); int c = fc(cfarray, i, j);
if (c == 1) { if (c == 1) {
red[i][j] = sum[0] / sum[3]; red[i][j] = sum[0] / sum[3];
@ -149,12 +157,12 @@ void RawImageSource::fast_demosaic()
for (int i1 = imin; i1 < imax; i1++) for (int i1 = imin; i1 < imax; i1++)
for (int j1 = j - 1; j1 < jmax; j1++) { for (int j1 = j - 1; j1 < jmax; j1++) {
int c = FC(i1, j1); int c = fc(cfarray, i1, j1);
sum[c] += rawData[i1][j1]; sum[c] += rawData[i1][j1];
sum[c + 3]++; sum[c + 3]++;
} }
int c = FC(i, j); int c = fc(cfarray, i, j);
if (c == 1) { if (c == 1) {
red[i][j] = sum[0] / sum[3]; red[i][j] = sum[0] / sum[3];
@ -192,12 +200,12 @@ void RawImageSource::fast_demosaic()
for (int i1 = max(0, i - 1); i1 < i + 2; i1++) for (int i1 = max(0, i - 1); i1 < i + 2; i1++)
for (int j1 = j - 1; j1 < j + 2; j1++) { for (int j1 = j - 1; j1 < j + 2; j1++) {
int c = FC(i1, j1); int c = fc(cfarray, i1, j1);
sum[c] += rawData[i1][j1]; sum[c] += rawData[i1][j1];
sum[c + 3]++; sum[c + 3]++;
} }
int c = FC(i, j); int c = fc(cfarray, i, j);
if (c == 1) { if (c == 1) {
red[i][j] = sum[0] / sum[3]; red[i][j] = sum[0] / sum[3];
@ -223,12 +231,12 @@ void RawImageSource::fast_demosaic()
for (int i1 = i - 1; i1 < min(i + 2, H); i1++) for (int i1 = i - 1; i1 < min(i + 2, H); i1++)
for (int j1 = j - 1; j1 < j + 2; j1++) { for (int j1 = j - 1; j1 < j + 2; j1++) {
int c = FC(i1, j1); int c = fc(cfarray, i1, j1);
sum[c] += rawData[i1][j1]; sum[c] += rawData[i1][j1];
sum[c + 3]++; sum[c + 3]++;
} }
int c = FC(i, j); int c = fc(cfarray, i, j);
if (c == 1) { if (c == 1) {
red[i][j] = sum[0] / sum[3]; red[i][j] = sum[0] / sum[3];
@ -281,7 +289,7 @@ void RawImageSource::fast_demosaic()
vmask selmask; vmask selmask;
vmask andmask = _mm_set_epi32( 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff ); vmask andmask = _mm_set_epi32( 0xffffffff, 0xffffffff, 0xffffffff, 0xffffffff );
if(FC(top, left) == 1) { if(fc(cfarray, top, left) == 1) {
selmask = _mm_set_epi32( 0, 0xffffffff, 0, 0xffffffff ); selmask = _mm_set_epi32( 0, 0xffffffff, 0, 0xffffffff );
} else { } else {
selmask = _mm_set_epi32( 0xffffffff, 0, 0xffffffff, 0 ); selmask = _mm_set_epi32( 0xffffffff, 0, 0xffffffff, 0 );
@ -311,7 +319,7 @@ void RawImageSource::fast_demosaic()
for (; j < right; j++, cc++) { for (; j < right; j++, cc++) {
if (FC(i, j) == 1) { if (fc(cfarray, i, j) == 1) {
greentile[rr * TS + cc] = rawData[i][j]; greentile[rr * TS + cc] = rawData[i][j];
} else { } else {
@ -332,7 +340,7 @@ void RawImageSource::fast_demosaic()
#else #else
for (int j = left, cc = 0; j < right; j++, cc++) { for (int j = left, cc = 0; j < right; j++, cc++) {
if (FC(i, j) == 1) { if (fc(cfarray, i, j) == 1) {
greentile[rr * TS + cc] = rawData[i][j]; greentile[rr * TS + cc] = rawData[i][j];
} else { } else {
//compute directional weights using image gradients //compute directional weights using image gradients
@ -358,7 +366,7 @@ void RawImageSource::fast_demosaic()
#endif #endif
for (int i = top + 1, rr = 1; i < bottom - 1; i++, rr++) { for (int i = top + 1, rr = 1; i < bottom - 1; i++, rr++) {
if (FC(i, left + (FC(i, 2) & 1) + 1) == 0) if (fc(cfarray, i, left + (fc(cfarray, i, 2) & 1) + 1) == 0)
#ifdef __SSE2__ #ifdef __SSE2__
for (int j = left + 1, cc = 1; j < right - 1; j += 4, cc += 4) { for (int j = left + 1, cc = 1; j < right - 1; j += 4, cc += 4) {
//interpolate B/R colors at R/B sites //interpolate B/R colors at R/B sites
@ -368,7 +376,7 @@ void RawImageSource::fast_demosaic()
#else #else
for (int cc = (FC(i, 2) & 1) + 1, j = left + cc; j < right - 1; j += 2, cc += 2) { for (int cc = (fc(cfarray, i, 2) & 1) + 1, j = left + cc; j < right - 1; j += 2, cc += 2) {
//interpolate B/R colors at R/B sites //interpolate B/R colors at R/B sites
bluetile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + (cc - 1)] + greentile[(rr - 1) * TS + (cc + 1)] + greentile[(rr + 1) * TS + cc + 1] + greentile[(rr + 1) * TS + cc - 1]) - bluetile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + (cc - 1)] + greentile[(rr - 1) * TS + (cc + 1)] + greentile[(rr + 1) * TS + cc + 1] + greentile[(rr + 1) * TS + cc - 1]) -
min(clip_pt, rawData[i - 1][j - 1] + rawData[i - 1][j + 1] + rawData[i + 1][j + 1] + rawData[i + 1][j - 1])); min(clip_pt, rawData[i - 1][j - 1] + rawData[i - 1][j + 1] + rawData[i + 1][j + 1] + rawData[i + 1][j - 1]));
@ -385,7 +393,7 @@ void RawImageSource::fast_demosaic()
#else #else
for (int cc = (FC(i, 2) & 1) + 1, j = left + cc; j < right - 1; j += 2, cc += 2) { for (int cc = (fc(cfarray, i, 2) & 1) + 1, j = left + cc; j < right - 1; j += 2, cc += 2) {
//interpolate B/R colors at R/B sites //interpolate B/R colors at R/B sites
redtile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + cc - 1] + greentile[(rr - 1) * TS + cc + 1] + greentile[(rr + 1) * TS + cc + 1] + greentile[(rr + 1) * TS + cc - 1]) - redtile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + cc - 1] + greentile[(rr - 1) * TS + cc + 1] + greentile[(rr + 1) * TS + cc + 1] + greentile[(rr + 1) * TS + cc - 1]) -
min(clip_pt, rawData[i - 1][j - 1] + rawData[i - 1][j + 1] + rawData[i + 1][j + 1] + rawData[i + 1][j - 1])); min(clip_pt, rawData[i - 1][j - 1] + rawData[i - 1][j + 1] + rawData[i + 1][j + 1] + rawData[i + 1][j - 1]));
@ -404,7 +412,7 @@ void RawImageSource::fast_demosaic()
for (int i = top + 2, rr = 2; i < bottom - 2; i++, rr++) { for (int i = top + 2, rr = 2; i < bottom - 2; i++, rr++) {
#ifdef __SSE2__ #ifdef __SSE2__
for (int cc = 2 + (FC(i, 2) & 1), j = left + cc; j < right - 2; j += 4, cc += 4) { for (int cc = 2 + (fc(cfarray, i, 2) & 1), j = left + cc; j < right - 2; j += 4, cc += 4) {
// no need to take care about the borders of the tile. There's enough free space. // no need to take care about the borders of the tile. There's enough free space.
//interpolate R and B colors at G sites //interpolate R and B colors at G sites
greenv = LVFU(greentile[rr * TS + cc]); greenv = LVFU(greentile[rr * TS + cc]);
@ -428,7 +436,7 @@ void RawImageSource::fast_demosaic()
#else #else
for (int cc = 2 + (FC(i, 2) & 1), j = left + cc; j < right - 2; j += 2, cc += 2) { for (int cc = 2 + (fc(cfarray, i, 2) & 1), j = left + cc; j < right - 2; j += 2, cc += 2) {
//interpolate R and B colors at G sites //interpolate R and B colors at G sites
redtile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + cc] - redtile[(rr - 1) * TS + cc]) + (greentile[(rr + 1) * TS + cc] - redtile[(rr + 1) * TS + cc]) + redtile[rr * TS + cc] = greentile[rr * TS + cc] - 0.25f * ((greentile[(rr - 1) * TS + cc] - redtile[(rr - 1) * TS + cc]) + (greentile[(rr + 1) * TS + cc] - redtile[(rr + 1) * TS + cc]) +
(greentile[rr * TS + cc - 1] - redtile[rr * TS + cc - 1]) + (greentile[rr * TS + cc + 1] - redtile[rr * TS + cc + 1])); (greentile[rr * TS + cc - 1] - redtile[rr * TS + cc - 1]) + (greentile[rr * TS + cc + 1] - redtile[rr * TS + cc + 1]));

1
rtengine/filmnegativeproc.cc
View File

@ -25,6 +25,7 @@
#include "coord.h" #include "coord.h"
#include "mytime.h" #include "mytime.h"
#include "opthelper.h" #include "opthelper.h"
#include "pixelsmap.h"
#include "procparams.h" #include "procparams.h"
#include "rt_algo.h" #include "rt_algo.h"
#include "rtengine.h" #include "rtengine.h"

2
rtengine/hphd_demosaic_RT.cc
View File

@ -352,7 +352,7 @@ void RawImageSource::hphd_demosaic ()
interpolate_row_rb_mul_pp(rawData, red[i], blue[i], green[i - 1], green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1); interpolate_row_rb_mul_pp(rawData, red[i], blue[i], green[i - 1], green[i], green[i + 1], i, 1.0, 1.0, 1.0, 0, W, 1);
} }
border_interpolate2(W, H, 4, rawData, red, green, blue); border_interpolate(W, H, 4, rawData, red, green, blue);
if (plistener) { if (plistener) {
plistener->setProgress(1.0); plistener->setProgress(1.0);

1
rtengine/improccoordinator.cc
View File

@ -38,7 +38,6 @@
#include "refreshmap.h" #include "refreshmap.h"
#include "../rtgui/options.h" #include "../rtgui/options.h"
#include "../rtgui/ppversion.h"
#ifdef _OPENMP #ifdef _OPENMP
#include <omp.h> #include <omp.h>

1
rtengine/improcfun.cc
View File

@ -45,7 +45,6 @@
#include "satandvalueblendingcurve.h" #include "satandvalueblendingcurve.h"
#include "StopWatch.h" #include "StopWatch.h"
#include "procparams.h" #include "procparams.h"
#include "../rtgui/ppversion.h"
#include "../rtgui/editcallbacks.h" #include "../rtgui/editcallbacks.h"
#ifdef _DEBUG #ifdef _DEBUG

1
rtengine/pdaflinesfilter.cc
View File

@ -22,6 +22,7 @@
#include "camconst.h" #include "camconst.h"
#include "pdaflinesfilter.h" #include "pdaflinesfilter.h"
#include "pixelsmap.h"
#include "rawimage.h" #include "rawimage.h"
#include "settings.h" #include "settings.h"

23
rtengine/pixelshift.cc
View File

@ -38,6 +38,10 @@
namespace namespace
{ {
unsigned fc(const unsigned int cfa[2][2], int r, int c) {
return cfa[r & 1][c & 1];
}
float greenDiff(float a, float b, float stddevFactor, float eperIso, float nreadIso, float prnu) float greenDiff(float a, float b, float stddevFactor, float eperIso, float nreadIso, float prnu)
{ {
// calculate the difference between two green samples // calculate the difference between two green samples
@ -323,6 +327,7 @@ BENCHFUN
bayerParams.pixelShiftShowMotion = false; bayerParams.pixelShiftShowMotion = false;
} }
const unsigned int cfarray[2][2] = {{FC(0,0), FC(0,1)}, {FC(1,0), FC(1,1)}};
const bool showMotion = bayerParams.pixelShiftShowMotion; const bool showMotion = bayerParams.pixelShiftShowMotion;
const bool showOnlyMask = bayerParams.pixelShiftShowMotionMaskOnly && showMotion; const bool showOnlyMask = bayerParams.pixelShiftShowMotionMaskOnly && showMotion;
const float smoothFactor = 1.0 - bayerParams.pixelShiftSmoothFactor; const float smoothFactor = 1.0 - bayerParams.pixelShiftSmoothFactor;
@ -641,11 +646,11 @@ BENCHFUN
for(int i = winy + 1; i < winh - 1; ++i) { for(int i = winy + 1; i < winh - 1; ++i) {
int j = winx + 1; int j = winx + 1;
int c = FC(i, j); int c = fc(cfarray, i, j);
bool bluerow = (c + FC(i, j + 1)) == 3; bool bluerow = (c + fc(cfarray, i, j + 1)) == 3;
for(int j = winx + 1, offset = FC(i, j) & 1; j < winw - 1; ++j, offset ^= 1) { for(int j = winx + 1, offset = fc(cfarray, i, j) & 1; j < winw - 1; ++j, offset ^= 1) {
(*histogreenThr[1 - offset])[(*rawDataFrames[1 - offset])[i - offset + 1][j]]++; (*histogreenThr[1 - offset])[(*rawDataFrames[1 - offset])[i - offset + 1][j]]++;
(*histogreenThr[3 - offset])[(*rawDataFrames[3 - offset])[i + offset][j + 1]]++; (*histogreenThr[3 - offset])[(*rawDataFrames[3 - offset])[i + offset][j + 1]]++;
@ -726,9 +731,9 @@ BENCHFUN
}; };
int ng = 0; int ng = 0;
int j = winx + 1; int j = winx + 1;
int c = FC(i, j); int c = fc(cfarray, i, j);
if((c + FC(i, j + 1)) == 3) { if((c + fc(cfarray, i, j + 1)) == 3) {
// row with blue pixels => swap destination pointers for non green pixels // row with blue pixels => swap destination pointers for non green pixels
std::swap(nonGreenDest0, nonGreenDest1); std::swap(nonGreenDest0, nonGreenDest1);
ng ^= 1; ng ^= 1;
@ -783,7 +788,7 @@ BENCHFUN
for(int i = winy + border - offsY; i < winh - (border + offsY); ++i) { for(int i = winy + border - offsY; i < winh - (border + offsY); ++i) {
// offset to keep the code short. It changes its value between 0 and 1 for each iteration of the loop // offset to keep the code short. It changes its value between 0 and 1 for each iteration of the loop
unsigned int offset = FC(i, winx + border - offsX) & 1; unsigned int offset = fc(cfarray, i, winx + border - offsX) & 1;
for(int j = winx + border - offsX; j < winw - (border + offsX); ++j, offset ^= 1) { for(int j = winx + border - offsX; j < winw - (border + offsX); ++j, offset ^= 1) {
psMask[i][j] = noMotion; psMask[i][j] = noMotion;
@ -919,7 +924,7 @@ BENCHFUN
float *blueDest = blue[i + offsY]; float *blueDest = blue[i + offsY];
// offset to keep the code short. It changes its value between 0 and 1 for each iteration of the loop // offset to keep the code short. It changes its value between 0 and 1 for each iteration of the loop
unsigned int offset = FC(i, winx + border - offsX) & 1; unsigned int offset = fc(cfarray, i, winx + border - offsX) & 1;
for(int j = winx + border - offsX; j < winw - (border + offsX); ++j, offset ^= 1) { for(int j = winx + border - offsX; j < winw - (border + offsX); ++j, offset ^= 1) {
if(showOnlyMask) { if(showOnlyMask) {
@ -969,9 +974,9 @@ BENCHFUN
float *nonGreenDest1 = blue[i]; float *nonGreenDest1 = blue[i];
int ng = 0; int ng = 0;
int j = winx + 1; int j = winx + 1;
int c = FC(i, j); int c = fc(cfarray, i, j);
if((c + FC(i, j + 1)) == 3) { if((c + fc(cfarray, i, j + 1)) == 3) {
// row with blue pixels => swap destination pointers for non green pixels // row with blue pixels => swap destination pointers for non green pixels
std::swap(nonGreenDest0, nonGreenDest1); std::swap(nonGreenDest0, nonGreenDest1);
ng ^= 1; ng ^= 1;

7
rtengine/rawimagesource.h
View File

@ -26,12 +26,12 @@
#include "colortemp.h" #include "colortemp.h"
#include "iimage.h" #include "iimage.h"
#include "imagesource.h" #include "imagesource.h"
#include "pixelsmap.h"
#define HR_SCALE 2 #define HR_SCALE 2
namespace rtengine namespace rtengine
{ {
class PixelsMap;
class RawImage; class RawImage;
class DiagonalCurve; class DiagonalCurve;
class RetinextransmissionCurve; class RetinextransmissionCurve;
@ -266,8 +266,6 @@ protected:
void eahd_demosaic(); void eahd_demosaic();
void hphd_demosaic(); void hphd_demosaic();
void vng4_demosaic(const array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue); void vng4_demosaic(const array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue);
void ppg_demosaic();
void jdl_interpolate_omp();
void igv_interpolate(int winw, int winh); void igv_interpolate(int winw, int winh);
void lmmse_interpolate_omp(int winw, int winh, array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue, int iterations); void lmmse_interpolate_omp(int winw, int winh, array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue, int iterations);
void amaze_demosaic_RT(int winx, int winy, int winw, int winh, const array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue, size_t chunkSize = 1, bool measure = false);//Emil's code for AMaZE void amaze_demosaic_RT(int winx, int winy, int winw, int winh, const array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue, size_t chunkSize = 1, bool measure = false);//Emil's code for AMaZE
@ -276,8 +274,7 @@ protected:
void dcb_demosaic(int iterations, bool dcb_enhance); void dcb_demosaic(int iterations, bool dcb_enhance);
void ahd_demosaic(); void ahd_demosaic();
void rcd_demosaic(size_t chunkSize = 1, bool measure = false); void rcd_demosaic(size_t chunkSize = 1, bool measure = false);
void border_interpolate(unsigned int border, float (*image)[4], unsigned int start = 0, unsigned int end = 0); void border_interpolate(int winw, int winh, int lborders, const array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue);
void border_interpolate2(int winw, int winh, int lborders, const array2D<float> &rawData, array2D<float> &red, array2D<float> &green, array2D<float> &blue);
void dcb_initTileLimits(int &colMin, int &rowMin, int &colMax, int &rowMax, int x0, int y0, int border); void dcb_initTileLimits(int &colMin, int &rowMin, int &colMax, int &rowMax, int x0, int y0, int border);
void fill_raw( float (*cache )[3], int x0, int y0, float** rawData); void fill_raw( float (*cache )[3], int x0, int y0, float** rawData);
void fill_border( float (*cache )[3], int border, int x0, int y0); void fill_border( float (*cache )[3], int border, int x0, int y0);

24
rtengine/rcd_demosaic.cc
View File

@ -26,6 +26,13 @@
using namespace std; using namespace std;
namespace
{
unsigned fc(const unsigned int cfa[2][2], int r, int c) {
return cfa[r & 1][c & 1];
}
}
namespace rtengine namespace rtengine
{ {
@ -55,6 +62,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
plistener->setProgress(progress); plistener->setProgress(progress);
} }
const unsigned int cfarray[2][2] = {{FC(0,0), FC(0,1)}, {FC(1,0), FC(1,1)}};
constexpr int rcdBorder = 9; constexpr int rcdBorder = 9;
constexpr int tileSize = 214; constexpr int tileSize = 214;
constexpr int tileSizeN = tileSize - 2 * rcdBorder; constexpr int tileSizeN = tileSize - 2 * rcdBorder;
@ -97,8 +105,8 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
for (int row = rowStart; row < rowEnd; row++) { for (int row = rowStart; row < rowEnd; row++) {
int indx = (row - rowStart) * tileSize; int indx = (row - rowStart) * tileSize;
int c0 = FC(row, colStart); int c0 = fc(cfarray, row, colStart);
int c1 = FC(row, colStart + 1); int c1 = fc(cfarray, row, colStart + 1);
int col = colStart; int col = colStart;
for (; col < colEnd - 1; col+=2, indx+=2) { for (; col < colEnd - 1; col+=2, indx+=2) {
@ -131,7 +139,7 @@ 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 // 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 row = 2; row < tileRows - 2; row++) {
for (int col = 2 + (FC(row, 0) & 1), indx = row * tileSize + col; col < tilecols - 2; col += 2, indx += 2) { 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]); 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]);
} }
} }
@ -141,7 +149,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
*/ */
// 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 < tileRows - 4; row++) { for (int row = 4; row < tileRows - 4; row++) {
for (int col = 4 + (FC(row, 0) & 1), indx = row * tileSize + col; col < tilecols - 4; col += 2, indx += 2) { for (int col = 4 + (fc(cfarray, row, 0) & 1), indx = row * tileSize + col; col < tilecols - 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];
@ -176,7 +184,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
// Step 4.1: Calculate P/Q diagonal local discrimination // Step 4.1: Calculate P/Q diagonal local discrimination
for (int row = rcdBorder - 4; row < tileRows - rcdBorder + 4; row++) { for (int row = rcdBorder - 4; row < tileRows - rcdBorder + 4; row++) {
for (int col = rcdBorder - 4 + (FC(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; col < tilecols - rcdBorder + 4; col += 2, indx += 2) {
const float cfai = cfa[indx]; 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 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]);
@ -189,7 +197,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
// 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 = rcdBorder - 3; row < tileRows - rcdBorder + 3; row++) { for (int row = rcdBorder - 3; row < tileRows - rcdBorder + 3; row++) {
for (int col = rcdBorder - 3 + (FC(row, rcdBorder - 1) & 1), indx = row * tileSize + col, c = 2 - FC(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); col < tilecols - rcdBorder + 3; 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];
@ -221,7 +229,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
// 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 = rcdBorder; row < tileRows - rcdBorder; row++) { for (int row = rcdBorder; row < tileRows - rcdBorder; row++) {
for (int col = rcdBorder + (FC(row, rcdBorder - 1) & 1), indx = row * tileSize + col; col < tilecols - rcdBorder; col += 2, indx += 2) { 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 // Refined vertical and horizontal local discrimination
float VH_Central_Value = VH_Dir[indx]; float VH_Central_Value = VH_Dir[indx];
@ -296,7 +304,7 @@ void RawImageSource::rcd_demosaic(size_t chunkSize, bool measure)
free(PQ_Dir); free(PQ_Dir);
} }
border_interpolate2(W, H, rcdBorder, rawData, red, green, blue); border_interpolate(W, H, rcdBorder, rawData, red, green, blue);
if (plistener) { if (plistener) {
plistener->setProgress(1); plistener->setProgress(1);

5
rtengine/rtthumbnail.cc
View File

@ -31,15 +31,12 @@
#include "colortemp.h" #include "colortemp.h"
#include "curves.h" #include "curves.h"
#include "dcp.h" #include "dcp.h"
#include "iccmatrices.h"
#include "iccstore.h" #include "iccstore.h"
#include "image8.h" #include "image8.h"
#include "improccoordinator.h"
#include "improcfun.h" #include "improcfun.h"
#include "jpeg.h" #include "jpeg.h"
#include "labimage.h" #include "labimage.h"
#include "median.h" #include "median.h"
#include "mytime.h"
#include "procparams.h" #include "procparams.h"
#include "rawimage.h" #include "rawimage.h"
#include "rawimagesource.h" #include "rawimagesource.h"
@ -50,8 +47,6 @@
#include "StopWatch.h" #include "StopWatch.h"
#include "utils.h" #include "utils.h"
#include "../rtgui/ppversion.h"
namespace namespace
{ {

2
rtengine/vng4_demosaic_RT.cc
View File

@ -384,7 +384,7 @@ void RawImageSource::vng4_demosaic (const array2D<float> &rawData, array2D<float
#endif #endif
{ {
// let the first thread, which is out of work, do the border interpolation // let the first thread, which is out of work, do the border interpolation
border_interpolate2(W, H, 3, rawData, red, green, blue); border_interpolate(W, H, 3, rawData, red, green, blue);
} }
} }