diff --git a/rtengine/amaze_demosaic_RT.cc b/rtengine/amaze_demosaic_RT.cc index 3b367ee2b..2720521aa 100644 --- a/rtengine/amaze_demosaic_RT.cc +++ b/rtengine/amaze_demosaic_RT.cc @@ -9,6 +9,7 @@ // incorporating ideas of Luis Sanz Rodrigues and Paul Lee // // code dated: May 27, 2010 +// latest modification: Ingo Weyrich, January 25, 2016 // // amaze_interpolate_RT.cc is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by @@ -52,9 +53,15 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, const float clip_pt = 1.0 / initialGain; const float clip_pt8 = 0.8 / initialGain; - -#define TS 160 // Tile size; the image is processed in square tiles to lower memory requirements and facilitate multi-threading -#define TSH 80 // half of Tile size +// this allows to pass AMAZETS to the code. On some machines larger AMAZETS is faster +// If AMAZETS is undefined it will be set to 160, which is the fastest on modern x86/64 machines +#ifndef AMAZETS +#define AMAZETS 160 +#endif + // Tile size; the image is processed in square tiles to lower memory requirements and facilitate multi-threading + // We assure that Tile size is a multiple of 32 in the range [96;992] + constexpr int ts = (AMAZETS & 992) < 96 ? 96 : (AMAZETS & 992); + constexpr int tsh = ts / 2; // half of Tile size //offset of R pixel within a Bayer quartet int ex, ey; @@ -79,27 +86,27 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, } //shifts of pointer value to access pixels in vertical and diagonal directions - static const int v1 = TS, v2 = 2 * TS, v3 = 3 * TS, p1 = -TS + 1, p2 = -2 * TS + 2, p3 = -3 * TS + 3, m1 = TS + 1, m2 = 2 * TS + 2, m3 = 3 * TS + 3; + constexpr int v1 = ts, v2 = 2 * ts, v3 = 3 * ts, p1 = -ts + 1, p2 = -2 * ts + 2, p3 = -3 * ts + 3, m1 = ts + 1, m2 = 2 * ts + 2, m3 = 3 * ts + 3; //tolerance to avoid dividing by zero - static const float eps = 1e-5, epssq = 1e-10; //tolerance to avoid dividing by zero + constexpr float eps = 1e-5, epssq = 1e-10; //tolerance to avoid dividing by zero //adaptive ratios threshold - static const float arthresh = 0.75; + constexpr float arthresh = 0.75; //gaussian on 5x5 quincunx, sigma=1.2 - static const float gaussodd[4] = {0.14659727707323927f, 0.103592713382435f, 0.0732036125103057f, 0.0365543548389495f}; + constexpr float gaussodd[4] = {0.14659727707323927f, 0.103592713382435f, 0.0732036125103057f, 0.0365543548389495f}; //nyquist texture test threshold - static const float nyqthresh = 0.5; + constexpr float nyqthresh = 0.5; //gaussian on 5x5, sigma=1.2, multiplied with nyqthresh to save some time later in loop // Is this really sigma=1.2????, seems more like sigma = 1.672 - static const float gaussgrad[6] = {nyqthresh * 0.07384411893421103f, nyqthresh * 0.06207511968171489f, nyqthresh * 0.0521818194747806f, - nyqthresh * 0.03687419286733595f, nyqthresh * 0.03099732204057846f, nyqthresh * 0.018413194161458882f - }; + constexpr float gaussgrad[6] = {nyqthresh * 0.07384411893421103f, nyqthresh * 0.06207511968171489f, nyqthresh * 0.0521818194747806f, + nyqthresh * 0.03687419286733595f, nyqthresh * 0.03099732204057846f, nyqthresh * 0.018413194161458882f + }; //gaussian on 5x5 alt quincunx, sigma=1.5 - static const float gausseven[2] = {0.13719494435797422f, 0.05640252782101291f}; + constexpr float gausseven[2] = {0.13719494435797422f, 0.05640252782101291f}; //guassian on quincunx grid - static const float gquinc[4] = {0.169917f, 0.108947f, 0.069855f, 0.0287182f}; + constexpr float gquinc[4] = {0.169917f, 0.108947f, 0.069855f, 0.0287182f}; typedef struct { float h; @@ -112,115 +119,87 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, { int progresscounter = 0; -#define CLF 1 + constexpr int cldf = 2; // factor to multiply cache line distance. 1 = 64 bytes, 2 = 128 bytes ... // assign working space - char *buffer = (char *) calloc(13 * sizeof(float) * TS * TS + sizeof(float) * TS * TSH + sizeof(char) * TS * TSH + 18 * CLF * 64 + 63, 1); + char *buffer = (char *) calloc(14 * sizeof(float) * ts * ts + sizeof(char) * ts * tsh + 18 * cldf * 64 + 63, 1); // aligned to 64 byte boundary char *data = (char*)( ( uintptr_t(buffer) + uintptr_t(63)) / 64 * 64); // green values float *rgbgreen = (float (*)) data; // sum of square of horizontal gradient and square of vertical gradient - float *delhvsqsum = (float (*)) ((char*)rgbgreen + sizeof(float) * TS * TS + CLF * 64); + float *delhvsqsum = (float (*)) ((char*)rgbgreen + sizeof(float) * ts * ts + cldf * 64); // 1 // gradient based directional weights for interpolation - float *dirwts0 = (float (*)) ((char*)delhvsqsum + sizeof(float) * TS * TS + CLF * 64); - float *dirwts1 = (float (*)) ((char*)dirwts0 + sizeof(float) * TS * TS + CLF * 64); - // vertically interpolated color differences G-R, G-B - float *vcd = (float (*)) ((char*)dirwts1 + sizeof(float) * TS * TS + CLF * 64); - // horizontally interpolated color differences - float *hcd = (float (*)) ((char*)vcd + sizeof(float) * TS * TS + CLF * 64); + float *dirwts0 = (float (*)) ((char*)delhvsqsum + sizeof(float) * ts * ts + cldf * 64); // 1 + float *dirwts1 = (float (*)) ((char*)dirwts0 + sizeof(float) * ts * ts + cldf * 64); // 1 + // vertically interpolated colour differences G-R, G-B + float *vcd = (float (*)) ((char*)dirwts1 + sizeof(float) * ts * ts + cldf * 64); // 1 + // horizontally interpolated colour differences + float *hcd = (float (*)) ((char*)vcd + sizeof(float) * ts * ts + cldf * 64); // 1 // alternative vertical interpolation - float *vcdalt = (float (*)) ((char*)hcd + sizeof(float) * TS * TS + CLF * 64); + float *vcdalt = (float (*)) ((char*)hcd + sizeof(float) * ts * ts + cldf * 64); // 1 // alternative horizontal interpolation - float *hcdalt = (float (*)) ((char*)vcdalt + sizeof(float) * TS * TS + CLF * 64); - // square of average color difference - float *cddiffsq = (float (*)) ((char*)hcdalt + sizeof(float) * TS * TS + CLF * 64); + float *hcdalt = (float (*)) ((char*)vcdalt + sizeof(float) * ts * ts + cldf * 64); // 1 + // square of average colour difference + float *cddiffsq = (float (*)) ((char*)hcdalt + sizeof(float) * ts * ts + cldf * 64); // 1 // weight to give horizontal vs vertical interpolation - float *hvwt = (float (*)) ((char*)cddiffsq + sizeof(float) * TS * TS + 2 * CLF * 64); - // final interpolated color difference - float (*Dgrb)[TS * TSH] = (float (*)[TS * TSH])vcdalt; // there is no overlap in buffer usage => share + float *hvwt = (float (*)) ((char*)cddiffsq + sizeof(float) * ts * ts + 2 * cldf * 64); // 1 + // final interpolated colour difference + float (*Dgrb)[ts * tsh] = (float (*)[ts * tsh])vcdalt; // there is no overlap in buffer usage => share // gradient in plus (NE/SW) direction float *delp = (float (*))cddiffsq; // there is no overlap in buffer usage => share // gradient in minus (NW/SE) direction - float *delm = (float (*)) ((char*)delp + sizeof(float) * TS * TSH + CLF * 64); + float *delm = (float (*)) ((char*)delp + sizeof(float) * ts * tsh + cldf * 64); // diagonal interpolation of R+B float *rbint = (float (*))delm; // there is no overlap in buffer usage => share // horizontal and vertical curvature of interpolated G (used to refine interpolation in Nyquist texture regions) - s_hv *Dgrb2 = (s_hv (*)) ((char*)hvwt + sizeof(float) * TS * TSH + CLF * 64); + s_hv *Dgrb2 = (s_hv (*)) ((char*)hvwt + sizeof(float) * ts * tsh + cldf * 64); // 1 // difference between up/down interpolations of G float *dgintv = (float (*))Dgrb2; // there is no overlap in buffer usage => share // difference between left/right interpolations of G - float *dginth = (float (*)) ((char*)dgintv + sizeof(float) * TS * TS + CLF * 64); + float *dginth = (float (*)) ((char*)dgintv + sizeof(float) * ts * ts + cldf * 64); // 1 // square of diagonal colour differences - float *Dgrbsq1m = (float (*)) ((char*)dginth + sizeof(float) * TS * TS + CLF * 64); - float *Dgrbsq1p = (float (*)) ((char*)Dgrbsq1m + sizeof(float) * TS * TSH + CLF * 64); + float *Dgrbsq1m = (float (*)) ((char*)dginth + sizeof(float) * ts * ts + cldf * 64); // 1 + float *Dgrbsq1p = (float (*)) ((char*)Dgrbsq1m + sizeof(float) * ts * tsh + cldf * 64); // 1 // tile raw data - float *cfa = (float (*)) ((char*)Dgrbsq1p + sizeof(float) * TS * TSH + CLF * 64); + float *cfa = (float (*)) ((char*)Dgrbsq1p + sizeof(float) * ts * tsh + cldf * 64); // 1 // relative weight for combining plus and minus diagonal interpolations float *pmwt = (float (*))delhvsqsum; // there is no overlap in buffer usage => share - // interpolated color difference R-B in minus and plus direction + // interpolated colour difference R-B in minus and plus direction float *rbm = (float (*))vcd; // there is no overlap in buffer usage => share - float *rbp = (float (*)) ((char*)rbm + sizeof(float) * TS * TSH + CLF * 64); - // nyquist texture flag 1=nyquist, 0=not nyquist - unsigned char *nyquist = (unsigned char (*)) ((char*)cfa + sizeof(float) * TS * TS + CLF * 64); - /* - rgbgreen = (float (*)) data; //pointers to array - delhvsqsum = (float (*)) ((char*)rgbgreen + sizeof(float) * TS * TS + CLF * 64); - dirwts0 = (float (*)) ((char*)delhvsqsum + sizeof(float) * TS * TS + CLF * 64); - dirwts1 = (float (*)) ((char*)dirwts0 + sizeof(float) * TS * TS + CLF * 64); - vcd = (float (*)) ((char*)dirwts1 + sizeof(float) * TS * TS + CLF * 64); - hcd = (float (*)) ((char*)vcd + sizeof(float) * TS * TS + CLF * 64); - vcdalt = (float (*)) ((char*)hcd + sizeof(float) * TS * TS + CLF * 64); - hcdalt = (float (*)) ((char*)vcdalt + sizeof(float) * TS * TS + CLF * 64); - cddiffsq = (float (*)) ((char*)hcdalt + sizeof(float) * TS * TS + CLF * 64); - hvwt = (float (*)) ((char*)cddiffsq + sizeof(float) * TS * TS + CLF * 64); - Dgrb = (float (*)[TS * TSH]) ((char*)hvwt + sizeof(float) * TS * TSH + CLF * 64); - delp = (float (*)) ((char*)Dgrb + sizeof(float) * TS * TS + CLF * 64); - delm = (float (*)) ((char*)delp + sizeof(float) * TS * TSH + CLF * 64); - rbint = (float (*)) ((char*)delm + sizeof(float) * TS * TSH + CLF * 64); - Dgrb2 = (s_hv (*)) ((char*)rbint + sizeof(float) * TS * TSH + CLF * 64); - dgintv = (float (*)) ((char*)Dgrb2 + sizeof(float) * TS * TS + CLF * 64); - dginth = (float (*)) ((char*)dgintv + sizeof(float) * TS * TS + CLF * 64); - Dgrbsq1m = (float (*)) ((char*)dginth + sizeof(float) * TS * TS + CLF * 64); - Dgrbsq1p = (float (*)) ((char*)Dgrbsq1m + sizeof(float) * TS * TSH + CLF * 64); - cfa = (float (*)) ((char*)Dgrbsq1p + sizeof(float) * TS * TSH + CLF * 64); - pmwt = (float (*)) ((char*)cfa + sizeof(float) * TS * TS + CLF * 64); - rbm = (float (*)) ((char*)pmwt + sizeof(float) * TS * TSH + CLF * 64); - rbp = (float (*)) ((char*)rbm + sizeof(float) * TS * TSH + CLF * 64); - - nyquist = (char (*)) ((char*)rbp + sizeof(float) * TS * TSH + CLF * 64); - */ -#undef CLF + float *rbp = (float (*)) ((char*)rbm + sizeof(float) * ts * tsh + cldf * 64); + // nyquist texture flags 1=nyquist, 0=not nyquist + unsigned char *nyquist = (unsigned char (*)) ((char*)cfa + sizeof(float) * ts * ts + cldf * 64); // 1 + unsigned char *nyquist2 = (unsigned char (*))cddiffsq; + float *nyqutest = (float(*)) ((char*)nyquist + sizeof(unsigned char) * ts * tsh + cldf * 64); // 1 // Main algorithm: Tile loop - - // Issue 1676 // use collapse(2) to collapse the 2 loops to one large loop, so there is better scaling #ifdef _OPENMP #pragma omp for schedule(dynamic) collapse(2) nowait #endif - for (int top = winy - 16; top < winy + height; top += TS - 32) - for (int left = winx - 16; left < winx + width; left += TS - 32) { + for (int top = winy - 16; top < winy + height; top += ts - 32) { + for (int left = winx - 16; left < winx + width; left += ts - 32) { #ifdef __SSE2__ // Using SSE2 we can zero the memory without cache pollution vfloat zerov = ZEROV; - for(int i = 3 * TSH; i < (TS - 6)*TSH; i += 16) { + for(int i = 3 * tsh; i < (ts - 6)*tsh; i += 16) { _mm_stream_ps((float*)&nyquist[i], zerov); } #else - memset(&nyquist[3 * TSH], 0, sizeof(unsigned char) * (TS - 6) * TSH); + memset(&nyquist[3 * tsh], 0, sizeof(unsigned char) * (ts - 6) * tsh); #endif //location of tile bottom edge - const int bottom = min(top + TS, winy + height + 16); + int bottom = min(top + ts, winy + height + 16); //location of tile right edge - const int right = min(left + TS, winx + width + 16); - //tile width (=TS except for right edge of image) - const int rr1 = bottom - top; - //tile height (=TS except for bottom edge of image) - const int cc1 = right - left; + int right = min(left + ts, winx + width + 16); + //tile width (=ts except for right edge of image) + int rr1 = bottom - top; + //tile height (=ts except for bottom edge of image) + int cc1 = right - left; // bookkeeping for borders // min and max row/column in the tile int rrmin = top < winy ? 16 : 0; @@ -232,54 +211,87 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, // rgb values should be floating point number between 0 and 1 // after white balance multipliers are applied // a 16 pixel border is added to each side of the image + + // begin of tile initialization #ifdef __SSE2__ - const vfloat c65535v = F2V( 65535.0f ); + vfloat c65535v = F2V( 65535.f ); //fill upper border if (rrmin > 0) { - for (int rr = 0; rr < 16; rr++) - for (int cc = ccmin, row = 32 - rr + top; cc < ccmax; cc++) { - cfa[rr * TS + cc] = (rawData[row][cc + left]) / 65535.0f; - rgbgreen[rr * TS + cc] = cfa[rr * TS + cc]; + for (int rr = 0; rr < 16; rr++) { + int row = 32 - rr + top; + + for (int cc = ccmin; cc < ccmax; cc += 4) { + int indx1 = rr * ts + cc; + vfloat tempv = LVFU(rawData[row][cc + left]) / c65535v; + STVF(cfa[indx1], tempv); + STVF(rgbgreen[indx1], tempv ); } + } } // fill inner part for (int rr = rrmin; rr < rrmax; rr++) { int row = rr + top; - int cc = ccmin; - for (; cc < ccmax - 3; cc += 4) { - int indx1 = rr * TS + cc; + for (int cc = ccmin; cc < ccmax; cc += 4) { + int indx1 = rr * ts + cc; vfloat tempv = LVFU(rawData[row][cc + left]) / c65535v; STVF(cfa[indx1], tempv ); STVF(rgbgreen[indx1], tempv ); } - - for (; cc < ccmax; cc++) { - int indx1 = rr * TS + cc; - cfa[indx1] = (rawData[row][cc + left]) / 65535.0f; - rgbgreen[indx1] = cfa[indx1]; - } } //fill lower border if (rrmax < rr1) { for (int rr = 0; rr < 16; rr++) for (int cc = ccmin; cc < ccmax; cc += 4) { - int indx1 = (rrmax + rr) * TS + cc; + int indx1 = (rrmax + rr) * ts + cc; vfloat tempv = LVFU(rawData[(winy + height - rr - 2)][left + cc]) / c65535v; STVF(cfa[indx1], tempv ); STVF(rgbgreen[indx1], tempv ); } } +#else + + //fill upper border + if (rrmin > 0) { + for (int rr = 0; rr < 16; rr++) + for (int cc = ccmin, row = 32 - rr + top; cc < ccmax; cc++) { + cfa[rr * ts + cc] = (rawData[row][cc + left]) / 65535.f; + rgbgreen[rr * ts + cc] = cfa[rr * ts + cc]; + } + } + + // fill inner part + for (int rr = rrmin; rr < rrmax; rr++) { + int row = rr + top; + + for (int cc = ccmin; cc < ccmax; cc++) { + int indx1 = rr * ts + cc; + cfa[indx1] = (rawData[row][cc + left]) / 65535.f; + rgbgreen[indx1] = cfa[indx1]; + } + } + + //fill lower border + if (rrmax < rr1) { + for (int rr = 0; rr < 16; rr++) + for (int cc = ccmin; cc < ccmax; cc++) { + cfa[(rrmax + rr)*ts + cc] = (rawData[(winy + height - rr - 2)][left + cc]) / 65535.f; + rgbgreen[(rrmax + rr)*ts + cc] = cfa[(rrmax + rr) * ts + cc]; + } + } + +#endif + //fill left border if (ccmin > 0) { for (int rr = rrmin; rr < rrmax; rr++) for (int cc = 0, row = rr + top; cc < 16; cc++) { - cfa[rr * TS + cc] = (rawData[row][32 - cc + left]) / 65535.0f; - rgbgreen[rr * TS + cc] = cfa[rr * TS + cc]; + cfa[rr * ts + cc] = (rawData[row][32 - cc + left]) / 65535.f; + rgbgreen[rr * ts + cc] = cfa[rr * ts + cc]; } } @@ -287,87 +299,8 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, if (ccmax < cc1) { for (int rr = rrmin; rr < rrmax; rr++) for (int cc = 0; cc < 16; cc++) { - cfa[rr * TS + ccmax + cc] = (rawData[(top + rr)][(winx + width - cc - 2)]) / 65535.0f; - rgbgreen[rr * TS + ccmax + cc] = cfa[rr * TS + ccmax + cc]; - } - } - - //also, fill the image corners - if (rrmin > 0 && ccmin > 0) { - for (int rr = 0; rr < 16; rr++) - for (int cc = 0; cc < 16; cc += 4) { - int indx1 = (rr) * TS + cc; - vfloat tempv = LVFU(rawData[winy + 32 - rr][winx + 32 - cc]) / c65535v; - STVF(cfa[indx1], tempv ); - STVF(rgbgreen[indx1], tempv ); - } - } - - if (rrmax < rr1 && ccmax < cc1) { - for (int rr = 0; rr < 16; rr++) - for (int cc = 0; cc < 16; cc += 4) { - int indx1 = (rrmax + rr) * TS + ccmax + cc; - vfloat tempv = LVFU(rawData[(winy + height - rr - 2)][(winx + width - cc - 2)]) / c65535v; - STVFU(cfa[indx1], tempv ); - STVFU(rgbgreen[indx1], tempv ); - } - } - - if (rrmin > 0 && ccmax < cc1) { - for (int rr = 0; rr < 16; rr++) - for (int cc = 0; cc < 16; cc++) { - cfa[(rr)*TS + ccmax + cc] = (rawData[(winy + 32 - rr)][(winx + width - cc - 2)]) / 65535.0f; - rgbgreen[(rr)*TS + ccmax + cc] = cfa[(rr) * TS + ccmax + cc]; - } - } - - if (rrmax < rr1 && ccmin > 0) { - for (int rr = 0; rr < 16; rr++) - for (int cc = 0; cc < 16; cc++) { - cfa[(rrmax + rr)*TS + cc] = (rawData[(winy + height - rr - 2)][(winx + 32 - cc)]) / 65535.0f; - rgbgreen[(rrmax + rr)*TS + cc] = cfa[(rrmax + rr) * TS + cc]; - } - } - -#else - - for (int rr = rrmin; rr < rrmax; rr++) - for (int row = rr + top, cc = ccmin; cc < ccmax; cc++) { - int indx1 = rr * TS + cc; - cfa[indx1] = (rawData[row][cc + left]) / 65535.0f; - rgbgreen[indx1] = cfa[indx1]; - } - - //fill borders - if (rrmin > 0) { - for (int rr = 0; rr < 16; rr++) - for (int cc = ccmin, row = 32 - rr + top; cc < ccmax; cc++) { - cfa[rr * TS + cc] = (rawData[row][cc + left]) / 65535.0f; - rgbgreen[rr * TS + cc] = cfa[rr * TS + cc]; - } - } - - if (rrmax < rr1) { - for (int rr = 0; rr < 16; rr++) - for (int cc = ccmin; cc < ccmax; cc++) { - cfa[(rrmax + rr)*TS + cc] = (rawData[(winy + height - rr - 2)][left + cc]) / 65535.0f; - rgbgreen[(rrmax + rr)*TS + cc] = cfa[(rrmax + rr) * TS + cc]; - } - } - - if (ccmin > 0) { - for (int rr = rrmin; rr < rrmax; rr++) - for (int cc = 0, row = rr + top; cc < 16; cc++) { - cfa[rr * TS + cc] = (rawData[row][32 - cc + left]) / 65535.0f; - rgbgreen[rr * TS + cc] = cfa[rr * TS + cc]; - } - } - - if (ccmax < cc1) { - for (int rr = rrmin; rr < rrmax; rr++) - for (int cc = 0; cc < 16; cc++) { - cfa[rr * TS + ccmax + cc] = (rawData[(top + rr)][(winx + width - cc - 2)]) / 65535.0f; - rgbgreen[rr * TS + ccmax + cc] = cfa[rr * TS + ccmax + cc]; + cfa[rr * ts + ccmax + cc] = (rawData[(top + rr)][(winx + width - cc - 2)]) / 65535.f; + rgbgreen[rr * ts + ccmax + cc] = cfa[rr * ts + ccmax + cc]; } } @@ -375,43 +308,43 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, if (rrmin > 0 && ccmin > 0) { for (int rr = 0; rr < 16; rr++) for (int cc = 0; cc < 16; cc++) { - cfa[(rr)*TS + cc] = (rawData[winy + 32 - rr][winx + 32 - cc]) / 65535.0f; - rgbgreen[(rr)*TS + cc] = cfa[(rr) * TS + cc]; + cfa[(rr)*ts + cc] = (rawData[winy + 32 - rr][winx + 32 - cc]) / 65535.f; + rgbgreen[(rr)*ts + cc] = cfa[(rr) * ts + cc]; } } if (rrmax < rr1 && ccmax < cc1) { for (int rr = 0; rr < 16; rr++) for (int cc = 0; cc < 16; cc++) { - cfa[(rrmax + rr)*TS + ccmax + cc] = (rawData[(winy + height - rr - 2)][(winx + width - cc - 2)]) / 65535.0f; - rgbgreen[(rrmax + rr)*TS + ccmax + cc] = cfa[(rrmax + rr) * TS + ccmax + cc]; + cfa[(rrmax + rr)*ts + ccmax + cc] = (rawData[(winy + height - rr - 2)][(winx + width - cc - 2)]) / 65535.f; + rgbgreen[(rrmax + rr)*ts + ccmax + cc] = cfa[(rrmax + rr) * ts + ccmax + cc]; } } if (rrmin > 0 && ccmax < cc1) { for (int rr = 0; rr < 16; rr++) for (int cc = 0; cc < 16; cc++) { - cfa[(rr)*TS + ccmax + cc] = (rawData[(winy + 32 - rr)][(winx + width - cc - 2)]) / 65535.0f; - rgbgreen[(rr)*TS + ccmax + cc] = cfa[(rr) * TS + ccmax + cc]; + cfa[(rr)*ts + ccmax + cc] = (rawData[(winy + 32 - rr)][(winx + width - cc - 2)]) / 65535.f; + rgbgreen[(rr)*ts + ccmax + cc] = cfa[(rr) * ts + ccmax + cc]; } } if (rrmax < rr1 && ccmin > 0) { for (int rr = 0; rr < 16; rr++) for (int cc = 0; cc < 16; cc++) { - cfa[(rrmax + rr)*TS + cc] = (rawData[(winy + height - rr - 2)][(winx + 32 - cc)]) / 65535.0f; - rgbgreen[(rrmax + rr)*TS + cc] = cfa[(rrmax + rr) * TS + cc]; + cfa[(rrmax + rr)*ts + cc] = (rawData[(winy + height - rr - 2)][(winx + 32 - cc)]) / 65535.f; + rgbgreen[(rrmax + rr)*ts + cc] = cfa[(rrmax + rr) * ts + cc]; } } -#endif + // end of tile initialization - //end of border fill + // horizontal and vertical gradients #ifdef __SSE2__ - const vfloat epsv = F2V( eps ); + vfloat epsv = F2V( eps ); for (int rr = 2; rr < rr1 - 2; rr++) { - for (int indx = rr * TS; indx < rr * TS + cc1; indx += 4) { + for (int indx = rr * ts; indx < rr * ts + cc1; indx += 4) { vfloat delhv = vabsf( LVFU( cfa[indx + 1] ) - LVFU( cfa[indx - 1] ) ); vfloat delvv = vabsf( LVF( cfa[indx + v1] ) - LVF( cfa[indx - v1] ) ); STVF(dirwts1[indx], epsv + vabsf( LVFU( cfa[indx + 2] ) - LVF( cfa[indx] )) + vabsf( LVF( cfa[indx] ) - LVFU( cfa[indx - 2] )) + delhv ); @@ -423,53 +356,55 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #else for (int rr = 2; rr < rr1 - 2; rr++) - for (int cc = 2, indx = (rr) * TS + cc; cc < cc1 - 2; cc++, indx++) { - // horizontal and vedrtical gradient + for (int cc = 2, indx = (rr) * ts + cc; cc < cc1 - 2; cc++, indx++) { float delh = fabsf(cfa[indx + 1] - cfa[indx - 1]); float delv = fabsf(cfa[indx + v1] - cfa[indx - v1]); dirwts0[indx] = eps + fabsf(cfa[indx + v2] - cfa[indx]) + fabsf(cfa[indx] - cfa[indx - v2]) + delv; - dirwts1[indx] = eps + fabsf(cfa[indx + 2] - cfa[indx]) + fabsf(cfa[indx] - cfa[indx - 2]) + delh; //+fabsf(cfa[indx+2]-cfa[indx-2]); + dirwts1[indx] = eps + fabsf(cfa[indx + 2] - cfa[indx]) + fabsf(cfa[indx] - cfa[indx - 2]) + delh; delhvsqsum[indx] = SQR(delh) + SQR(delv); } #endif - //interpolate vertical and horizontal color differences + //interpolate vertical and horizontal colour differences #ifdef __SSE2__ vfloat sgnv; if( !(FC(4, 4) & 1) ) { - sgnv = _mm_set_ps( 1.0f, -1.0f, 1.0f, -1.0f ); + sgnv = _mm_set_ps( 1.f, -1.f, 1.f, -1.f ); } else { - sgnv = _mm_set_ps( -1.0f, 1.0f, -1.0f, 1.0f ); + sgnv = _mm_set_ps( -1.f, 1.f, -1.f, 1.f ); } - vfloat zd5v = F2V( 0.5f ); - vfloat onev = F2V( 1.0f ); - vfloat arthreshv = F2V( arthresh ); - vfloat clip_pt8v = F2V( clip_pt8 ); + vfloat zd5v = F2V( 0.5f ); + vfloat onev = F2V( 1.f ); + vfloat arthreshv = F2V( arthresh ); + vfloat clip_pt8v = F2V( clip_pt8 ); for (int rr = 4; rr < rr1 - 4; rr++) { sgnv = -sgnv; - for (int indx = rr * TS + 4; indx < rr * TS + cc1 - 7; indx += 4) { - //color ratios in each cardinal direction + for (int indx = rr * ts + 4; indx < rr * ts + cc1 - 7; indx += 4) { + //colour ratios in each cardinal direction vfloat cfav = LVF(cfa[indx]); vfloat cruv = LVF(cfa[indx - v1]) * (LVF(dirwts0[indx - v2]) + LVF(dirwts0[indx])) / (LVF(dirwts0[indx - v2]) * (epsv + cfav) + LVF(dirwts0[indx]) * (epsv + LVF(cfa[indx - v2]))); vfloat crdv = LVF(cfa[indx + v1]) * (LVF(dirwts0[indx + v2]) + LVF(dirwts0[indx])) / (LVF(dirwts0[indx + v2]) * (epsv + cfav) + LVF(dirwts0[indx]) * (epsv + LVF(cfa[indx + v2]))); vfloat crlv = LVFU(cfa[indx - 1]) * (LVFU(dirwts1[indx - 2]) + LVF(dirwts1[indx])) / (LVFU(dirwts1[indx - 2]) * (epsv + cfav) + LVF(dirwts1[indx]) * (epsv + LVFU(cfa[indx - 2]))); vfloat crrv = LVFU(cfa[indx + 1]) * (LVFU(dirwts1[indx + 2]) + LVF(dirwts1[indx])) / (LVFU(dirwts1[indx + 2]) * (epsv + cfav) + LVF(dirwts1[indx]) * (epsv + LVFU(cfa[indx + 2]))); + //G interpolated in vert/hor directions using Hamilton-Adams method vfloat guhav = LVF(cfa[indx - v1]) + zd5v * (cfav - LVF(cfa[indx - v2])); vfloat gdhav = LVF(cfa[indx + v1]) + zd5v * (cfav - LVF(cfa[indx + v2])); vfloat glhav = LVFU(cfa[indx - 1]) + zd5v * (cfav - LVFU(cfa[indx - 2])); vfloat grhav = LVFU(cfa[indx + 1]) + zd5v * (cfav - LVFU(cfa[indx + 2])); + //G interpolated in vert/hor directions using adaptive ratios vfloat guarv = vself(vmaskf_lt(vabsf(onev - cruv), arthreshv), cfav * cruv, guhav); vfloat gdarv = vself(vmaskf_lt(vabsf(onev - crdv), arthreshv), cfav * crdv, gdhav); vfloat glarv = vself(vmaskf_lt(vabsf(onev - crlv), arthreshv), cfav * crlv, glhav); vfloat grarv = vself(vmaskf_lt(vabsf(onev - crrv), arthreshv), cfav * crrv, grhav); + //adaptive weights for vertical/horizontal directions vfloat hwtv = LVFU(dirwts1[indx - 1]) / (LVFU(dirwts1[indx - 1]) + LVFU(dirwts1[indx + 1])); vfloat vwtv = LVF(dirwts0[indx - v1]) / (LVF(dirwts0[indx + v1]) + LVF(dirwts0[indx - v1])); @@ -477,7 +412,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, vfloat Ginthhav = vintpf(hwtv, grhav, glhav); vfloat Gintvhav = vintpf(vwtv, gdhav, guhav); - //interpolated color differences + //interpolated colour differences vfloat hcdaltv = sgnv * (Ginthhav - cfav); vfloat vcdaltv = sgnv * (Gintvhav - cfav); STVF(hcdalt[indx], hcdaltv); @@ -488,13 +423,14 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, gdarv = vself( clipmask, gdhav, gdarv); glarv = vself( clipmask, glhav, glarv); grarv = vself( clipmask, grhav, grarv); + + //use HA if highlights are (nearly) clipped STVF(vcd[indx], vself( clipmask, vcdaltv, sgnv * (vintpf(vwtv, gdarv, guarv) - cfav))); STVF(hcd[indx], vself( clipmask, hcdaltv, sgnv * (vintpf(hwtv, grarv, glarv) - cfav))); - //differences of interpolations in opposite directions + //differences of interpolations in opposite directions STVF(dgintv[indx], vminf(SQRV(guhav - gdhav), SQRV(guarv - gdarv))); STVF(dginth[indx], vminf(SQRV(glhav - grhav), SQRV(glarv - grarv))); - } } @@ -503,9 +439,9 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, for (int rr = 4; rr < rr1 - 4; rr++) { bool fcswitch = FC(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++) { - //color ratios in each cardinal direction + //colour ratios in each cardinal direction float cru = cfa[indx - v1] * (dirwts0[indx - v2] + dirwts0[indx]) / (dirwts0[indx - v2] * (eps + cfa[indx]) + dirwts0[indx] * (eps + cfa[indx - v2])); float crd = cfa[indx + v1] * (dirwts0[indx + v2] + dirwts0[indx]) / (dirwts0[indx + v2] * (eps + cfa[indx]) + dirwts0[indx] * (eps + cfa[indx + v2])); float crl = cfa[indx - 1] * (dirwts1[indx - 2] + dirwts1[indx]) / (dirwts1[indx - 2] * (eps + cfa[indx]) + dirwts1[indx] * (eps + cfa[indx - 2])); @@ -520,25 +456,25 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, //G interpolated in vert/hor directions using adaptive ratios float guar, gdar, glar, grar; - if (fabsf(1.0f - cru) < arthresh) { + if (fabsf(1.f - cru) < arthresh) { guar = cfa[indx] * cru; } else { guar = guha; } - if (fabsf(1.0f - crd) < arthresh) { + if (fabsf(1.f - crd) < arthresh) { gdar = cfa[indx] * crd; } else { gdar = gdha; } - if (fabsf(1.0f - crl) < arthresh) { + if (fabsf(1.f - crl) < arthresh) { glar = cfa[indx] * crl; } else { glar = glha; } - if (fabsf(1.0f - crr) < arthresh) { + if (fabsf(1.f - crr) < arthresh) { grar = cfa[indx] * crr; } else { grar = grha; @@ -549,19 +485,19 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, float vwt = dirwts0[indx - v1] / (dirwts0[indx + v1] + dirwts0[indx - v1]); //interpolated G via adaptive weights of cardinal evaluations - float Gintvha = vwt * gdha + (1.0f - vwt) * guha; - float Ginthha = hwt * grha + (1.0f - hwt) * glha; + float Gintvha = vwt * gdha + (1.f - vwt) * guha; + float Ginthha = hwt * grha + (1.f - hwt) * glha; - //interpolated color differences + //interpolated colour differences if (fcswitch) { - vcd[indx] = cfa[indx] - (vwt * gdar + (1.0f - vwt) * guar); - hcd[indx] = cfa[indx] - (hwt * grar + (1.0f - hwt) * glar); + vcd[indx] = cfa[indx] - (vwt * gdar + (1.f - vwt) * guar); + hcd[indx] = cfa[indx] - (hwt * grar + (1.f - hwt) * glar); vcdalt[indx] = cfa[indx] - Gintvha; hcdalt[indx] = cfa[indx] - Ginthha; } else { - //interpolated color differences - vcd[indx] = (vwt * gdar + (1.0f - vwt) * guar) - cfa[indx]; - hcd[indx] = (hwt * grar + (1.0f - hwt) * glar) - cfa[indx]; + //interpolated colour differences + vcd[indx] = (vwt * gdar + (1.f - vwt) * guar) - cfa[indx]; + hcd[indx] = (hwt * grar + (1.f - hwt) * glar) - cfa[indx]; vcdalt[indx] = Gintvha - cfa[indx]; hcdalt[indx] = Ginthha - cfa[indx]; } @@ -583,8 +519,6 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, dginth[indx] = min(SQR(glha - grha), SQR(glar - grar)); } - - } #endif @@ -596,9 +530,9 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, vfloat sgn3v; if( !(FC(4, 4) & 1) ) { - sgnv = _mm_set_ps( 1.0f, -1.0f, 1.0f, -1.0f ); + sgnv = _mm_set_ps( 1.f, -1.f, 1.f, -1.f ); } else { - sgnv = _mm_set_ps( -1.0f, 1.0f, -1.0f, 1.0f ); + sgnv = _mm_set_ps( -1.f, 1.f, -1.f, 1.f ); } sgn3v = sgnv + sgnv + sgnv; @@ -608,7 +542,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, sgnv = -sgnv; sgn3v = -sgn3v; - for (int indx = rr * TS + 4; indx < rr * TS + cc1 - 4; indx += 4) { + for (int indx = rr * ts + 4; indx < rr * ts + cc1 - 4; indx += 4) { vfloat hcdv = LVF( hcd[indx] ); vfloat hcdvarv = SQRV(LVFU(hcd[indx - 2]) - hcdv) + SQRV(LVFU(hcd[indx - 2]) - LVFU(hcd[indx + 2])) + SQRV(hcdv - LVFU(hcd[indx + 2])); vfloat hcdaltv = LVF( hcdalt[indx] ); @@ -622,6 +556,8 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, hcdv = vself( vmaskf_lt( hcdaltvarv, hcdvarv ), hcdaltv, hcdv); vcdv = vself( vmaskf_lt( vcdaltvarv, vcdvarv ), vcdaltv, vcdv); + //bound the interpolation in regions of high saturation + //vertical and horizontal G interpolations vfloat Ginthv = sgnv * hcdv + LVF( cfa[indx] ); vfloat temp2v = sgn3v * hcdv; vfloat hwtv = onev + temp2v / ( epsv + Ginthv + LVF( cfa[indx])); @@ -651,11 +587,11 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #else 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++) { - float hcdvar = 3.0f * (SQR(hcd[indx - 2]) + SQR(hcd[indx]) + SQR(hcd[indx + 2])) - SQR(hcd[indx - 2] + hcd[indx] + hcd[indx + 2]); - float hcdaltvar = 3.0f * (SQR(hcdalt[indx - 2]) + SQR(hcdalt[indx]) + SQR(hcdalt[indx + 2])) - SQR(hcdalt[indx - 2] + hcdalt[indx] + hcdalt[indx + 2]); - float vcdvar = 3.0f * (SQR(vcd[indx - v2]) + SQR(vcd[indx]) + SQR(vcd[indx + v2])) - SQR(vcd[indx - v2] + vcd[indx] + vcd[indx + v2]); - float vcdaltvar = 3.0f * (SQR(vcdalt[indx - v2]) + SQR(vcdalt[indx]) + SQR(vcdalt[indx + v2])) - SQR(vcdalt[indx - v2] + vcdalt[indx] + vcdalt[indx + v2]); + for (int cc = 4, indx = rr * ts + cc, c = FC(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 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 vcdaltvar = 3.f * (SQR(vcdalt[indx - v2]) + SQR(vcdalt[indx]) + SQR(vcdalt[indx + v2])) - SQR(vcdalt[indx - v2] + vcdalt[indx] + vcdalt[indx + v2]); //choose the smallest variance; this yields a smoother interpolation if (hcdaltvar < hcdvar) { @@ -667,7 +603,6 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, } //bound the interpolation in regions of high saturation - //vertical and horizontal G interpolations float Gintv, Ginth; @@ -676,33 +611,31 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, Gintv = -vcd[indx] + cfa[indx]; //B or R if (hcd[indx] > 0) { - if (3.0f * hcd[indx] > (Ginth + cfa[indx])) { + if (3.f * hcd[indx] > (Ginth + cfa[indx])) { hcd[indx] = -ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) + cfa[indx]; } else { - float hwt = 1.0f - 3.0f * hcd[indx] / (eps + Ginth + cfa[indx]); - hcd[indx] = hwt * hcd[indx] + (1.0f - hwt) * (-ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) + cfa[indx]); + float hwt = 1.f - 3.f * hcd[indx] / (eps + Ginth + cfa[indx]); + hcd[indx] = hwt * hcd[indx] + (1.f - hwt) * (-ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) + cfa[indx]); } } if (vcd[indx] > 0) { - if (3.0f * vcd[indx] > (Gintv + cfa[indx])) { + if (3.f * vcd[indx] > (Gintv + cfa[indx])) { vcd[indx] = -ULIM(Gintv, cfa[indx - v1], cfa[indx + v1]) + cfa[indx]; } else { - float vwt = 1.0f - 3.0f * vcd[indx] / (eps + Gintv + cfa[indx]); - vcd[indx] = vwt * vcd[indx] + (1.0f - vwt) * (-ULIM(Gintv, cfa[indx - v1], cfa[indx + v1]) + cfa[indx]); + float vwt = 1.f - 3.f * vcd[indx] / (eps + Gintv + cfa[indx]); + vcd[indx] = vwt * vcd[indx] + (1.f - vwt) * (-ULIM(Gintv, cfa[indx - v1], cfa[indx + v1]) + cfa[indx]); } } if (Ginth > clip_pt) { - hcd[indx] = -ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) + cfa[indx]; //for RT implementation + hcd[indx] = -ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) + cfa[indx]; } if (Gintv > clip_pt) { vcd[indx] = -ULIM(Gintv, cfa[indx - v1], cfa[indx + v1]) + cfa[indx]; } - //if (Ginth > pre_mul[c]) hcd[indx]=-ULIM(Ginth,cfa[indx-1],cfa[indx+1])+cfa[indx];//for dcraw implementation - //if (Gintv > pre_mul[c]) vcd[indx]=-ULIM(Gintv,cfa[indx-v1],cfa[indx+v1])+cfa[indx]; } else {//R or B site @@ -710,25 +643,25 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, Gintv = vcd[indx] + cfa[indx]; if (hcd[indx] < 0) { - if (3.0f * hcd[indx] < -(Ginth + cfa[indx])) { + if (3.f * hcd[indx] < -(Ginth + cfa[indx])) { hcd[indx] = ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) - cfa[indx]; } else { - float hwt = 1.0f + 3.0f * hcd[indx] / (eps + Ginth + cfa[indx]); - hcd[indx] = hwt * hcd[indx] + (1.0f - hwt) * (ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) - cfa[indx]); + float hwt = 1.f + 3.f * hcd[indx] / (eps + Ginth + cfa[indx]); + hcd[indx] = hwt * hcd[indx] + (1.f - hwt) * (ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) - cfa[indx]); } } if (vcd[indx] < 0) { - if (3.0f * vcd[indx] < -(Gintv + cfa[indx])) { + if (3.f * vcd[indx] < -(Gintv + cfa[indx])) { vcd[indx] = ULIM(Gintv, cfa[indx - v1], cfa[indx + v1]) - cfa[indx]; } else { - float vwt = 1.0f + 3.0f * vcd[indx] / (eps + Gintv + cfa[indx]); - vcd[indx] = vwt * vcd[indx] + (1.0f - vwt) * (ULIM(Gintv, cfa[indx - v1], cfa[indx + v1]) - cfa[indx]); + float vwt = 1.f + 3.f * vcd[indx] / (eps + Gintv + cfa[indx]); + vcd[indx] = vwt * vcd[indx] + (1.f - vwt) * (ULIM(Gintv, cfa[indx - v1], cfa[indx + v1]) - cfa[indx]); } } if (Ginth > clip_pt) { - hcd[indx] = ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) - cfa[indx]; //for RT implementation + hcd[indx] = ULIM(Ginth, cfa[indx - 1], cfa[indx + 1]) - cfa[indx]; } if (Gintv > clip_pt) { @@ -750,7 +683,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, vfloat epssqv = F2V( epssq ); 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(rr, 2) & 1); indx < rr * ts + cc1 - 6; indx += 8) { //compute colour difference variances in cardinal directions vfloat tempv = LC2VFU(vcd[indx]); vfloat uavev = tempv + LC2VFU(vcd[indx - v1]) + LC2VFU(vcd[indx - v2]) + LC2VFU(vcd[indx - v3]); @@ -758,8 +691,8 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, vfloat Dgrbvvaruv = SQRV(tempv - uavev) + SQRV(LC2VFU(vcd[indx - v1]) - uavev) + SQRV(LC2VFU(vcd[indx - v2]) - uavev) + SQRV(LC2VFU(vcd[indx - v3]) - uavev); vfloat Dgrbvvardv = SQRV(tempv - davev) + SQRV(LC2VFU(vcd[indx + v1]) - davev) + SQRV(LC2VFU(vcd[indx + v2]) - davev) + SQRV(LC2VFU(vcd[indx + v3]) - davev); - vfloat hwtv = LC2VFU(dirwts1[indx - 1]) / (LC2VFU(dirwts1[indx - 1]) + LC2VFU(dirwts1[indx + 1])); - vfloat vwtv = LC2VFU(dirwts0[indx - v1]) / (LC2VFU(dirwts0[indx + v1]) + LC2VFU(dirwts0[indx - v1])); + vfloat hwtv = vadivapb(LC2VFU(dirwts1[indx - 1]), LC2VFU(dirwts1[indx + 1])); + vfloat vwtv = vadivapb(LC2VFU(dirwts0[indx - v1]), LC2VFU(dirwts0[indx + v1])); tempv = LC2VFU(hcd[indx]); vfloat lavev = tempv + vaddc2vfu(hcd[indx - 3]) + LC2VFU(hcd[indx - 1]); @@ -772,7 +705,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, vfloat vcdvarv = epssqv + vintpf(vwtv, Dgrbvvardv, Dgrbvvaruv); vfloat hcdvarv = epssqv + vintpf(hwtv, Dgrbhvarrv, Dgrbhvarlv); - //compute fluctuations in up/down and left/right interpolations of colors + //compute fluctuations in up/down and left/right interpolations of colours Dgrbvvaruv = LC2VFU(dgintv[indx - v1]) + LC2VFU(dgintv[indx - v2]); Dgrbvvardv = LC2VFU(dgintv[indx + v1]) + LC2VFU(dgintv[indx + v2]); @@ -796,16 +729,16 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #else 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(rr, 2) & 1), indx = rr * ts + cc; cc < cc1 - 6; cc += 2, indx += 2) { - //compute color difference variances in cardinal directions + //compute colour difference variances in cardinal directions float uave = vcd[indx] + vcd[indx - v1] + vcd[indx - v2] + vcd[indx - v3]; float dave = vcd[indx] + vcd[indx + v1] + vcd[indx + v2] + vcd[indx + v3]; float lave = hcd[indx] + hcd[indx - 1] + hcd[indx - 2] + hcd[indx - 3]; float rave = hcd[indx] + hcd[indx + 1] + hcd[indx + 2] + hcd[indx + 3]; - //color difference (G-R or G-B) variance in up/down/left/right directions + //colour difference (G-R or G-B) variance in up/down/left/right directions float Dgrbvvaru = SQR(vcd[indx] - uave) + SQR(vcd[indx - v1] - uave) + SQR(vcd[indx - v2] - uave) + SQR(vcd[indx - v3] - uave); float Dgrbvvard = SQR(vcd[indx] - dave) + SQR(vcd[indx + v1] - dave) + SQR(vcd[indx + v2] - dave) + SQR(vcd[indx + v3] - dave); float Dgrbhvarl = SQR(hcd[indx] - lave) + SQR(hcd[indx - 1] - lave) + SQR(hcd[indx - 2] - lave) + SQR(hcd[indx - 3] - lave); @@ -814,17 +747,17 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, float hwt = dirwts1[indx - 1] / (dirwts1[indx - 1] + dirwts1[indx + 1]); float vwt = dirwts0[indx - v1] / (dirwts0[indx + v1] + dirwts0[indx - v1]); - float vcdvar = epssq + vwt * Dgrbvvard + (1.0f - vwt) * Dgrbvvaru; - float hcdvar = epssq + hwt * Dgrbhvarr + (1.0f - hwt) * Dgrbhvarl; + float vcdvar = epssq + vwt * Dgrbvvard + (1.f - vwt) * Dgrbvvaru; + float hcdvar = epssq + hwt * Dgrbhvarr + (1.f - hwt) * Dgrbhvarl; - //compute fluctuations in up/down and left/right interpolations of colors + //compute fluctuations in up/down and left/right interpolations of colours Dgrbvvaru = (dgintv[indx]) + (dgintv[indx - v1]) + (dgintv[indx - v2]); Dgrbvvard = (dgintv[indx]) + (dgintv[indx + v1]) + (dgintv[indx + v2]); Dgrbhvarl = (dginth[indx]) + (dginth[indx - 1]) + (dginth[indx - 2]); Dgrbhvarr = (dginth[indx]) + (dginth[indx + 1]) + (dginth[indx + 2]); - float vcdvar1 = epssq + vwt * Dgrbvvard + (1.0f - vwt) * Dgrbvvaru; - float hcdvar1 = epssq + hwt * Dgrbhvarr + (1.0f - hwt) * Dgrbhvarl; + float vcdvar1 = epssq + vwt * Dgrbvvard + (1.f - vwt) * Dgrbvvaru; + float hcdvar1 = epssq + hwt * Dgrbhvarr + (1.f - hwt) * Dgrbhvarl; //determine adaptive weights for G interpolation float varwt = hcdvar / (vcdvar + hcdvar); @@ -837,52 +770,94 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, } else { hvwt[indx >> 1] = diffwt; } - } } #endif +#ifdef __SSE2__ + vfloat gaussg0 = F2V(gaussgrad[0]); + vfloat gaussg1 = F2V(gaussgrad[1]); + vfloat gaussg2 = F2V(gaussgrad[2]); + vfloat gaussg3 = F2V(gaussgrad[3]); + vfloat gaussg4 = F2V(gaussgrad[4]); + vfloat gaussg5 = F2V(gaussgrad[5]); + vfloat gausso0 = F2V(gaussodd[0]); + vfloat gausso1 = F2V(gaussodd[1]); + vfloat gausso2 = F2V(gaussodd[2]); + vfloat gausso3 = F2V(gaussodd[3]); + +#endif + + // precompute nyquist + for (int rr = 6; rr < rr1 - 6; rr++) { + int cc = 6 + (FC(rr, 2) & 1); + int indx = rr * ts + cc; + +#ifdef __SSE2__ + + for (; cc < cc1 - 7; cc += 8, indx += 8) { + vfloat valv = (gausso0 * LC2VFU(cddiffsq[indx]) + + gausso1 * (LC2VFU(cddiffsq[(indx - m1)]) + LC2VFU(cddiffsq[(indx + p1)]) + + LC2VFU(cddiffsq[(indx - p1)]) + LC2VFU(cddiffsq[(indx + m1)])) + + gausso2 * (LC2VFU(cddiffsq[(indx - v2)]) + LC2VFU(cddiffsq[(indx - 2)]) + + LC2VFU(cddiffsq[(indx + 2)]) + LC2VFU(cddiffsq[(indx + v2)])) + + gausso3 * (LC2VFU(cddiffsq[(indx - m2)]) + LC2VFU(cddiffsq[(indx + p2)]) + + LC2VFU(cddiffsq[(indx - p2)]) + LC2VFU(cddiffsq[(indx + m2)]))) - + (gaussg0 * LC2VFU(delhvsqsum[indx]) + + gaussg1 * (LC2VFU(delhvsqsum[indx - v1]) + LC2VFU(delhvsqsum[indx - 1]) + + LC2VFU(delhvsqsum[indx + 1]) + LC2VFU(delhvsqsum[indx + v1])) + + gaussg2 * (LC2VFU(delhvsqsum[indx - m1]) + LC2VFU(delhvsqsum[indx + p1]) + + LC2VFU(delhvsqsum[indx - p1]) + LC2VFU(delhvsqsum[indx + m1])) + + gaussg3 * (LC2VFU(delhvsqsum[indx - v2]) + LC2VFU(delhvsqsum[indx - 2]) + + LC2VFU(delhvsqsum[indx + 2]) + LC2VFU(delhvsqsum[indx + v2])) + + gaussg4 * (LC2VFU(delhvsqsum[indx - v2 - 1]) + LC2VFU(delhvsqsum[indx - v2 + 1]) + + LC2VFU(delhvsqsum[indx - ts - 2]) + LC2VFU(delhvsqsum[indx - ts + 2]) + + LC2VFU(delhvsqsum[indx + ts - 2]) + LC2VFU(delhvsqsum[indx + ts + 2]) + + LC2VFU(delhvsqsum[indx + v2 - 1]) + LC2VFU(delhvsqsum[indx + v2 + 1])) + + gaussg5 * (LC2VFU(delhvsqsum[indx - m2]) + LC2VFU(delhvsqsum[indx + p2]) + + LC2VFU(delhvsqsum[indx - p2]) + LC2VFU(delhvsqsum[indx + m2]))); + STVFU(nyqutest[indx >> 1], valv); + + } + +#endif + + for (; cc < cc1 - 6; cc += 2, indx += 2) { + nyqutest[indx >> 1] = (gaussodd[0] * cddiffsq[indx] + + gaussodd[1] * (cddiffsq[(indx - m1)] + cddiffsq[(indx + p1)] + + cddiffsq[(indx - p1)] + cddiffsq[(indx + m1)]) + + gaussodd[2] * (cddiffsq[(indx - v2)] + cddiffsq[(indx - 2)] + + cddiffsq[(indx + 2)] + cddiffsq[(indx + v2)]) + + gaussodd[3] * (cddiffsq[(indx - m2)] + cddiffsq[(indx + p2)] + + cddiffsq[(indx - p2)] + cddiffsq[(indx + m2)])) - + (gaussgrad[0] * delhvsqsum[indx] + + gaussgrad[1] * (delhvsqsum[indx - v1] + delhvsqsum[indx + 1] + + delhvsqsum[indx - 1] + delhvsqsum[indx + v1]) + + gaussgrad[2] * (delhvsqsum[indx - m1] + delhvsqsum[indx + p1] + + delhvsqsum[indx - p1] + delhvsqsum[indx + m1]) + + gaussgrad[3] * (delhvsqsum[indx - v2] + delhvsqsum[indx - 2] + + delhvsqsum[indx + 2] + delhvsqsum[indx + v2]) + + gaussgrad[4] * (delhvsqsum[indx - v2 - 1] + delhvsqsum[indx - v2 + 1] + + delhvsqsum[indx - ts - 2] + delhvsqsum[indx - ts + 2] + + delhvsqsum[indx + ts - 2] + delhvsqsum[indx + ts + 2] + + delhvsqsum[indx + v2 - 1] + delhvsqsum[indx + v2 + 1]) + + gaussgrad[5] * (delhvsqsum[indx - m2] + delhvsqsum[indx + p2] + + delhvsqsum[indx - p2] + delhvsqsum[indx + m2])); + } + } // Nyquist test int nystartrow = 0; int nyendrow = 0; - int nystartcol = TS + 1; + int nystartcol = ts + 1; int nyendcol = 0; 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(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 - // TODO_INGO: currently this part needs 10 float mults, 36 float adds, 4 int mults and 44 int adds for every second pixel - // it reads 304 bytes for every second pixel and writes <= 1 byte for every second pixel - // a precalculated vectorized version could do this with 1/4 of the operations - // but it would read 304 bytes for every second pixel and write 8 bytes for every second pixel for the precalculation - // (though the vectorized read should be faster than the scalar version) - // and read 8 bytes for every second pixel and write 1 byte for every second pixel for final calculation (maybe this last step can be avoided too) - float nyqtest1 = gaussodd[0] * cddiffsq[indx] + - gaussodd[1] * (cddiffsq[(indx - m1)] + cddiffsq[(indx + p1)] + - cddiffsq[(indx - p1)] + cddiffsq[(indx + m1)]) + - gaussodd[2] * (cddiffsq[(indx - v2)] + cddiffsq[(indx - 2)] + - cddiffsq[(indx + 2)] + cddiffsq[(indx + v2)]) + - gaussodd[3] * (cddiffsq[(indx - m2)] + cddiffsq[(indx + p2)] + - cddiffsq[(indx - p2)] + cddiffsq[(indx + m2)]); - float nyqtest2 = gaussgrad[0] * delhvsqsum[indx] + - gaussgrad[1] * (delhvsqsum[indx - v1] + delhvsqsum[indx + 1] + - delhvsqsum[indx - 1] + delhvsqsum[indx + v1]) + - gaussgrad[2] * (delhvsqsum[indx - m1] + delhvsqsum[indx + p1] + - delhvsqsum[indx - p1] + delhvsqsum[indx + m1]) + - gaussgrad[3] * (delhvsqsum[indx - v2] + delhvsqsum[indx - 2] + - delhvsqsum[indx + 2] + delhvsqsum[indx + v2]) + - gaussgrad[4] * (delhvsqsum[indx - 2 * TS - 1] + delhvsqsum[indx - 2 * TS + 1] + - delhvsqsum[indx - TS - 2] + delhvsqsum[indx - TS + 2] + - delhvsqsum[indx + TS - 2] + delhvsqsum[indx + TS + 2] + - delhvsqsum[indx + 2 * TS - 1] + delhvsqsum[indx + 2 * TS + 1]) + - gaussgrad[5] * (delhvsqsum[indx - m2] + delhvsqsum[indx + p2] + - delhvsqsum[indx - p2] + delhvsqsum[indx + m2]); - - - if(nyqtest1 > nyqtest2) { + if(nyqutest[indx >> 1] > 0.f) { nyquist[indx >> 1] = 1; //nyquist=1 for nyquist region nystartrow = nystartrow ? nystartrow : rr; nyendrow = rr; @@ -903,46 +878,60 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, nyendrow = std::min(rr1 - 8, nyendrow); nystartcol = std::max(8, nystartcol); nyendcol = std::min(cc1 - 8, nyendcol); + memset(&nyquist2[4 * tsh], 0, sizeof(char) * (ts - 8) * tsh); + +#ifdef __SSE2__ + vint fourvb = _mm_set1_epi8(4); + vint onevb = _mm_set1_epi8(1); + +#endif for (int rr = nystartrow; rr < nyendrow; rr++) { - for (int indx = rr * TS + nystartcol + (FC(rr, 2) & 1); indx < rr * TS + nyendcol; indx += 2) { - // TODO_INGO: if you look at the comments below, it does not seem to be correct to include nyquist[indx >> 1] into the summation - // Also this implementation has loop dependencies, which are not correct IMHO - // An implementation which uses a second buffer could avoid this dependencies and could be vectorized by factor 16 too (we're working with single bytes here) - // That would lead to differences in output compared to current code, but also would lead to more consistent output when changing TS - unsigned int nyquistneighbours = (nyquist[(indx - v2) >> 1] + nyquist[(indx - m1) >> 1] + nyquist[(indx + p1) >> 1] + - nyquist[(indx - 2) >> 1] + nyquist[indx >> 1] + nyquist[(indx + 2) >> 1] + - nyquist[(indx - p1) >> 1] + nyquist[(indx + m1) >> 1] + nyquist[(indx + v2) >> 1]); +#ifdef __SSE2__ - //if most of your neighbours are named Nyquist, it's likely that you're one too - if (nyquistneighbours > 4) { - nyquist[indx >> 1] = 1; - } - - //or not - if (nyquistneighbours < 4) { - nyquist[indx >> 1] = 0; - } + for (int indx = rr * ts; indx < rr * ts + cc1; indx += 32) { + vint nyquisttemp1v = _mm_adds_epi8(_mm_load_si128((vint*)&nyquist[(indx - v2) >> 1]), _mm_load_si128((vint*)&nyquist[(indx - m1) >> 1])); + vint nyquisttemp2v = _mm_adds_epi8(_mm_load_si128((vint*)&nyquist[(indx + p1) >> 1]), _mm_load_si128((vint*)&nyquist[(indx - 2) >> 1])); + vint nyquisttemp3v = _mm_adds_epi8(_mm_load_si128((vint*)&nyquist[(indx + 2) >> 1]), _mm_load_si128((vint*)&nyquist[(indx - p1) >> 1])); + vint valv = _mm_load_si128((vint*)&nyquist[indx >> 1]); + vint nyquisttemp4v = _mm_adds_epi8(_mm_load_si128((vint*)&nyquist[(indx + m1) >> 1]), _mm_load_si128((vint*)&nyquist[(indx + v2) >> 1])); + nyquisttemp1v = _mm_adds_epi8(nyquisttemp1v, nyquisttemp3v); + nyquisttemp2v = _mm_adds_epi8(nyquisttemp2v, nyquisttemp4v); + nyquisttemp1v = _mm_adds_epi8(nyquisttemp1v, nyquisttemp2v); + valv = vselc(_mm_cmpgt_epi8(nyquisttemp1v, fourvb), onevb, valv); + valv = vselinotzero(_mm_cmplt_epi8(nyquisttemp1v, fourvb), valv); + _mm_store_si128((vint*)&nyquist2[indx >> 1], valv); } + +#else + + for (int indx = rr * ts + nystartcol + (FC(rr, 2) & 1); indx < rr * ts + nyendcol; indx += 2) { + 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 - p1) >> 1] + nyquist[(indx + m1) >> 1] + nyquist[(indx + v2) >> 1]); + //if most of your neighbours are named Nyquist, it's likely that you're one too, or not + nyquist2[indx >> 1] = nyquisttemp > 4 ? 1 : (nyquisttemp < 4 ? 0 : nyquist[indx >> 1]); + } + +#endif } // end of Nyquist test // in areas of Nyquist texture, do area interpolation 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(rr, 2) & 1); indx < rr * ts + nyendcol; indx += 2) { - if (nyquist[indx >> 1]) { + if (nyquist2[indx >> 1]) { // area interpolation float sumcfa = 0.f, sumh = 0.f, sumv = 0.f, sumsqh = 0.f, sumsqv = 0.f, areawt = 0.f; for (int i = -6; i < 7; i += 2) { - int indx1 = indx + (i * TS) - 6; + int indx1 = indx + (i * ts) - 6; for (int j = -6; j < 7; j += 2, indx1 += 2) { - - if (nyquist[indx1 >> 1]) { + if (nyquist2[indx1 >> 1]) { float cfatemp = cfa[indx1]; sumcfa += cfatemp; sumh += (cfa[indx1 - 1] + cfa[indx1 + 1]); @@ -954,11 +943,10 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, } } - //horizontal and vertical color differences, and adaptive weight + //horizontal and vertical colour differences, and adaptive weight sumh = sumcfa - xdiv2f(sumh); sumv = sumcfa - xdiv2f(sumv); - sumsqh = xdiv2f(sumsqh); - sumsqv = xdiv2f(sumsqv); + areawt = xdiv2f(areawt); float hcdvar = epssq + fabsf(areawt * sumsqh - sumh * sumh); float vcdvar = epssq + fabsf(areawt * sumsqv - sumv * sumv); hvwt[indx >> 1] = hcdvar / (vcdvar + hcdvar); @@ -972,7 +960,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, //populate G at R/B sites 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(rr, 2) & 1); indx < rr * ts + cc1 - 8; indx += 2) { //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); @@ -980,13 +968,13 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, hvwt[indx >> 1] = fabsf(0.5f - hvwt[indx >> 1]) < fabsf(0.5f - hvwtalt) ? hvwtalt : hvwt[indx >> 1]; //a better result was obtained from the neighbours - Dgrb[0][indx >> 1] = intp(hvwt[indx >> 1], vcd[indx], hcd[indx]); //evaluate color differences + Dgrb[0][indx >> 1] = intp(hvwt[indx >> 1], vcd[indx], hcd[indx]); //evaluate colour differences rgbgreen[indx] = cfa[indx] + Dgrb[0][indx >> 1]; //evaluate G (finally!) //local curvature in G (preparation for nyquist refinement step) - Dgrb2[indx >> 1].h = nyquist[indx >> 1] ? SQR(rgbgreen[indx] - xdiv2f(rgbgreen[indx - 1] + rgbgreen[indx + 1])) : 0.f; - Dgrb2[indx >> 1].v = nyquist[indx >> 1] ? SQR(rgbgreen[indx] - xdiv2f(rgbgreen[indx - v1] + rgbgreen[indx + v1])) : 0.f; + Dgrb2[indx >> 1].h = nyquist2[indx >> 1] ? SQR(rgbgreen[indx] - xdiv2f(rgbgreen[indx - 1] + rgbgreen[indx + 1])) : 0.f; + Dgrb2[indx >> 1].v = nyquist2[indx >> 1] ? SQR(rgbgreen[indx] - xdiv2f(rgbgreen[indx - v1] + rgbgreen[indx + v1])) : 0.f; } @@ -995,10 +983,11 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, // refine Nyquist areas using G curvatures if(doNyquist) { for (int rr = nystartrow; rr < nyendrow; rr++) - // TODO_INGO: maybe this part is also worth vectorizing using _mm_movemask_ps - for (int indx = rr * TS + nystartcol + (FC(rr, 2) & 1); indx < rr * TS + nyendcol; indx += 2) { - if (nyquist[indx >> 1]) { + // TODO_INGO: maybe this part is also worth vectorizing using _mm_movemask_ps + for (int indx = rr * ts + nystartcol + (FC(rr, 2) & 1); indx < rr * ts + nyendcol; indx += 2) { + + if (nyquist2[indx >> 1]) { //local averages (over Nyquist pixels only) of G curvature squared float gvarh = epssq + (gquinc[0] * Dgrb2[indx >> 1].h + gquinc[1] * (Dgrb2[(indx - m1) >> 1].h + Dgrb2[(indx + p1) >> 1].h + Dgrb2[(indx - p1) >> 1].h + Dgrb2[(indx + m1) >> 1].h) + @@ -1020,7 +1009,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, for (int rr = 6; rr < rr1 - 6; rr++) { if((FC(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 Dgrbsq1pv = (SQRV(tempv - LC2VFU(cfa[indx + 1 - p1])) + SQRV(tempv - LC2VFU(cfa[indx + 1 + p1]))); STVFU(delp[indx >> 1], vabsf(LC2VFU(cfa[indx + p1]) - LC2VFU(cfa[indx - p1]))); @@ -1030,7 +1019,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, STVFU(Dgrbsq1p[indx >> 1], Dgrbsq1pv ); } } else { - 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]); vfloat Dgrbsq1pv = (SQRV(tempv - LC2VFU(cfa[indx - p1])) + SQRV(tempv - LC2VFU(cfa[indx + p1]))); STVFU(delp[indx >> 1], vabsf(LC2VFU(cfa[indx + 1 + p1]) - LC2VFU(cfa[indx + 1 - p1]))); @@ -1046,14 +1035,14 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, for (int rr = 6; rr < rr1 - 6; rr++) { if((FC(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]); delm[indx >> 1] = fabsf(cfa[indx + m1] - cfa[indx - m1]); Dgrbsq1p[indx >> 1] = (SQR(cfa[indx + 1] - cfa[indx + 1 - p1]) + SQR(cfa[indx + 1] - cfa[indx + 1 + p1])); Dgrbsq1m[indx >> 1] = (SQR(cfa[indx + 1] - cfa[indx + 1 - m1]) + SQR(cfa[indx + 1] - cfa[indx + 1 + m1])); } } else { - 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) { Dgrbsq1p[indx >> 1] = (SQR(cfa[indx] - cfa[indx - p1]) + SQR(cfa[indx] - cfa[indx + p1])); Dgrbsq1m[indx >> 1] = (SQR(cfa[indx] - cfa[indx - m1]) + SQR(cfa[indx] - cfa[indx + m1])); delp[indx >> 1] = fabsf(cfa[indx + 1 + p1] - cfa[indx + 1 - p1]); @@ -1074,9 +1063,9 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, for (int rr = 8; rr < rr1 - 8; rr++) { #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(rr, 2) & 1), indx1 = indx >> 1; indx < rr * ts + cc1 - 8; indx += 8, indx1 += 4) { - //diagonal color ratios + //diagonal colour ratios vfloat cfav = LC2VFU(cfa[indx]); vfloat temp1v = LC2VFU(cfa[indx + m1]); @@ -1139,36 +1128,36 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #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(rr, 2) & 1), indx = rr * ts + cc, indx1 = indx >> 1; cc < cc1 - 8; cc += 2, indx += 2, indx1++) { - //diagonal color ratios + //diagonal colour ratios float crse = xmul2f(cfa[indx + m1]) / (eps + cfa[indx] + (cfa[indx + m2])); float crnw = xmul2f(cfa[indx - m1]) / (eps + cfa[indx] + (cfa[indx - m2])); float crne = xmul2f(cfa[indx + p1]) / (eps + cfa[indx] + (cfa[indx + p2])); float crsw = xmul2f(cfa[indx - p1]) / (eps + cfa[indx] + (cfa[indx - p2])); - //color differences in diagonal directions + //colour differences in diagonal directions float rbse, rbnw, rbne, rbsw; //assign B/R at R/B sites - if (fabsf(1.0f - crse) < arthresh) { + if (fabsf(1.f - crse) < arthresh) { rbse = cfa[indx] * crse; //use this if more precise diag interp is necessary } else { rbse = (cfa[indx + m1]) + xdiv2f(cfa[indx] - cfa[indx + m2]); } - if (fabsf(1.0f - crnw) < arthresh) { + if (fabsf(1.f - crnw) < arthresh) { rbnw = cfa[indx] * crnw; } else { rbnw = (cfa[indx - m1]) + xdiv2f(cfa[indx] - cfa[indx - m2]); } - if (fabsf(1.0f - crne) < arthresh) { + if (fabsf(1.f - crne) < arthresh) { rbne = cfa[indx] * crne; } else { rbne = (cfa[indx + p1]) + xdiv2f(cfa[indx] - cfa[indx + p2]); } - if (fabsf(1.0f - crsw) < arthresh) { + if (fabsf(1.f - crsw) < arthresh) { rbsw = cfa[indx] * crsw; } else { rbsw = (cfa[indx - p1]) + xdiv2f(cfa[indx] - cfa[indx - p2]); @@ -1198,7 +1187,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, rbp[indx1] = ULIM(rbp[indx1] , cfa[indx - p1], cfa[indx + p1]); } else { float pwt = xmul2f(cfa[indx] - rbp[indx1]) / (eps + rbp[indx1] + cfa[indx]); - rbp[indx1] = pwt * rbp[indx1] + (1.0f - pwt) * ULIM(rbp[indx1], cfa[indx - p1], cfa[indx + p1]); + rbp[indx1] = pwt * rbp[indx1] + (1.f - pwt) * ULIM(rbp[indx1], cfa[indx - p1], cfa[indx + p1]); } } @@ -1207,7 +1196,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, rbm[indx1] = ULIM(rbm[indx1] , cfa[indx - m1], cfa[indx + m1]); } else { float mwt = xmul2f(cfa[indx] - rbm[indx1]) / (eps + rbm[indx1] + cfa[indx]); - rbm[indx1] = mwt * rbm[indx1] + (1.0f - mwt) * ULIM(rbm[indx1], cfa[indx - m1], cfa[indx + m1]); + rbm[indx1] = mwt * rbm[indx1] + (1.f - mwt) * ULIM(rbm[indx1], cfa[indx - m1], cfa[indx + m1]); } } @@ -1229,7 +1218,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, for (int rr = 10; rr < rr1 - 10; rr++) #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(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 vfloat pmwtaltv = zd25v * (LVFU(pmwt[(indx - m1) >> 1]) + LVFU(pmwt[(indx + p1) >> 1]) + LVFU(pmwt[(indx - p1) >> 1]) + LVFU(pmwt[(indx + m1) >> 1])); @@ -1241,7 +1230,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #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(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 float pmwtalt = xdivf(pmwt[(indx - m1) >> 1] + pmwt[(indx + p1) >> 1] + pmwt[(indx - p1) >> 1] + pmwt[(indx + m1) >> 1], 2); @@ -1250,20 +1239,20 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, pmwt[indx1] = pmwtalt; //a better result was obtained from the neighbours } - rbint[indx1] = xdiv2f(cfa[indx] + rbm[indx1] * (1.0f - pmwt[indx1]) + rbp[indx1] * pmwt[indx1]); //this is R+B, interpolated + rbint[indx1] = xdiv2f(cfa[indx] + rbm[indx1] * (1.f - pmwt[indx1]) + rbp[indx1] * pmwt[indx1]); //this is R+B, interpolated } #endif for (int rr = 12; rr < rr1 - 12; rr++) #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(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]))); - if(_mm_movemask_ps((vfloat)copymask)) { // if for any of the 4 pixels the condition is true, do the math 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 //now interpolate G vertically/horizontally using R+B values - //unfortunately, since G interpolation cannot be done diagonally this may lead to color shifts - //color ratios for G interpolation + //unfortunately, since G interpolation cannot be done diagonally this may lead to colour shifts + //colour ratios for G interpolation vfloat rbintv = LVFU(rbint[indx1]); //interpolated G via adaptive ratios or Hamilton-Adams in each cardinal direction @@ -1312,16 +1301,16 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #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(rr, 2) & 1), indx = rr * ts + cc, indx1 = indx >> 1; cc < cc1 - 12; cc += 2, indx += 2, indx1++) { if (fabsf(0.5 - pmwt[indx >> 1]) < fabsf(0.5 - hvwt[indx >> 1]) ) { continue; } //now interpolate G vertically/horizontally using R+B values - //unfortunately, since G interpolation cannot be done diagonally this may lead to color shifts + //unfortunately, since G interpolation cannot be done diagonally this may lead to colour shifts - //color ratios for G interpolation + //colour ratios for G interpolation float cru = cfa[indx - v1] * 2.0 / (eps + rbint[indx1] + rbint[(indx1 - v1)]); float crd = cfa[indx + v1] * 2.0 / (eps + rbint[indx1] + rbint[(indx1 + v1)]); float crl = cfa[indx - 1] * 2.0 / (eps + rbint[indx1] + rbint[(indx1 - 1)]); @@ -1397,7 +1386,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, //fancy chrominance interpolation //(ey,ex) is location of R site for (int rr = 13 - ey; rr < rr1 - 12; rr += 2) - for (int indx1 = (rr * TS + 13 - ex) >> 1; indx1 < (rr * TS + cc1 - 12) >> 1; indx1++) { //B coset + for (int indx1 = (rr * ts + 13 - ex) >> 1; indx1 < (rr * ts + cc1 - 12) >> 1; indx1++) { //B coset Dgrb[1][indx1] = Dgrb[0][indx1]; //split out G-B from G-R Dgrb[0][indx1] = 0; } @@ -1410,7 +1399,7 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, for (int rr = 14; rr < rr1 - 14; rr++) #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(rr, 2) & 1), indx = rr * ts + cc, c = 1 - FC(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 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]))); @@ -1426,11 +1415,11 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #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) { - float wtnw = 1.0f / (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.0f / (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.0f / (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 wtse = 1.0f / (eps + fabsf(Dgrb[c][(indx + m1) >> 1] - Dgrb[c][(indx - m1) >> 1]) + fabsf(Dgrb[c][(indx + m1) >> 1] - Dgrb[c][(indx - p3) >> 1]) + fabsf(Dgrb[c][(indx - m1) >> 1] - Dgrb[c][(indx + m3) >> 1])); + 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) { + 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 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 wtse = 1.f / (eps + fabsf(Dgrb[c][(indx + m1) >> 1] - Dgrb[c][(indx - m1) >> 1]) + fabsf(Dgrb[c][(indx + m1) >> 1] - Dgrb[c][(indx - p3) >> 1]) + fabsf(Dgrb[c][(indx - m1) >> 1] - Dgrb[c][(indx + m3) >> 1])); Dgrb[c][indx >> 1] = (wtnw * (1.325f * Dgrb[c][(indx - m1) >> 1] - 0.175f * Dgrb[c][(indx - m3) >> 1] - 0.075f * Dgrb[c][(indx - m1 - 2) >> 1] - 0.075f * Dgrb[c][(indx - m1 - v2) >> 1] ) + wtne * (1.325f * Dgrb[c][(indx + p1) >> 1] - 0.175f * Dgrb[c][(indx + p3) >> 1] - 0.075f * Dgrb[c][(indx + p1 + 2) >> 1] - 0.075f * Dgrb[c][(indx + p1 + v2) >> 1] ) + @@ -1439,15 +1428,6 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, } #endif - //tile vars - //counters for pixel location in the image - int row, col; - //counters for pixel location within the tile - int cc; - //pointer counters within the tile - int indx; - - // end of tile initialization #ifdef __SSE2__ int offset; @@ -1465,12 +1445,14 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #endif for (int rr = 16; rr < rr1 - 16; rr++) { + int row = rr + top; + int col = left + 16; + int indx = rr * ts + 16; #ifdef __SSE2__ offset = 1 - offset; selmask = vnotm(selmask); - for (cc = 16, indx = rr * TS + cc, row = rr + top; cc < cc1 - 18 - (cc1 & 1); cc += 4, indx += 4) { - col = cc + left; + for (; indx < rr * ts + cc1 - 18 - (cc1 & 1); indx += 4, col += 4) { vfloat greenv = LVF(rgbgreen[indx]); vfloat temp00v = vdup(LVF(hvwt[(indx - v1) >> 1])); vfloat temp01v = vdup(LVF(hvwt[(indx + v1) >> 1])); @@ -1485,94 +1467,86 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, } if(offset == 0) { - for (indx = rr * TS + cc; cc < cc1 - 16 - (cc1 & 1); cc += 2, indx++) { - col = cc + left; - float temp = 1.0f / (hvwt[(indx - v1) >> 1] + 2.0f - hvwt[(indx + 1) >> 1] - hvwt[(indx - 1) >> 1] + hvwt[(indx + v1) >> 1]); - red[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[0][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[0][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[0][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[0][(indx + v1) >> 1]) * + 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]); + 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]) * temp); - blue[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * + blue[row][col] = 65535.f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * temp); indx++; col++; - red[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); - blue[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); + red[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); + blue[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); } if(cc1 & 1) { // width of tile is odd - col = cc + left; - float temp = 1.0f / (hvwt[(indx - v1) >> 1] + 2.0f - hvwt[(indx + 1) >> 1] - hvwt[(indx - 1) >> 1] + hvwt[(indx + v1) >> 1]); - red[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[0][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[0][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[0][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[0][(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]) * temp); - blue[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * + blue[row][col] = 65535.f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * temp); } } else { - for (indx = rr * TS + cc; cc < cc1 - 16 - (cc1 & 1); cc += 2, indx++) { - col = cc + left; - red[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); - blue[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); + for (; indx < rr * ts + cc1 - 16 - (cc1 & 1); indx++, col++) { + red[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); + blue[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); indx++; col++; - float temp = 1.0f / (hvwt[(indx - v1) >> 1] + 2.0f - hvwt[(indx + 1) >> 1] - hvwt[(indx - 1) >> 1] + hvwt[(indx + v1) >> 1]); - red[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[0][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[0][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[0][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[0][(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]) * temp); - blue[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * + blue[row][col] = 65535.f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * temp); } if(cc1 & 1) { // width of tile is odd - col = cc + left; - red[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); - blue[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); + red[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); + blue[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); } } #else if((FC(rr, 2) & 1) == 1) { - for (cc = 16, indx = rr * TS + cc, row = rr + top; cc < cc1 - 16 - (cc1 & 1); cc += 2, indx++) { - col = cc + left; - float temp = 1.0f / (hvwt[(indx - v1) >> 1] + 2.0f - hvwt[(indx + 1) >> 1] - hvwt[(indx - 1) >> 1] + hvwt[(indx + v1) >> 1]); - red[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[0][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[0][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[0][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[0][(indx + v1) >> 1]) * + 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]); + 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]) * temp); - blue[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * + blue[row][col] = 65535.f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * temp); indx++; col++; - red[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); - blue[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); + red[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); + blue[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); } if(cc1 & 1) { // width of tile is odd - col = cc + left; - float temp = 1.0f / (hvwt[(indx - v1) >> 1] + 2.0f - hvwt[(indx + 1) >> 1] - hvwt[(indx - 1) >> 1] + hvwt[(indx + v1) >> 1]); - red[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[0][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[0][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[0][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[0][(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]) * temp); - blue[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * + blue[row][col] = 65535.f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * temp); } } else { - for (cc = 16, indx = rr * TS + cc, row = rr + top; cc < cc1 - 16 - (cc1 & 1); cc += 2, indx++) { - col = cc + left; - red[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); - blue[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); + for (; indx < rr * ts + cc1 - 16 - (cc1 & 1); indx++, col++) { + red[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); + blue[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); indx++; col++; - float temp = 1.0f / (hvwt[(indx - v1) >> 1] + 2.0f - hvwt[(indx + 1) >> 1] - hvwt[(indx - 1) >> 1] + hvwt[(indx + v1) >> 1]); - red[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[0][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[0][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[0][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[0][(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]) * temp); - blue[row][col] = 65535.0f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.0f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.0f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * + blue[row][col] = 65535.f * (rgbgreen[indx] - ((hvwt[(indx - v1) >> 1]) * Dgrb[1][(indx - v1) >> 1] + (1.f - hvwt[(indx + 1) >> 1]) * Dgrb[1][(indx + 1) >> 1] + (1.f - hvwt[(indx - 1) >> 1]) * Dgrb[1][(indx - 1) >> 1] + (hvwt[(indx + v1) >> 1]) * Dgrb[1][(indx + v1) >> 1]) * temp); } if(cc1 & 1) { // width of tile is odd - col = cc + left; - red[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); - blue[row][col] = 65535.0f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); + red[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[0][indx >> 1]); + blue[row][col] = 65535.f * (rgbgreen[indx] - Dgrb[1][indx >> 1]); } } @@ -1586,18 +1560,16 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #ifdef __SSE2__ for (; cc < cc1 - 19; cc += 4) { - STVFU(green[row][cc + left], LVF(rgbgreen[rr * TS + cc]) * c65535v); + STVFU(green[row][cc + left], LVF(rgbgreen[rr * ts + cc]) * c65535v); } #endif for (; cc < cc1 - 16; cc++) { - green[row][cc + left] = 65535.0f * rgbgreen[rr * TS + cc]; + green[row][cc + left] = 65535.f * rgbgreen[rr * ts + cc]; } } - //end of main loop - if(plistener) { progresscounter++; @@ -1606,13 +1578,14 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, #pragma omp critical (amazeprogress) #endif { - progress += (double)32 * ((TS - 32) * (TS - 32)) / (height * width); + progress += (double)32 * ((ts - 32) * (ts - 32)) / (height * width); progress = progress > 1.0 ? 1.0 : progress; plistener->setProgress(progress); } } } } + } //end of main loop // clean up free(buffer); @@ -1622,9 +1595,5 @@ SSEFUNCTION void RawImageSource::amaze_demosaic_RT(int winx, int winy, int winw, plistener->setProgress(1.0); } - // done - -#undef TS - } } diff --git a/rtengine/helpersse2.h b/rtengine/helpersse2.h index 7bc480861..3f2bf6299 100644 --- a/rtengine/helpersse2.h +++ b/rtengine/helpersse2.h @@ -39,8 +39,15 @@ typedef __m128i vint2; #define STVFU(x,y) _mm_storeu_ps(&x,y) #endif -// Load 8 floats from a and combine a[0],a[2],a[4] and a[6] into a vector of 4 floats -#define LC2VFU(a) _mm_shuffle_ps( LVFU(a), _mm_loadu_ps( (&a) + 4 ), _MM_SHUFFLE( 2,0,2,0 ) ) + +static INLINE vfloat LC2VFU(float &a) +{ + // Load 8 floats from a and combine a[0],a[2],a[4] and a[6] into a vector of 4 floats + vfloat a1 = _mm_loadu_ps( &a ); + vfloat a2 = _mm_loadu_ps( (&a) + 4 ); + return _mm_shuffle_ps(a1,a2,_MM_SHUFFLE( 2,0,2,0 )); +} + // Store a vector of 4 floats in a[0],a[2],a[4] and a[6] #if defined(__x86_64__) && defined(__SSE4_1__) diff --git a/rtengine/rt_math.h b/rtengine/rt_math.h index 951d397bf..44c29fd97 100644 --- a/rtengine/rt_math.h +++ b/rtengine/rt_math.h @@ -83,8 +83,8 @@ template inline const _Tp intp(const _Tp a, const _Tp b, const _Tp c) { // calculate a * b + (1 - a) * c // following is valid: - // intp(a, b+x, c+x) = vintpf(a, b, c) + x - // intp(a, b*x, c*x) = vintpf(a, b, c) * x + // intp(a, b+x, c+x) = intp(a, b, c) + x + // intp(a, b*x, c*x) = intp(a, b, c) * x return a * (b-c) + c; } diff --git a/rtengine/sleefsseavx.c b/rtengine/sleefsseavx.c index a0300a1cc..6fed6d3d1 100644 --- a/rtengine/sleefsseavx.c +++ b/rtengine/sleefsseavx.c @@ -910,11 +910,20 @@ static INLINE vfloat vnegf(vfloat f) { return (vfloat)vxorm((vmask)f, (vmask)vca static INLINE vfloat vself(vmask mask, vfloat x, vfloat y) { return _mm_blendv_ps(y,x,(vfloat)mask); } + + static INLINE vint vselc(vmask mask, vint x, vint y) { + return _mm_blendv_epi8(y,x,mask); + } + #else // three instructions when using SSE2 static INLINE vfloat vself(vmask mask, vfloat x, vfloat y) { return (vfloat)vorm(vandm(mask, (vmask)x), vandnotm(mask, (vmask)y)); } + + static INLINE vint vselc(vmask mask, vint x, vint y) { + return vorm(vandm(mask, (vmask)x), vandnotm(mask, (vmask)y)); + } #endif static INLINE vfloat vselfzero(vmask mask, vfloat x) { @@ -928,6 +937,16 @@ static INLINE vfloat vselfnotzero(vmask mask, vfloat x) { return _mm_andnot_ps((vfloat)mask, x); } +static INLINE vint vselizero(vmask mask, vint x) { + // returns value of x if corresponding mask bits are 1, else returns 0 + // faster than vselc(mask, x, ZEROV) + return _mm_and_si128(mask, x); +} +static INLINE vint vselinotzero(vmask mask, vint x) { + // returns value of x if corresponding mask bits are 0, else returns 0 + // faster than vselc(mask, ZEROV, x) + return _mm_andnot_si128(mask, x); +} static INLINE vint2 vseli2_lt(vfloat f0, vfloat f1, vint2 x, vint2 y) { vint2 m2 = vcast_vi2_vm(vmaskf_lt(f0, f1)); @@ -1362,9 +1381,12 @@ static INLINE vfloat vaddc2vfu(float &a) // loads a[0]..a[7] and returns { a[0]+a[1], a[2]+a[3], a[4]+a[5], a[6]+a[7] } vfloat a1 = _mm_loadu_ps( &a ); vfloat a2 = _mm_loadu_ps( (&a) + 4 ); - return _mm_shuffle_ps(a1,a2,_MM_SHUFFLE( 2,0,2,0 )) + _mm_shuffle_ps(a1,a2,_MM_SHUFFLE( 3,1,3,1 )); + return _mm_shuffle_ps(a1,a2,_MM_SHUFFLE( 2,0,2,0 )) + _mm_shuffle_ps(a1,a2,_MM_SHUFFLE( 3,1,3,1 )); } +static INLINE vfloat vadivapb (vfloat a, vfloat b) { + return a / (a+b); +} #endif // __SSE2__ #endif // SLEEFSSEAVX