From 92f8d987296a16c1bd094a0c731f8eed5eb4f653 Mon Sep 17 00:00:00 2001 From: Franco Comida Date: Wed, 21 Aug 2019 20:22:12 +0200 Subject: [PATCH] hdr fusion: work in linear colorspace #185 --- librtprocess/CMakeLists.txt | 3 +- src/HdrCreation/debevec.cpp | 58 +-- src/HdrWizard/HdrCreationManager.cpp | 30 +- src/Libpfs/colorspace/gamma.h | 4 + src/Libpfs/io/jpegreader.cpp | 19 +- src/Libpfs/io/rawreader.cpp | 110 +++-- src/gauss.h | 676 --------------------------- 7 files changed, 116 insertions(+), 784 deletions(-) delete mode 100644 src/gauss.h diff --git a/librtprocess/CMakeLists.txt b/librtprocess/CMakeLists.txt index b8f783f2a..b1687a06d 100644 --- a/librtprocess/CMakeLists.txt +++ b/librtprocess/CMakeLists.txt @@ -25,7 +25,8 @@ ADD_LIBRARY(librtprocess STATIC ${CMAKE_CURRENT_SOURCE_DIR}/src/demosaic/ahd.cc ${CMAKE_CURRENT_SOURCE_DIR}/src/demosaic/rcd.cc ${CMAKE_CURRENT_SOURCE_DIR}/src/demosaic/vng4.cc ${CMAKE_CURRENT_SOURCE_DIR}/src/demosaic/xtransfast.cc - ${CMAKE_CURRENT_SOURCE_DIR}/src/postprocess/hilite_recon.cc) + ${CMAKE_CURRENT_SOURCE_DIR}/src/postprocess/hilite_recon.cc + ${CMAKE_CURRENT_SOURCE_DIR}/src/preprocess/CA_correct.cc) SET(LUMINANCE_MODULES_CLI ${LUMINANCE_MODULES_CLI} librtprocess PARENT_SCOPE) SET(LUMINANCE_MODULES_GUI ${LUMINANCE_MODULES_GUI} librtprocess PARENT_SCOPE) diff --git a/src/HdrCreation/debevec.cpp b/src/HdrCreation/debevec.cpp index 4dd32640a..f8c9a40e6 100644 --- a/src/HdrCreation/debevec.cpp +++ b/src/HdrCreation/debevec.cpp @@ -25,6 +25,7 @@ //! \author Giuseppe Rota //! \author Davide Anastasia +//! \author Franco Comida #include "HdrCreation/debevec.h" #include @@ -118,33 +119,6 @@ void DebevecOperator::computeFusion(ResponseCurve &response, images[i].frame()->getXYZChannels(Ch[0], Ch[1], Ch[2]); Array2Df *imagesCh[channels] = {Ch[0], Ch[1], Ch[2]}; - float cmax[3]; - float cmin[3]; -#ifdef _OPENMP - #pragma omp parallel for num_threads(nestedthreads) if (nestedthreads>1) -#endif - for (int c = 0; c < channels; c++) { - float minval = numeric_limits::max(); - float maxval = numeric_limits::min(); - for (size_t k = 0; k < size; k++) { - minval = std::min(minval, (*imagesCh[c])(k)); - maxval = std::max(maxval, (*imagesCh[c])(k)); - } - cmax[c] = maxval; - cmin[c] = minval; - } - - float Max = max(cmax[0], max(cmax[1], cmax[2])); - float Min = min(cmin[0], min(cmin[1], cmin[2])); - -#ifdef _OPENMP - #pragma omp parallel for num_threads(nestedthreads) if (nestedthreads>1) -#endif - for (int c = 0; c < channels; c++) { - transform(Ch[c]->begin(), Ch[c]->end(), Ch[c]->begin(), - Normalizer(Min, Max)); - } - Array2Df response_img(W, H); Array2Df w(W, H); @@ -212,22 +186,6 @@ void DebevecOperator::computeFusion(ResponseCurve &response, } } } - float cmax[3]; - for (int c = 0; c < channels; c++) { - float max = numeric_limits::min(); -#ifdef _OPENMP - #pragma omp parallel for reduction(max:max) -#endif - for (size_t k = 0; k < size; k++) { - float val = (*resultCh[c])(k); - if(std::isnormal(val)) { - max = std::max(max, val); - } - } - cmax[c] = max; - } - - float Max = max(cmax[0], max(cmax[1], cmax[2])); for (int c = 0; c < channels; c++) { #ifdef _OPENMP @@ -236,18 +194,10 @@ void DebevecOperator::computeFusion(ResponseCurve &response, for (size_t k = 0; k < size; k++) { float val = (*resultCh[c])(k); if(!std::isnormal(val)) { - (*resultCh[c])(k) = Max; + (*resultCh[c])(k) = 0.f; } - } - } - - // TODO: Investigate why scaling hdr yields better result - for (int c = 0; c < channels; c++) { -#ifdef _OPENMP - #pragma omp parallel for -#endif - for (size_t k = 0; k < size; k++) { - (*resultCh[c])(k) *= 0.1f; + (*resultCh[c])(k) = (val < 0.f) ? 0.f : val; + (*resultCh[c])(k) = (val > 1.f) ? 1.f : val; } } diff --git a/src/HdrWizard/HdrCreationManager.cpp b/src/HdrWizard/HdrCreationManager.cpp index 840db5a2f..5b2789ec8 100644 --- a/src/HdrWizard/HdrCreationManager.cpp +++ b/src/HdrWizard/HdrCreationManager.cpp @@ -561,21 +561,21 @@ pfs::Frame *HdrCreationManager::doAntiGhosting(bool patches[][agGridSize], std::unique_ptr ghosted(createHdr()); ghosted->getXYZChannels(Ch[0], Ch[1], Ch[2]); - for (int c = 0; c < 3; c++) { - cmax[c] = *max_element(Ch[c]->begin(), Ch[c]->end()); - cmin[c] = *min_element(Ch[c]->begin(), Ch[c]->end()); - } - Max = std::max(cmax[0], std::max(cmax[1], cmax[2])); - Min = std::min(cmin[0], std::min(cmin[1], cmin[2])); - - for (int c = 0; c < 3; c++) { - replace_if(Ch[c]->begin(), Ch[c]->end(), - [](float f) { return !isnormal(f); }, Max); - replace_if(Ch[c]->begin(), Ch[c]->end(), - [](float f) { return !isfinite(f); }, Max); - transform(Ch[c]->begin(), Ch[c]->end(), Ch[c]->begin(), - Normalizer(Min, Max)); - } + /* for (int c = 0; c < 3; c++) { */ + /* cmax[c] = *max_element(Ch[c]->begin(), Ch[c]->end()); */ + /* cmin[c] = *min_element(Ch[c]->begin(), Ch[c]->end()); */ + /* } */ + /* Max = std::max(cmax[0], std::max(cmax[1], cmax[2])); */ + /* Min = std::min(cmin[0], std::min(cmin[1], cmin[2])); */ + + /* for (int c = 0; c < 3; c++) { */ + /* replace_if(Ch[c]->begin(), Ch[c]->end(), */ + /* [](float f) { return !isnormal(f); }, Max); */ + /* replace_if(Ch[c]->begin(), Ch[c]->end(), */ + /* [](float f) { return !isfinite(f); }, Max); */ + /* transform(Ch[c]->begin(), Ch[c]->end(), Ch[c]->begin(), */ + /* Normalizer(Min, Max)); */ + /* } */ Rc = Ch[0]; Gc = Ch[1]; diff --git a/src/Libpfs/colorspace/gamma.h b/src/Libpfs/colorspace/gamma.h index 266d5a494..7083b5bf1 100644 --- a/src/Libpfs/colorspace/gamma.h +++ b/src/Libpfs/colorspace/gamma.h @@ -29,6 +29,10 @@ namespace pfs { namespace colorspace { +struct Inv_Gamma2_2 { + static float gamma() { return 2.2f; } +}; + struct Gamma2_2 { static float gamma() { return 1.f / 2.2f; } }; diff --git a/src/Libpfs/io/jpegreader.cpp b/src/Libpfs/io/jpegreader.cpp index a8f13e6c9..101ea18fe 100644 --- a/src/Libpfs/io/jpegreader.cpp +++ b/src/Libpfs/io/jpegreader.cpp @@ -25,6 +25,7 @@ #include #include #include +#include #include #include #include @@ -35,6 +36,8 @@ #include #include +#undef NDEBUG + using namespace pfs; #ifndef NDEBUG @@ -392,14 +395,16 @@ void JpegReader::read(Frame &frame, const Params ¶ms) { switch (m_data->cinfo()->jpeg_color_space) { case JCS_RGB: case JCS_YCbCr: { - if (xform) { - PRINT_DEBUG("Use LCMS RGB"); + /* TODO */ + /* if (xform) { */ + /* PRINT_DEBUG("Use LCMS RGB"); */ + /* read3Components(m_data->cinfo(), tempFrame, */ + /* colorspace::Convert3LCMS3(xform.data())); */ + /* } else { */ + PRINT_DEBUG("NO COLOR PROFILE"); read3Components(m_data->cinfo(), tempFrame, - colorspace::Convert3LCMS3(xform.data())); - } else { - read3Components(m_data->cinfo(), tempFrame, - colorspace::Copy()); - } + colorspace::Gamma()); + /* } */ } break; case JCS_CMYK: case JCS_YCCK: { diff --git a/src/Libpfs/io/rawreader.cpp b/src/Libpfs/io/rawreader.cpp index f5b522606..fe4eecd6b 100644 --- a/src/Libpfs/io/rawreader.cpp +++ b/src/Libpfs/io/rawreader.cpp @@ -31,6 +31,7 @@ #include #include +#include #include #include #include @@ -38,9 +39,6 @@ #include "sleef.c" #include "opthelper.h" - -#define ABS(a) ((a)<0?-(a):(a)) - bool callback(double a) { return false; } using namespace pfs; @@ -109,8 +107,8 @@ static void temperatureToRGB(double T, double RGB[3]) { struct RAWReaderParams { RAWReaderParams() - : gamma0_(1 / 2.4), - gamma1_(12.92), + : gamma0_(1), + gamma1_(1), fourColorRGB_(0), useFujiRotate_(-1), userQuality_(10), @@ -470,6 +468,8 @@ void RAWReader::read(Frame &frame, const Params ¶ms) { OUT.no_interpolation = 0; } + p.wbMethod_ = 3; + if (m_processor.dcraw_process() != LIBRAW_SUCCESS) { m_processor.recycle(); throw pfs::io::ReadException("Error Processing RAW File"); @@ -517,6 +517,21 @@ void RAWReader::read(Frame &frame, const Params ¶ms) { PRINT_DEBUG("Data size: " << image->data_size << " " << W * H * 3 * sizeof(uint16_t)); PRINT_DEBUG("W: " << W << " H: " << H); + float rCamMul = m_processor.imgdata.color.cam_mul[0]; + float gCamMul = m_processor.imgdata.color.cam_mul[1]; + float bCamMul = m_processor.imgdata.color.cam_mul[2]; + float minMult = min(min(rCamMul, gCamMul), bCamMul); + rCamMul /= minMult; + gCamMul /= minMult; + bCamMul /= minMult; + float rPreMul = m_processor.imgdata.color.pre_mul[0]; + float gPreMul = m_processor.imgdata.color.pre_mul[1]; + float bPreMul = m_processor.imgdata.color.pre_mul[2]; + minMult = min(min(rPreMul, gPreMul), bPreMul); + rPreMul /= minMult; + gPreMul /= minMult; + bPreMul /= minMult; + unsigned cf_array[2][2]; if ( isBayer() ) { PRINT_DEBUG("Bayer"); @@ -546,6 +561,7 @@ void RAWReader::read(Frame &frame, const Params ¶ms) { } if ( isBayer() ) { + const float mult = 1.0f / 65535.f; #ifdef _OPENMP #pragma omp parallel for #endif @@ -554,12 +570,13 @@ void RAWReader::read(Frame &frame, const Params ¶ms) { for(unsigned j = 0; j < W; j++) { unsigned c = cf_array[i%2][j%2]; unsigned pos = k + i*W*3; - rawdata[i][j] = raw_data[pos + c] / 65535.f; + rawdata[i][j] = raw_data[pos + c] * mult; k += 3; } } } else if ( isXtrans() ) { + const float mult = 1.0f / 65535.f; #ifdef _OPENMP #pragma omp parallel for #endif @@ -568,7 +585,7 @@ void RAWReader::read(Frame &frame, const Params ¶ms) { for(unsigned j = 0; j < W; j++) { unsigned c = xtrans[(i) % 6][(j) % 6]; unsigned pos = k + i*W*3; - rawdata[i][j] = raw_data[pos + c] / 65535.f; + rawdata[i][j] = raw_data[pos + c] * mult; k += 3; } } @@ -584,15 +601,21 @@ void RAWReader::read(Frame &frame, const Params ¶ms) { r = (float **) malloc(H * sizeof(float *)); g = (float **) malloc(H * sizeof(float *)); b = (float **) malloc(H * sizeof(float *)); - r[0] = (float *)malloc(W*H * sizeof(float)); - g[0] = (float *)malloc(W*H * sizeof(float)); - b[0] = (float *)malloc(W*H * sizeof(float)); + r[0] = Xc->data(); + g[0] = Yc->data(); + b[0] = Zc->data(); for(unsigned i = 1; i < H; i++) { r[i] = r[i-1] + W; g[i] = g[i-1] + W; b[i] = b[i-1] + W; } + if (p.chromaAberation_) { + PRINT_DEBUG("CA_correct"); + double fitparams[2][2][16]; + CA_correct(0, 0, W, H, true, 1, 0.0, 0.0, true, rawdata, rawdata, cf_array, callback, fitparams, false); + } + try { if ( isBayer() ) { switch (p.userQuality_) { @@ -637,39 +660,64 @@ void RAWReader::read(Frame &frame, const Params ¶ms) { } catch(...) { PRINT_DEBUG("DEMOSAICING FAILED"); + free ( b ); + free ( g ); + free ( r ); + free ( rawdata[0] ); + free ( rawdata ); throw pfs::io::ReadException("DEMOSICING FAILED"); } - //HLRecovery_inpaint(W, H, r, g, b, chmax, clmax, callback); + if (p.highlightsMethod_ >= 3) { + PRINT_DEBUG("HL_recovery"); + auto minmaxR = std::minmax_element(Xc->begin(), Xc->end()); + auto minmaxG = std::minmax_element(Yc->begin(), Yc->end()); + auto minmaxB = std::minmax_element(Zc->begin(), Zc->end()); -#ifdef _OPENMP - #pragma omp parallel for -#endif - for (unsigned i = 0; i < H; ++i) { - unsigned j = 0; -#ifdef __SSE2__ - for (; j < W - 3; j += 4) { - STVFU((*Xc)(j, i), LVFU(r[i][j])); - STVFU((*Yc)(j, i), LVFU(g[i][j])); - STVFU((*Zc)(j, i), LVFU(b[i][j])); - } -#endif - for (; j < W; ++j) { - (*Xc)(j, i) = r[i][j]; - (*Yc)(j, i) = g[i][j]; - (*Zc)(j, i) = b[i][j]; - } + PRINT_DEBUG("maxR: " << *minmaxR.second); + PRINT_DEBUG("minR: " << *minmaxR.first); + PRINT_DEBUG("maxG: " << *minmaxG.second); + PRINT_DEBUG("minG: " << *minmaxG.first); + PRINT_DEBUG("maxB: " << *minmaxB.second); + PRINT_DEBUG("minB: " << *minmaxB.first); + + const float chmax[3] = {*minmaxR.second, *minmaxG.second, *minmaxB.second}; + //Max clip point: + const float clmax[3] = {rCamMul, gCamMul, bCamMul}; + HLRecovery_inpaint(W, H, r, g, b, chmax, clmax, callback); } - free ( b[0] ); free ( b ); - free ( g[0] ); free ( g ); - free ( r[0] ); free ( r ); free ( rawdata[0] ); free ( rawdata ); + Array2Df *Ch[3] = {Xc, Yc, Zc}; +#ifdef _OPENMP + #pragma omp parallel for +#endif + for (size_t k = 0; k < W*H; k++) { + float r = (*Ch[0])(k); + float g = (*Ch[1])(k); + float b = (*Ch[2])(k); + if(!std::isnormal(r) || !std::isnormal(g) || !std::isnormal(b)) { + (*Ch[0])(k) = 0.f; + (*Ch[1])(k) = 0.f; + (*Ch[2])(k) = 0.f; + } + if ( (r < 0.f) || (g < 0.f) || (b < 0.f)) { + (*Ch[0])(k) = 0.f; + (*Ch[1])(k) = 0.f; + (*Ch[2])(k) = 0.f; + } + if ( (r > 1.f) || (g > 1.f) || (b > 1.f)) { + (*Ch[0])(k) = 1.f; + (*Ch[1])(k) = 1.f; + (*Ch[2])(k) = 1.f; + } + } + FrameReader::read(tempFrame, params); frame.swap(tempFrame); } diff --git a/src/gauss.h b/src/gauss.h deleted file mode 100644 index bbd3295a2..000000000 --- a/src/gauss.h +++ /dev/null @@ -1,676 +0,0 @@ -/* - * This file is part of Luminance HDR. - * - * Luminance HDR is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * Luminance HDR is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with Luminance HDR. If not, see . - * - * This file was copied from RawTherapee on 23 Nov 2017, commit d61df9d. - * - * Copyright (c) 2004-2010 Gabor Horvath - * Copyright (c) Ingo Weyrich -*/ - -#ifndef _GAUSS_H_ -#define _GAUSS_H_ - -#include -#include -#include -#include "opthelper.h" - - -namespace -{ - -template void calculateYvVFactors( const T sigma, T &b1, T &b2, T &b3, T &B, T M[3][3]) -{ - // coefficient calculation - T q; - - if (sigma < 2.5) { - q = 3.97156 - 4.14554 * sqrt (1.0 - 0.26891 * sigma); - } else { - q = 0.98711 * sigma - 0.96330; - } - - T b0 = 1.57825 + 2.44413 * q + 1.4281 * q * q + 0.422205 * q * q * q; - b1 = 2.44413 * q + 2.85619 * q * q + 1.26661 * q * q * q; - b2 = -1.4281 * q * q - 1.26661 * q * q * q; - b3 = 0.422205 * q * q * q; - 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 - 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); - -} - -// classical filtering if the support window is small and src != dst -template void gauss3x3 (T** RESTRICT src, T** RESTRICT dst, const int W, const int H, const T c0, const T c1, const T c2, const T b0, const T b1) -{ - - // first row -#ifdef _OPENMP - #pragma omp single nowait -#endif - { - dst[0][0] = src[0][0]; - - for (int j = 1; j < W - 1; j++) - { - dst[0][j] = b1 * (src[0][j - 1] + src[0][j + 1]) + b0 * src[0][j]; - } - - dst[0][W - 1] = src[0][W - 1]; - } - -#ifdef _OPENMP - #pragma omp for nowait -#endif - - for (int i = 1; i < H - 1; i++) { - dst[i][0] = b1 * (src[i - 1][0] + src[i + 1][0]) + b0 * src[i][0]; - - for (int j = 1; j < W - 1; j++) { - dst[i][j] = c2 * (src[i - 1][j - 1] + src[i - 1][j + 1] + src[i + 1][j - 1] + src[i + 1][j + 1]) + c1 * (src[i - 1][j] + src[i][j - 1] + src[i][j + 1] + src[i + 1][j]) + c0 * src[i][j]; - } - - dst[i][W - 1] = b1 * (src[i - 1][W - 1] + src[i + 1][W - 1]) + b0 * src[i][W - 1]; - } - - // last row -#ifdef _OPENMP - #pragma omp single -#endif - { - dst[H - 1][0] = src[H - 1][0]; - - for (int j = 1; j < W - 1; j++) { - dst[H - 1][j] = b1 * (src[H - 1][j - 1] + src[H - 1][j + 1]) + b0 * src[H - 1][j]; - } - - dst[H - 1][W - 1] = src[H - 1][W - 1]; - } -} - - -// use separated filter if the support window is small and src == dst -template void gaussHorizontal3 (T** src, T** dst, int W, int H, const float c0, const float c1) -{ - T temp[W] ALIGNED16; -#ifdef _OPENMP - #pragma omp for -#endif - - for (int i = 0; i < H; i++) { - for (int j = 1; j < W - 1; j++) { - temp[j] = (T)(c1 * (src[i][j - 1] + src[i][j + 1]) + c0 * src[i][j]); - } - - dst[i][0] = src[i][0]; - memcpy (dst[i] + 1, temp + 1, (W - 2)*sizeof(T)); - - dst[i][W - 1] = src[i][W - 1]; - } -} - -#ifdef __SSE2__ -template SSEFUNCTION void gaussVertical3 (T** src, T** dst, int W, int H, const float c0, const float c1) -{ - vfloat Tv = F2V(0.f), Tm1v, Tp1v; - vfloat Tv1 = F2V(0.f), Tm1v1, Tp1v1; - vfloat c0v, c1v; - c0v = F2V(c0); - c1v = F2V(c1); - -#ifdef _OPENMP - #pragma omp for nowait -#endif - - // process 8 columns per iteration for better usage of cpu cache - for (int i = 0; i < W - 7; i += 8) { - Tm1v = LVFU( src[0][i] ); - Tm1v1 = LVFU( src[0][i + 4] ); - STVFU( dst[0][i], Tm1v); - STVFU( dst[0][i + 4], Tm1v1); - - if (H > 1) { - Tv = LVFU( src[1][i]); - Tv1 = LVFU( src[1][i + 4]); - } - - for (int j = 1; j < H - 1; j++) { - Tp1v = LVFU( src[j + 1][i]); - Tp1v1 = LVFU( src[j + 1][i + 4]); - STVFU( dst[j][i], c1v * (Tp1v + Tm1v) + Tv * c0v); - STVFU( dst[j][i + 4], c1v * (Tp1v1 + Tm1v1) + Tv1 * c0v); - Tm1v = Tv; - Tm1v1 = Tv1; - Tv = Tp1v; - Tv1 = Tp1v1; - } - - STVFU( dst[H - 1][i], LVFU( src[H - 1][i])); - STVFU( dst[H - 1][i + 4], LVFU( src[H - 1][i + 4])); - } - -// Borders are done without SSE - float temp[H] ALIGNED16; -#ifdef _OPENMP - #pragma omp single -#endif - - for (int i = W - (W % 8); i < W; i++) { - for (int j = 1; j < H - 1; j++) { - temp[j] = c1 * (src[j - 1][i] + src[j + 1][i]) + c0 * src[j][i]; - } - - dst[0][i] = src[0][i]; - - for (int j = 1; j < H - 1; j++) { - dst[j][i] = temp[j]; - } - - dst[H - 1][i] = src[H - 1][i]; - } -} -#else -template void gaussVertical3 (T** src, T** dst, int W, int H, const float c0, const float c1) -{ - T temp[H] ALIGNED16; -#ifdef _OPENMP - #pragma omp for -#endif - - for (int i = 0; i < W; i++) { - for (int j = 1; j < H - 1; j++) { - temp[j] = (T)(c1 * (src[j - 1][i] + src[j + 1][i]) + c0 * src[j][i]); - } - - dst[0][i] = src[0][i]; - - for (int j = 1; j < H - 1; j++) { - dst[j][i] = temp[j]; - } - - dst[H - 1][i] = src[H - 1][i]; - } -} -#endif - -#ifdef __SSE2__ -// fast gaussian approximation if the support window is large -template SSEFUNCTION void gaussHorizontalSse (T** src, T** dst, const int W, const int H, const float sigma) -{ - double b1, b2, b3, B, M[3][3]; - calculateYvVFactors(sigma, b1, b2, b3, B, M); - - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) { - M[i][j] *= (1.0 + b2 + (b1 - b3) * b3); - M[i][j] /= (1.0 + b1 - b2 + b3) * (1.0 - b1 - b2 - b3); - } - - vfloat Rv; - vfloat Tv, Tm2v, Tm3v; - vfloat Bv, b1v, b2v, b3v; - vfloat temp2W, temp2Wp1; - float tmp[W][4] ALIGNED16; - Bv = F2V(B); - b1v = F2V(b1); - b2v = F2V(b2); - b3v = F2V(b3); - -#ifdef _OPENMP - #pragma omp for nowait -#endif - - for (int i = 0; i < H - 3; i += 4) { - Tv = _mm_set_ps(src[i][0], src[i + 1][0], src[i + 2][0], src[i + 3][0]); - Tm3v = Tv * (Bv + b1v + b2v + b3v); - STVF( tmp[0][0], Tm3v ); - - Tm2v = _mm_set_ps(src[i][1], src[i + 1][1], src[i + 2][1], src[i + 3][1]) * Bv + Tm3v * b1v + Tv * (b2v + b3v); - STVF( tmp[1][0], Tm2v ); - - Rv = _mm_set_ps(src[i][2], src[i + 1][2], src[i + 2][2], src[i + 3][2]) * Bv + Tm2v * b1v + Tm3v * b2v + Tv * b3v; - STVF( tmp[2][0], Rv ); - - for (int j = 3; j < W; j++) { - Tv = Rv; - Rv = _mm_set_ps(src[i][j], src[i + 1][j], src[i + 2][j], src[i + 3][j]) * Bv + Tv * b1v + Tm2v * b2v + Tm3v * b3v; - STVF( tmp[j][0], Rv ); - Tm3v = Tm2v; - Tm2v = Tv; - } - - Tv = _mm_set_ps(src[i][W - 1], src[i + 1][W - 1], src[i + 2][W - 1], src[i + 3][W - 1]); - - temp2Wp1 = Tv + F2V(M[2][0]) * (Rv - Tv) + F2V(M[2][1]) * ( Tm2v - Tv ) + F2V(M[2][2]) * (Tm3v - Tv); - temp2W = Tv + F2V(M[1][0]) * (Rv - Tv) + F2V(M[1][1]) * (Tm2v - Tv) + F2V(M[1][2]) * (Tm3v - Tv); - - Rv = Tv + F2V(M[0][0]) * (Rv - Tv) + F2V(M[0][1]) * (Tm2v - Tv) + F2V(M[0][2]) * (Tm3v - Tv); - STVF(tmp[W - 1][0], Rv); - - Tm2v = Bv * Tm2v + b1v * Rv + b2v * temp2W + b3v * temp2Wp1; - STVF(tmp[W - 2][0], Tm2v); - - Tm3v = Bv * Tm3v + b1v * Tm2v + b2v * Rv + b3v * temp2W; - STVF(tmp[W - 3][0], Tm3v); - - Tv = Rv; - Rv = Tm3v; - Tm3v = Tv; - - for (int j = W - 4; j >= 0; j--) { - Tv = Rv; - Rv = LVF(tmp[j][0]) * Bv + Tv * b1v + Tm2v * b2v + Tm3v * b3v; - STVF(tmp[j][0], Rv); - Tm3v = Tm2v; - Tm2v = Tv; - } - - for (int j = 0; j < W; j++) { - dst[i + 3][j] = tmp[j][0]; - dst[i + 2][j] = tmp[j][1]; - dst[i + 1][j] = tmp[j][2]; - dst[i + 0][j] = tmp[j][3]; - } - - - } - -// Borders are done without SSE -#ifdef _OPENMP - #pragma omp single -#endif - - for (int i = H - (H % 4); i < H; i++) { - tmp[0][0] = src[i][0] * (B + b1 + b2 + b3); - tmp[1][0] = B * src[i][1] + b1 * tmp[0][0] + src[i][0] * (b2 + b3); - tmp[2][0] = B * src[i][2] + b1 * tmp[1][0] + b2 * tmp[0][0] + b3 * src[i][0]; - - for (int j = 3; j < W; j++) { - tmp[j][0] = B * src[i][j] + b1 * tmp[j - 1][0] + b2 * tmp[j - 2][0] + b3 * tmp[j - 3][0]; - } - - float temp2Wm1 = src[i][W - 1] + M[0][0] * (tmp[W - 1][0] - src[i][W - 1]) + M[0][1] * (tmp[W - 2][0] - src[i][W - 1]) + M[0][2] * (tmp[W - 3][0] - src[i][W - 1]); - float temp2W = src[i][W - 1] + M[1][0] * (tmp[W - 1][0] - src[i][W - 1]) + M[1][1] * (tmp[W - 2][0] - src[i][W - 1]) + M[1][2] * (tmp[W - 3][0] - src[i][W - 1]); - float temp2Wp1 = src[i][W - 1] + M[2][0] * (tmp[W - 1][0] - src[i][W - 1]) + M[2][1] * (tmp[W - 2][0] - src[i][W - 1]) + M[2][2] * (tmp[W - 3][0] - src[i][W - 1]); - - tmp[W - 1][0] = temp2Wm1; - tmp[W - 2][0] = B * tmp[W - 2][0] + b1 * tmp[W - 1][0] + b2 * temp2W + b3 * temp2Wp1; - tmp[W - 3][0] = B * tmp[W - 3][0] + b1 * tmp[W - 2][0] + b2 * tmp[W - 1][0] + b3 * temp2W; - - for (int j = W - 4; 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 < W; j++) { - dst[i][j] = tmp[j][0]; - } - } -} -#endif - -// fast gaussian approximation if the support window is large -template void gaussHorizontal (T** src, T** dst, const int W, const int H, const double sigma) -{ - double b1, b2, b3, B, M[3][3]; - calculateYvVFactors(sigma, b1, b2, b3, B, M); - - 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); - } - - double temp2[W] ALIGNED16; - -#ifdef _OPENMP - #pragma omp for -#endif - - for (int i = 0; i < H; i++) { - - temp2[0] = B * src[i][0] + b1 * src[i][0] + b2 * src[i][0] + b3 * src[i][0]; - temp2[1] = B * src[i][1] + b1 * temp2[0] + b2 * src[i][0] + b3 * src[i][0]; - temp2[2] = B * src[i][2] + b1 * temp2[1] + b2 * temp2[0] + b3 * src[i][0]; - - for (int j = 3; j < W; j++) { - temp2[j] = B * src[i][j] + b1 * temp2[j - 1] + b2 * temp2[j - 2] + b3 * temp2[j - 3]; - } - - double temp2Wm1 = src[i][W - 1] + M[0][0] * (temp2[W - 1] - src[i][W - 1]) + M[0][1] * (temp2[W - 2] - src[i][W - 1]) + M[0][2] * (temp2[W - 3] - src[i][W - 1]); - double temp2W = src[i][W - 1] + M[1][0] * (temp2[W - 1] - src[i][W - 1]) + M[1][1] * (temp2[W - 2] - src[i][W - 1]) + M[1][2] * (temp2[W - 3] - src[i][W - 1]); - double temp2Wp1 = src[i][W - 1] + M[2][0] * (temp2[W - 1] - src[i][W - 1]) + M[2][1] * (temp2[W - 2] - src[i][W - 1]) + M[2][2] * (temp2[W - 3] - src[i][W - 1]); - - temp2[W - 1] = temp2Wm1; - temp2[W - 2] = B * temp2[W - 2] + b1 * temp2[W - 1] + b2 * temp2W + b3 * temp2Wp1; - temp2[W - 3] = B * temp2[W - 3] + b1 * temp2[W - 2] + b2 * temp2[W - 1] + b3 * temp2W; - - for (int j = W - 4; j >= 0; j--) { - temp2[j] = B * temp2[j] + b1 * temp2[j + 1] + b2 * temp2[j + 2] + b3 * temp2[j + 3]; - } - - for (int j = 0; j < W; j++) { - dst[i][j] = (T)temp2[j]; - } - - } -} - -#ifdef __SSE2__ -template SSEFUNCTION void gaussVerticalSse (T** src, T** dst, const int W, const int H, const float sigma) -{ - double b1, b2, b3, B, M[3][3]; - calculateYvVFactors(sigma, b1, b2, b3, B, M); - - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) { - M[i][j] *= (1.0 + b2 + (b1 - b3) * b3); - M[i][j] /= (1.0 + b1 - b2 + b3) * (1.0 - b1 - b2 - b3); - } - - float tmp[H][8] ALIGNED16; - vfloat Rv; - vfloat Tv, Tm2v, Tm3v; - vfloat Rv1; - vfloat Tv1, Tm2v1, Tm3v1; - vfloat Bv, b1v, b2v, b3v; - vfloat temp2W, temp2Wp1; - vfloat temp2W1, temp2Wp11; - Bv = F2V(B); - b1v = F2V(b1); - b2v = F2V(b2); - b3v = F2V(b3); - -#ifdef _OPENMP - #pragma omp for nowait -#endif - - // process 8 columns per iteration for better usage of cpu cache - for (int i = 0; i < W - 7; i += 8) { - Tv = LVFU( src[0][i]); - Tv1 = LVFU( src[0][i + 4]); - Rv = Tv * (Bv + b1v + b2v + b3v); - Rv1 = Tv1 * (Bv + b1v + b2v + b3v); - Tm3v = Rv; - Tm3v1 = Rv1; - STVF( tmp[0][0], Rv ); - STVF( tmp[0][4], Rv1 ); - - Rv = LVFU(src[1][i]) * Bv + Rv * b1v + Tv * (b2v + b3v); - Rv1 = LVFU(src[1][i + 4]) * Bv + Rv1 * b1v + Tv1 * (b2v + b3v); - Tm2v = Rv; - Tm2v1 = Rv1; - STVF( tmp[1][0], Rv ); - STVF( tmp[1][4], Rv1 ); - - Rv = LVFU(src[2][i]) * Bv + Rv * b1v + Tm3v * b2v + Tv * b3v; - Rv1 = LVFU(src[2][i + 4]) * Bv + Rv1 * b1v + Tm3v1 * b2v + Tv1 * b3v; - STVF( tmp[2][0], Rv ); - STVF( tmp[2][4], Rv1 ); - - for (int j = 3; j < H; j++) { - Tv = Rv; - Tv1 = Rv1; - Rv = LVFU(src[j][i]) * Bv + Tv * b1v + Tm2v * b2v + Tm3v * b3v; - Rv1 = LVFU(src[j][i + 4]) * Bv + Tv1 * b1v + Tm2v1 * b2v + Tm3v1 * b3v; - STVF( tmp[j][0], Rv ); - STVF( tmp[j][4], Rv1 ); - Tm3v = Tm2v; - Tm3v1 = Tm2v1; - Tm2v = Tv; - Tm2v1 = Tv1; - } - - Tv = LVFU(src[H - 1][i]); - Tv1 = LVFU(src[H - 1][i + 4]); - - temp2Wp1 = Tv + F2V(M[2][0]) * (Rv - Tv) + F2V(M[2][1]) * (Tm2v - Tv) + F2V(M[2][2]) * (Tm3v - Tv); - temp2Wp11 = Tv1 + F2V(M[2][0]) * (Rv1 - Tv1) + F2V(M[2][1]) * (Tm2v1 - Tv1) + F2V(M[2][2]) * (Tm3v1 - Tv1); - temp2W = Tv + F2V(M[1][0]) * (Rv - Tv) + F2V(M[1][1]) * (Tm2v - Tv) + F2V(M[1][2]) * (Tm3v - Tv); - temp2W1 = Tv1 + F2V(M[1][0]) * (Rv1 - Tv1) + F2V(M[1][1]) * (Tm2v1 - Tv1) + F2V(M[1][2]) * (Tm3v1 - Tv1); - - Rv = Tv + F2V(M[0][0]) * (Rv - Tv) + F2V(M[0][1]) * (Tm2v - Tv) + F2V(M[0][2]) * (Tm3v - Tv); - Rv1 = Tv1 + F2V(M[0][0]) * (Rv1 - Tv1) + F2V(M[0][1]) * (Tm2v1 - Tv1) + F2V(M[0][2]) * (Tm3v1 - Tv1); - STVFU( dst[H - 1][i], Rv ); - STVFU( dst[H - 1][i + 4], Rv1 ); - - Tm2v = Bv * Tm2v + b1v * Rv + b2v * temp2W + b3v * temp2Wp1; - Tm2v1 = Bv * Tm2v1 + b1v * Rv1 + b2v * temp2W1 + b3v * temp2Wp11; - STVFU( dst[H - 2][i], Tm2v ); - STVFU( dst[H - 2][i + 4], Tm2v1 ); - - Tm3v = Bv * Tm3v + b1v * Tm2v + b2v * Rv + b3v * temp2W; - Tm3v1 = Bv * Tm3v1 + b1v * Tm2v1 + b2v * Rv1 + b3v * temp2W1; - STVFU( dst[H - 3][i], Tm3v ); - STVFU( dst[H - 3][i + 4], Tm3v1 ); - - Tv = Rv; - Tv1 = Rv1; - Rv = Tm3v; - Rv1 = Tm3v1; - Tm3v = Tv; - Tm3v1 = Tv1; - - for (int j = H - 4; j >= 0; j--) { - Tv = Rv; - Tv1 = Rv1; - Rv = LVF(tmp[j][0]) * Bv + Tv * b1v + Tm2v * b2v + Tm3v * b3v; - Rv1 = LVF(tmp[j][4]) * Bv + Tv1 * b1v + Tm2v1 * b2v + Tm3v1 * b3v; - STVFU( dst[j][i], Rv ); - STVFU( dst[j][i + 4], Rv1 ); - Tm3v = Tm2v; - Tm3v1 = Tm2v1; - Tm2v = Tv; - Tm2v1 = Tv1; - } - } - -// Borders are done without SSE -#ifdef _OPENMP - #pragma omp single -#endif - - for (int i = W - (W % 8); i < W; i++) { - tmp[0][0] = src[0][i] * (B + b1 + b2 + b3); - tmp[1][0] = B * src[1][i] + b1 * tmp[0][0] + src[0][i] * (b2 + b3); - tmp[2][0] = B * src[2][i] + b1 * tmp[1][0] + b2 * tmp[0][0] + b3 * src[0][i]; - - for (int j = 3; j < H; j++) { - tmp[j][0] = B * src[j][i] + b1 * tmp[j - 1][0] + b2 * tmp[j - 2][0] + b3 * tmp[j - 3][0]; - } - - float temp2Hm1 = src[H - 1][i] + M[0][0] * (tmp[H - 1][0] - src[H - 1][i]) + M[0][1] * (tmp[H - 2][0] - src[H - 1][i]) + M[0][2] * (tmp[H - 3][0] - src[H - 1][i]); - float temp2H = src[H - 1][i] + M[1][0] * (tmp[H - 1][0] - src[H - 1][i]) + M[1][1] * (tmp[H - 2][0] - src[H - 1][i]) + M[1][2] * (tmp[H - 3][0] - src[H - 1][i]); - float temp2Hp1 = src[H - 1][i] + M[2][0] * (tmp[H - 1][0] - src[H - 1][i]) + M[2][1] * (tmp[H - 2][0] - src[H - 1][i]) + M[2][2] * (tmp[H - 3][0] - src[H - 1][i]); - - tmp[H - 1][0] = temp2Hm1; - tmp[H - 2][0] = B * tmp[H - 2][0] + b1 * tmp[H - 1][0] + b2 * temp2H + b3 * temp2Hp1; - tmp[H - 3][0] = B * tmp[H - 3][0] + b1 * tmp[H - 2][0] + b2 * tmp[H - 1][0] + b3 * temp2H; - - for (int j = H - 4; 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 < H; j++) { - dst[j][i] = tmp[j][0]; - } - - } -} -#endif - -template void gaussVertical (T** src, T** dst, const int W, const int H, const double sigma) -{ - double b1, b2, b3, B, M[3][3]; - calculateYvVFactors(sigma, b1, b2, b3, B, M); - - 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); - } - - // process 'numcols' columns for better usage of L1 cpu cache (especially faster for large values of H) - static const int numcols = 8; - double temp2[H][numcols] ALIGNED16; - double temp2Hm1[numcols], temp2H[numcols], temp2Hp1[numcols]; -#ifdef _OPENMP - #pragma omp for nowait -#endif - - for (unsigned int i = 0; i < static_cast(std::max(0, W - numcols + 1)); i += numcols) { - for (int k = 0; k < numcols; k++) { - temp2[0][k] = B * src[0][i + k] + b1 * src[0][i + k] + b2 * src[0][i + k] + b3 * src[0][i + k]; - temp2[1][k] = B * src[1][i + k] + b1 * temp2[0][k] + b2 * src[0][i + k] + b3 * src[0][i + k]; - temp2[2][k] = B * src[2][i + k] + b1 * temp2[1][k] + b2 * temp2[0][k] + b3 * src[0][i + k]; - } - - for (int j = 3; j < H; j++) { - for (int k = 0; k < numcols; k++) { - temp2[j][k] = B * src[j][i + k] + b1 * temp2[j - 1][k] + b2 * temp2[j - 2][k] + b3 * temp2[j - 3][k]; - } - } - - for (int k = 0; k < numcols; k++) { - temp2Hm1[k] = src[H - 1][i + k] + M[0][0] * (temp2[H - 1][k] - src[H - 1][i + k]) + M[0][1] * (temp2[H - 2][k] - src[H - 1][i + k]) + M[0][2] * (temp2[H - 3][k] - src[H - 1][i + k]); - temp2H[k] = src[H - 1][i + k] + M[1][0] * (temp2[H - 1][k] - src[H - 1][i + k]) + M[1][1] * (temp2[H - 2][k] - src[H - 1][i + k]) + M[1][2] * (temp2[H - 3][k] - src[H - 1][i + k]); - temp2Hp1[k] = src[H - 1][i + k] + M[2][0] * (temp2[H - 1][k] - src[H - 1][i + k]) + M[2][1] * (temp2[H - 2][k] - src[H - 1][i + k]) + M[2][2] * (temp2[H - 3][k] - src[H - 1][i + k]); - } - - for (int k = 0; k < numcols; k++) { - dst[H - 1][i + k] = temp2[H - 1][k] = temp2Hm1[k]; - dst[H - 2][i + k] = temp2[H - 2][k] = B * temp2[H - 2][k] + b1 * temp2[H - 1][k] + b2 * temp2H[k] + b3 * temp2Hp1[k]; - dst[H - 3][i + k] = temp2[H - 3][k] = B * temp2[H - 3][k] + b1 * temp2[H - 2][k] + b2 * temp2[H - 1][k] + b3 * temp2H[k]; - } - - for (int j = H - 4; j >= 0; j--) { - for (int k = 0; k < numcols; k++) { - dst[j][i + k] = temp2[j][k] = B * temp2[j][k] + b1 * temp2[j + 1][k] + b2 * temp2[j + 2][k] + b3 * temp2[j + 3][k]; - } - } - } - -#ifdef _OPENMP - #pragma omp single -#endif - - // process remaining columns - for (int i = W - (W % numcols); i < W; i++) { - temp2[0][0] = B * src[0][i] + b1 * src[0][i] + b2 * src[0][i] + b3 * src[0][i]; - temp2[1][0] = B * src[1][i] + b1 * temp2[0][0] + b2 * src[0][i] + b3 * src[0][i]; - temp2[2][0] = B * src[2][i] + b1 * temp2[1][0] + b2 * temp2[0][0] + b3 * src[0][i]; - - for (int j = 3; j < H; j++) { - temp2[j][0] = B * src[j][i] + b1 * temp2[j - 1][0] + b2 * temp2[j - 2][0] + b3 * temp2[j - 3][0]; - } - - double temp2Hm1 = src[H - 1][i] + M[0][0] * (temp2[H - 1][0] - src[H - 1][i]) + M[0][1] * (temp2[H - 2][0] - src[H - 1][i]) + M[0][2] * (temp2[H - 3][0] - src[H - 1][i]); - double temp2H = src[H - 1][i] + M[1][0] * (temp2[H - 1][0] - src[H - 1][i]) + M[1][1] * (temp2[H - 2][0] - src[H - 1][i]) + M[1][2] * (temp2[H - 3][0] - src[H - 1][i]); - double temp2Hp1 = src[H - 1][i] + M[2][0] * (temp2[H - 1][0] - src[H - 1][i]) + M[2][1] * (temp2[H - 2][0] - src[H - 1][i]) + M[2][2] * (temp2[H - 3][0] - src[H - 1][i]); - - dst[H - 1][i] = temp2[H - 1][0] = temp2Hm1; - dst[H - 2][i] = temp2[H - 2][0] = B * temp2[H - 2][0] + b1 * temp2[H - 1][0] + b2 * temp2H + b3 * temp2Hp1; - dst[H - 3][i] = temp2[H - 3][0] = B * temp2[H - 3][0] + b1 * temp2[H - 2][0] + b2 * temp2[H - 1][0] + b3 * temp2H; - - for (int j = H - 4; j >= 0; j--) { - dst[j][i] = temp2[j][0] = B * temp2[j][0] + b1 * temp2[j + 1][0] + b2 * temp2[j + 2][0] + b3 * temp2[j + 3][0]; - } - } -} - -template void gaussianBlurImpl(T** src, T** dst, const int W, const int H, const double sigma) -{ - static constexpr auto GAUSS_SKIP = 0.25; - static constexpr auto GAUSS_3X3_LIMIT = 0.6; - static constexpr auto GAUSS_DOUBLE = 70.0; - - if (sigma < GAUSS_SKIP) { - // don't perform filtering - if (src != dst) { - memcpy (dst[0], src[0], static_cast(W * H) * sizeof(T)); - } - } else if (sigma < GAUSS_3X3_LIMIT) { - if(src != dst) { - // If src != dst we can take the fast way - // compute 3x3 kernel values - double c0 = 1.0; - double c1 = exp( -0.5 * (lhdrengine::SQR(1.0 / sigma)) ); - double c2 = exp( -lhdrengine::SQR(1.0 / sigma) ); - - // normalize kernel values - double sum = c0 + 4.0 * (c1 + c2); - c0 /= sum; - c1 /= sum; - c2 /= sum; - // compute kernel values for border pixels - double b1 = exp (-1.0 / (2.0 * sigma * sigma)); - double bsum = 2.0 * b1 + 1.0; - b1 /= bsum; - double b0 = 1.0 / bsum; - gauss3x3 (src, dst, W, H, c0, c1, c2, b0, b1); - } else { - // compute kernel values for separated 3x3 gaussian blur - double c1 = exp (-1.0 / (2.0 * sigma * sigma)); - double csum = 2.0 * c1 + 1.0; - c1 /= csum; - double c0 = 1.0 / csum; - gaussHorizontal3 (src, dst, W, H, c0, c1); - gaussVertical3 (dst, dst, W, H, c0, c1); - } - } else { -#ifdef __SSE2__ - - if (sigma < GAUSS_DOUBLE) { - gaussHorizontalSse (src, dst, W, H, sigma); - gaussVerticalSse (dst, dst, W, H, sigma); - } else { // large sigma only with double precision - gaussHorizontal (src, dst, W, H, sigma); - gaussVertical (dst, dst, W, H, sigma); - } - -#else - - if (sigma < GAUSS_DOUBLE) { - gaussHorizontal (src, dst, W, H, sigma); - gaussVertical (dst, dst, W, H, sigma); - } else { // large sigma only with double precision - gaussHorizontal (src, dst, W, H, sigma); - gaussVertical (dst, dst, W, H, sigma); - } - -#endif - } -} -} - -void gaussianBlur(float** src, float** dst, const int W, const int H, const double sigma) -{ - gaussianBlurImpl(src, dst, W, H, sigma); -} - -void gaussianBlur(float** src, float** dst, const int W, const int H, const double sigma); - -#endif