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

run tmo_fattal02.cc through astyle

Browse Source
This commit is contained in:
Alberto Griggio
2017-11-18 17:03:57 +01:00
parent e384edba55
commit e9d3053265
1 changed files with 683 additions and 628 deletions
Download Patch File Download Diff File Expand all files Collapse all files

231
rtengine/tmo_fattal02.cc
View File

@@ -73,7 +73,8 @@
#include "StopWatch.h" #include "StopWatch.h"
#include "sleef.c" #include "sleef.c"
#include "opthelper.h" #include "opthelper.h"
namespace rtengine { namespace rtengine
{
/****************************************************************************** /******************************************************************************
* RT code * RT code
@@ -84,9 +85,11 @@ extern MyMutex *fftwMutex;
using namespace std; using namespace std;
namespace { namespace
{
class Array2Df: public array2D<float> { class Array2Df: public array2D<float>
{
typedef array2D<float> Super; typedef array2D<float> Super;
public: public:
Array2Df(): Super() {} Array2Df(): Super() {}
@@ -153,10 +156,8 @@ void downSample(const Array2Df& A, Array2Df& B)
// speed improvements. The main issue is the pde solver and in case of the // speed improvements. The main issue is the pde solver and in case of the
// fft solver uses optimised threaded fftw routines. // fft solver uses optimised threaded fftw routines.
//#pragma omp parallel for //#pragma omp parallel for
for ( int y=0 ; y<height ; y++ ) for ( int y = 0 ; y < height ; y++ ) {
{ for ( int x = 0 ; x < width ; x++ ) {
for ( int x=0 ; x<width ; x++ )
{
float p = A (2 * x, 2 * y); float p = A (2 * x, 2 * y);
p += A (2 * x + 1, 2 * y); p += A (2 * x + 1, 2 * y);
p += A (2 * x, 2 * y + 1); p += A (2 * x, 2 * y + 1);
@@ -177,6 +178,7 @@ void gaussianBlur(const Array2Df& I, Array2Df& L)
L (i) = I (i); L (i) = I (i);
} }
} }
return; return;
} }
@@ -184,25 +186,24 @@ void gaussianBlur(const Array2Df& I, Array2Df& L)
//--- X blur //--- X blur
#pragma omp parallel for shared(I, T) #pragma omp parallel for shared(I, T)
for ( int y=0 ; y<height ; y++ )
{ for ( int y = 0 ; y < height ; y++ ) {
for ( int x=1 ; x<width-1 ; x++ ) for ( int x = 1 ; x < width - 1 ; x++ ) {
{
float t = 2.f * I (x, y); float t = 2.f * I (x, y);
t += I (x - 1, y); t += I (x - 1, y);
t += I (x + 1, y); t += I (x + 1, y);
T (x, y) = t * 0.25f; // t / 4.f; T (x, y) = t * 0.25f; // t / 4.f;
} }
T (0, y) = ( 3.f * I (0, y) + I (1, y) ) * 0.25f; // / 4.f; T (0, y) = ( 3.f * I (0, y) + I (1, y) ) * 0.25f; // / 4.f;
T (width - 1, y) = ( 3.f * I (width - 1, y) + I (width - 2, y) ) * 0.25f; // / 4.f; T (width - 1, y) = ( 3.f * I (width - 1, y) + I (width - 2, y) ) * 0.25f; // / 4.f;
} }
//--- Y blur //--- Y blur
#pragma omp parallel for #pragma omp parallel for
for ( int x=0 ; x<width-7 ; x+=8 )
{ for ( int x = 0 ; x < width - 7 ; x += 8 ) {
for ( int y=1 ; y<height-1 ; y++ ) for ( int y = 1 ; y < height - 1 ; y++ ) {
{
for (int xx = 0; xx < 8; ++xx) { for (int xx = 0; xx < 8; ++xx) {
float t = 2.f * T (x + xx, y); float t = 2.f * T (x + xx, y);
t += T (x + xx, y - 1); t += T (x + xx, y - 1);
@@ -210,20 +211,21 @@ void gaussianBlur(const Array2Df& I, Array2Df& L)
L (x + xx, y) = t * 0.25f; // t/4.0f; L (x + xx, y) = t * 0.25f; // t/4.0f;
} }
} }
for (int xx = 0; xx < 8; ++xx) { for (int xx = 0; xx < 8; ++xx) {
L (x + xx, 0) = ( 3.f * T (x + xx, 0) + T (x + xx, 1) ) * 0.25f; // / 4.0f; L (x + xx, 0) = ( 3.f * T (x + xx, 0) + T (x + xx, 1) ) * 0.25f; // / 4.0f;
L (x + xx, height - 1) = ( 3.f * T (x + xx, height - 1) + T (x + xx, height - 2) ) * 0.25f; // / 4.0f; L (x + xx, height - 1) = ( 3.f * T (x + xx, height - 1) + T (x + xx, height - 2) ) * 0.25f; // / 4.0f;
} }
} }
for ( int x = width - (width % 8) ; x<width ; x++ )
{ for ( int x = width - (width % 8) ; x < width ; x++ ) {
for ( int y=1 ; y<height-1 ; y++ ) for ( int y = 1 ; y < height - 1 ; y++ ) {
{
float t = 2.f * T (x, y); float t = 2.f * T (x, y);
t += T (x, y - 1); t += T (x, y - 1);
t += T (x, y + 1); t += T (x, y + 1);
L (x, y) = t * 0.25f; // t/4.0f; L (x, y) = t * 0.25f; // t/4.0f;
} }
L (x, 0) = ( 3.f * T (x, 0) + T (x, 1) ) * 0.25f; // / 4.0f; L (x, 0) = ( 3.f * T (x, 0) + T (x, 1) ) * 0.25f; // / 4.0f;
L (x, height - 1) = ( 3.f * T (x, height - 1) + T (x, height - 2) ) * 0.25f; // / 4.0f; L (x, height - 1) = ( 3.f * T (x, height - 1) + T (x, height - 2) ) * 0.25f; // / 4.0f;
} }
@@ -236,15 +238,16 @@ void createGaussianPyramids( Array2Df* H, Array2Df** pyramids, int nlevels)
const int size = width * height; const int size = width * height;
pyramids[0] = new Array2Df (width, height); pyramids[0] = new Array2Df (width, height);
//#pragma omp parallel for shared(pyramids, H) //#pragma omp parallel for shared(pyramids, H)
for( int i=0 ; i<size ; i++ ) for ( int i = 0 ; i < size ; i++ ) {
(*pyramids[0]) (i) = (*H) (i); (*pyramids[0]) (i) = (*H) (i);
}
Array2Df* L = new Array2Df (width, height); Array2Df* L = new Array2Df (width, height);
gaussianBlur ( *pyramids[0], *L ); gaussianBlur ( *pyramids[0], *L );
for ( int k=1 ; k<nlevels ; k++ ) for ( int k = 1 ; k < nlevels ; k++ ) {
{
if (width > 2 && height > 2) { if (width > 2 && height > 2) {
width /= 2; width /= 2;
height /= 2; height /= 2;
@@ -255,6 +258,7 @@ void createGaussianPyramids( Array2Df* H, Array2Df** pyramids, int nlevels)
// there), so it might happen that we have to add some padding to // there), so it might happen that we have to add some padding to
// the gaussian pyramids // the gaussian pyramids
pyramids[k] = new Array2Df (width, height); pyramids[k] = new Array2Df (width, height);
for (int j = 0, n = width * height; j < n; ++j) { for (int j = 0, n = width * height; j < n; ++j) {
(*pyramids[k]) (j) = (*L) (j); (*pyramids[k]) (j) = (*L) (j);
} }
@@ -280,12 +284,12 @@ float calculateGradients(Array2Df* H, Array2Df* G, int k)
float avgGrad = 0.0f; float avgGrad = 0.0f;
#pragma omp parallel for reduction(+:avgGrad) #pragma omp parallel for reduction(+:avgGrad)
for( int y=0 ; y<height ; y++ )
{ for ( int y = 0 ; y < height ; y++ ) {
int n = (y == 0 ? 0 : y - 1); int n = (y == 0 ? 0 : y - 1);
int s = (y + 1 == height ? y : y + 1); int s = (y + 1 == height ? y : y + 1);
for( int x=0 ; x<width ; x++ )
{ for ( int x = 0 ; x < width ; x++ ) {
float gx, gy; float gx, gy;
int w, e; int w, e;
w = (x == 0 ? 0 : x - 1); w = (x == 0 ? 0 : x - 1);
@@ -317,10 +321,8 @@ void upSample(const Array2Df& A, Array2Df& B)
const int aheight = A.getRows(); const int aheight = A.getRows();
//#pragma omp parallel for shared(A, B) //#pragma omp parallel for shared(A, B)
for ( int y=0 ; y<height ; y++ ) for ( int y = 0 ; y < height ; y++ ) {
{ for ( int x = 0 ; x < width ; x++ ) {
for ( int x=0 ; x<width ; x++ )
{
int ax = static_cast<int> (x * 0.5f); //x / 2.f; int ax = static_cast<int> (x * 0.5f); //x / 2.f;
int ay = static_cast<int> (y * 0.5f); //y / 2.f; int ay = static_cast<int> (y * 0.5f); //y / 2.f;
ax = (ax < awidth) ? ax : awidth - 1; ax = (ax < awidth) ? ax : awidth - 1;
@@ -329,6 +331,7 @@ void upSample(const Array2Df& A, Array2Df& B)
B (x, y) = A (ax, ay); B (x, y) = A (ax, ay);
} }
} }
//--- this code below produces 'use of uninitialized value error' //--- this code below produces 'use of uninitialized value error'
// int width = A->getCols(); // int width = A->getCols();
// int height = A->getRows(); // int height = A->getRows();
@@ -355,65 +358,59 @@ void calculateFiMatrix(Array2Df* FI, Array2Df* gradients[],
Array2Df** fi = new Array2Df*[nlevels]; Array2Df** fi = new Array2Df*[nlevels];
fi[nlevels - 1] = new Array2Df (width, height); fi[nlevels - 1] = new Array2Df (width, height);
if (newfattal)
{ if (newfattal) {
//#pragma omp parallel for shared(fi) //#pragma omp parallel for shared(fi)
for ( int k = 0 ; k < width*height ; k++ ) for ( int k = 0 ; k < width * height ; k++ ) {
{
(*fi[nlevels - 1]) (k) = 1.0f; (*fi[nlevels - 1]) (k) = 1.0f;
} }
} }
for ( int k = nlevels-1; k >= 0 ; k-- ) for ( int k = nlevels - 1; k >= 0 ; k-- ) {
{
width = gradients[k]->getCols(); width = gradients[k]->getCols();
height = gradients[k]->getRows(); height = gradients[k]->getRows();
// only apply gradients to levels>=detail_level but at least to the coarsest // only apply gradients to levels>=detail_level but at least to the coarsest
if ( k >= detail_level if ( k >= detail_level
|| k == nlevels - 1 || k == nlevels - 1
|| newfattal == false) || newfattal == false) {
{
//DEBUG_STR << "calculateFiMatrix: apply gradient to level " << k << endl; //DEBUG_STR << "calculateFiMatrix: apply gradient to level " << k << endl;
#pragma omp parallel for shared(fi,avgGrad) #pragma omp parallel for shared(fi,avgGrad)
for ( int y = 0; y < height; y++ ) for ( int y = 0; y < height; y++ ) {
{ for ( int x = 0; x < width; x++ ) {
for ( int x = 0; x < width; x++ )
{
float grad = ((*gradients[k]) (x, y) < 1e-4f) ? 1e-4 : (*gradients[k]) (x, y); float grad = ((*gradients[k]) (x, y) < 1e-4f) ? 1e-4 : (*gradients[k]) (x, y);
float a = alfa * avgGrad[k]; float a = alfa * avgGrad[k];
float value = pow ((grad + noise) / a, beta - 1.0f); float value = pow ((grad + noise) / a, beta - 1.0f);
if (newfattal) if (newfattal) {
(*fi[k]) (x, y) *= value; (*fi[k]) (x, y) *= value;
else } else {
(*fi[k]) (x, y) = value; (*fi[k]) (x, y) = value;
} }
} }
} }
}
// create next level // create next level
if ( k>1 ) if ( k > 1 ) {
{
width = gradients[k - 1]->getCols(); width = gradients[k - 1]->getCols();
height = gradients[k - 1]->getRows(); height = gradients[k - 1]->getRows();
fi[k - 1] = new Array2Df (width, height); fi[k - 1] = new Array2Df (width, height);
} } else {
else
fi[0] = FI; // highest level -> result fi[0] = FI; // highest level -> result
}
if ( k>0 && newfattal ) if ( k > 0 && newfattal ) {
{
upSample (*fi[k], *fi[k - 1]); // upsample to next level upSample (*fi[k], *fi[k - 1]); // upsample to next level
gaussianBlur (*fi[k - 1], *fi[k - 1]); gaussianBlur (*fi[k - 1], *fi[k - 1]);
} }
} }
for ( int k=1 ; k<nlevels ; k++ ) for ( int k = 1 ; k < nlevels ; k++ ) {
{
delete fi[k]; delete fi[k];
} }
delete[] fi; delete[] fi;
} }
@@ -445,6 +442,7 @@ void findMaxMinPercentile(const Array2Df& I,
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp for nowait #pragma omp for nowait
#endif #endif
for (int i = 0; i < size; ++i) { for (int i = 0; i < size; ++i) {
// values are in [0;1] range, so we have to multiply with 65535 to get the histogram index // values are in [0;1] range, so we have to multiply with 65535 to get the histogram index
histothr[ (unsigned int) (65535.f * data[i])]++; histothr[ (unsigned int) (65535.f * data[i])]++;
@@ -464,6 +462,7 @@ void findMaxMinPercentile(const Array2Df& I,
while (count < minPrct * size) { while (count < minPrct * size) {
count += histo[k++]; count += histo[k++];
} }
if (k > 0) { // interpolate if (k > 0) { // interpolate
int count_ = count - histo[k - 1]; int count_ = count - histo[k - 1];
float c0 = count - minPrct * size; float c0 = count - minPrct * size;
@@ -477,6 +476,7 @@ void findMaxMinPercentile(const Array2Df& I,
while (count < maxPrct * size) { while (count < maxPrct * size) {
count += histo[k++]; count += histo[k++];
} }
if (k > 0) { // interpolate if (k > 0) { // interpolate
int count_ = count - histo[k - 1]; int count_ = count - histo[k - 1];
float c0 = count - maxPrct * size; float c0 = count - maxPrct * size;
@@ -507,9 +507,15 @@ void tmo_fattal02(size_t width,
static const float black_point = 0.1f; static const float black_point = 0.1f;
static const float white_point = 0.5f; static const float white_point = 0.5f;
static const float gamma = 1.0f; // 0.8f; static const float gamma = 1.0f; // 0.8f;
// static const int detail_level = 3; // static const int detail_level = 3;
if ( detail_level < 0 ) detail_level = 0; if ( detail_level < 0 ) {
if ( detail_level > 3 ) detail_level = 3; detail_level = 0;
}
if ( detail_level > 3 ) {
detail_level = 3;
}
// ph.setValue(2); // ph.setValue(2);
// if (ph.canceled()) return; // if (ph.canceled()) return;
@@ -538,8 +544,8 @@ void tmo_fattal02(size_t width,
float maxLum = Y (0, 0); float maxLum = Y (0, 0);
#pragma omp parallel for reduction(max:maxLum) #pragma omp parallel for reduction(max:maxLum)
for ( int i=0 ; i<size ; i++ )
{ for ( int i = 0 ; i < size ; i++ ) {
maxLum = std::max (maxLum, Y (i)); maxLum = std::max (maxLum, Y (i));
} }
@@ -553,13 +559,17 @@ void tmo_fattal02(size_t width,
vfloat tempv = F2V (temp); vfloat tempv = F2V (temp);
#endif #endif
#pragma omp for schedule(dynamic,16) #pragma omp for schedule(dynamic,16)
for (size_t i = 0 ; i < height ; ++i) { for (size_t i = 0 ; i < height ; ++i) {
size_t j = 0; size_t j = 0;
#ifdef __SSE2__ #ifdef __SSE2__
for (; j < width - 3; j += 4) { for (; j < width - 3; j += 4) {
STVFU ((*H)[i][j], xlogf (tempv * LVFU (Y[i][j]) + epsv)); STVFU ((*H)[i][j], xlogf (tempv * LVFU (Y[i][j]) + epsv));
} }
#endif #endif
for (; j < width; ++j) { for (; j < width; ++j) {
(*H)[i][j] = xlogf (temp * Y[i][j] + eps); (*H)[i][j] = xlogf (temp * Y[i][j] + eps);
} }
@@ -590,6 +600,7 @@ void tmo_fattal02(size_t width,
int fullheight = height; int fullheight = height;
int dim = std::max (width, height); int dim = std::max (width, height);
Array2Df *fullH = nullptr; Array2Df *fullH = nullptr;
if (dim > RT_dimension_cap) { if (dim > RT_dimension_cap) {
float s = float (RT_dimension_cap) / float (dim); float s = float (RT_dimension_cap) / float (dim);
Array2Df *HH = new Array2Df (width * s, height * s); Array2Df *HH = new Array2Df (width * s, height * s);
@@ -599,6 +610,7 @@ void tmo_fattal02(size_t width,
width = H->getCols(); width = H->getCols();
height = H->getRows(); height = H->getRows();
} }
/** RT */ /** RT */
// create gaussian pyramids // create gaussian pyramids
@@ -621,23 +633,25 @@ void tmo_fattal02(size_t width,
// calculate gradients and its average values on pyramid levels // calculate gradients and its average values on pyramid levels
Array2Df** gradients = new Array2Df*[nlevels]; Array2Df** gradients = new Array2Df*[nlevels];
float* avgGrad = new float[nlevels]; float* avgGrad = new float[nlevels];
for ( int k=0 ; k<nlevels ; k++ )
{ for ( int k = 0 ; k < nlevels ; k++ ) {
gradients[k] = new Array2Df (pyramids[k]->getCols(), pyramids[k]->getRows()); gradients[k] = new Array2Df (pyramids[k]->getCols(), pyramids[k]->getRows());
avgGrad[k] = calculateGradients (pyramids[k], gradients[k], k); avgGrad[k] = calculateGradients (pyramids[k], gradients[k], k);
delete pyramids[k]; delete pyramids[k];
} }
delete[] pyramids; delete[] pyramids;
// ph.setValue(12); // ph.setValue(12);
// calculate fi matrix // calculate fi matrix
Array2Df* FI = new Array2Df (width, height); Array2Df* FI = new Array2Df (width, height);
calculateFiMatrix (FI, gradients, avgGrad, nlevels, detail_level, alfa, beta, noise); calculateFiMatrix (FI, gradients, avgGrad, nlevels, detail_level, alfa, beta, noise);
// dumpPFS( "FI.pfs", FI, "Y" ); // dumpPFS( "FI.pfs", FI, "Y" );
for ( int i=0 ; i<nlevels ; i++ ) for ( int i = 0 ; i < nlevels ; i++ ) {
{
delete gradients[i]; delete gradients[i];
} }
delete[] gradients; delete[] gradients;
delete[] avgGrad; delete[] avgGrad;
// ph.setValue(16); // ph.setValue(16);
@@ -658,6 +672,7 @@ void tmo_fattal02(size_t width,
delete H; delete H;
H = fullH; H = fullH;
} }
/** RT */ /** RT */
// attenuate gradients // attenuate gradients
@@ -670,11 +685,12 @@ void tmo_fattal02(size_t width,
// Neumann conditions assume U(-1)=U(0)), see also divergence calculation // Neumann conditions assume U(-1)=U(0)), see also divergence calculation
// if (fftsolver) // if (fftsolver)
#pragma omp parallel for #pragma omp parallel for
for ( size_t y = 0 ; y < height ; y++ ) { for ( size_t y = 0 ; y < height ; y++ ) {
// sets index+1 based on the boundary assumption H(N+1)=H(N-1) // sets index+1 based on the boundary assumption H(N+1)=H(N-1)
unsigned int yp1 = (y + 1 >= height ? height - 2 : y + 1); unsigned int yp1 = (y + 1 >= height ? height - 2 : y + 1);
for ( size_t x=0 ; x<width ; x++ )
{ for ( size_t x = 0 ; x < width ; x++ ) {
// sets index+1 based on the boundary assumption H(N+1)=H(N-1) // sets index+1 based on the boundary assumption H(N+1)=H(N-1)
unsigned int xp1 = (x + 1 >= width ? width - 2 : x + 1); unsigned int xp1 = (x + 1 >= width ? width - 2 : x + 1);
// forward differences in H, so need to use between-points approx of FI // forward differences in H, so need to use between-points approx of FI
@@ -682,26 +698,38 @@ void tmo_fattal02(size_t width,
(*Gy) (x, y) = ((*H) (x, yp1) - (*H) (x, y)) * 0.5 * ((*FI) (x, yp1) + (*FI) (x, y)); (*Gy) (x, y) = ((*H) (x, yp1) - (*H) (x, y)) * 0.5 * ((*FI) (x, yp1) + (*FI) (x, y));
} }
} }
delete H; delete H;
// calculate divergence // calculate divergence
#pragma omp parallel for #pragma omp parallel for
for ( size_t y = 0; y < height; ++y )
{ for ( size_t y = 0; y < height; ++y ) {
for ( size_t x = 0; x < width; ++x ) for ( size_t x = 0; x < width; ++x ) {
{
(*FI) (x, y) = (*Gx) (x, y) + (*Gy) (x, y); (*FI) (x, y) = (*Gx) (x, y) + (*Gy) (x, y);
if ( x > 0 ) (*FI)(x,y) -= (*Gx)(x-1,y);
if ( y > 0 ) (*FI)(x,y) -= (*Gy)(x,y-1); if ( x > 0 ) {
(*FI) (x, y) -= (*Gx) (x - 1, y);
}
if ( y > 0 ) {
(*FI) (x, y) -= (*Gy) (x, y - 1);
}
// if (fftsolver) // if (fftsolver)
{ {
if (x==0) (*FI)(x,y) += (*Gx)(x,y); if (x == 0) {
if (y==0) (*FI)(x,y) += (*Gy)(x,y); (*FI) (x, y) += (*Gx) (x, y);
}
if (y == 0) {
(*FI) (x, y) += (*Gy) (x, y);
}
} }
} }
} }
//delete Gx; // RT - reused as temp buffer in solve_pde_fft, deleted later //delete Gx; // RT - reused as temp buffer in solve_pde_fft, deleted later
delete Gy; delete Gy;
@@ -732,13 +760,17 @@ void tmo_fattal02(size_t width,
vfloat gammav = F2V (gamma); vfloat gammav = F2V (gamma);
#endif #endif
#pragma omp for schedule(dynamic,16) #pragma omp for schedule(dynamic,16)
for (size_t i = 0 ; i < height ; i++) { for (size_t i = 0 ; i < height ; i++) {
size_t j = 0; size_t j = 0;
#ifdef __SSE2__ #ifdef __SSE2__
for (; j < width - 3; j += 4) { for (; j < width - 3; j += 4) {
STVFU (L[i][j], xexpf (gammav * LVFU (L[i][j]))); STVFU (L[i][j], xexpf (gammav * LVFU (L[i][j])));
} }
#endif #endif
for (; j < width; j++) { for (; j < width; j++) {
L[i][j] = xexpf ( gamma * L[i][j]); L[i][j] = xexpf ( gamma * L[i][j]);
} }
@@ -755,6 +787,7 @@ void tmo_fattal02(size_t width,
float dividor = (maxLum - minLum); float dividor = (maxLum - minLum);
#pragma omp parallel for #pragma omp parallel for
for (size_t i = 0; i < height; ++i) { for (size_t i = 0; i < height; ++i) {
for (size_t j = 0; j < width; ++j) { for (size_t j = 0; j < width; ++j) {
L[i][j] = std::max ((L[i][j] - minLum) / dividor, 0.f); L[i][j] = std::max ((L[i][j] - minLum) / dividor, 0.f);
@@ -845,17 +878,18 @@ void transform_ev2normal(Array2Df *A, Array2Df *T)
// the discrete cosine transform is not exactly the transform needed // the discrete cosine transform is not exactly the transform needed
// need to scale input values to get the right transformation // need to scale input values to get the right transformation
#pragma omp parallel for #pragma omp parallel for
for(int y=1 ; y<height-1 ; y++ )
for(int x=1 ; x<width-1 ; x++ )
(*A)(x,y)*=0.25f;
for(int x=1 ; x<width-1 ; x++ ) for (int y = 1 ; y < height - 1 ; y++ )
{ for (int x = 1 ; x < width - 1 ; x++ ) {
(*A) (x, y) *= 0.25f;
}
for (int x = 1 ; x < width - 1 ; x++ ) {
(*A) (x, 0) *= 0.5f; (*A) (x, 0) *= 0.5f;
(*A) (x, height - 1) *= 0.5f; (*A) (x, height - 1) *= 0.5f;
} }
for(int y=1 ; y<height-1 ; y++ )
{ for (int y = 1 ; y < height - 1 ; y++ ) {
(*A) (0, y) *= 0.5; (*A) (0, y) *= 0.5;
(*A) (width - 1, y) *= 0.5f; (*A) (width - 1, y) *= 0.5f;
} }
@@ -895,16 +929,18 @@ void transform_normal2ev(Array2Df *A, Array2Df *T)
// need to scale the output matrix to get the right transform // need to scale the output matrix to get the right transform
float factor = (1.0f / ((height - 1) * (width - 1))); float factor = (1.0f / ((height - 1) * (width - 1)));
#pragma omp parallel for #pragma omp parallel for
for (int y = 0 ; y < height ; y++ ) for (int y = 0 ; y < height ; y++ )
for(int x=0 ; x<width ; x++ ) for (int x = 0 ; x < width ; x++ ) {
(*T) (x, y) *= factor; (*T) (x, y) *= factor;
for(int x=0 ; x<width ; x++ ) }
{
for (int x = 0 ; x < width ; x++ ) {
(*T) (x, 0) *= 0.5f; (*T) (x, 0) *= 0.5f;
(*T) (x, height - 1) *= 0.5f; (*T) (x, height - 1) *= 0.5f;
} }
for(int y=0 ; y<height ; y++ )
{ for (int y = 0 ; y < height ; y++ ) {
(*T) (0, y) *= 0.5f; (*T) (0, y) *= 0.5f;
(*T) (width - 1, y) *= 0.5f; (*T) (width - 1, y) *= 0.5f;
} }
@@ -915,8 +951,8 @@ std::vector<double> get_lambda(int n)
{ {
assert (n > 1); assert (n > 1);
std::vector<double> v (n); std::vector<double> v (n);
for (int i=0; i<n; i++)
{ for (int i = 0; i < n; i++) {
v[i] = -4.0 * SQR (sin ((double)i / (2 * (n - 1)) * RT_PI)); v[i] = -4.0 * SQR (sin ((double)i / (2 * (n - 1)) * RT_PI));
} }
@@ -979,10 +1015,12 @@ void solve_pde_fft(Array2Df *F, Array2Df *U, Array2Df *buf, bool multithread)/*,
// activate parallel execution of fft routines // activate parallel execution of fft routines
#ifdef RT_FFTW3F_OMP #ifdef RT_FFTW3F_OMP
if (multithread) { if (multithread) {
fftwf_init_threads(); fftwf_init_threads();
fftwf_plan_with_nthreads ( omp_get_max_threads() ); fftwf_plan_with_nthreads ( omp_get_max_threads() );
} }
// #else // #else
// fftwf_plan_with_nthreads( 2 ); // fftwf_plan_with_nthreads( 2 );
#endif #endif
@@ -1020,13 +1058,13 @@ void solve_pde_fft(Array2Df *F, Array2Df *U, Array2Df *buf, bool multithread)/*,
std::vector<double> l2 = get_lambda (width); std::vector<double> l2 = get_lambda (width);
#pragma omp parallel for #pragma omp parallel for
for(int y=0 ; y<height ; y++ )
{ for (int y = 0 ; y < height ; y++ ) {
for(int x=0 ; x<width ; x++ ) for (int x = 0 ; x < width ; x++ ) {
{
(*F_tr) (x, y) = (*F_tr) (x, y) / (l1[y] + l2[x]); (*F_tr) (x, y) = (*F_tr) (x, y) / (l1[y] + l2[x]);
} }
} }
(*F_tr) (0, 0) = 0.f; // any value ok, only adds a const to the solution (*F_tr) (0, 0) = 0.f; // any value ok, only adds a const to the solution
// transforms U_tr back to the normal space // transforms U_tr back to the normal space
@@ -1042,20 +1080,24 @@ void solve_pde_fft(Array2Df *F, Array2Df *U, Array2Df *buf, bool multithread)/*,
//DEBUG_STR << "solve_pde_fft: removing constant from solution" << std::endl; //DEBUG_STR << "solve_pde_fft: removing constant from solution" << std::endl;
float max = 0.f; float max = 0.f;
#pragma omp parallel for reduction(max:max) #pragma omp parallel for reduction(max:max)
for (int i = 0; i < width * height; i++) { for (int i = 0; i < width * height; i++) {
max = std::max (max, (*U) (i)); max = std::max (max, (*U) (i));
} }
#pragma omp parallel for #pragma omp parallel for
for (int i = 0; i < width * height; i++) { for (int i = 0; i < width * height; i++) {
(*U) (i) -= max; (*U) (i) -= max;
} }
// fft parallel threads cleanup, better handled outside this function? // fft parallel threads cleanup, better handled outside this function?
#ifdef RT_FFTW3F_OMP #ifdef RT_FFTW3F_OMP
if (multithread) { if (multithread) {
fftwf_cleanup_threads(); fftwf_cleanup_threads();
} }
#endif #endif
// ph.setValue(90); // ph.setValue(90);
@@ -1123,8 +1165,10 @@ void rescale_bilinear(const Array2Df &src, Array2Df &dst, bool multithread)
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp parallel for if (multithread) #pragma omp parallel for if (multithread)
#endif #endif
for (int y = 0; y < dst.getRows(); ++y) { for (int y = 0; y < dst.getRows(); ++y) {
float ymrs = y * row_scale; float ymrs = y * row_scale;
for (int x = 0; x < dst.getCols(); ++x) { for (int x = 0; x < dst.getCols(); ++x) {
dst (x, y) = get_bilinear_value (src, x * col_scale, ymrs); dst (x, y) = get_bilinear_value (src, x * col_scale, ymrs);
} }
@@ -1141,8 +1185,10 @@ void rescale_nearest(const Array2Df &src, Array2Df &dst, bool multithread)
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp parallel for if (multithread) #pragma omp parallel for if (multithread)
#endif #endif
for (int y = 0; y < nh; ++y) { for (int y = 0; y < nh; ++y) {
int sy = y * height / nh; int sy = y * height / nh;
for (int x = 0; x < nw; ++x) { for (int x = 0; x < nw; ++x) {
int sx = x * width / nw; int sx = x * width / nw;
dst (x, y) = src (sx, sy); dst (x, y) = src (sx, sy);
@@ -1198,11 +1244,13 @@ inline int find_fast_dim(int dim)
d1 / 8 * 7, d1 / 8 * 7,
d1 d1
}; };
for (size_t i = 0; i < d.size(); ++i) { for (size_t i = 0; i < d.size(); ++i) {
if (d[i] >= dim) { if (d[i] >= dim) {
return d[i]; return d[i];
} }
} }
assert (false); assert (false);
return dim; return dim;
} }
@@ -1216,6 +1264,7 @@ void ImProcFunctions::ToneMapFattal02(Imagefloat *rgb)
const int detail_level = 3; const int detail_level = 3;
float alpha = 1.f; float alpha = 1.f;
if (params->fattal.threshold < 0) { if (params->fattal.threshold < 0) {
alpha += (params->fattal.threshold * 0.9f) / 100.f; alpha += (params->fattal.threshold * 0.9f) / 100.f;
} else if (params->fattal.threshold > 0) { } else if (params->fattal.threshold > 0) {
@@ -1242,11 +1291,13 @@ void ImProcFunctions::ToneMapFattal02(Imagefloat *rgb)
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp parallel for if (multiThread) #pragma omp parallel for if (multiThread)
#endif #endif
for (int y = 0; y < h; y++) { for (int y = 0; y < h; y++) {
for (int x = 0; x < w; x++) { for (int x = 0; x < w; x++) {
Yr (x, y) = std::max (luminance (rgb->r (y, x), rgb->g (y, x), rgb->b (y, x), ws), min_luminance); // clip really black pixels Yr (x, y) = std::max (luminance (rgb->r (y, x), rgb->g (y, x), rgb->b (y, x), ws), min_luminance); // clip really black pixels
} }
} }
// median filter on the deep shadows, to avoid boosting noise // median filter on the deep shadows, to avoid boosting noise
// because w2 >= w and h2 >= h, we can use the L buffer as temporary buffer for Median_Denoise() // because w2 >= w and h2 >= h, we can use the L buffer as temporary buffer for Median_Denoise()
int w2 = find_fast_dim (w) + 1; int w2 = find_fast_dim (w) + 1;
@@ -1260,6 +1311,7 @@ void ImProcFunctions::ToneMapFattal02(Imagefloat *rgb)
#endif #endif
float r = float (std::max (w, h)) / float (RT_dimension_cap); float r = float (std::max (w, h)) / float (RT_dimension_cap);
Median med; Median med;
if (r >= 3) { if (r >= 3) {
med = Median::TYPE_7X7; med = Median::TYPE_7X7;
} else if (r >= 2) { } else if (r >= 2) {
@@ -1269,6 +1321,7 @@ void ImProcFunctions::ToneMapFattal02(Imagefloat *rgb)
} else { } else {
med = Median::TYPE_3X3_STRONG; med = Median::TYPE_3X3_STRONG;
} }
Median_Denoise (Yr, Yr, luminance_noise_floor, w, h, med, 1, num_threads, L); Median_Denoise (Yr, Yr, luminance_noise_floor, w, h, med, 1, num_threads, L);
} }
@@ -1287,8 +1340,10 @@ void ImProcFunctions::ToneMapFattal02(Imagefloat *rgb)
#ifdef _OPENMP #ifdef _OPENMP
#pragma omp parallel for if(multiThread) #pragma omp parallel for if(multiThread)
#endif #endif
for (int y = 0; y < h; y++) { for (int y = 0; y < h; y++) {
int yy = y * h2 / h; int yy = y * h2 / h;
for (int x = 0; x < w; x++) { for (int x = 0; x < w; x++) {
int xx = x * w2 / w; int xx = x * w2 / w;
float Y = Yr (x, y); float Y = Yr (x, y);