From cde922284ee5032eca1e01be80da8837ea75d25f Mon Sep 17 00:00:00 2001 From: Ingo Date: Wed, 27 Feb 2013 23:02:42 +0100 Subject: [PATCH] Fixes compilation problem with SSE-includes and mingw32 and Linux --- rtengine/gauss.h | 797 +++++++++++++++++++++--------------------- rtengine/ipsharpen.cc | 48 +-- 2 files changed, 416 insertions(+), 429 deletions(-) diff --git a/rtengine/gauss.h b/rtengine/gauss.h index 895f2ac97..c719b6082 100644 --- a/rtengine/gauss.h +++ b/rtengine/gauss.h @@ -26,7 +26,13 @@ #ifdef _OPENMP #include #endif -#include +#ifdef __SSE__ +#if defined( WIN32 ) && defined(__x86_64__) + #include +#else + #include +#endif +#endif // classical filtering if the support window is small: @@ -70,11 +76,257 @@ template void gaussVertical3 (T** src, T** dst, AlignedBufferMP dst[H-1][i] = src[H-1][i]; } } + +#ifdef __SSE__ +template void gaussVertical3Sse (T** src, T** dst, int W, int H, const float c0, const float c1) { + + __m128 Tv,Tm1v,Tp1v; + __m128 c0v,c1v; + c0v = _mm_set1_ps(c0); + c1v = _mm_set1_ps(c1); +#ifdef _OPENMP +#pragma omp for +#endif + for (int i=0; i1) + Tv = _mm_loadu_ps( &src[1][i]); + for (int j=1; j void gaussHorizontal3Sse (T** src, T** dst, int W, int H, const float c0, const float c1) { + + float tmp[W][4] __attribute__ ((aligned (16))); + + __m128 Tv,Tm1v,Tp1v; + __m128 c0v,c1v; + c0v = _mm_set1_ps(c0); + c1v = _mm_set1_ps(c1); +#ifdef _OPENMP +#pragma omp for +#endif + for (int i=0; i1) + Tv = _mm_set_ps( src[i][1], src[i+1][1], src[i+2][1], src[i+3][1] ); + for (int j=1; j void gaussHorizontalSse (T** src, T** dst, int W, int H, float sigma) { + + if (sigma<0.25) { + // dont perform filtering + if (src!=dst) +#pragma omp for + for (int i = 0; i (src, dst, W, H, c0, c1); + return; + } + + // coefficient calculation + float q = 0.98711 * sigma - 0.96330; + if (sigma<2.5) + q = 3.97156 - 4.14554 * sqrt (1.0 - 0.26891 * sigma); + float b0 = 1.57825 + 2.44413*q + 1.4281*q*q + 0.422205*q*q*q; + float b1 = 2.44413*q + 2.85619*q*q + 1.26661*q*q*q; + float b2 = -1.4281*q*q - 1.26661*q*q*q; + float b3 = 0.422205*q*q*q; + float B = 1.0 - (b1+b2+b3) / b0; + + b1 /= b0; + b2 /= b0; + b3 /= b0; + + // From: Bill Triggs, Michael Sdika: Boundary Conditions for Young-van Vliet Recursive Filtering + float M[3][3]; + M[0][0] = -b3*b1+1.0-b3*b3-b2; + M[0][1] = (b3+b1)*(b2+b3*b1); + M[0][2] = b3*(b1+b3*b2); + M[1][0] = b1+b3*b2; + M[1][1] = -(b2-1.0)*(b2+b3*b1); + M[1][2] = -(b3*b1+b3*b3+b2-1.0)*b3; + M[2][0] = b3*b1+b2+b1*b1-b2*b2; + M[2][1] = b1*b2+b3*b2*b2-b1*b3*b3-b3*b3*b3-b3*b2+b3; + M[2][2] = b3*(b1+b3*b2); + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + M[i][j] /= (1.0+b1-b2+b3)*(1.0+b2+(b1-b3)*b3); + + float tmp[W][4] __attribute__ ((aligned (16))); + float tmpV[4] __attribute__ ((aligned (16))); + __m128 Rv; + __m128 Tv,Tm2v,Tm3v; + __m128 Bv,b1v,b2v,b3v; + __m128 temp2W,temp2Wp1; + Bv = _mm_set1_ps(B); + b1v = _mm_set1_ps(b1); + b2v = _mm_set1_ps(b2); + b3v = _mm_set1_ps(b3); +#pragma omp for + for (int i=0; i=0; j--) { + Tv = Rv; + Rv = _mm_load_ps(&tmp[j][0]) * Bv + Tv * b1v + Tm2v * b2v + Tm3v * b3v; + _mm_store_ps( &tmp[j][0], Rv ); + Tm3v = Tm2v; + Tm2v = Tv; + } + + for (int j=0; j=0; j--) + tmp[j][0] = B * tmp[j][0] + b1*tmp[j+1][0] + b2*tmp[j+2][0] + b3*tmp[j+3][0]; + + for (int j=0; j void gaussHorizontal (T** src, T** dst, AlignedBufferMP &buffer, int W, int H, double sigma) { - + +#ifdef __SSE__ + gaussHorizontalSse (src, dst, W, H, sigma); + return; +#else if (sigma<0.25) { // dont perform filtering if (src!=dst) @@ -149,11 +401,157 @@ template void gaussHorizontal (T** src, T** dst, AlignedBufferMP void gaussVerticalSse (T** src, T** dst, int W, int H, float sigma) { + if (sigma<0.25) { + // dont perform filtering + if (src!=dst) +#pragma omp for + for (int i = 0; i (src, dst, W, H, c0, c1); + return; + } + + // coefficient calculation + double q = 0.98711 * sigma - 0.96330; + if (sigma<2.5) + q = 3.97156 - 4.14554 * sqrt (1.0 - 0.26891 * sigma); + double b0 = 1.57825 + 2.44413*q + 1.4281*q*q + 0.422205*q*q*q; + double b1 = 2.44413*q + 2.85619*q*q + 1.26661*q*q*q; + double b2 = -1.4281*q*q - 1.26661*q*q*q; + double b3 = 0.422205*q*q*q; + double B = 1.0 - (b1+b2+b3) / b0; + + b1 /= b0; + b2 /= b0; + b3 /= b0; + + // From: Bill Triggs, Michael Sdika: Boundary Conditions for Young-van Vliet Recursive Filtering + double M[3][3]; + M[0][0] = -b3*b1+1.0-b3*b3-b2; + M[0][1] = (b3+b1)*(b2+b3*b1); + M[0][2] = b3*(b1+b3*b2); + M[1][0] = b1+b3*b2; + M[1][1] = -(b2-1.0)*(b2+b3*b1); + M[1][2] = -(b3*b1+b3*b3+b2-1.0)*b3; + M[2][0] = b3*b1+b2+b1*b1-b2*b2; + M[2][1] = b1*b2+b3*b2*b2-b1*b3*b3-b3*b3*b3-b3*b2+b3; + M[2][2] = b3*(b1+b3*b2); + for (int i=0; i<3; i++) + for (int j=0; j<3; j++) + M[i][j] /= (1.0+b1-b2+b3)*(1.0+b2+(b1-b3)*b3); + + float tmp[H][4] __attribute__ ((aligned (16))); + __m128 Rv; + __m128 Tv,Tm2v,Tm3v; + __m128 Bv,b1v,b2v,b3v; + __m128 temp2W,temp2Wp1; + Bv = _mm_set1_ps(B); + b1v = _mm_set1_ps(b1); + b2v = _mm_set1_ps(b2); + b3v = _mm_set1_ps(b3); + + +#ifdef _OPENMP +#pragma omp for +#endif + for (int i=0; i=0; j--) { + Tv = Rv; + Rv = _mm_load_ps(&tmp[j][0]) * Bv + Tv * b1v + Tm2v * b2v + Tm3v * b3v; + _mm_storeu_ps( &dst[j][i], Rv ); + Tm3v = Tm2v; + Tm2v = Tv; + } + } +// Borders are done without SSE +#pragma omp for + for(int i=W-(W%4);i=0; j--) + tmp[j][0] = B * tmp[j][0] + b1*tmp[j+1][0] + b2*tmp[j+2][0] + b3*tmp[j+3][0]; + + for (int j=0; j void gaussVertical (T** src, T** dst, AlignedBufferMP &buffer, int W, int H, double sigma) { + +#ifdef __SSE__ + gaussVerticalSse (src, dst, W, H, sigma); + return; +#else if (sigma<0.25) { // dont perform filtering @@ -230,7 +628,8 @@ template void gaussVertical (T** src, T** dst, AlignedBufferMP dst[j][i] = (T)temp2[j]; buffer.release(pBuf); - } + } +#endif } @@ -421,396 +820,4 @@ template void gaussDerivV (T** src, T** dst, AlignedBufferMP &b } } -template void gaussVertical3Sse (T** src, T** dst, int W, int H, const float c0, const float c1) { - - __m128 Tv,Tm1v,Tp1v; - __m128 c0v,c1v; - c0v = _mm_set1_ps(c0); - c1v = _mm_set1_ps(c1); -#ifdef _OPENMP -#pragma omp for -#endif - for (int i=0; i1) - Tv = _mm_loadu_ps( &src[1][i]); - for (int j=1; j void gaussHorizontal3Sse (T** src, T** dst, int W, int H, const float c0, const float c1) { - - float tmp[W][4] __attribute__ ((aligned (16))); - - __m128 Tv,Tm1v,Tp1v; - __m128 c0v,c1v; - c0v = _mm_set1_ps(c0); - c1v = _mm_set1_ps(c1); -#ifdef _OPENMP -#pragma omp for -#endif - for (int i=0; i1) - Tv = _mm_set_ps( src[i][1], src[i+1][1], src[i+2][1], src[i+3][1] ); - for (int j=1; j void gaussHorizontalSse (T** src, T** dst, int W, int H, double sigma) { - - if (sigma<0.25) { - // dont perform filtering - if (src!=dst) -#pragma omp for - for (int i = 0; i (src, dst, W, H, c0, c1); - return; - } - - // coefficient calculation - float q = 0.98711 * sigma - 0.96330; - if (sigma<2.5) - q = 3.97156 - 4.14554 * sqrt (1.0 - 0.26891 * sigma); - float b0 = 1.57825 + 2.44413*q + 1.4281*q*q + 0.422205*q*q*q; - float b1 = 2.44413*q + 2.85619*q*q + 1.26661*q*q*q; - float b2 = -1.4281*q*q - 1.26661*q*q*q; - float b3 = 0.422205*q*q*q; - float B = 1.0 - (b1+b2+b3) / b0; - - b1 /= b0; - b2 /= b0; - b3 /= b0; - - // From: Bill Triggs, Michael Sdika: Boundary Conditions for Young-van Vliet Recursive Filtering - float M[3][3]; - M[0][0] = -b3*b1+1.0-b3*b3-b2; - M[0][1] = (b3+b1)*(b2+b3*b1); - M[0][2] = b3*(b1+b3*b2); - M[1][0] = b1+b3*b2; - M[1][1] = -(b2-1.0)*(b2+b3*b1); - M[1][2] = -(b3*b1+b3*b3+b2-1.0)*b3; - M[2][0] = b3*b1+b2+b1*b1-b2*b2; - M[2][1] = b1*b2+b3*b2*b2-b1*b3*b3-b3*b3*b3-b3*b2+b3; - M[2][2] = b3*(b1+b3*b2); - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - M[i][j] /= (1.0+b1-b2+b3)*(1.0+b2+(b1-b3)*b3); - - float tmp[W][4] __attribute__ ((aligned (16))); - float tmpV[4] __attribute__ ((aligned (16))); - __m128 Rv; - __m128 Tv,Tm2v,Tm3v; - __m128 Bv,b1v,b2v,b3v; - __m128 temp2W,temp2Wp1; - Bv = _mm_set1_ps(B); - b1v = _mm_set1_ps(b1); - b2v = _mm_set1_ps(b2); - b3v = _mm_set1_ps(b3); -#pragma omp for - for (int i=0; i=0; j--) { - Tv = Rv; - Rv = _mm_add_ps( _mm_mul_ps( _mm_load_ps(&tmp[j][0]), Bv ), _mm_add_ps( _mm_mul_ps( Tv, b1v), _mm_add_ps( _mm_mul_ps( Tm2v, b2v), _mm_mul_ps( Tm3v, b3v) ))); - _mm_store_ps( &tmp[j][0], Rv ); - Tm3v = Tm2v; - Tm2v = Tv; - } - - for (int j=0; j=0; j--) - tmp[j][0] = B * tmp[j][0] + b1*tmp[j+1][0] + b2*tmp[j+2][0] + b3*tmp[j+3][0]; - - for (int j=0; j void gaussVerticalSse (T** src, T** dst, int W, int H, double sigma) { - - if (sigma<0.25) { - // dont perform filtering - if (src!=dst) -#pragma omp for - for (int i = 0; i (src, dst, W, H, c0, c1); - return; - } - - // coefficient calculation - double q = 0.98711 * sigma - 0.96330; - if (sigma<2.5) - q = 3.97156 - 4.14554 * sqrt (1.0 - 0.26891 * sigma); - double b0 = 1.57825 + 2.44413*q + 1.4281*q*q + 0.422205*q*q*q; - double b1 = 2.44413*q + 2.85619*q*q + 1.26661*q*q*q; - double b2 = -1.4281*q*q - 1.26661*q*q*q; - double b3 = 0.422205*q*q*q; - double B = 1.0 - (b1+b2+b3) / b0; - - b1 /= b0; - b2 /= b0; - b3 /= b0; - - // From: Bill Triggs, Michael Sdika: Boundary Conditions for Young-van Vliet Recursive Filtering - double M[3][3]; - M[0][0] = -b3*b1+1.0-b3*b3-b2; - M[0][1] = (b3+b1)*(b2+b3*b1); - M[0][2] = b3*(b1+b3*b2); - M[1][0] = b1+b3*b2; - M[1][1] = -(b2-1.0)*(b2+b3*b1); - M[1][2] = -(b3*b1+b3*b3+b2-1.0)*b3; - M[2][0] = b3*b1+b2+b1*b1-b2*b2; - M[2][1] = b1*b2+b3*b2*b2-b1*b3*b3-b3*b3*b3-b3*b2+b3; - M[2][2] = b3*(b1+b3*b2); - for (int i=0; i<3; i++) - for (int j=0; j<3; j++) - M[i][j] /= (1.0+b1-b2+b3)*(1.0+b2+(b1-b3)*b3); - - float tmp[H][4] __attribute__ ((aligned (16))); - float tmpV[4] __attribute__ ((aligned (16))); - __m128 Rv; - __m128 Tv,Tm2v,Tm3v; - __m128 Bv,b1v,b2v,b3v; - __m128 temp2W,temp2Wp1; - Bv = _mm_set1_ps(B); - b1v = _mm_set1_ps(b1); - b2v = _mm_set1_ps(b2); - b3v = _mm_set1_ps(b3); - - -#ifdef _OPENMP -#pragma omp for -#endif - for (int i=0; i=0; j--) { - Tv = Rv; - Rv = _mm_add_ps( _mm_mul_ps( _mm_load_ps(&tmp[j][0]), Bv ), _mm_add_ps( _mm_mul_ps( Tv, b1v), _mm_add_ps( _mm_mul_ps( Tm2v, b2v), _mm_mul_ps( Tm3v, b3v) ))); - _mm_storeu_ps( &dst[j][i], Rv ); - Tm3v = Tm2v; - Tm2v = Tv; - } - } -// Borders are done without SSE - int start = W-(W%4); -#pragma omp for - for(int i=start;i=0; j--) - tmp[j][0] = B * tmp[j][0] + b1*tmp[j+1][0] + b2*tmp[j+2][0] + b3*tmp[j+3][0]; - - for (int j=0; j buffer(max(W,H)); -#endif float damping = params->sharpening.deconvdamping / 5.0; bool needdamp = params->sharpening.deconvdamping > 0; for (int k=0; ksharpening.deconviter; k++) { // apply blur function (gaussian blur) -#ifdef __SSE__ - gaussHorizontalSse (tmpI, tmp, W, H, params->sharpening.deconvradius / scale); - gaussVerticalSse (tmp, tmp, W, H, params->sharpening.deconvradius / scale); -#else - gaussHorizontal (tmpI, tmp, buffer, W, H, params->sharpening.deconvradius / scale); - gaussVertical (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); -#endif + gaussHorizontal (tmpI, tmp, buffer, W, H, params->sharpening.deconvradius / scale); + gaussVertical (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); if (!needdamp) { #ifdef _OPENMP @@ -109,13 +102,8 @@ void ImProcFunctions::deconvsharpening (LabImage* lab, float** b2) { else dcdamping (tmp, lab->L, damping, W, H); -#ifdef __SSE__ - gaussHorizontalSse (tmp, tmp, W, H, params->sharpening.deconvradius / scale); - gaussVerticalSse (tmp, tmp, W, H, params->sharpening.deconvradius / scale); -#else - gaussHorizontal (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); - gaussVertical (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); -#endif + gaussHorizontal (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); + gaussVertical (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); #ifdef _OPENMP #pragma omp for @@ -139,7 +127,8 @@ void ImProcFunctions::deconvsharpening (LabImage* lab, float** b2) { for (int i=0; i buffer(max(W,H)); -#endif + AlignedBufferMP buffer(max(W,H)); + float damping = params->sharpening.deconvdamping / 5.0; bool needdamp = params->sharpening.deconvdamping > 0; for (int k=0; ksharpening.deconviter; k++) { // apply blur function (gaussian blur) -#ifdef __SSE__ - gaussHorizontalSse (tmpI, tmp, W, H, params->sharpening.deconvradius / scale); - gaussVerticalSse (tmp, tmp, W, H, params->sharpening.deconvradius / scale); -#else - gaussHorizontal (tmpI, tmp, buffer, W, H, params->sharpening.deconvradius / scale); - gaussVertical (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); -#endif + gaussHorizontal (tmpI, tmp, buffer, W, H, params->sharpening.deconvradius / scale); + gaussVertical (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); if (!needdamp) { #ifdef _OPENMP @@ -877,13 +860,9 @@ void ImProcFunctions::deconvsharpeningcam (CieImage* ncie, float** b2) { else dcdamping (tmp, ncie->sh_p, damping, W, H); -#ifdef __SSE__ - gaussHorizontalSse (tmp, tmp, W, H, params->sharpening.deconvradius / scale); - gaussVerticalSse (tmp, tmp, W, H, params->sharpening.deconvradius / scale); -#else - gaussHorizontal (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); - gaussVertical (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); -#endif + gaussHorizontal (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); + gaussVertical (tmp, tmp, buffer, W, H, params->sharpening.deconvradius / scale); + #ifdef _OPENMP #pragma omp for @@ -909,6 +888,7 @@ void ImProcFunctions::deconvsharpeningcam (CieImage* ncie, float** b2) { for (int i=0; i