diff --git a/CMakeLists.txt b/CMakeLists.txt
index a62d68d53..662a7ff3d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -183,7 +183,9 @@ endif (WIN32)
 # you may need lcms v1.xx for older version : pkg_check_modules (LCMS REQUIRED lcms<=1.99)
 pkg_check_modules (LCMS REQUIRED lcms2)
 find_package (EXPAT REQUIRED expat>=2.0)
+pkg_check_modules (FFTW3F REQUIRED fftw3f)
 pkg_check_modules (IPTCDATA REQUIRED libiptcdata)
+pkg_check_modules(FFTW3 fftw3)
 find_package (JPEG REQUIRED)
 find_package (PNG REQUIRED)
 find_package (TIFF REQUIRED)
diff --git a/rtdata/languages/default b/rtdata/languages/default
index 8e86989a2..cadc7e1ee 100644
--- a/rtdata/languages/default
+++ b/rtdata/languages/default
@@ -395,6 +395,7 @@ HISTORY_MSG_169;'CH' curve
 HISTORY_MSG_170;Vibrance - Curve
 HISTORY_MSG_171;'LC' curve
 HISTORY_MSG_172;Restrict 'LC' to Red and Skin Tones
+HISTORY_MSG_173;NR - Luminance Detail
 HISTORY_NEWSNAPSHOTAS;As...
 HISTORY_NEWSNAPSHOT;Add
 HISTORY_NEWSSDIALOGLABEL;Label of the snapshot:
@@ -886,7 +887,8 @@ TP_DEFRINGE_THRESHOLD;Threshold
 TP_DETAIL_AMOUNT;Amount
 TP_DIRPYRDENOISE_CHROMA;Chrominance
 TP_DIRPYRDENOISE_GAMMA;Gamma
-TP_DIRPYRDENOISE_LABEL;Noise reduction
+TP_DIRPYRDENOISE_LABEL;Noise reduction (RAW images ONLY)
+TP_DIRPYRDENOISE_LDETAIL;Luminance Detail
 TP_DIRPYRDENOISE_LUMA;Luminance
 TP_DIRPYREQUALIZER_LABEL;Contrast by detail levels
 TP_DIRPYREQUALIZER_LUMACOARSEST;Coarsest
diff --git a/rtdata/profiles/BW-1.pp3 b/rtdata/profiles/BW-1.pp3
index a082a2a71..2cbe7629a 100644
--- a/rtdata/profiles/BW-1.pp3
+++ b/rtdata/profiles/BW-1.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/BW-2.pp3 b/rtdata/profiles/BW-2.pp3
index 36ae302a6..1c6b8bc65 100644
--- a/rtdata/profiles/BW-2.pp3
+++ b/rtdata/profiles/BW-2.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/BW-3.pp3 b/rtdata/profiles/BW-3.pp3
index 0cbf52d8a..d531f7e49 100644
--- a/rtdata/profiles/BW-3.pp3
+++ b/rtdata/profiles/BW-3.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/BW-4.pp3 b/rtdata/profiles/BW-4.pp3
index 19447157c..fe91696e7 100644
--- a/rtdata/profiles/BW-4.pp3
+++ b/rtdata/profiles/BW-4.pp3
@@ -94,11 +94,13 @@ Enabled=false
 Radius=2.0
 Threshold=25
 
+
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/Default-ISO-High.pp3 b/rtdata/profiles/Default-ISO-High.pp3
index 2c99a83b7..399b1566a 100644
--- a/rtdata/profiles/Default-ISO-High.pp3
+++ b/rtdata/profiles/Default-ISO-High.pp3
@@ -86,7 +86,7 @@ Temperature=5745
 Green=1.0
 
 [Impulse Denoising]
-Enabled=true
+Enabled=false
 Threshold=80
 
 [Defringing]
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=true
-Luma=0
-Chroma=15
-Gamma=1.2
+Luma=50
+Ldetail=80
+Chroma=50
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/Default-ISO-Medium.pp3 b/rtdata/profiles/Default-ISO-Medium.pp3
index 62eacf3d9..7d0474b37 100644
--- a/rtdata/profiles/Default-ISO-Medium.pp3
+++ b/rtdata/profiles/Default-ISO-Medium.pp3
@@ -86,7 +86,7 @@ Temperature=5745
 Green=1.0
 
 [Impulse Denoising]
-Enabled=true
+Enabled=false
 Threshold=50
 
 [Defringing]
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=true
-Luma=0
-Chroma=5
-Gamma=1.2
+Luma=20
+Ldetail=80
+Chroma=30
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/Default.pp3 b/rtdata/profiles/Default.pp3
index bb1bca21b..c29e3256b 100644
--- a/rtdata/profiles/Default.pp3
+++ b/rtdata/profiles/Default.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/Highkey-1.pp3 b/rtdata/profiles/Highkey-1.pp3
index 6098b6be1..d6b9ac28b 100644
--- a/rtdata/profiles/Highkey-1.pp3
+++ b/rtdata/profiles/Highkey-1.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/Natural-1.pp3 b/rtdata/profiles/Natural-1.pp3
index cba04b30e..153f5a6b3 100644
--- a/rtdata/profiles/Natural-1.pp3
+++ b/rtdata/profiles/Natural-1.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/Natural-2.pp3 b/rtdata/profiles/Natural-2.pp3
index 954d54ff4..1f65f5f67 100644
--- a/rtdata/profiles/Natural-2.pp3
+++ b/rtdata/profiles/Natural-2.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/Neutral.pp3 b/rtdata/profiles/Neutral.pp3
index 9b71f7579..f985bad8d 100644
--- a/rtdata/profiles/Neutral.pp3
+++ b/rtdata/profiles/Neutral.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/Punchy-1.pp3 b/rtdata/profiles/Punchy-1.pp3
index 2fbdc0fbd..5c4b9c235 100644
--- a/rtdata/profiles/Punchy-1.pp3
+++ b/rtdata/profiles/Punchy-1.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtdata/profiles/Punchy-2.pp3 b/rtdata/profiles/Punchy-2.pp3
index 097ce2ebd..c1eb8c263 100644
--- a/rtdata/profiles/Punchy-2.pp3
+++ b/rtdata/profiles/Punchy-2.pp3
@@ -96,9 +96,10 @@ Threshold=25
 
 [Directional Pyramid Denoising]
 Enabled=false
-Luma=5
-Chroma=5
-Gamma=2
+Luma=0
+Ldetail=80
+Chroma=15
+Gamma=1.7
 
 [EPD]
 Enabled=false
diff --git a/rtengine/CMakeLists.txt b/rtengine/CMakeLists.txt
index fad3c3cbc..a58894da7 100644
--- a/rtengine/CMakeLists.txt
+++ b/rtengine/CMakeLists.txt
@@ -1,9 +1,9 @@
 include_directories (${EXTRA_INCDIR} ${GTHREAD_INCLUDE_DIRS} ${GOBJECT_INCLUDE_DIRS} ${GLIB2_INCLUDE_DIRS}
-    ${GLIBMM_INCLUDE_DIRS} ${IPTCDATA_INCLUDE_DIRS} ${LCMS_INCLUDE_DIRS} ${EXPAT_INCLUDE_DIRS}
+    ${GLIBMM_INCLUDE_DIRS} ${IPTCDATA_INCLUDE_DIRS} ${LCMS_INCLUDE_DIRS} ${EXPAT_INCLUDE_DIRS} ${FFTW3F_INCLUDE_DIRS}
     ${GTKMM_INCLUDE_DIRS} ${GTK_INCLUDE_DIRS})
 link_directories ("${PROJECT_SOURCE_DIR}/rtexif" ${EXTRA_LIBDIR} ${GTHREAD_LIBRARY_DIRS}
     ${GOBJECT_LIBRARY_DIRS} ${GLIB2_LIBRARY_DIRS} ${GLIBMM_LIBRARY_DIRS}
-    ${IPTCDATA_LIBRARY_DIRS} ${LCMS_LIBRARY_DIRS} ${EXPAT_LIBRARY_DIRS})
+    ${IPTCDATA_LIBRARY_DIRS} ${LCMS_LIBRARY_DIRS} ${EXPAT_LIBRARY_DIRS} ${FFTW3F_LIBRARY_DIRS})
 
 set (RTENGINESOURCEFILES safegtk.cc colortemp.cc curves.cc flatcurves.cc diagonalcurves.cc dcraw.cc iccstore.cc color.cc
     dfmanager.cc ffmanager.cc rawimage.cc image8.cc image16.cc imagefloat.cc imagedata.cc imageio.cc improcfun.cc init.cc dcrop.cc
@@ -12,8 +12,9 @@ set (RTENGINESOURCEFILES safegtk.cc colortemp.cc curves.cc flatcurves.cc diagona
     processingjob.cc rtthumbnail.cc utils.cc labimage.cc slicer.cc
     iplab2rgb.cc ipsharpen.cc iptransform.cc ipresize.cc ipvibrance.cc
 	jpeg_memsrc.cc jdatasrc.cc
-	PF_correct_RT.cc
-    dirpyrLab_denoise.cc dirpyrLab_equalizer.cc dirpyr_equalizer.cc
+	EdgePreserveLab.cc EdgePreservingDecomposition.cc cplx_wavelet_dec.cc FTblockDN.cc
+	PF_correct_RT.cc 
+    dirpyr_equalizer.cc
     calc_distort.cc lcp.cc dcp.cc
     klt/convolve.cc klt/error.cc klt/klt.cc klt/klt_util.cc klt/pnmio.cc klt/pyramid.cc klt/selectGoodFeatures.cc
     klt/storeFeatures.cc klt/trackFeatures.cc klt/writeFeatures.cc
@@ -30,5 +31,5 @@ ENDIF (BUILD_SHARED_LIBS)
 set_target_properties (rtengine PROPERTIES COMPILE_FLAGS "${RTENGINE_CXX_FLAGS}")
 
 target_link_libraries (rtengine rtexif ${EXTRA_LIB} ${GOBJECT_LIBRARIES} ${GTHREAD_LIBRARIES}
-    ${GLIB2_LIBRARIES} ${GLIBMM_LIBRARIES} ${LCMS_LIBRARIES} ${EXPAT_LIBRARIES} ${IPTCDATA_LIBRARIES}
+    ${GLIB2_LIBRARIES} ${GLIBMM_LIBRARIES} ${LCMS_LIBRARIES} ${EXPAT_LIBRARIES} ${FFTW3F_LIBRARIES} ${IPTCDATA_LIBRARIES}
     ${JPEG_LIBRARIES} ${PNG_LIBRARIES} ${TIFF_LIBRARIES} ${ZLIB_LIBRARIES})
diff --git a/rtengine/EdgePreserveLab.cc b/rtengine/EdgePreserveLab.cc
new file mode 100644
index 000000000..967555a94
--- /dev/null
+++ b/rtengine/EdgePreserveLab.cc
@@ -0,0 +1,229 @@
+#include "EdgePreserveLab.h"
+#include "boxblur.h"
+#include <cstdlib>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+//#define MAX(a,b) ((a)<(b)?(b):(a))
+//#define MIN(a,b) ((a)>(b)?(b):(a))
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+EdgePreserveLab::EdgePreserveLab(unsigned int width, unsigned int height){
+	w = width;
+	h = height;
+	n = w*h;
+
+	//Initialize the matrix just once at construction.
+	A = new MultiDiagonalSymmetricMatrix(n, 5);
+	if(!(
+		A->CreateDiagonal(0, 0) &&
+		A->CreateDiagonal(1, 1) &&
+		A->CreateDiagonal(2, w - 1) &&
+		A->CreateDiagonal(3, w) &&
+		A->CreateDiagonal(4, w + 1))){
+		delete A;
+		A = NULL;
+		printf("Error in EdgePreserveLab construction: out of memory.\n");
+	}else{
+		a0    = A->Diagonals[0];
+		a_1   = A->Diagonals[1];
+		a_w1  = A->Diagonals[2];
+		a_w   = A->Diagonals[3];
+		a_w_1 = A->Diagonals[4];
+	}
+}
+
+EdgePreserveLab::~EdgePreserveLab(){
+	delete A;
+}
+
+float *EdgePreserveLab::CreateBlur(float *Source, float LScale, float abScale, float EdgeStoppingLuma, float EdgeStoppingChroma, unsigned int Iterates, float *Blur, bool UseBlurForEdgeStop){
+	if(Blur == NULL)
+		UseBlurForEdgeStop = false,	//Use source if there's no supplied Blur.
+		Blur = new float[3*n];
+	if(LScale == 0.0f){
+		memcpy(Blur, Source, 3*n*sizeof(float));
+		return Blur;
+	}
+
+	//Create the edge stopping function a, rotationally symmetric and just one instead of (ax, ay). Maybe don't need Blur yet, so use its memory.
+	float *a, *b, *g;
+	if(UseBlurForEdgeStop) a = new float[n], g = Blur;
+	else a = Blur, g = Source;
+	//b = new float[n];
+	
+	unsigned int x, y, i;
+	unsigned int w1 = w - 1, h1 = h - 1;
+	float eps = 0.0001f;
+	float scL = powf(100.0f,LScale);
+	float scab = powf(200.0f,abScale);
+	
+	float * var = new float[w*h];
+	rtengine::boxvar(g, var, 1, 1, w, h);
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for(y = 0; y < h1; y++){
+		float *rg = &g[w*y];
+		for(x = 0; x < w1; x++){
+			//Estimate the central difference gradient in the center of a four pixel square. (gx, gy) is actually 2*gradient.
+			/*float gx = (fabs((rg[x + 1] - rg[x]) + (rg[x + w + 1] - rg[x + w])));
+			float gy = (fabs((rg[x + w] - rg[x]) + (rg[x + w + 1] - rg[x + 1])));
+						
+			//TODO: combine this with gx, gy if not needing separate quantities
+			float hx =  (fabs((rg[x + 1 + n] - rg[x + n]) + (rg[x + w + 1 + n] - rg[x + w + n])) + \
+						 fabs((rg[x + 1 + 2*n] - rg[x + 2*n]) + (rg[x + w + 1 + 2*n] - rg[x + w + 2*n]))); 
+			float hy = (fabs((rg[x + w + n] - rg[x + n]) + (rg[x + w + 1 + n] - rg[x + 1 + n])) + \
+						fabs((rg[x + w + 2*n] - rg[x + 2*n]) + (rg[x + w + 1 + 2*n] - rg[x + 1 + 2*n])));
+			*/
+			//float gradtot = (gx+gy+hx+hy);
+			//gradhisto[MAX(0,MIN(32767,(int)gradtot))] ++;
+			
+			//Apply power to the magnitude of the gradient to get the edge stopping function.
+			//a[x + w*y] = scL*expf(-100.0f*(gx + gy + hx + hy)/(EdgeStoppingLuma));
+			//a[x + w*y] = scL*expf(-var[y*w+x]/SQR(0.02*EdgeStoppingLuma));///(0.1+rg[x]);
+			a[x + w*y] = scL*expf(-50.0f*sqrt(var[y*w+x])/(EdgeStoppingLuma+eps));///(0.1+rg[x]);
+
+			//b[x + w*y] = (scab)*expf(-20.0f*(gx + gy + Lave*(hx + hy))/(EdgeStoppingChroma));
+			//b[x + w*y] = (scab)*expf(-400.0f*SQR(gx + gy + Lave*(hx + hy))/SQR(EdgeStoppingChroma));;
+
+		}
+	}
+
+	/* Now setup the linear problem. I use the Maxima CAS, here's code for making an FEM formulation for the smoothness term:
+		p(x, y) := (1 - x)*(1 - y);
+		P(m, n) := A[m][n]*p(x, y) + A[m + 1][n]*p(1 - x, y) + A[m + 1][n + 1]*p(1 - x, 1 - y) + A[m][n + 1]*p(x, 1 - y);
+		Integrate(f) := integrate(integrate(f, x, 0, 1), y, 0, 1);
+
+		Integrate(diff(P(u, v), x)*diff(p(x, y), x) + diff(P(u, v), y)*diff(p(x, y), y));
+		Integrate(diff(P(u - 1, v), x)*diff(p(1 - x, y), x) + diff(P(u - 1, v), y)*diff(p(1 - x, y), y));
+		Integrate(diff(P(u - 1, v - 1), x)*diff(p(1 - x, 1 - y), x) + diff(P(u - 1, v - 1), y)*diff(p(1 - x, 1 - y), y));
+		Integrate(diff(P(u, v - 1), x)*diff(p(x, 1 - y), x) + diff(P(u, v - 1), y)*diff(p(x, 1 - y), y));
+	So yeah. Use the numeric results of that to fill the matrix A.*/
+	memset(a_1, 0, A->DiagonalLength(1)*sizeof(float));
+	memset(a_w1, 0, A->DiagonalLength(w - 1)*sizeof(float));
+	memset(a_w, 0, A->DiagonalLength(w)*sizeof(float));
+	memset(a_w_1, 0, A->DiagonalLength(w + 1)*sizeof(float));
+
+//TODO: OMP here?
+	for(i = y = 0; y != h; y++){
+		for(x = 0; x != w; x++, i++){
+			float ac;
+			a0[i] = 1.0;
+
+			//Remember, only fill the lower triangle. Memory for upper is never made. It's symmetric. Trust.
+			if(x > 0 && y > 0)
+				ac = a[i - w - 1]/6.0f,
+				a_w_1[i - w - 1] -= 2.0f*ac, a_w[i - w] -= ac,
+				a_1[i - 1]       -=      ac, a0[i] += 4.0f*ac;
+
+			if(x < w1 && y > 0)
+				ac = a[i - w]/6.0f,
+				a_w[i - w] -= ac, a_w1[i - w + 1] -= 2.0f*ac,
+				a0[i] += 4.0f*ac;
+
+			if(x > 0 && y < h1)
+				ac = a[i - 1]/6.0f,
+				a_1[i - 1] -= ac, a0[i] += 4.0f*ac;
+
+			if(x < w1 && y < h1)
+				a0[i] += 4.0f*a[i]/6.0f;
+		}
+	}
+	if(UseBlurForEdgeStop) delete[] a;
+	
+
+	//Solve & return.
+	A->CreateIncompleteCholeskyFactorization(1);	//Fill-in of 1 seems to work really good. More doesn't really help and less hurts (slightly).
+	if(!UseBlurForEdgeStop) memcpy(Blur, Source, 3*n*sizeof(float));
+	// blur L channel
+	SparseConjugateGradient(A->PassThroughVectorProduct, Source, n, false, Blur, 0.0f, (void *)A, Iterates, A->PassThroughCholeskyBackSolve);
+	
+	//reset A for ab channels
+	/*memset(a_1, 0, A->DiagonalLength(1)*sizeof(float));
+	memset(a_w1, 0, A->DiagonalLength(w - 1)*sizeof(float));
+	memset(a_w, 0, A->DiagonalLength(w)*sizeof(float));
+	memset(a_w_1, 0, A->DiagonalLength(w + 1)*sizeof(float));
+	for(i = y = 0; y != h; y++){
+		for(x = 0; x != w; x++, i++){
+			float ac;
+			a0[i] = 1.0;
+			
+			//Remember, only fill the lower triangle. Memory for upper is never made. It's symmetric. Trust.
+			if(x > 0 && y > 0)
+				ac = b[i - w - 1]/6.0f,
+				a_w_1[i - w - 1] -= 2.0f*ac, a_w[i - w] -= ac,
+				a_1[i - 1]       -=      ac, a0[i] += 4.0f*ac;
+			
+			if(x < w1 && y > 0)
+				ac = b[i - w]/6.0f,
+				a_w[i - w] -= ac, a_w1[i - w + 1] -= 2.0f*ac,
+				a0[i] += 4.0f*ac;
+			
+			if(x > 0 && y < h1)
+				ac = b[i - 1]/6.0f,
+				a_1[i - 1] -= ac, a0[i] += 4.0f*ac;
+			
+			if(x < w1 && y < h1)
+				a0[i] += 4.0f*b[i]/6.0f;
+		}
+	}*/
+	/*if(UseBlurForEdgeStop)*/ //delete[] b;
+	
+	// blur ab channels
+	//A->CreateIncompleteCholeskyFactorization(1);	//Fill-in of 1 seems to work really good. More doesn't really help and less hurts (slightly).
+	//SparseConjugateGradient(A->PassThroughVectorProduct, Source+n, n, false, Blur+n, 0.0f, (void *)A, Iterates, A->PassThroughCholeskyBackSolve);
+	//SparseConjugateGradient(A->PassThroughVectorProduct, Source+2*n, n, false, Blur+2*n, 0.0f, (void *)A, Iterates, A->PassThroughCholeskyBackSolve);
+	
+	A->KillIncompleteCholeskyFactorization();
+	return Blur;
+}
+
+float *EdgePreserveLab::CreateIteratedBlur(float *Source, float LScale, float abScale, float EdgeStoppingLuma, float EdgeStoppingChroma, unsigned int Iterates, unsigned int Reweightings, float *Blur){
+	//Simpler outcome?
+	if(Reweightings == 0) return CreateBlur(Source, LScale, abScale, EdgeStoppingLuma, EdgeStoppingChroma, Iterates, Blur);
+
+	//Create a blur here, initialize.
+	if(Blur == NULL) Blur = new float[3*n];
+	memcpy(Blur, Source, 3*n*sizeof(float));
+
+	//Iteratively improve the blur.
+	Reweightings++;
+	for(unsigned int i = 0; i != Reweightings; i++)
+		CreateBlur(Source, LScale, abScale, EdgeStoppingLuma, EdgeStoppingChroma, Iterates, Blur, true);
+
+	return Blur;
+}
+
+float *EdgePreserveLab::CompressDynamicRange(float *Source, float LScale, float abScale, float EdgeStoppingLuma, float EdgeStoppingChroma, float CompressionExponent, float DetailBoost, unsigned int Iterates, unsigned int Reweightings, float *Compressed){
+	//We're working with luminance, which does better logarithmic.
+	unsigned int i;
+	//for(i = 0; i != n; i++)
+	//	Source[i] = logf(Source[i] + 0.0001f);
+	
+	//Blur. Also setup memory for Compressed (we can just use u since each element of u is used in one calculation).
+	float *u = CreateIteratedBlur(Source, LScale, abScale, EdgeStoppingLuma, EdgeStoppingChroma, Iterates, Reweightings);
+	if(Compressed == NULL) Compressed = u;
+	
+	//Apply compression, detail boost, unlogging. Compression is done on the logged data and detail boost on unlogged.
+	for(i = 0; i != n; i++){
+		//float ce = expf(Source[i] + u[i]*(CompressionExponent - 1.0f)) - 0.0001f;
+		//float ue = expf(u[i]) - 0.0001f;
+		//Source[i] = expf(Source[i]) - 0.0001f;
+		//Compressed[i] = ce + DetailBoost*(Source[i] - ue);
+		Compressed[i] = u[i];//ue;//for testing, to display blur
+	}
+	
+	if(Compressed != u) delete[] u;
+	return Compressed;
+}
+
+
+
+
diff --git a/rtengine/EdgePreserveLab.h b/rtengine/EdgePreserveLab.h
new file mode 100644
index 000000000..7da795247
--- /dev/null
+++ b/rtengine/EdgePreserveLab.h
@@ -0,0 +1,76 @@
+#pragma once
+/*
+The EdgePreserveLab files contain standard C++ (standard except the first line) code for creating and, to a
+limited extent (create your own uses!), messing with multi scale edge preserving decompositions of a 32 bit single channel
+image. As a byproduct it contains a lot of linear algebra which can be useful for optimization problems that
+you want to solve in rectangles on rectangular grids.
+
+Anyway. Basically, this is an implementation of what's presented in the following papers:
+	Edge-Preserving Decompositions for Multi-Scale Tone and Detail Manipulation
+	An Iterative Solution Method for Linear Systems of Which the Coefficient Matrix is a Symetric M-Matrix
+	Color correction for tone mapping
+	Wikipedia, the free encyclopedia
+
+First one is most of what matters, next two are details, last everything else. I did a few things differently, especially:
+	Reformulated the minimization with finite elements instead of finite differences. This results in better conditioning,
+	slightly better accuracy (less artifacts), the possibility of a better picked edge stopping function, but more memory consumption.
+
+	A single rotationally invariant edge stopping function is used instead of two non-invariant ones.
+
+	Incomplete Cholseky factorization instead of Szeliski's LAHBF. Slower, but not subject to any patents.
+
+	For tone mapping, original images are decomposed instead of their logarithms, and just one decomposition is made;
+	I find that this way works plenty good (theirs isn't better or worse... just different) and is simpler.
+
+Written by ben_pcc in Portland, Oregon, USA; code modified by Emil Martinec.
+
+ // EdgePreserveLab.h and EdgePreserveLab.cpp are 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.
+ //
+ //	This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+
+
+
+#include <cmath>
+#include <stdio.h>
+#include <string.h>
+#include "EdgePreservingDecomposition.h"
+
+class EdgePreserveLab{
+public:
+	EdgePreserveLab(unsigned int width, unsigned int height);
+	~EdgePreserveLab();
+
+	//Create an edge preserving blur of Source. Will create and return, or fill into Blur if not NULL. In place not ok.
+	//If UseBlurForEdgeStop is true, supplied not NULL Blur is used to calculate the edge stopping function instead of Source.
+	float *CreateBlur(float *Source, float LScale, float abScale, float EdgeStoppingLuma, float EdgeStoppingChroma, unsigned int Iterates, float *Blur = NULL, bool UseBlurForEdgeStop = false);
+
+	//Iterates CreateBlur such that the smoothness term approaches a specific norm via iteratively reweighted least squares. In place not ok.
+	float *CreateIteratedBlur(float *Source, float LScale, float abScale, float EdgeStoppingLuma, float EdgeStoppingChroma, unsigned int Iterates, unsigned int Reweightings, float *Blur = NULL);
+
+	/*Lowers global contrast while preserving or boosting local contrast. Can fill into Compressed. The smaller Compression
+	the more compression is applied, with Compression = 1 giving no effect and above 1 the opposite effect. You can totally
+	use Compression = 1 and play with DetailBoost for some really sweet unsharp masking. If working on luma/grey, consider giving it a logarithm.	
+	In place calculation to save memory (Source == Compressed) is totally ok. Reweightings > 0 invokes CreateIteratedBlur instead of CreateBlur. */
+	float *CompressDynamicRange(float *Source, float LScale = 1.0f, float abScale = 5.0f, float EdgeStoppingLuma = 1.0f, float EdgeStoppingChroma = 1.0f, 
+								float CompressionExponent = 0.8f, float DetailBoost = 0.1f, unsigned int Iterates = 20, 
+								unsigned int Reweightings = 0, float *Compressed = NULL);
+
+private:
+	MultiDiagonalSymmetricMatrix *A;	//The equations are simple enough to not mandate a matrix class, but fast solution NEEDS a complicated preconditioner.
+	unsigned int w, h, n;
+
+	//Convenient access to the data in A.
+	float *a0, *a_1, *a_w, *a_w_1, *a_w1;
+};
diff --git a/rtengine/EdgePreservingDecomposition.cc b/rtengine/EdgePreservingDecomposition.cc
index 339ec4b76..5a69daac1 100644
--- a/rtengine/EdgePreservingDecomposition.cc
+++ b/rtengine/EdgePreservingDecomposition.cc
@@ -185,14 +185,14 @@ bool MultiDiagonalSymmetricMatrix::CreateIncompleteCholeskyFactorization(unsigne
 	MaxFillAbove++;	//Conceptually, now "fill" includes an existing diagonal. Simpler in the math that follows.
 	unsigned int i, j, mic;
 	for(mic = i = 1; i != m; i++)
-		mic += min(StartRows[i] - StartRows[i - 1], MaxFillAbove);	//Guarunteed positive since StartRows must be created in increasing order.
+		mic += rtengine::min(StartRows[i] - StartRows[i - 1], MaxFillAbove);	//Guarunteed positive since StartRows must be created in increasing order.
 
 	//Initialize the decomposition - setup memory, start rows, etc.
 	MultiDiagonalSymmetricMatrix *ic = new MultiDiagonalSymmetricMatrix(n, mic);
 	ic->CreateDiagonal(0, 0);	//There's always a main diagonal in this type of decomposition.
 	for(mic = i = 1; i != m; i++){
 		//Set j to the number of diagonals to be created corresponding to a diagonal on this source matrix...
-		j = min(StartRows[i] - StartRows[i - 1], MaxFillAbove);
+		j = rtengine::min(StartRows[i] - StartRows[i - 1], MaxFillAbove);
 
 		//...and create those diagonals. I want to take a moment to tell you about how much I love minimalistic loops: very much.
 		while(j-- != 0)
diff --git a/rtengine/FTblockDN.cc b/rtengine/FTblockDN.cc
new file mode 100644
index 000000000..4c417fd44
--- /dev/null
+++ b/rtengine/FTblockDN.cc
@@ -0,0 +1,912 @@
+////////////////////////////////////////////////////////////////
+//
+//			CFA denoise by wavelet transform, FT filtering
+//
+//	copyright (c) 2008-2012  Emil Martinec <ejmartin@uchicago.edu>
+//
+//
+//  code dated: March 9, 2012
+//
+//	FTblockDN.cc 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.
+//
+//	This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.
+//
+////////////////////////////////////////////////////////////////
+
+
+
+
+#include <math.h>
+#include <fftw3.h>
+
+//#include "bilateral2.h"
+#include "gauss.h"
+
+#include "rtengine.h"
+#include "improcfun.h"
+#include "LUT.h"
+#include "array2D.h"
+#include "iccmatrices.h"
+#include "boxblur.h"
+#include "rt_math.h"
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#include "cplx_wavelet_dec.h"
+
+//#define MIN(a,b) ((a) < (b) ? (a) : (b))
+//#define MAX(a,b) ((a) > (b) ? (a) : (b))
+//#define LIM(x,min,max) MAX(min,MIN(x,max))
+//#define ULIM(x,y,z) ((y) < (z) ? LIM(x,y,z) : LIM(x,z,y))
+//#define CLIP(x) LIM(x,0,65535)
+
+#define TS 64		// Tile size
+#define offset 25	// shift between tiles
+#define fTS ((TS/2+1))	// second dimension of Fourier tiles
+#define blkrad 1	// radius of block averaging
+
+#define epsilon 0.001f/(TS*TS) //tolerance
+
+namespace rtengine {
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	/*
+	 Structure of the algorithm: 
+	 
+	 1. Compute an initial denoise of the image via undecimated wavelet transform 
+	 and universal thresholding modulated by user input.
+	 2. Decompose the residual image into TSxTS size tiles, shifting by 'offset' each step 
+	 (so roughly each pixel is in (TS/offset)^2 tiles); Discrete Cosine transform the tiles.
+	 3. Filter the DCT data to pick out patterns missed by the wavelet denoise
+	 4. Inverse DCT the denoised tile data and combine the tiles into a denoised output image.
+	 
+	 */
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	
+	
+	
+	
+	void ImProcFunctions::RGB_denoise(Imagefloat * src, Imagefloat * dst, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe)
+	{	
+		
+		//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+		
+		/*if (plistener) {
+		 plistener->setProgressStr ("Denoise...");
+		 plistener->setProgress (0.0);
+		 }*/
+		
+		volatile double progress = 0.0;
+		
+		//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+		
+		const short int imheight=src->height, imwidth=src->width;
+		
+		if (dnparams.luma==0 && dnparams.chroma==0) {
+            //nothing to do; copy src to dst
+            memcpy(dst->data,src->data,dst->width*dst->height*3*sizeof(float));
+			return;
+		}
+		
+		//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+		// gamma transform for input data
+		float gam = dnparams.gamma;
+		float gamthresh = 0.001;
+		float gamslope = exp(log((double)gamthresh)/gam)/gamthresh;
+		
+		LUTf gamcurve(65536,0);
+		
+		for (int i=0; i<65536; i++) {
+			gamcurve[i] = (Color::gamma((double)i/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)) * 32768.0f;
+		}
+		
+		// inverse gamma transform for output data
+		float igam = 1/gam;
+		float igamthresh = gamthresh*gamslope;
+		float igamslope = 1/gamslope;
+		
+		LUTf igamcurve(65536,0);
+		
+		for (int i=0; i<65536; i++) {
+			igamcurve[i] = (Color::gamma((float)i/32768.0f, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0f);
+		}
+		//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+		
+		
+		//srand((unsigned)time(0));//test with random data
+		
+		const float gain = pow (2.0, dnparams.expcomp);
+		
+		float noisevar_Ldetail = SQR((SQR(100-dnparams.Ldetail) + 50*(100-dnparams.Ldetail)) * TS * 0.5f);
+		
+		
+		array2D<float> tilemask_in(TS,TS);
+		array2D<float> tilemask_out(TS,TS);
+		
+		const int border = MAX(2,TS/16);
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int i=0; i<TS; i++) {
+			float i1 = abs((i>TS/2 ? i-TS+1 : i));
+			float vmask = (i1<border ? SQR(sin((M_PI*i1)/(2*border))) : 1.0f);
+			float vmask2 = (i1<2*border ? SQR(sin((M_PI*i1)/(2*border))) : 1.0f);
+			for (int j=0; j<TS; j++) {
+				float j1 = abs((j>TS/2 ? j-TS+1 : j));
+				tilemask_in[i][j] = (vmask * (j1<border ? SQR(sin((M_PI*j1)/(2*border))) : 1.0f)) + epsilon;
+				
+				tilemask_out[i][j] = (vmask2 * (j1<2*border ? SQR(sin((M_PI*j1)/(2*border))) : 1.0f)) + epsilon;
+				
+			}
+		}
+		
+		//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+		//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+		// begin tile processing of image
+		
+		//output buffer
+		Imagefloat * dsttmp = new Imagefloat(imwidth,imheight);
+		for (int n=0; n<3*imwidth*imheight; n++) dsttmp->data[n] = 0;
+		
+		const int tilesize = 1024;
+		const int overlap = 128;
+
+		int numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip;
+
+		if (imwidth<tilesize) {
+			numtiles_W = 1;
+			tileWskip = imwidth;
+			tilewidth = imwidth;
+		} else {
+			numtiles_W = ceil(((float)(imwidth))/(tilesize-overlap));
+			tilewidth  = ceil(((float)(imwidth))/(numtiles_W))+overlap;
+			tilewidth += (tilewidth&1);
+			tileWskip = tilewidth-overlap;
+		}
+		if (imheight<tilesize) {
+			numtiles_H = 1;
+			tileHskip = imheight;
+			tileheight = imheight;
+		} else {
+			numtiles_H = ceil(((float)(imheight))/(tilesize-overlap));
+			tileheight = ceil(((float)(imheight))/(numtiles_H))+overlap;
+			tileheight += (tileheight&1);
+			tileHskip = tileheight-overlap;
+		}
+		//now we have tile dimensions, overlaps
+		//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+//adding omp here slows it down
+		for (int tiletop=0; tiletop<imheight; tiletop+=tileHskip) {
+			for (int tileleft=0; tileleft<imwidth; tileleft+=tileWskip) {
+				
+				int tileright = MIN(imwidth,tileleft+tilewidth);
+				int tilebottom = MIN(imheight,tiletop+tileheight);
+				int width  = tileright-tileleft;
+				int height = tilebottom-tiletop;
+				
+				//input L channel
+				array2D<float> Lin(width,height,ARRAY2D_CLEAR_DATA);
+				//wavelet denoised image
+				LabImage * labdn = new LabImage(width,height);
+				//residual between input and denoised L channel
+				array2D<float> Ldetail(width,height,ARRAY2D_CLEAR_DATA);
+				//pixel weight
+				array2D<float> totwt(width,height,ARRAY2D_CLEAR_DATA);//weight for combining DCT blocks
+
+//#ifdef _OPENMP
+//#pragma omp parallel for
+//#endif
+//TODO: implement using AlignedBufferMP
+				//fill tile from image; convert RGB to "luma/chroma"
+				for (int i=tiletop/*, i1=0*/; i<tilebottom; i++/*, i1++*/) {
+					int i1 = i - tiletop;
+					for (int j=tileleft/*, j1=0*/; j<tileright; j++/*, j1++*/) {
+						int j1 = j - tileleft;
+						
+						float X = gain*src->r[i][j];//xyz_prophoto[0][0]*src->r[i][j] + xyz_prophoto[0][1]*src->g[i][j] + xyz_prophoto[0][2]*src->b[i][j];
+						float Y = gain*src->g[i][j];//xyz_prophoto[1][0]*src->r[i][j] + xyz_prophoto[1][1]*src->g[i][j] + xyz_prophoto[1][2]*src->b[i][j];
+						float Z = gain*src->b[i][j];//xyz_prophoto[2][0]*src->r[i][j] + xyz_prophoto[2][1]*src->g[i][j] + xyz_prophoto[2][2]*src->b[i][j];
+						
+						X = X<65535.0f ? gamcurve[X] : (Color::gamma((double)X/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)*32768.0f);
+						Y = Y<65535.0f ? gamcurve[Y] : (Color::gamma((double)Y/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)*32768.0f);
+						Z = Z<65535.0f ? gamcurve[Z] : (Color::gamma((double)Z/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)*32768.0f);
+						
+						labdn->L[i1][j1] = Y;
+						labdn->a[i1][j1] = (X-Y);
+						labdn->b[i1][j1] = (Y-Z);
+						
+						Ldetail[i1][j1] = 0;
+						Lin[i1][j1] = Y;
+						totwt[i1][j1] = 0;
+					}
+				}
+				
+				//initial impulse denoise
+				if (dnparams.luma>0.01) {
+					impulse_nr (labdn, MIN(50.0f,dnparams.luma)/20.0f);
+				}
+				
+				int datalen = labdn->W * labdn->H;
+				
+				//now perform basic wavelet denoise
+				//last two arguments of wavelet decomposition are max number of wavelet decomposition levels; 
+				//and whether to subsample the image after wavelet filtering.  Subsampling is coded as 
+				//binary 1 or 0 for each level, eg subsampling = 0 means no subsampling, 1 means subsample 
+				//the first level only, 7 means subsample the first three levels, etc.
+				wavelet_decomposition Ldecomp(labdn->data, labdn->W, labdn->H, 5/*maxlevels*/, 0/*subsampling*/ );
+				wavelet_decomposition adecomp(labdn->data+datalen, labdn->W, labdn->H, 5, 1 );
+				wavelet_decomposition bdecomp(labdn->data+2*datalen, labdn->W, labdn->H, 5, 1 );
+				
+				float noisevarL	 = SQR((dnparams.luma/125.0f)*(1+ dnparams.luma/25.0f));
+				float noisevarab = SQR(dnparams.chroma/10.0f);
+				
+				//WaveletDenoiseAll_BiShrink(Ldecomp, adecomp, bdecomp, noisevarL, noisevarab);
+				WaveletDenoiseAll(Ldecomp, adecomp, bdecomp, noisevarL, noisevarab);
+
+				Ldecomp.reconstruct(labdn->data);
+				adecomp.reconstruct(labdn->data+datalen);
+				bdecomp.reconstruct(labdn->data+2*datalen);
+				
+				//TODO: at this point wavelet coefficients storage can be freed
+				
+				//second impulse denoise
+				if (dnparams.luma>0.01) {
+					impulse_nr (labdn, MIN(50.0f,dnparams.luma)/20.0f);
+				}				
+				//PF_correct_RT(dst, dst, defringe.radius, defringe.threshold);
+				
+				//wavelet denoised L channel
+				array2D<float> Lwavdn(width,height);
+				float * Lwavdnptr = Lwavdn;
+				memcpy (Lwavdnptr, labdn->data, width*height*sizeof(float));
+				
+				//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+				//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+				// now do detail recovery using block DCT to detect  
+				// patterns missed by wavelet denoise
+				// blocks are not the same thing as tiles!
+				
+				
+				// allocate DCT data structures
+				
+				// calculation for detail recovery blocks
+				const int numblox_W = ceil(((float)(width))/(offset))+2*blkrad;
+				const int numblox_H = ceil(((float)(height))/(offset))+2*blkrad;
+				//const int nrtiles = numblox_W*numblox_H;
+				// end of tiling calc
+				
+				//DCT block data storage
+				float * Lblox  = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float));
+				float * fLblox = (float *) fftwf_malloc (numblox_W*TS*TS * sizeof (float));
+				
+				
+				//make a plan for FFTW
+				fftwf_plan plan_forward_blox, plan_backward_blox;
+				
+				int nfwd[2]={TS,TS};
+				
+				//for DCT:
+				const fftw_r2r_kind fwdkind[2] = {FFTW_REDFT10, FFTW_REDFT10};
+				const fftw_r2r_kind bwdkind[2] = {FFTW_REDFT01, FFTW_REDFT01};
+				
+				plan_forward_blox  = fftwf_plan_many_r2r(2, nfwd, numblox_W, Lblox, NULL, 1, TS*TS, fLblox, NULL, 1, TS*TS, fwdkind, FFTW_ESTIMATE );
+				plan_backward_blox = fftwf_plan_many_r2r(2, nfwd, numblox_W, fLblox, NULL, 1, TS*TS, Lblox, NULL, 1, TS*TS, bwdkind, FFTW_ESTIMATE );
+				
+				//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+				//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+				// Main detail recovery algorithm: Block loop
+//OpenMP here			
+//adding omp here leads to artifacts	
+				for (int vblk=0; vblk<numblox_H; vblk++) {
+					//printf("vblock=%d",vblk);
+					int vblkmod = vblk%8;
+					
+					int top = (vblk-blkrad)*offset;
+					
+					float * buffer = new float [width + TS + 2*blkrad*offset];
+					float * datarow = buffer+blkrad*offset;
+
+//#ifdef _OPENMP
+//#pragma omp parallel for
+//#endif
+//TODO: implement using AlignedBufferMP
+					for (int i=0/*, row=top*/; i<TS; i++/*, row++*/) {
+						int row = top + i;
+						int rr = row;
+						if (row<0) {
+							rr = MIN(-row,height-1);
+						} else if (row>=height) {
+							rr = MAX(0,2*height-2-row);
+						} 
+						
+						for (int j=0; j<labdn->W; j++) {
+							datarow[j] = (Lin[rr][j]-Lwavdn[rr][j]);
+						}
+						
+						for (int j=-blkrad*offset; j<0; j++) {
+							datarow[j] = datarow[MIN(-j,width-1)];
+						}
+						for (int j=width; j<width+TS+blkrad*offset; j++) {
+							datarow[j] = datarow[MAX(0,2*width-2-j)];
+						}//now we have a padded data row
+						
+						//now fill this row of the blocks with Lab high pass data
+//OMP here does not add speed, better handled on the outside loop
+						for (int hblk=0; hblk<numblox_W; hblk++) {
+							int left = (hblk-blkrad)*offset;
+							int indx = (hblk)*TS;//index of block in malloc
+							
+							for (int j=0; j<TS; j++) {
+								Lblox[(indx + i)*TS+j]  = tilemask_in[i][j]*datarow[left+j];// luma data
+								if (top+i>=0 && top+i<height && left+j>=0 && left+j<width) {
+									totwt[top+i][left+j] += tilemask_in[i][j]*tilemask_out[i][j];
+								}
+							}
+						}
+					}//end of filling block row
+					delete[] buffer;
+					
+					//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+					
+					//fftwf_print_plan (plan_forward_blox);
+					fftwf_execute_r2r(plan_forward_blox,Lblox,fLblox);		// DCT an entire row of tiles
+					
+					//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+					// now process the vblk row of blocks for noise reduction
+					for (int hblk=0; hblk<numblox_W; hblk++) {
+						
+						RGBtile_denoise (fLblox, vblk, hblk, numblox_H, numblox_W, noisevar_Ldetail );
+						
+					}//end of horizontal block loop
+					
+					//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+					
+					//now perform inverse FT of an entire row of blocks
+					fftwf_execute_r2r(plan_backward_blox,fLblox,Lblox);	//for DCT
+					
+					int topproc = (vblk-blkrad)*offset;
+					
+					//add row of blocks to output image tile
+					RGBoutput_tile_row (Lblox, Ldetail, tilemask_out, height, width, topproc );
+					
+					//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+					
+				}//end of vertical block loop
+				
+				//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+				
+				// clean up
+				//#pragma omp single nowait
+				fftwf_destroy_plan( plan_forward_blox );
+				//#pragma omp single nowait
+				fftwf_destroy_plan( plan_backward_blox );
+								
+				fftwf_free ( Lblox);
+				fftwf_free ( fLblox);
+				
+				fftwf_cleanup();
+				
+				//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+				for (int i=0; i<height; i++) {
+					for (int j=0; j<width; j++) {
+						//may want to include masking threshold for large hipass data to preserve edges/detail
+						float hpdn = Ldetail[i][j]/totwt[i][j];//note that labdn initially stores the denoised hipass data
+						
+						labdn->L[i][j] = Lwavdn[i][j] + hpdn;
+						
+					}
+				}
+				
+				//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+				// transform denoised "Lab" to output RGB
+				
+				//calculate mask for feathering output tile overlaps
+				float * Vmask = new float [height];
+				float * Hmask = new float [width];
+
+				for (int i=0; i<height; i++) {
+					Vmask[i] = 1;
+				}
+				for (int j=0; j<width; j++) {
+					Hmask[j] = 1;
+				}
+				for (int i=0; i<overlap; i++) {
+					float mask = SQR(sin((M_PI*i)/(2*overlap)));
+					if (tiletop>0) Vmask[i] = mask;
+					if (tilebottom<imheight) Vmask[height-1-i] = mask;
+					if (tileleft>0) Hmask[i] = mask;
+					if (tileright<imwidth) Hmask[width-1-i] = mask;
+				}
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+				//convert back to RGB and write to destination array
+				for (int i=tiletop/* i1=0*/; i<tilebottom; i++/*, i1++*/){
+					int i1 = i-tiletop;
+					float X,Y,Z;
+					for (int j=tileleft/*, j1=0*/; j<tileright; j++/*, j1++*/) {
+						int j1=j-tileleft;
+						
+						Y = labdn->L[i1][j1];
+						X = (labdn->a[i1][j1]) + Y;
+						Z = Y - (labdn->b[i1][j1]);
+						
+						X = X<32768.0f ? igamcurve[X] : (Color::gamma((float)X/32768.0f, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0f);
+						Y = Y<32768.0f ? igamcurve[Y] : (Color::gamma((float)Y/32768.0f, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0f);
+						Z = Z<32768.0f ? igamcurve[Z] : (Color::gamma((float)Z/32768.0f, igam, igamthresh, igamslope, 1.0, 0.0) * 65535.0f);
+						
+						//Y = 65535.0f*(0.05+0.1*((float)rand()/(float)RAND_MAX));//test with random data
+						
+						float factor = Vmask[i1]*Hmask[j1]/gain;
+						
+						dsttmp->r[i][j] += factor*X;//prophoto_xyz[0][0]*X + prophoto_xyz[0][1]*Y + prophoto_xyz[0][2]*Z;
+						dsttmp->g[i][j] += factor*Y;//prophoto_xyz[1][0]*X + prophoto_xyz[1][1]*Y + prophoto_xyz[1][2]*Z;
+						dsttmp->b[i][j] += factor*Z;//prophoto_xyz[2][0]*X + prophoto_xyz[2][1]*Y + prophoto_xyz[2][2]*Z;
+						
+					}
+				}
+				
+				//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+				
+				delete labdn;
+				delete[] Vmask;
+				delete[] Hmask;
+				
+			}//end of tile row
+		}//end of tile loop
+		
+	//copy denoised image to output
+	memcpy (dst->data, dsttmp->data, 3*dst->width*dst->height*sizeof(float));
+	
+	delete dsttmp;
+	
+}//end of main RGB_denoise
+
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	
+	
+	void ImProcFunctions::RGBtile_denoise (float * fLblox, int vblproc, int hblproc, int numblox_H, int numblox_W, float noisevar_Ldetail )	//for DCT
+	{			
+		float * nbrwt  = new float[TS*TS];	//for DCT
+		
+		int blkstart = hblproc*TS*TS;
+		
+		boxabsblur(fLblox+blkstart, nbrwt, 3, 3, TS, TS);//blur neighbor weights for more robust estimation	//for DCT
+
+		for (int n=0; n<TS*TS; n++) {		//for DCT
+			fLblox[blkstart+n] *= (1-expf(-SQR(nbrwt[n])/noisevar_Ldetail));
+		}//output neighbor averaged result
+		
+		delete[] nbrwt;
+		
+		//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+		//printf("vblk=%d  hlk=%d  wsqave=%f   ||   ",vblproc,hblproc,wsqave);
+		
+	}//end of function tile_denoise
+	
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	void ImProcFunctions::RGBoutput_tile_row (float *bloxrow_L, float ** Ldetail, float ** tilemask_out, int height, int width, int top )
+	{
+		const int numblox_W = ceil(((float)(width))/(offset));
+		const float DCTnorm = 1.0f/(4*TS*TS); //for DCT
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif			
+		//add row of tiles to output image
+		for (int hblk=0; hblk < numblox_W; hblk++) {
+			int left = (hblk-blkrad)*offset;
+			int bottom = MIN( top+TS,height);
+			int right  = MIN(left+TS, width);
+			int imin = MAX(0,-top);
+			int jmin = MAX(0,-left);
+			int imax = bottom - top;
+			int jmax = right - left;
+			
+			int indx = hblk*TS;
+			
+			for (int i=imin; i<imax; i++)
+				for (int j=jmin; j<jmax; j++) {
+					
+					Ldetail[top+i][left+j] += tilemask_out[i][j]*bloxrow_L[(indx + i)*TS+j]*DCTnorm; //for DCT
+					
+				}
+		}
+		
+	}
+	
+#undef TS
+#undef fTS
+#undef offset
+#undef epsilon
+	
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	float ImProcFunctions::MadMax(float * DataList, int & max, int datalen) {
+		
+		//computes Median Absolute Deviation and Maximum of DataList
+		//DataList values should mostly have abs val < 65535
+		
+		int * histo = new int[65536];
+		//memset(histo, 0, 65536*sizeof(histo));
+		for (int i=0; i<65536; i++) histo[i]=0;
+		
+		//calculate histogram of absolute values of HH wavelet coeffs
+		for (int i=0; i<datalen; i++) {
+			histo[MAX(0,MIN(65535,abs((int)DataList[i])))]++;
+		}
+		
+		//find median of histogram
+		int median=0, count=0;
+		while (count<datalen/2) {
+			count += histo[median];
+			median++;
+		}
+		
+		//find max of histogram
+		max=65535;
+		while (histo[max]==0) {
+			max--;
+		}
+		
+		int count_ = count - histo[median-1];
+		
+		delete[] histo;
+		
+		// interpolate 
+		return (( (median-1) + (datalen/2-count_)/((float)(count-count_)) )/0.6745);
+		
+	}
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	
+	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decomposition &WaveletCoeffs_L, wavelet_decomposition &WaveletCoeffs_a, 
+													 wavelet_decomposition &WaveletCoeffs_b, float noisevar_L, float noisevar_ab ) 
+	{
+		int maxlvl = WaveletCoeffs_L.maxlevel();
+		const float eps = 0.01f;
+		int max;
+		float parfrac = 0.05;
+		
+		float madL[8][3], mada[8][3], madb[8][3];
+
+//OpenMP here		
+		for (int lvl=0; lvl<maxlvl; lvl++) {
+			// compute median absolute deviation (MAD) of detail coefficients as robust noise estimator
+			
+			int Wlvl_L = WaveletCoeffs_L.level_W(lvl);
+			int Hlvl_L = WaveletCoeffs_L.level_H(lvl);
+			
+			int Wlvl_ab = WaveletCoeffs_a.level_W(lvl);
+			int Hlvl_ab = WaveletCoeffs_a.level_H(lvl);
+			
+			float ** WavCoeffs_L = WaveletCoeffs_L.level_coeffs(lvl);
+			float ** WavCoeffs_a = WaveletCoeffs_a.level_coeffs(lvl);
+			float ** WavCoeffs_b = WaveletCoeffs_b.level_coeffs(lvl);
+
+			for (int dir=1; dir<4; dir++) {
+				madL[lvl][dir-1] = SQR(MadMax(WavCoeffs_L[dir], max, Wlvl_L*Hlvl_L));
+				mada[lvl][dir-1] = SQR(MadMax(WavCoeffs_a[dir], max, Wlvl_ab*Hlvl_ab));
+				madb[lvl][dir-1] = SQR(MadMax(WavCoeffs_b[dir], max, Wlvl_ab*Hlvl_ab));
+			}
+		}
+		
+		for (int lvl=maxlvl-1; lvl>=0; lvl--) {//for levels less than max, use level diff to make edge mask
+			//for (int lvl=0; lvl<maxlvl; lvl++) {
+			
+			int Wlvl_L = WaveletCoeffs_L.level_W(lvl);
+			int Hlvl_L = WaveletCoeffs_L.level_H(lvl);
+			
+			int Wlvl_ab = WaveletCoeffs_a.level_W(lvl);
+			int Hlvl_ab = WaveletCoeffs_a.level_H(lvl);
+			
+			float skip_L = WaveletCoeffs_L.level_stride(lvl);
+			float skip_ab = WaveletCoeffs_a.level_stride(lvl);
+			
+			float ** WavCoeffs_L = WaveletCoeffs_L.level_coeffs(lvl);
+			float ** WavCoeffs_a = WaveletCoeffs_a.level_coeffs(lvl);
+			float ** WavCoeffs_b = WaveletCoeffs_b.level_coeffs(lvl);
+			
+			if (lvl==maxlvl-1) {
+				ShrinkAll(WavCoeffs_L, WavCoeffs_a, WavCoeffs_b, lvl, Wlvl_L, Hlvl_L, Wlvl_ab, Hlvl_ab,
+						  skip_L, skip_ab, noisevar_L, noisevar_ab);//TODO: this implies redundant evaluation of MAD
+			} else {
+				
+				float ** WavPars_L = WaveletCoeffs_L.level_coeffs(lvl+1);
+				//float ** WavPars_a = WaveletCoeffs_a.level_coeffs(lvl+1);
+				//float ** WavPars_b = WaveletCoeffs_b.level_coeffs(lvl+1);
+				
+				//simple wavelet shrinkage
+				float * sfave = new float[Wlvl_L*Hlvl_L]; 
+				array2D<float> edge(Wlvl_L,Hlvl_L);
+				AlignedBuffer<double>* buffer = new AlignedBuffer<double> (MAX(Wlvl_L,Hlvl_L));
+				
+				printf("\n level=%d  \n",lvl);
+				
+				for (int dir=1; dir<4; dir++) {
+					float mad_L = madL[lvl][dir-1];
+					float mad_a = noisevar_ab*mada[lvl][dir-1];
+					float mad_b = noisevar_ab*madb[lvl][dir-1];
+					//float mad_Lpar = madL[lvl+1][dir-1];
+					//float mad_apar = mada[lvl+1][dir-1];
+					//float mad_bpar = mada[lvl+1][dir-1];
+					
+					//float skip_ab_ratio = WaveletCoeffs_a.level_stride(lvl+1)/skip_ab;
+					float skip_L_ratio =  WaveletCoeffs_L.level_stride(lvl+1)/skip_L;
+					
+					if (noisevar_ab>0.01) {
+						
+						printf("  dir=%d  mad_L=%f		mad_a=%f		mad_b=%f	\n",dir,sqrt(mad_L),sqrt(mad_a),sqrt(mad_b));
+						
+//OpenMP here						
+						for (int i=0; i<Hlvl_ab; i++) {
+							for (int j=0; j<Wlvl_ab; j++) {
+								
+								int coeffloc_ab = i*Wlvl_ab+j;
+								//int coeffloc_abpar = (MAX(0,i-skip_ab)*Wlvl_ab+MAX(0,j-skip_ab))/skip_ab_ratio;
+								
+								int coeffloc_L	= ((i*skip_L)/skip_ab)*Wlvl_L + ((j*skip_L)/skip_ab);
+								
+								float mag_L = SQR(WavCoeffs_L[dir][coeffloc_L ])+eps;
+								float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab])+eps;
+								float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab])+eps;
+								
+								//float edgefactor = 1-exp(-mag_L/(9*mad_L));// * exp(-mag_a/(4*mad_a)) * exp(-mag_b/(4*mad_b));
+								
+								//float coeff_a = sqrt(SQR(WavCoeffs_a[dir][coeffloc_ab])/mad_a+SQR(parfrac*WavPars_a[dir][coeffloc_abpar])/mad_apar);
+								//float coeff_b = sqrt(SQR(WavCoeffs_b[dir][coeffloc_ab])/mad_b+SQR(parfrac*WavPars_b[dir][coeffloc_abpar])/mad_bpar);
+								
+								// 'firm' threshold of chroma coefficients
+								//WavCoeffs_a[dir][coeffloc_ab] *= edgefactor*(coeff_a>2 ? 1 : (coeff_a<1 ? 0 : (coeff_a - 1)));
+								//WavCoeffs_b[dir][coeffloc_ab] *= edgefactor*(coeff_b>2 ? 1 : (coeff_b<1 ? 0 : (coeff_b - 1)));
+								
+								//float satfactor_a = mad_a/(mad_a+0.5*SQR(WavCoeffs_a[0][coeffloc_ab]));
+								//float satfactor_b = mad_b/(mad_b+0.5*SQR(WavCoeffs_b[0][coeffloc_ab]));
+								
+								WavCoeffs_a[dir][coeffloc_ab] *= SQR(1-exp(-(mag_a/mad_a)-(mag_L/(9*mad_L)))/*satfactor_a*/);
+								WavCoeffs_b[dir][coeffloc_ab] *= SQR(1-exp(-(mag_b/mad_b)-(mag_L/(9*mad_L)))/*satfactor_b*/);
+								
+							}
+						}//now chrominance coefficients are denoised
+					}
+					
+					if (noisevar_L>0.01) {
+						mad_L *= noisevar_L*5/(lvl+1);
+//OpenMP here						
+						for (int i=0; i<Hlvl_L; i++) 
+							for (int j=0; j<Wlvl_L; j++) {
+								
+								int coeffloc_L = i*Wlvl_L+j;
+								int coeffloc_Lpar = (MAX(0,i-skip_L)*Wlvl_L+MAX(0,j-skip_L))/skip_L_ratio;
+								
+								float mag_L = SQR(WavCoeffs_L[dir][coeffloc_L]);
+								//float mag_Lpar = SQR(parfrac*WavPars_L[dir][coeffloc_Lpar]);
+								//float sf_L = SQR(1-expf(-(mag_L/mad_L)-(mag_Lpar/mad_L)));
+								float sf_L = mag_L/(mag_L+mad_L*exp(-mag_L/(9*mad_L))+eps);
+								
+								sfave[coeffloc_L] = sf_L;
+								
+								//edge[i][j] = (WavCoeffs_L[dir][coeffloc_L] - WavPars_L[dir][coeffloc_Lpar]);
+							}
+						
+						//blur edge measure
+						//gaussHorizontal<float> (edge, edge, buffer, Wlvl_L, Hlvl_L, 1<<(lvl+1), false /*multiThread*/);
+						//gaussVertical<float>   (edge, edge, buffer, Wlvl_L, Hlvl_L, 1<<(lvl+1), false);
+						
+						boxblur(sfave, sfave, lvl+2, lvl+2, Wlvl_L, Hlvl_L);//increase smoothness by locally averaging shrinkage
+//OpenMP here						
+						for (int i=0; i<Hlvl_L; i++) 
+							for (int j=0; j<Wlvl_L; j++) {
+								
+								int coeffloc_L = i*Wlvl_L+j;
+								
+								float mag_L = SQR(WavCoeffs_L[dir][coeffloc_L]);
+								//float sf_L = SQR(1-expf(-(mag_L/mad_L)-(mag_Lpar/mad_L)));
+								
+								float edgefactor = 1;//expf(-SQR(edge[i][j])/mad_L);
+								
+								float sf_L = mag_L/(mag_L + edgefactor*mad_L*exp(-mag_L/(9*mad_L))+eps);
+								
+								//use smoothed shrinkage unless local shrinkage is much less
+								WavCoeffs_L[dir][coeffloc_L] *= (SQR(edgefactor*sfave[coeffloc_L])+SQR(sf_L))/(edgefactor*sfave[coeffloc_L]+sf_L+eps);
+								
+							}//now luminance coeffs are denoised
+						
+					}
+					
+				}
+				delete[] sfave;
+				delete buffer;
+				
+			}
+		}
+	}
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	void ImProcFunctions::WaveletDenoiseAll(wavelet_decomposition &WaveletCoeffs_L, wavelet_decomposition &WaveletCoeffs_a, 
+											wavelet_decomposition &WaveletCoeffs_b, float noisevar_L, float noisevar_ab ) 
+	{
+		int maxlvl = WaveletCoeffs_L.maxlevel();
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int lvl=0; lvl<maxlvl; lvl++) {
+			
+			int Wlvl_L = WaveletCoeffs_L.level_W(lvl);
+			int Hlvl_L = WaveletCoeffs_L.level_H(lvl);
+			
+			int Wlvl_ab = WaveletCoeffs_a.level_W(lvl);
+			int Hlvl_ab = WaveletCoeffs_a.level_H(lvl);
+			
+			float skip_L = WaveletCoeffs_L.level_stride(lvl);
+			float skip_ab = WaveletCoeffs_a.level_stride(lvl);
+			
+			float ** WavCoeffs_L = WaveletCoeffs_L.level_coeffs(lvl);
+			float ** WavCoeffs_a = WaveletCoeffs_a.level_coeffs(lvl);
+			float ** WavCoeffs_b = WaveletCoeffs_b.level_coeffs(lvl);
+			
+			ShrinkAll(WavCoeffs_L, WavCoeffs_a, WavCoeffs_b, lvl, Wlvl_L, Hlvl_L, Wlvl_ab, Hlvl_ab,
+					  skip_L, skip_ab, noisevar_L, noisevar_ab);
+		}
+		
+	}
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCoeffs_a, float ** WavCoeffs_b, int level, 
+									int W_L, int H_L, int W_ab, int H_ab, int skip_L, int skip_ab, float noisevar_L, float noisevar_ab) 
+	{	
+		//simple wavelet shrinkage
+		const float eps = 0.01f;
+		float * sfave = new float[W_L*H_L];
+		float * sfavea = new float[W_L*H_L];
+		float * sfaveb = new float[W_L*H_L];
+
+		int max;
+		
+		printf("\n level=%d  \n",level);
+	
+		for (int dir=1; dir<4; dir++) {
+			float madL = SQR(MadMax(WavCoeffs_L[dir], max, W_L*H_L));	
+			float mada = SQR(MadMax(WavCoeffs_a[dir], max, W_ab*H_ab));
+			float madb = SQR(MadMax(WavCoeffs_b[dir], max, W_ab*H_ab));
+			
+			
+			printf("  dir=%d  mad_L=%f	mad_a=%f	mad_b=%f	\n",dir,sqrt(madL),sqrt(mada),sqrt(madb));
+
+			float mad_L = madL*noisevar_L*5/(level+1);
+			float mad_a = mada*noisevar_ab;
+			float mad_b = madb*noisevar_ab;
+			
+			if (noisevar_ab>0.01) {
+//OpenMP here
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif				
+				for (int i=0; i<H_ab; i++) {
+					for (int j=0; j<W_ab; j++) {
+						
+						int coeffloc_ab = i*W_ab+j;
+						int coeffloc_L	= ((i*skip_L)/skip_ab)*W_L + ((j*skip_L)/skip_ab);
+						
+						float mag_L = SQR(WavCoeffs_L[dir][coeffloc_L ])+eps;
+						float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab])+eps;
+						float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab])+eps;
+						
+						sfavea[coeffloc_ab] = (1-exp(-(mag_a/mad_a)-(mag_L/(9*madL))));
+						sfaveb[coeffloc_ab] = (1-exp(-(mag_b/mad_b)-(mag_L/(9*madL))));
+						
+						// 'firm' threshold of chroma coefficients
+						//WavCoeffs_a[dir][coeffloc_ab] *= (1-exp(-(mag_a/mad_a)-(mag_L/(9*madL))));//(coeff_a>2*thresh_a ? 1 : (coeff_a<thresh_a ? 0 : (coeff_a/thresh_a - 1)));
+						//WavCoeffs_b[dir][coeffloc_ab] *= (1-exp(-(mag_b/mad_b)-(mag_L/(9*madL))));//(coeff_b>2*thresh_b ? 1 : (coeff_b<thresh_b ? 0 : (coeff_b/thresh_b - 1)));
+												
+						
+					}
+				}//now chrominance coefficients are denoised
+				
+				boxblur(sfavea, sfavea, level+2, level+2, W_ab, H_ab);//increase smoothness by locally averaging shrinkage
+				boxblur(sfaveb, sfaveb, level+2, level+2, W_ab, H_ab);//increase smoothness by locally averaging shrinkage
+
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+				for (int i=0; i<H_ab; i++) 
+					for (int j=0; j<W_ab; j++) {
+						
+						int coeffloc_ab = i*W_ab+j;
+						int coeffloc_L	= ((i*skip_L)/skip_ab)*W_L + ((j*skip_L)/skip_ab);
+						
+						float mag_L = SQR(WavCoeffs_L[dir][coeffloc_L ])+eps;
+						float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab])+eps;
+						float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab])+eps;
+						
+						float sfa = (1-exp(-(mag_a/mad_a)-(mag_L/(9*madL))));
+						float sfb = (1-exp(-(mag_b/mad_b)-(mag_L/(9*madL))));
+						
+						//use smoothed shrinkage unless local shrinkage is much less
+						WavCoeffs_a[dir][coeffloc_ab] *= (SQR(sfavea[coeffloc_ab])+SQR(sfa))/(sfavea[coeffloc_ab]+sfa+eps);
+						WavCoeffs_b[dir][coeffloc_ab] *= (SQR(sfaveb[coeffloc_ab])+SQR(sfb))/(sfaveb[coeffloc_ab]+sfb+eps);
+						
+					}//now chrominance coefficients are denoised
+			}
+			
+			if (noisevar_L>0.01) {
+//OpenMP here				
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+				for (int i=0; i<W_L*H_L; i++) {
+
+					float mag = SQR(WavCoeffs_L[dir][i]);
+					float shrinkfactor = mag/(mag+mad_L*exp(-mag/(9*mad_L))+eps);
+										
+					//WavCoeffs_L[dir][i] *= shrinkfactor;
+					sfave[i] = shrinkfactor;
+				}
+//OpenMP here				
+				boxblur(sfave, sfave, level+2, level+2, W_L, H_L);//increase smoothness by locally averaging shrinkage
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+				for (int i=0; i<W_L*H_L; i++) {
+					
+					
+					float mag = SQR(WavCoeffs_L[dir][i]);
+					float sf = mag/(mag+mad_L*exp(-mag/(9*mad_L))+eps);
+					
+					//use smoothed shrinkage unless local shrinkage is much less
+					WavCoeffs_L[dir][i] *= (SQR(sfave[i])+SQR(sf))/(sfave[i]+sf+eps);
+					
+				}//now luminance coefficients are denoised
+			}
+			
+			
+		}
+		delete[] sfave;
+		delete[] sfavea;
+		delete[] sfaveb;
+
+	}
+
+	
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	
+}
+
diff --git a/rtengine/bilateral2.h b/rtengine/bilateral2.h
index 4154387a6..d80dfcd21 100644
--- a/rtengine/bilateral2.h
+++ b/rtengine/bilateral2.h
@@ -30,6 +30,7 @@
 #include "mytime.h"
 #include "gauss.h"
 
+#include "array2D.h"
 #ifdef _OPENMP
 #include <omp.h>
 #endif
@@ -437,34 +438,39 @@ template<class T, class A> void bilateral (T** src, T** dst, T** buffer, int W,
 // START OF EXPERIMENTAL CODE: O(1) bilateral box filter
 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-#define BINBIT 12
-#define TRANSBIT 4
+#define MAXVAL  65535.0
+#define CLIP(a) ((a)>0.0?((a)<MAXVAL?(a):MAXVAL):0.0)
+
+#define BINBIT 7 //bit depth of histogram -- there are 2^BINBIT discrete levels
+#define TRANSBIT 9 //bit shift = 16-BINBIT taking short ints to bit depth BINBIT
 
 template<class T> void bilateral (T** src, T** dst, int W, int H, int sigmar, double sigmas, int row_from, int row_to) {
 
     // range weights
-    LUTf ec(0x20000);
+    /*LUTf ec(0x20000);
     for (int i=0; i<0x20000; i++) 
-        ec[i] = exp(-(double)(i-0x10000)*(double)(i-0x10000) / (2.0*sigmar*sigmar)); 
+        ec[i] = exp(-(double)(i-0x10000)*(double)(i-0x10000) / (2.0*sigmar*sigmar)); */
 
     // histogram
     LUTu hist (1<<BINBIT);
     LUTu rhist(1<<BINBIT);
 
     // buffer for the final image
-    float** buff_final = new float*[H];
+    /*float** buff_final = new float*[H];
     float * real_buff_final = new float[W*H];
     for (int i=0; i<H; i++) {
         buff_final[i] = real_buff_final + i*W;
         memset (buff_final[i], 0, W*sizeof(float));
-    }
+    }*/
+	
+	array2D<float> buff_final(W,H,ARRAY2D_CLEAR_DATA);
     
     int r = sigmas;
    
     // calculate histogram at the beginning of the row
     rhist.clear();
-    for (int x = std::max(0,row_from-r-1); x<row_from+r; x++)
-        for (int y = 0; y<r; y++)
+    for (int x = MAX(0,row_from-r); x<=MIN(H,row_from+r); x++)
+        for (int y = 0; y<r+1; y++)
             rhist[((int)src[x][y])>>TRANSBIT]++;
 
     sigmar*=2;
@@ -473,52 +479,52 @@ template<class T> void bilateral (T** src, T** dst, int W, int H, int sigmar, do
         
         // calculate histogram at the beginning of the row
         if (i>r)
-            for (int x = 0; x<r; x++) 
+            for (int x = 0; x<=MIN(H,r); x++) 
                 rhist[((int)src[i-r-1][x])>>TRANSBIT]--;
         if (i<H-r)
-            for (int x = 0; x<r; x++) 
+            for (int x = 0; x<=MIN(H,r); x++) 
                 rhist[((int)src[i+r][x])>>TRANSBIT]++;
             
         hist=rhist;
         for (int j=0; j<W; j++) {
         
-            // substract pixels at the left and add pixels at the right
+            // subtract pixels at the left and add pixels at the right
             if (j>r)
-                for (int x=std::max(0,i-r); x<=std::min(i+r,H-1); x++)
+                for (int x=MAX(0,i-r); x<=MIN(i+r,H-1); x++) 
                     hist[(int)(src[x][j-r-1])>>TRANSBIT]--;
             if (j<W-r)
-                for (int x=std::max(0,i-r); x<=std::min(i+r,H-1); x++)
+                for (int x=MAX(0,i-r); x<=MIN(i+r,H-1); x++) 
                     hist[((int)src[x][j+r])>>TRANSBIT]++;
 
             // calculate pixel value
-            float weight = 0.0;
+            float totwt = 0.0, weight;
             for (int k=0; k<=(sigmar>>TRANSBIT); k++) {               
                 float w = 1.0 - (double)k/(sigmar>>TRANSBIT);
                 int v = (int)(src[i][j])>>TRANSBIT;
+				//float frac = ((float)(src[i][j]-(v<<TRANSBIT)))/(1<<TRANSBIT);
                 if (v-k >= 0) {
-                    weight += hist [v-k] * w;
-                    buff_final[i][j] += hist [v-k] * w * (src[i][j]-(k<<TRANSBIT));
+                    weight = hist [v-k/*+frac*/] * w;
+					totwt += weight;
+                    buff_final[i][j] += weight * (src[i][j]-(k<<TRANSBIT));
                 }
                 if (v+k < (1<<BINBIT)) {
-                    weight += hist [v+k] * w;
-                    buff_final[i][j] += hist [v+k] * w * (src[i][j]+(k<<TRANSBIT));
+                    weight = hist [v+k/*+frac*/] * w;
+					totwt += weight;
+                    buff_final[i][j] += weight * (src[i][j]+(k<<TRANSBIT));
                 }
             }
-            buff_final[i][j] /= weight;
+            buff_final[i][j] /= totwt;
         }
     }
     for (int i=row_from; i<row_to; i++) 
         for (int j=0; j<W; j++) 
             dst[i][j] = (T)CLIP(buff_final[i][j]);
 
-    //delete [] ec;
-    //delete [] hist;
-    //delete [] rhist;
-    //for (int i=0; i<H; i++)
-    //    delete [] buff_final[i];
-    delete [] real_buff_final;
-    delete [] buff_final;
+    //delete [] real_buff_final;
+    //delete [] buff_final;
 }
+#undef MAXVAL
+#undef CLIP
 
 #undef BINBIT
 #undef TRANSBIT
diff --git a/rtengine/boxblur.h b/rtengine/boxblur.h
new file mode 100644
index 000000000..40ab39755
--- /dev/null
+++ b/rtengine/boxblur.h
@@ -0,0 +1,699 @@
+/*
+ *  This file is part of RawTherapee.
+ *
+ *  Copyright © 2010 Emil Martinec <ejmartin@uchicago.edu>
+ *
+ *  RawTherapee 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.
+ * 
+ *  RawTherapee 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 RawTherapee.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#ifndef _BOXBLUR_H_
+#define _BOXBLUR_H_
+
+#include <stdlib.h>
+#include <string.h> 
+#include <math.h>
+#include "alignedbuffer.h"
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+#include "rt_math.h"
+//using namespace rtengine;
+
+namespace rtengine {
+
+// classical filtering if the support window is small:
+
+template<class T, class A> void boxblur (T** src, A** dst, int radx, int rady, int W, int H) {
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//box blur image; box range = (radx,rady)
+	
+	AlignedBuffer<float>* buffer = new AlignedBuffer<float> (W*H);
+	float* temp = buffer->data;
+	
+	if (radx==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				temp[row*W+col] = (float)src[row][col];
+			}
+	} else {
+		//horizontal blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row = 0; row < H; row++) {
+			int len = radx + 1;
+			temp[row*W+0] = (float)src[row][0]/len;
+			for (int j=1; j<=radx; j++) {
+				temp[row*W+0] += (float)src[row][j]/len;
+			}
+			for (int col=1; col<=radx; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len + (float)src[row][col+radx])/(len+1);
+				len ++;
+			}
+			for (int col = radx+1; col < W-radx; col++) {
+				temp[row*W+col] = temp[row*W+col-1] + ((float)(src[row][col+radx] - src[row][col-radx-1]))/len;
+			}
+			for (int col=W-radx; col<W; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len - src[row][col-radx-1])/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	if (rady==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				dst[row][col] = temp[row*W+col];
+			}
+	} else {
+		//vertical blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int col = 0; col < W; col++) {
+			int len = rady + 1;
+			dst[0][col] = temp[0*W+col]/len;
+			for (int i=1; i<=rady; i++) {
+				dst[0][col] += temp[i*W+col]/len;
+			}
+			for (int row=1; row<=rady; row++) {
+				dst[row][col] = (dst[(row-1)][col]*len + temp[(row+rady)*W+col])/(len+1);
+				len ++;
+			}
+			for (int row = rady+1; row < H-rady; row++) {
+				dst[row][col] = dst[(row-1)][col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
+			}
+			for (int row=H-rady; row<H; row++) {
+				dst[row][col] = (dst[(row-1)][col]*len - temp[(row-rady-1)*W+col])/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	delete buffer;
+	
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+template<class T, class A> void boxblur (T* src, A* dst, int radx, int rady, int W, int H) {
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//box blur image; box range = (radx,rady) i.e. box size is (2*radx+1)x(2*rady+1)
+	
+	AlignedBuffer<float>* buffer = new AlignedBuffer<float> (W*H);
+	float* temp = buffer->data;
+		
+	if (radx==0) {
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				temp[row*W+col] = src[row*W+col];
+			}
+	} else {
+		//horizontal blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row = 0; row < H; row++) {
+			int len = radx + 1;
+			temp[row*W+0] = (float)src[row*W+0]/len;
+			for (int j=1; j<=radx; j++) {
+				temp[row*W+0] += (float)src[row*W+j]/len;
+			}
+			for (int col=1; col<=radx; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len + src[row*W+col+radx])/(len+1);
+				len ++;
+			}
+			for (int col = radx+1; col < W-radx; col++) {
+				temp[row*W+col] = temp[row*W+col-1] + ((float)(src[row*W+col+radx] - src[row*W+col-radx-1]))/len;
+			}
+			for (int col=W-radx; col<W; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len - src[row*W+col-radx-1])/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	if (rady==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				dst[row*W+col] = temp[row*W+col];
+			}
+	} else {
+		//vertical blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int col = 0; col < W; col++) {
+			int len = rady + 1;
+			dst[0*W+col] = temp[0*W+col]/len;
+			for (int i=1; i<=rady; i++) {
+				dst[0*W+col] += temp[i*W+col]/len;
+			}
+			for (int row=1; row<=rady; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
+				len ++;
+			}
+			for (int row = rady+1; row < H-rady; row++) {
+				dst[row*W+col] = dst[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
+			}
+			for (int row=H-rady; row<H; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	delete buffer;
+	
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+template<typename T> void boxvar (T* src, T* dst, int radx, int rady, int W, int H) {
+	
+	AlignedBuffer<float> buffer1(W*H);
+	AlignedBuffer<float> buffer2(W*H);
+	float* tempave = buffer1.data;
+	float* tempsqave = buffer2.data;
+	
+	AlignedBufferMP<float> buffer3(H);
+
+	//float image_ave = 0;
+	
+	//box blur image channel; box size = 2*box+1
+	//horizontal blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int row = 0; row < H; row++) {
+		int len = radx + 1;
+		tempave[row*W+0] = src[row*W+0]/len;
+		tempsqave[row*W+0] = SQR(src[row*W+0])/len;
+		for (int j=1; j<=radx; j++) {
+			tempave[row*W+0] += src[row*W+j]/len;
+			tempsqave[row*W+0] += SQR(src[row*W+j])/len;
+		}
+		for (int col=1; col<=radx; col++) {
+			tempave[row*W+col] = (tempave[row*W+col-1]*len + src[row*W+col+radx])/(len+1);
+			tempsqave[row*W+col] = (tempsqave[row*W+col-1]*len + SQR(src[row*W+col+radx]))/(len+1);
+			len ++;
+		}
+		for (int col = radx+1; col < W-radx; col++) {
+			tempave[row*W+col] = tempave[row*W+col-1] + (src[row*W+col+radx] - src[row*W+col-radx-1])/len;
+			tempsqave[row*W+col] = tempsqave[row*W+col-1] + (SQR(src[row*W+col+radx]) - SQR(src[row*W+col-radx-1]))/len;
+		}
+		for (int col=W-radx; col<W; col++) {
+			tempave[row*W+col] = (tempave[row*W+col-1]*len - src[row*W+col-radx-1])/(len-1);
+			tempsqave[row*W+col] = (tempsqave[row*W+col-1]*len - SQR(src[row*W+col-radx-1]))/(len-1);
+			len --;
+		}
+	}
+	
+	//vertical blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+	for (int col = 0; col < W; col++) {
+		AlignedBuffer<float>* pBuf3 = buffer3.acquire();
+		T* tempave2=(T*)pBuf3->data;
+			
+		int len = rady + 1;
+		tempave2[0] = tempave[0*W+col]/len;
+		dst[0*W+col] = tempsqave[0*W+col]/len;
+		for (int i=1; i<=rady; i++) {
+			tempave2[0] += tempave[i*W+col]/len;
+			dst[0*W+col] += tempsqave[i*W+col]/len;
+		}
+		for (int row=1; row<=rady; row++) {
+			tempave2[row] = (tempave2[(row-1)]*len + tempave[(row+rady)*W+col])/(len+1);
+			dst[row*W+col] = (dst[(row-1)*W+col]*len + tempsqave[(row+rady)*W+col])/(len+1);
+			len ++;
+		}
+		for (int row = rady+1; row < H-rady; row++) {
+			tempave2[row] = tempave2[(row-1)] + (tempave[(row+rady)*W+col] - tempave[(row-rady-1)*W+col])/len;
+			dst[row*W+col] = dst[(row-1)*W+col] + (tempsqave[(row+rady)*W+col] - tempsqave[(row-rady-1)*W+col])/len;
+		}
+		for (int row=H-rady; row<H; row++) {
+			tempave2[row] = (tempave2[(row-1)]*len - tempave[(row-rady-1)*W+col])/(len-1);
+			dst[row*W+col] = (dst[(row-1)*W+col]*len - tempsqave[(row-rady-1)*W+col])/(len-1);
+			len --;
+		}
+		//now finish off
+		for (int row=0; row<H; row++) {
+			dst[row*W+col] = fabs(dst[row*W+col] - SQR(tempave2[row]));
+			//image_ave += src[row*W+col];
+		}
+		buffer3.release(pBuf3);
+	}
+	
+	//image_ave /= (W*H);	
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+template<typename T> void boxdev (T* src, T* dst, int radx, int rady, int W, int H) {
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//box blur image; box range = (radx,rady) i.e. box size is (2*radx+1)x(2*rady+1)
+	
+	AlignedBuffer<float>* buffer1 = new AlignedBuffer<float> (W*H);
+	float* temp = buffer1->data;
+	
+	AlignedBuffer<float>* buffer2 = new AlignedBuffer<float> (W*H);
+	float* tempave = buffer2->data;
+	
+	if (radx==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif	
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				temp[row*W+col] = src[row*W+col];
+			}
+	} else {
+		//horizontal blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row = 0; row < H; row++) {
+			int len = radx + 1;
+			temp[row*W+0] = (float)src[row*W+0]/len;
+			for (int j=1; j<=radx; j++) {
+				temp[row*W+0] += (float)src[row*W+j]/len;
+			}
+			for (int col=1; col<=radx; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len + src[row*W+col+radx])/(len+1);
+				len ++;
+			}
+			for (int col = radx+1; col < W-radx; col++) {
+				temp[row*W+col] = temp[row*W+col-1] + ((float)(src[row*W+col+radx] - src[row*W+col-radx-1]))/len;
+			}
+			for (int col=W-radx; col<W; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len - src[row*W+col-radx-1])/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	if (rady==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++){
+			for (int col=0; col<W; col++) {
+				tempave[row*W+col] = temp[row*W+col];
+			}
+		}
+	} else {
+		//vertical blur
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int col = 0; col < W; col++) {
+			int len = rady + 1;
+			tempave[0*W+col] = temp[0*W+col]/len;
+			for (int i=1; i<=rady; i++) {
+				tempave[0*W+col] += temp[i*W+col]/len;
+			}
+			for (int row=1; row<=rady; row++) {
+				tempave[row*W+col] = (tempave[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
+				len ++;
+			}
+			for (int row = rady+1; row < H-rady; row++) {
+				tempave[row*W+col] = tempave[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
+			}
+			for (int row=H-rady; row<H; row++) {
+				tempave[row*W+col] = (tempave[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//box blur absolute deviation
+
+	
+	if (radx==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				temp[row*W+col] = fabs(src[row*W+col]-tempave[row*W+col]);
+			}
+	} else {
+		//horizontal blur
+//OpenMP here
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif		
+		for (int row = 0; row < H; row++) {
+			int len = radx + 1;
+			temp[row*W+0] = fabs(src[row*W+0]-tempave[row*W+0])/len;
+			for (int j=1; j<=radx; j++) {
+				temp[row*W+0] += fabs(src[row*W+j]-tempave[row*W+j])/len;
+			}
+			for (int col=1; col<=radx; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len + fabs(src[row*W+col+radx]-tempave[row*W+col+radx]))/(len+1);
+				len ++;
+			}
+			for (int col = radx+1; col < W-radx; col++) {
+				temp[row*W+col] = temp[row*W+col-1] + (fabs(src[row*W+col+radx]-tempave[row*W+col+radx]) - \
+													   fabs(src[row*W+col-radx-1]-tempave[row*W+col-radx-1]))/len;
+			}
+			for (int col=W-radx; col<W; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len - fabs(src[row*W+col-radx-1]-tempave[row*W+col-radx-1]))/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	if (rady==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				dst[row*W+col] = temp[row*W+col];
+			}
+	} else {
+		//vertical blur
+//OpenMP here		
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int col = 0; col < W; col++) {
+			int len = rady + 1;
+			dst[0*W+col] = temp[0*W+col]/len;
+			for (int i=1; i<=rady; i++) {
+				dst[0*W+col] += temp[i*W+col]/len;
+			}
+			for (int row=1; row<=rady; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
+				len ++;
+			}
+			for (int row = rady+1; row < H-rady; row++) {
+				dst[row*W+col] = dst[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
+			}
+			for (int row=H-rady; row<H; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
+				len --;
+			}
+		}
+	}
+			
+	delete buffer1;
+	delete buffer2;
+
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+template<class T, class A> void boxsqblur (T* src, A* dst, int radx, int rady, int W, int H) {
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//box blur image; box range = (radx,rady) i.e. box size is (2*radx+1)x(2*rady+1)
+	
+	AlignedBuffer<float>* buffer = new AlignedBuffer<float> (W*H);
+	float* temp = buffer->data;
+	
+	if (radx==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				temp[row*W+col] = SQR(src[row*W+col]);
+			}
+	} else {
+		//horizontal blur
+//OpenMP here		
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row = 0; row < H; row++) {
+			int len = radx + 1;
+			temp[row*W+0] = SQR((float)src[row*W+0])/len;
+			for (int j=1; j<=radx; j++) {
+				temp[row*W+0] += SQR((float)src[row*W+j])/len;
+			}
+			for (int col=1; col<=radx; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len + SQR(src[row*W+col+radx]))/(len+1);
+				len ++;
+			}
+			for (int col = radx+1; col < W-radx; col++) {
+				temp[row*W+col] = temp[row*W+col-1] + ((float)(SQR(src[row*W+col+radx]) - SQR(src[row*W+col-radx-1])))/len;
+			}
+			for (int col=W-radx; col<W; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len - SQR(src[row*W+col-radx-1]))/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	if (rady==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				dst[row*W+col] = temp[row*W+col];
+			}
+	} else {
+		//vertical blur
+//OpenMP here		
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int col = 0; col < W; col++) {
+			int len = rady + 1;
+			dst[0*W+col] = temp[0*W+col]/len;
+			for (int i=1; i<=rady; i++) {
+				dst[0*W+col] += temp[i*W+col]/len;
+			}
+			for (int row=1; row<=rady; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
+				len ++;
+			}
+			for (int row = rady+1; row < H-rady; row++) {
+				dst[row*W+col] = dst[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
+			}
+			for (int row=H-rady; row<H; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
+				len --;
+			}
+		}
+	}	
+			
+	delete buffer;
+	
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+template<class T, class A> void boxcorrelate (T* src, A* dst, int dx, int dy, int radx, int rady, int W, int H) {
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//box blur image; box range = (radx,rady) i.e. box size is (2*radx+1)x(2*rady+1)
+	
+	AlignedBuffer<float>* buffer = new AlignedBuffer<float> (W*H);
+	float* temp = buffer->data;
+	
+	if (radx==0) {
+		for (int row=0; row<H; row++) {
+			int rr = min(H-1,max(0,row+dy));
+			for (int col=0; col<W; col++) {
+				int cc = min(W-1,max(0,col+dx));
+				temp[row*W+col] = dy>0 ? (src[row*W+col])*(src[rr*W+cc]) : 0;
+			}
+		}
+	} else {
+		//horizontal blur
+		for (int row = 0; row < H; row++) {
+			int len = radx + 1;
+			int rr = min(H-1,max(0,row+dy));
+			int cc = min(W-1,max(0,0+dx));
+			temp[row*W+0] = ((float)src[row*W+0])*(src[rr*W+cc])/len;
+			for (int j=1; j<=radx; j++) {
+				int cc = min(W-1,max(0,j+dx));
+				temp[row*W+0] += ((float)src[row*W+j])*(src[rr*W+cc])/len;
+			}
+			for (int col=1; col<=radx; col++) {
+				int cc = min(W-1,max(0,col+dx+radx));
+				temp[row*W+col] = (temp[row*W+col-1]*len + (src[row*W+col+radx])*(src[rr*W+cc]))/(len+1);
+				len ++;
+			}
+			for (int col = radx+1; col < W-radx; col++) {
+				int cc = min(W-1,max(0,col+dx+radx));
+				int cc1 = min(W-1,max(0,col+dx-radx-1));
+				temp[row*W+col] = temp[row*W+col-1] + ((float)((src[row*W+col+radx])*(src[rr*W+cc]) - 
+															   (src[row*W+col-radx-1])*(src[rr*W+cc1])))/len;
+			}
+			for (int col=W-radx; col<W; col++) {
+				int cc1 = min(W-1,max(0,col+dx-radx-1));
+				temp[row*W+col] = (temp[row*W+col-1]*len - (src[row*W+col-radx-1])*(src[rr*W+cc1]))/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	if (rady==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				dst[row*W+col] = temp[row*W+col];
+			}
+	} else {
+		//vertical blur
+//OpenMP here		
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int col = 0; col < W; col++) {
+			int len = rady + 1;
+			dst[0*W+col] = temp[0*W+col]/len;
+			for (int i=1; i<=rady; i++) {
+				dst[0*W+col] += temp[i*W+col]/len;
+			}
+			for (int row=1; row<=rady; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
+				len ++;
+			}
+			for (int row = rady+1; row < H-rady; row++) {
+				dst[row*W+col] = dst[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
+			}
+			for (int row=H-rady; row<H; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	delete buffer;
+	
+}
+
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+template<class T, class A> void boxabsblur (T* src, A* dst, int radx, int rady, int W, int H) {
+	
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//box blur image; box range = (radx,rady) i.e. box size is (2*radx+1)x(2*rady+1)
+	
+	AlignedBuffer<float>* buffer = new AlignedBuffer<float> (W*H);
+	float* temp = buffer->data;
+	
+	if (radx==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				temp[row*W+col] = fabs(src[row*W+col]);
+			}
+	} else {
+		//horizontal blur
+//OpenMP here		
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row = 0; row < H; row++) {
+			int len = radx + 1;
+			temp[row*W+0] = fabs((float)src[row*W+0])/len;
+			for (int j=1; j<=radx; j++) {
+				temp[row*W+0] += fabs((float)src[row*W+j])/len;
+			}
+			for (int col=1; col<=radx; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len + fabs(src[row*W+col+radx]))/(len+1);
+				len ++;
+			}
+			for (int col = radx+1; col < W-radx; col++) {
+				temp[row*W+col] = temp[row*W+col-1] + ((float)(fabs(src[row*W+col+radx]) - fabs(src[row*W+col-radx-1])))/len;
+			}
+			for (int col=W-radx; col<W; col++) {
+				temp[row*W+col] = (temp[row*W+col-1]*len - fabs(src[row*W+col-radx-1]))/(len-1);
+				len --;
+			}
+		}
+	}
+	
+	if (rady==0) {
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int row=0; row<H; row++) 
+			for (int col=0; col<W; col++) {
+				dst[row*W+col] = temp[row*W+col];
+			}
+	} else {
+		//vertical blur
+//OpenMP here		
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+		for (int col = 0; col < W; col++) {
+			int len = rady + 1;
+			dst[0*W+col] = temp[0*W+col]/len;
+			for (int i=1; i<=rady; i++) {
+				dst[0*W+col] += temp[i*W+col]/len;
+			}
+			for (int row=1; row<=rady; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len + temp[(row+rady)*W+col])/(len+1);
+				len ++;
+			}
+			for (int row = rady+1; row < H-rady; row++) {
+				dst[row*W+col] = dst[(row-1)*W+col] + (temp[(row+rady)*W+col] - temp[(row-rady-1)*W+col])/len;
+			}
+			for (int row=H-rady; row<H; row++) {
+				dst[row*W+col] = (dst[(row-1)*W+col]*len - temp[(row-rady-1)*W+col])/(len-1);
+				len --;
+			}
+	}
+}
+
+	delete buffer;
+	
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+}
+#endif /* _BOXBLUR_H_ */
diff --git a/rtengine/cplx_wavelet_dec.cc b/rtengine/cplx_wavelet_dec.cc
new file mode 100644
index 000000000..27b885827
--- /dev/null
+++ b/rtengine/cplx_wavelet_dec.cc
@@ -0,0 +1,48 @@
+/*
+ *  This file is part of RawTherapee.
+ *
+ *  RawTherapee 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.
+ * 
+ *  RawTherapee 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 RawTherapee.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *  2010 Ilya Popov <ilia_popov@rambler.ru>
+ *  2012 Emil Martinec <ejmartin@uchicago.edu>
+ */
+
+#include "cplx_wavelet_dec.h"
+
+namespace rtengine {
+	
+	cplx_wavelet_decomposition::~cplx_wavelet_decomposition()
+	{
+		for(int i = 0; i <= lvltot; i++) {
+			for (int j=0; j<4; j++) {
+				delete dual_tree[i][j];
+			}
+		}
+		delete[] first_lev_anal;
+		delete[] first_lev_synth;
+		delete[] wavfilt_anal;
+		delete[] wavfilt_synth;
+	}
+	
+	wavelet_decomposition::~wavelet_decomposition()
+	{
+		for(int i = 0; i <= lvltot; i++) {
+			delete wavelet_decomp[i];
+		}
+		delete[] wavfilt_anal;
+		delete[] wavfilt_synth;
+	}
+
+};
+
diff --git a/rtengine/cplx_wavelet_dec.h b/rtengine/cplx_wavelet_dec.h
new file mode 100644
index 000000000..0df25fa72
--- /dev/null
+++ b/rtengine/cplx_wavelet_dec.h
@@ -0,0 +1,454 @@
+/*
+ *  This file is part of RawTherapee.
+ *
+ *  RawTherapee 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.
+ * 
+ *  RawTherapee 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 RawTherapee.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *  2010 Ilya Popov <ilia_popov@rambler.ru>
+ *  2012 Emil Martinec <ejmartin@uchicago.edu>
+ */
+
+#ifndef CPLX_WAVELET_DEC_H_INCLUDED
+#define CPLX_WAVELET_DEC_H_INCLUDED
+
+#include <cstddef>
+#include <math.h>
+
+
+#include "cplx_wavelet_level.h"
+#include "cplx_wavelet_filter_coeffs.h"
+
+namespace rtengine {
+
+
+// %%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+template <typename A>
+void copy_out(A * a, A * b, size_t datalen)
+{// for standard wavelet decomposition
+	memcpy(b, a, datalen*sizeof(A));
+}
+
+template <typename A, typename B>
+void copy_out(A ** a, B * b, size_t datalen)
+{// for complex wavelet decomposition
+	for (size_t j=0; j<datalen; j++) {
+		b[j] = static_cast<B> (0.25*(a[0][j]+a[1][j]+a[2][j]+a[3][j]));
+	}
+}
+	
+template <typename A, typename B>
+void copy_out(A * a, B * b, size_t datalen)
+{// for standard wavelet decomposition
+	for (size_t j=0; j<datalen; j++) {
+		b[j] = static_cast<B> (a[j]);
+	}
+}
+
+// %%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+class cplx_wavelet_decomposition
+{
+public:
+
+    typedef float internal_type;
+
+private:
+
+    static const int maxlevels = 8;//should be greater than any conceivable order of decimation
+    
+    int lvltot, subsamp;
+    size_t m_w, m_h;//dimensions
+	
+	int first_lev_len, first_lev_offset;
+	float *first_lev_anal;
+	float *first_lev_synth;
+	
+	int wavfilt_len, wavfilt_offset;
+	float *wavfilt_anal;
+	float *wavfilt_synth;
+	
+	int testfilt_len, testfilt_offset;
+	float *testfilt_anal;
+	float *testfilt_synth;
+
+    wavelet_level<internal_type> * dual_tree[maxlevels][4];
+    
+public:
+
+    template<typename E>
+    cplx_wavelet_decomposition(E * src, int width, int height, int maxlvl, int subsampling);
+    
+    ~cplx_wavelet_decomposition();
+	
+	internal_type ** level_coeffs(int level, int branch) const
+	{
+		return dual_tree[level][branch]->subbands();
+	}
+	
+	int level_W(int level, int branch) const
+	{
+		return dual_tree[level][branch]->width();
+	}
+    
+	int level_H(int level, int branch) const
+	{
+		return dual_tree[level][branch]->height();
+	}
+	
+	int level_pad(int level, int branch) const
+	{
+		return dual_tree[level][branch]->padding();
+	}
+	
+	int maxlevel() const
+	{
+		return lvltot;
+	}
+	
+    template<typename E>
+    void reconstruct(E * dst);
+	
+};
+	
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	template<typename E>
+	cplx_wavelet_decomposition::cplx_wavelet_decomposition(E * src, int width, int height, int maxlvl, int subsampling)
+	: lvltot(0), subsamp(subsampling), m_w(width), m_h(height)
+	{
+		
+		//initialize wavelet filters
+		
+		first_lev_len = FSFarras_len;
+		first_lev_offset = FSFarras_offset;
+		first_lev_anal = new float[4*first_lev_len];
+		first_lev_synth = new float[4*first_lev_len];
+
+		for (int n=0; n<2; n++) {
+			for (int m=0; m<2; m++) {
+				for (int i=0; i<first_lev_len; i++) {
+					first_lev_anal[first_lev_len*(2*n+m)+i]  = FSFarras_anal[n][m][i]/sqrt(2);
+					first_lev_synth[first_lev_len*(2*n+m)+i] = FSFarras_anal[n][m][first_lev_len-1-i]/sqrt(2);
+				}
+			}
+		}
+		
+		/*first_lev_len = AntonB_len;
+		first_lev_offset = AntonB_offset;
+		first_lev_anal = new float[4*first_lev_len];
+		first_lev_synth = new float[4*first_lev_len];
+		
+		for (int n=0; n<2; n++) {
+			for (int m=0; m<2; m++) {
+				for (int i=0; i<first_lev_len; i++) {
+					first_lev_anal[first_lev_len*(2*n+m)+i]  = AntonB_anal[n][m][i];
+					first_lev_synth[first_lev_len*(2*n+m)+i] = 2*AntonB_synth[n][m][i];
+				}
+			}
+		}*/
+		
+		wavfilt_len = Kingsbury_len;
+		wavfilt_offset = Kingsbury_offset;
+		wavfilt_anal = new float[4*wavfilt_len];
+		wavfilt_synth = new float[4*wavfilt_len];
+		
+		for (int n=0; n<2; n++) {
+			for (int m=0; m<2; m++) {
+				for (int i=0; i<wavfilt_len; i++) {
+					wavfilt_anal[wavfilt_len*(2*n+m)+i]  = Kingsbury_anal[n][m][i]/sqrt(2);
+					wavfilt_synth[wavfilt_len*(2*n+m)+i] = Kingsbury_anal[n][m][first_lev_len-1-i]/sqrt(2);
+				}
+			}
+		}
+		
+		
+		// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+
+		
+		// after coefficient rotation, data structure is:
+		// dual_tree[scale][2*n+m=2*(Re/Im)+dir][channel={lo,hi1,hi2,hi3}][pixel_array]
+		
+		//srand((unsigned)time(0));
+		//for (int i=0; i<m_w*m_h; i++ )
+		//	src[i] = (float)rand()/(float)RAND_MAX;
+
+		for (int n=0; n<2; n++) {
+			for (int m=0; m<2; m++) {
+				lvltot=0;
+				float padding = 0;//1<<(maxlvl-1);
+				dual_tree[0][2*n+m] = new wavelet_level<internal_type>(src, lvltot, subsamp, padding, m_w, m_h, first_lev_anal+first_lev_len*2*n, \
+																				   first_lev_anal+first_lev_len*2*m, first_lev_len, first_lev_offset);
+				while(lvltot < maxlvl) {
+					lvltot++;
+					dual_tree[lvltot][2*n+m] = new wavelet_level<internal_type>(dual_tree[lvltot-1][2*n+m]->lopass()/*lopass*/, lvltot, subsamp, 0/*no padding*/, \
+																							dual_tree[lvltot-1][2*n+m]->width(), \
+																							dual_tree[lvltot-1][2*n+m]->height(), \
+																							wavfilt_anal+wavfilt_len*2*n, wavfilt_anal+wavfilt_len*2*m, \
+																							wavfilt_len, wavfilt_offset);
+				}
+			}
+		}
+		
+		
+		//rotate detail coefficients
+		float coeffave[5][4][3];
+
+		float root2 = sqrt(2);
+		for (int lvl=0; lvl<lvltot; lvl++) {
+			int Wlvl = dual_tree[lvl][0]->width();
+			int Hlvl = dual_tree[lvl][0]->height();
+			for (int n=0; n<4; n++) 
+				for (int m=1; m<4; m++)
+					coeffave[lvl][n][m-1]=0;
+
+			for (int m=1; m<4; m++) {//detail coefficients only
+				for (int i=0; i<Wlvl*Hlvl; i++) {//pixel
+					
+					float wavtmp = (dual_tree[lvl][0]->wavcoeffs[m][i] + dual_tree[lvl][3]->wavcoeffs[m][i])/root2;
+					dual_tree[lvl][3]->wavcoeffs[m][i] = (dual_tree[lvl][0]->wavcoeffs[m][i] - dual_tree[lvl][3]->wavcoeffs[m][i])/root2;
+					dual_tree[lvl][0]->wavcoeffs[m][i] = wavtmp;
+					
+					wavtmp = (dual_tree[lvl][1]->wavcoeffs[m][i] + dual_tree[lvl][2]->wavcoeffs[m][i])/root2;
+					dual_tree[lvl][2]->wavcoeffs[m][i] = (dual_tree[lvl][1]->wavcoeffs[m][i] - dual_tree[lvl][2]->wavcoeffs[m][i])/root2;
+					dual_tree[lvl][1]->wavcoeffs[m][i] = wavtmp;
+					
+					for (int n=0; n<4; n++) coeffave[lvl][n][m-1] += fabs(dual_tree[lvl][n]->wavcoeffs[m][i]);
+				}
+			}
+			for (int n=0; n<4; n++) 
+				for (int i=0; i<3; i++) 
+					coeffave[lvl][n][i] /= Wlvl*Hlvl;
+		}
+
+	}
+	
+	/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
+	/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
+	
+		
+	template<typename E>
+	void cplx_wavelet_decomposition::reconstruct(E * dst) { 
+		
+		// data structure is wavcoeffs[scale][2*n+m=2*(Re/Im)+dir][channel={lo,hi1,hi2,hi3}][pixel_array]
+		
+		//rotate detail coefficients
+		float root2 = sqrt(2);
+		for (int lvl=0; lvl<lvltot; lvl++) {
+			int Wlvl = dual_tree[lvl][0]->width();
+			int Hlvl = dual_tree[lvl][0]->height();
+			for (int i=0; i<Wlvl*Hlvl; i++) {//pixel
+				for (int m=1; m<4; m++) {//detail coefficients only
+					float wavtmp = (dual_tree[lvl][0]->wavcoeffs[m][i] + dual_tree[lvl][3]->wavcoeffs[m][i])/root2;
+					dual_tree[lvl][3]->wavcoeffs[m][i] = (dual_tree[lvl][0]->wavcoeffs[m][i] - dual_tree[lvl][3]->wavcoeffs[m][i])/root2;
+					dual_tree[lvl][0]->wavcoeffs[m][i] = wavtmp;
+					
+					wavtmp = (dual_tree[lvl][1]->wavcoeffs[m][i] + dual_tree[lvl][2]->wavcoeffs[m][i])/root2;
+					dual_tree[lvl][2]->wavcoeffs[m][i] = (dual_tree[lvl][1]->wavcoeffs[m][i] - dual_tree[lvl][2]->wavcoeffs[m][i])/root2;
+					dual_tree[lvl][1]->wavcoeffs[m][i] = wavtmp;
+				}
+			}
+		}
+		
+		internal_type ** tmp = new internal_type *[4];
+		for (int i=0; i<4; i++) {
+			tmp[i] = new internal_type[m_w*m_h];
+		}
+		
+		for (int n=0; n<2; n++) {
+			for (int m=0; m<2; m++) {
+				int skip=1<<(lvltot-1);
+				for (int lvl=lvltot-1; lvl>0; lvl--) {
+					dual_tree[lvl][2*n+m]->reconstruct_level(dual_tree[lvl-1][2*n+m]->wavcoeffs[0], wavfilt_synth+wavfilt_len*2*n, \
+																	wavfilt_synth+wavfilt_len*2*m, wavfilt_len, wavfilt_offset);
+					skip /=2;
+				}
+				dual_tree[0][2*n+m]->reconstruct_level(tmp[2*n+m], first_lev_synth+first_lev_len*2*n, 
+															  first_lev_synth+first_lev_len*2*m, first_lev_len, first_lev_offset);
+			}
+		}
+				
+		copy_out(tmp,dst,m_w*m_h);
+		
+		for (int i=0; i<4; i++) {
+			delete[] tmp[i];
+		}
+		delete[] tmp;
+		
+		
+		
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	
+	
+	class wavelet_decomposition
+	{
+	public:
+		
+		typedef float internal_type;
+		
+	private:
+		
+		static const int maxlevels = 8;//should be greater than any conceivable order of decimation
+		
+		int lvltot, subsamp;
+		size_t m_w, m_h;//dimensions
+		
+		int wavfilt_len, wavfilt_offset;
+		float *wavfilt_anal;
+		float *wavfilt_synth;
+		
+		int testfilt_len, testfilt_offset;
+		float *testfilt_anal;
+		float *testfilt_synth;
+		
+		wavelet_level<internal_type> * wavelet_decomp[maxlevels];
+		
+	public:
+		
+		template<typename E>
+		wavelet_decomposition(E * src, int width, int height, int maxlvl, int subsampling);
+		
+		~wavelet_decomposition();
+		
+		internal_type ** level_coeffs(int level) const
+		{
+			return wavelet_decomp[level]->subbands();
+		}
+		
+		int level_W(int level) const
+		{
+			return wavelet_decomp[level]->width();
+		}
+		
+		int level_H(int level) const
+		{
+			return wavelet_decomp[level]->height();
+		}
+		
+		int level_pad(int level) const
+		{
+			return wavelet_decomp[level]->padding();
+		}
+		
+		int level_stride(int level) const
+		{
+			return wavelet_decomp[level]->stride();
+		}
+		
+		int maxlevel() const
+		{
+			return lvltot;
+		}
+		
+		int subsample() const
+		{
+			return subsamp;
+		}
+		
+		template<typename E>
+		void reconstruct(E * dst);
+		
+	};
+	
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	template<typename E>
+	wavelet_decomposition::wavelet_decomposition(E * src, int width, int height, int maxlvl, int subsampling)
+	: lvltot(0), subsamp(subsampling), m_w(width), m_h(height)
+	{
+		
+		//initialize wavelet filters
+		
+		wavfilt_len = Daub4_len;
+		wavfilt_offset = Daub4_offset;
+		wavfilt_anal = new float[2*wavfilt_len];
+		wavfilt_synth = new float[2*wavfilt_len];
+		
+		for (int n=0; n<2; n++) {
+			for (int i=0; i<wavfilt_len; i++) {
+				wavfilt_anal[wavfilt_len*(n)+i]  = Daub4_anal[n][i];
+				wavfilt_synth[wavfilt_len*(n)+i] = Daub4_anal[n][wavfilt_len-1-i];
+				//n=0 lopass, n=1 hipass
+			}
+		}
+		
+		
+		// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+		
+		
+		// after coefficient rotation, data structure is:
+		// wavelet_decomp[scale][channel={lo,hi1,hi2,hi3}][pixel_array]
+		
+		//srand((unsigned)time(0));
+		//for (int i=0; i<m_w*m_h; i++ )
+		//	src[i] = (float)rand()/(float)RAND_MAX;
+		
+		int padding = 0;//pow(2, maxlvl);//must be a multiple of two
+		lvltot=0;
+		wavelet_decomp[lvltot] = new wavelet_level<internal_type>(src, lvltot/*level*/, subsamp, padding/*padding*/, m_w, m_h, \
+																  wavfilt_anal, wavfilt_anal, wavfilt_len, wavfilt_offset);
+		while(lvltot < maxlvl) {
+			lvltot++;
+			wavelet_decomp[lvltot] = new wavelet_level<internal_type>(wavelet_decomp[lvltot-1]->lopass()/*lopass*/, lvltot/*level*/, subsamp, 0/*no padding*/, \
+																	  wavelet_decomp[lvltot-1]->width(), wavelet_decomp[lvltot-1]->height(), \
+																	  wavfilt_anal, wavfilt_anal, wavfilt_len, wavfilt_offset);
+		}
+		
+		
+	}
+	
+	/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
+	/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
+	
+	
+	template<typename E>
+	void wavelet_decomposition::reconstruct(E * dst) { 
+		
+		// data structure is wavcoeffs[scale][channel={lo,hi1,hi2,hi3}][pixel_array]
+		
+		//int skip=1<<(lvltot-1);
+		for (int lvl=lvltot-1; lvl>0; lvl--) {
+			wavelet_decomp[lvl]->reconstruct_level(wavelet_decomp[lvl-1]->wavcoeffs[0], wavfilt_synth, wavfilt_synth, wavfilt_len, wavfilt_offset);
+			//skip /=2;
+		}
+		
+		internal_type * tmp = new internal_type[m_w*m_h];
+
+		wavelet_decomp[0]->reconstruct_level(tmp, wavfilt_synth, wavfilt_synth, wavfilt_len, wavfilt_offset);
+		
+		copy_out(tmp,dst,m_w*m_h);
+		
+		delete[] tmp;
+		
+		
+		
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	
+	
+
+};
+
+#endif
diff --git a/rtengine/cplx_wavelet_filter_coeffs.h b/rtengine/cplx_wavelet_filter_coeffs.h
new file mode 100644
index 000000000..ac537d0a1
--- /dev/null
+++ b/rtengine/cplx_wavelet_filter_coeffs.h
@@ -0,0 +1,129 @@
+/*
+ *  This file is part of RawTherapee.
+ *
+ *  RawTherapee 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.
+ * 
+ *  RawTherapee 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 RawTherapee.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *  2012 Emil Martinec <ejmartin@uchicago.edu>
+ */
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+//#include "improcfun.h"
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+
+namespace rtengine {
+
+
+// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
+	
+/*const int AntonB_len = 12;//length of filter
+const int AntonB_offset = 6;//offset
+	
+const float AntonB_anal[2][2][12] = {//analysis filter
+		{{0, -0.08838834764832, 0.08838834764832, 0.69587998903400, 0.69587998903400, 
+			0.08838834764832, -0.08838834764832, 0.01122679215254, 0.01122679215254, 0},
+		{0, 0, 0, 0.04563588155712, -0.02877176311425, -0.29563588155712 , 
+			0.55754352622850, -0.29563588155713, -0.02877176311425, 0.04563588155712, 0, 0}},
+		{{0 , 0 , 0.02674875741081, -0.01686411844287, -0.07822326652899, 0.26686411844288, 
+			0.60294901823636, 0.26686411844287, -0.07822326652899, -0.01686411844287, 0.02674875741081, 0},
+		{0 , 0 , 0, 0 , 0.04563588155712, -0.02877176311425, 
+			-0.29563588155712 , 0.55754352622850, -0.29563588155713, -0.02877176311425, 0.04563588155712 , 0}} };
+	
+const float AntonB_synth[2][2][12] = {//synthesis filter
+		{{0 , 0 , 0, -0.04563588155712, -0.02877176311425, 0.29563588155712, 
+			0.55754352622850, 0.29563588155713, -0.02877176311425, -0.04563588155712, 0, 0},
+		{0, 0.02674875741081, 0.01686411844287, -0.07822326652899, -0.26686411844288 , 0.60294901823636, 
+			-0.26686411844287, -0.07822326652899, 0.01686411844287, 0.02674875741081, 0, 0}},
+		{{0 , 0, -0.04563588155712, -0.02877176311425, 0.29563588155712 , 0.55754352622850 , 
+			0.29563588155713, -0.02877176311425, -0.04563588155712, 0, 0 , 0},
+		{0.02674875741081 , 0.01686411844287, -0.07822326652899, -0.26686411844288 , 0.60294901823636, -0.26686411844287, 
+			-0.07822326652899, 0.01686411844287 , 0.02674875741081 , 0 , 0, 0}} };*/
+
+	
+/*	for (int i=0; i<4; i++)
+		for (int n=0; n<12; n++) {
+			AntonB.synth[i][n] *= 2;
+		}*/
+
+// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
+
+
+const int FSFarras_len=10;//length of filter
+const int FSFarras_offset=4;//offset
+
+const float FSFarras_anal[2][2][10] = {//analysis filter
+	{{0, -0.08838834764832, 0.08838834764832, 0.69587998903400, 0.69587998903400, 0.08838834764832, -0.08838834764832, 0.01122679215254 , 0.01122679215254, 0},
+		{ 0, -0.01122679215254, 0.01122679215254, 0.08838834764832, 0.08838834764832, -0.69587998903400, 0.69587998903400, -0.08838834764832, -0.08838834764832, 0}},
+	{{0.01122679215254, 0.01122679215254, -0.08838834764832, 0.08838834764832, 0.69587998903400, 0.69587998903400, 0.08838834764832, -0.08838834764832, 0, 0},
+		{0, 0, -0.08838834764832, -0.08838834764832, 0.69587998903400, -0.69587998903400, 0.08838834764832, 0.08838834764832, 0.01122679215254, -0.01122679215254}} };
+
+//synthesis filter is the reverse (see cplx_wavelet_dec.h)
+
+
+// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%	
+
+/*
+ % Kingsbury Q-filters for the dual-tree complex DWT
+ %
+ % af{i},i=1,2-analysis filters for tree i
+ % sf{i},i=1,2-synthesis filters for tree i
+ % note:af{2} is the reverse of af{1}
+ % ordering is {af[1],af[2],sf[1],sf[2]}
+ % REFERENCE:% N.G.Kingsbury,"A dual-tree complex wavelet
+ %    transform with improved orthogonality and symmetry
+ %    properties",Proceedings of the IEEE Int.Conf.on
+ % Image Proc.(ICIP),2000  */
+
+const int Kingsbury_len=10;//length of filter
+const int Kingsbury_offset=4;//offset
+
+const float Kingsbury_anal[2][2][10] =	{//analysis filter
+	{{0.03516384000000, 0, -0.08832942000000, 0.23389032000000, 0.76027237000000, 0.58751830000000, 0, -0.11430184000000 , 0, 0},
+		{ 0, 0, -0.11430184000000, 0, 0.58751830000000, -0.76027237000000, 0.23389032000000, 0.08832942000000, 0, -0.03516384000000}},
+	{{0, 0, -0.11430184000000, 0, 0.58751830000000, 0.76027237000000, 0.23389032000000, -0.08832942000000, 0, 0.03516384000000},
+		{-0.03516384000000, 0, 0.08832942000000, 0.23389032000000, -0.76027237000000, 0.58751830000000, 0, -0.11430184000000, 0, 0}} };
+
+//synthesis filter is the reverse (see cplx_wavelet_dec.h)
+
+/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
+
+const int Haar_len=2;//length of filter
+const int Haar_offset=1;//offset
+
+const float Haar_anal[2][2] =	{{0.5,0.5}, {0.5,-0.5}};//analysis filter
+
+//synthesis filter is the reverse (see cplx_wavelet_dec.h)
+
+/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
+
+const int LeGall_len=6;
+const int LeGall_offset=2;
+const float LeGall_anal[2][6] = {{0, 0.25, 0.5, 0.25, 0, 0}, {0, -0.125, -0.25, 0.75, -0.25, -0.125}};
+const float LeGall_synth[2][6] = {{-0.125, 0.25, 0.75, 0.25, -0.125, 0}, {0, 0, -0.25, 0.5, -0.25, 0}};
+	
+/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
+
+const int Daub4_len=6;
+const int Daub4_offset=2;
+const float Daub4_anal[2][6] = {//analysis filter
+	{0, 0, 0.34150635, 0.59150635, 0.15849365, -0.091506351}, 
+	{-0.091506351, -0.15849365, 0.59150635, -0.34150635, 0, 0}};
+
+};
+
diff --git a/rtengine/cplx_wavelet_level.h b/rtengine/cplx_wavelet_level.h
new file mode 100644
index 000000000..756b5ce8f
--- /dev/null
+++ b/rtengine/cplx_wavelet_level.h
@@ -0,0 +1,603 @@
+/*
+ *  This file is part of RawTherapee.
+ *
+ *  RawTherapee 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.
+ * 
+ *  RawTherapee 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 RawTherapee.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ *  2010 Ilya Popov <ilia_popov@rambler.ru>
+ *  2012 Emil Martinec <ejmartin@uchicago.edu>
+ */
+
+#ifndef CPLX_WAVELET_LEVEL_H_INCLUDED
+#define CPLX_WAVELET_LEVEL_H_INCLUDED
+
+#include <cstddef>
+#include <algorithm>
+
+#include "gauss.h"
+#include "rt_math.h"
+
+namespace rtengine {
+
+	
+	//////////////////////////////////////////////////////////////////////////////
+	
+	template<typename T>
+	class wavelet_level
+	{
+		// full size
+		size_t m_w, m_h;
+		
+		// size of low frequency part
+		size_t m_w2, m_h2;
+		
+		// size of padded border
+		size_t m_pad;
+		
+		// level of decomposition
+		int lvl;
+		
+		// whether to subsample the output
+		bool subsamp_out;
+		
+		// spacing of filter taps
+		size_t skip;
+		
+		// allocation and destruction of data storage
+		T ** create(size_t n);
+		void destroy(T ** subbands);
+		
+		// load a row/column of input data, possibly with padding
+		template<typename E>
+		void loadbuffer(E * src, E * dst, int srclen, int pitch);
+		
+		void AnalysisFilter (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, 
+							 int taps, int offset, int pitch, int srclen);
+		void SynthesisFilter (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, 
+							  float *filterLo, float *filterHi, int taps, int offset, int pitch, int dstlen);
+		
+		void AnalysisFilterHaar (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen);
+		void SynthesisFilterHaar (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, int pitch, int dstlen);
+		
+		void AnalysisFilterSubsamp (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, 
+							 int taps, int offset, int pitch, int srclen);
+		void SynthesisFilterSubsamp (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, 
+							  float *filterLo, float *filterHi, int taps, int offset, int pitch, int dstlen);
+		
+		void AnalysisFilterSubsampHaar (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen);
+		void SynthesisFilterSubsampHaar (T * srcLo, T * srcHi, T * dst, int pitch, int dstlen);
+		
+		void imp_nr (T* src, int width, int height, double thresh);
+
+		
+	public:
+		
+		T ** wavcoeffs;
+		
+		template<typename E>
+		wavelet_level(E * src, int level, int subsamp, int padding, size_t w, size_t h, float *filterV, float *filterH, int len, int offset)
+		: m_w(w), m_h(h), m_w2(w), m_h2(h), m_pad(padding), wavcoeffs(NULL), lvl(level), skip(1<<level), subsamp_out((subsamp>>level)&1)
+		{
+			if (subsamp) {
+				skip = 1;
+				for (int n=0; n<level; n++) {
+					skip *= 2-((subsamp>>n)&1);
+				}
+			}
+			m_w2 = (subsamp_out ? ((w+1+2*skip*padding)/2) : (w+2*skip*padding));
+			m_h2 = (subsamp_out ? ((h+1+2*skip*padding)/2) : (h+2*skip*padding));
+			m_pad= skip*padding;
+			
+			wavcoeffs = create((m_w2)*(m_h2));
+			decompose_level(src, filterV, filterH, len, offset);
+			
+		}
+		
+		~wavelet_level()
+		{
+			destroy(wavcoeffs);
+		}
+		
+		T ** subbands() const
+		{
+			return wavcoeffs;
+		}
+		
+		T * lopass() const
+		{
+			return wavcoeffs[0];
+		}
+		
+		size_t width() const
+		{
+			return m_w2;
+		}
+		
+		size_t height() const
+		{
+			return m_h2;
+		}
+		
+		size_t padding() const
+		{
+			return m_pad/skip;
+		}
+		
+		size_t stride() const
+		{
+			return skip;
+		}
+		
+		template<typename E>
+		void decompose_level(E *src, float *filterV, float *filterH, int len, int offset);
+		
+		template<typename E>
+		void reconstruct_level(E *dst, float *filterV, float *filterH, int len, int offset);
+		
+	};
+	
+	//////////////////////////////////////////////////////////////////////////////
+	
+	
+	
+	template<typename T>
+	T ** wavelet_level<T>::create(size_t n)
+	{
+		T * data = new T[4*n];
+		T ** subbands = new T*[4];
+		for(size_t j = 0; j < 4; j++)
+		{
+			subbands[j] = data + n * j;
+		}
+		return subbands;
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	
+	
+	template<typename T>
+	void wavelet_level<T>::destroy(T ** subbands)
+	{
+		if(subbands)
+		{
+			delete[] subbands[0];
+			delete[] subbands;
+		}
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+
+	template<typename T> template<typename E>
+	void wavelet_level<T>::loadbuffer(E * src, E * dst, int pitch, int srclen)
+	{
+		E * tmp = dst + m_pad;
+		memset(dst, 0, (MAX(m_w2,m_h2))*sizeof(E));
+		
+		//create padded buffer from src data
+		for (size_t i = 0, j = 0; i<srclen; i++, j += pitch)
+		{
+			tmp[i] = src[j];
+		}	
+		
+		// extend each coset mirror-like by padding amount 'm_pad' and to a multiple of 'skip'
+		for (size_t i=1; i<=MIN(srclen-1,m_pad); i++) {
+			tmp[-i] = tmp[i];
+			tmp[srclen+i-1] = tmp[srclen-i-1];
+		}
+		for (size_t i=0; i<srclen%skip; i++) {
+			tmp[srclen+m_pad+i] = tmp[srclen+m_pad-i-2];
+		}
+		
+	}
+	
+	
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	template<typename T>
+	void wavelet_level<T>::AnalysisFilter (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, 
+												int taps, int offset, int pitch, int srclen) {
+		
+		/* Basic convolution code
+		 * Applies an FIR filter 'filter' with filter length 'taps', 
+		 * aligning the 'offset' element of the filter with
+		 * the input pixel, and skipping 'skip' pixels between taps 
+		 *
+		 */
+		
+				
+		for (size_t i = 0; i < (srclen); i++) {
+			float lo=0,hi=0;
+			if (i>skip*taps && i<srclen-skip*taps) {//bulk
+				for (int j=0, l=-skip*offset; j<taps; j++, l+=skip) {
+					lo += filterLo[j] * srcbuffer[i-l];//lopass channel
+					hi += filterHi[j] * srcbuffer[i-l];//hipass channel
+				}
+			} else {//boundary
+				for (int j=0; j<taps; j++) {
+					int arg = MAX(0,MIN(i+skip*(offset-j),srclen-1));//clamped BC's
+					lo += filterLo[j] * srcbuffer[arg];//lopass channel
+					hi += filterHi[j] * srcbuffer[arg];//hipass channel
+				}
+			}
+			
+			dstLo[(pitch*(i))] = lo;
+			dstHi[(pitch*(i))] = hi;
+		}
+		
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+
+	template<typename T>
+	void wavelet_level<T>::SynthesisFilter (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, float *filterLo, 
+												 float *filterHi, int taps, int offset, int pitch, int dstlen) {
+		
+		/* Basic convolution code
+		 * Applies an FIR filter 'filter' with filter length 'taps', 
+		 * aligning the 'offset' element of the filter with
+		 * the input pixel, and skipping 'skip' pixels between taps 
+		 *
+		 */
+		
+		
+		// load into buffer
+		
+		int srclen = (dstlen==m_w ? m_w2 : m_h2);//length of row/col in src (coarser level)
+		
+		for (size_t i=0, j=0; i<srclen; i++, j+=pitch) {
+			bufferLo[i]=srcLo[j];
+			bufferHi[i]=srcHi[j];
+		}
+		
+		int shift=skip*(taps-offset-1);
+		for(size_t i = m_pad; i < (dstlen+m_pad); i++) {
+			float tot=0;
+			if (i>skip*taps && i<(srclen-skip*taps)) {//bulk
+				for (int j=0, l=-shift; j<taps; j++, l+=skip) {
+					tot += (filterLo[j] * bufferLo[i-l] + filterHi[j] * bufferHi[i-l]);
+				}
+			} else {//boundary
+				for (int j=0, l=-shift; j<taps; j++, l+=skip) {
+					int arg = MAX(0,MIN((i-l),srclen-1));//clamped BC's
+					tot += (filterLo[j] * bufferLo[arg] + filterHi[j] * bufferHi[arg]);
+				}
+			}
+			
+			dst[pitch*(i-m_pad)] = tot;
+
+		}
+		
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	template<typename T>
+	void wavelet_level<T>::AnalysisFilterHaar (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen) {
+		
+		/* Basic convolution code
+		 * Applies a Haar filter 
+		 *
+		 */										
+		
+		for(size_t i = 0; i < (srclen - skip); i++) {
+			dstLo[(pitch*(i))] = 0.5*(srcbuffer[i] + srcbuffer[i+skip]);
+			dstHi[(pitch*(i))] = 0.5*(srcbuffer[i] - srcbuffer[i+skip]);
+		}
+		
+		for(size_t i = (srclen-skip); i < (srclen); i++) {
+			dstLo[(pitch*(i))] = 0.5*(srcbuffer[i] + srcbuffer[i-skip]);
+			dstHi[(pitch*(i))] = 0.5*(srcbuffer[i] - srcbuffer[i-skip]);
+		}
+																		 
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	
+	template<typename T>
+	void wavelet_level<T>::SynthesisFilterHaar (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, int pitch, int dstlen) {
+		
+		/* Basic convolution code
+		 * Applies a Haar filter 
+		 *
+		 */
+		
+		int srclen = (dstlen==m_w ? m_w2 : m_h2);//length of row/col in src (coarser level)
+		
+		for (size_t i=0, j=0; i<srclen; i++, j+=pitch) {
+			bufferLo[i]=srcLo[j];
+			bufferHi[i]=srcHi[j];
+		}
+		
+		for(size_t i = m_pad+skip; i < (dstlen+m_pad); i++) {
+			dst[pitch*(i-m_pad)] = 0.5*(bufferLo[i] + bufferHi[i] + bufferLo[i-skip] - bufferHi[i-skip]);			
+		}
+		
+		for(size_t i = (m_pad); i < (m_pad+skip); i++) {
+			dst[pitch*(i-m_pad)] = (bufferLo[i] + bufferHi[i]);			
+		}
+	}
+
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+
+
+	template<typename T>
+	void wavelet_level<T>::AnalysisFilterSubsamp (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, 
+												int taps, int offset, int pitch, int srclen) {
+		
+		/* Basic convolution code
+		 * Applies an FIR filter 'filter' with filter length 'taps', 
+		 * aligning the 'offset' element of the filter with
+		 * the input pixel, and skipping 'skip' pixels between taps 
+		 * Output is subsampled by two
+		 */
+		
+		// calculate coefficients
+		
+		for(int i = 0; i < (srclen); i+=2) {
+			float lo=0,hi=0;
+			if (i>skip*taps && i<srclen-skip*taps) {//bulk
+				for (int j=0, l=-skip*offset; j<taps; j++, l+=skip) {
+					lo += filterLo[j] * srcbuffer[i-l];//lopass channel
+					hi += filterHi[j] * srcbuffer[i-l];//hipass channel
+				}
+			} else {//boundary
+				for (int j=0; j<taps; j++) {
+					int arg = MAX(0,MIN(i+skip*(offset-j),srclen-1));//clamped BC's
+					lo += filterLo[j] * srcbuffer[arg];//lopass channel
+					hi += filterHi[j] * srcbuffer[arg];//hipass channel
+				}
+			}
+			
+			dstLo[(pitch*(i/2))] = lo;
+			dstHi[(pitch*(i/2))] = hi;
+		}
+		
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	template<typename T>
+	void wavelet_level<T>::SynthesisFilterSubsamp (T * srcLo, T * srcHi, T * dst, T *bufferLo, T *bufferHi, 
+												 float *filterLo, float *filterHi, int taps, int offset, int pitch, int dstlen) {
+		
+		/* Basic convolution code
+		 * Applies an FIR filter 'filter' with filter length 'taps', 
+		 * aligning the 'offset' element of the filter with
+		 * the input pixel, and skipping 'skip' pixels between taps 
+		 * Output is subsampled by two
+		 */
+		
+		
+		
+		// calculate coefficients
+		
+		int srclen = (dstlen==m_w ? m_w2 : m_h2);//length of row/col in src (coarser level)
+		
+		//fill a buffer with a given row/column of data
+		for (size_t i=0, j=0; i<srclen; i++, j+=pitch) {
+			bufferLo[i]=srcLo[j];
+			bufferHi[i]=srcHi[j];
+		}
+				
+		int shift=skip*(taps-offset-1);//align filter with data
+		for(size_t i = m_pad; i < (dstlen+m_pad); i++) {
+			
+			float tot=0;
+			//TODO: this is correct only if skip=1; otherwise, want to work with cosets of length 'skip'
+			int i_src = (i+shift)/2;
+			int begin = (i+shift)%2;
+			if (i>skip*taps && i<(srclen-skip*taps)) {//bulk
+				for (int j=begin, l=0; j<taps; j+=2, l+=skip) {
+					tot += 2*((filterLo[j] * bufferLo[i_src-l] + filterHi[j] * bufferHi[i_src-l]));
+				}
+			} else {//boundary
+				for (int j=begin, l=0; j<taps; j+=2, l+=skip) {
+					int arg = MAX(0,MIN((i_src-l),srclen-1));//clamped BC's
+					tot += 2*((filterLo[j] * bufferLo[arg] + filterHi[j] * bufferHi[arg]));
+				}
+			}
+			
+			dst[pitch*(i-m_pad)] = tot;
+			
+		}
+		
+	}
+	
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	
+	
+	template<typename T>
+	void wavelet_level<T>::AnalysisFilterSubsampHaar (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen) {
+		
+		/* Basic convolution code
+		 * Applies a Haar filter
+		 * Output is subsampled by two
+		 */
+		
+		// calculate coefficients
+		
+		for(size_t i = 0; i < (srclen - skip); i+=2) {
+			dstLo[(pitch*(i/2))] = 0.5*(srcbuffer[i] + srcbuffer[i+skip]);
+			dstHi[(pitch*(i/2))] = 0.5*(srcbuffer[i] - srcbuffer[i+skip]);
+		}
+		
+		for(size_t i = (srclen-skip)-((srclen-skip)&1); i < (srclen); i+=2) {
+			dstLo[(pitch*(i/2))] = 0.5*(srcbuffer[i] + srcbuffer[i-skip]);
+			dstHi[(pitch*(i/2))] = 0.5*(srcbuffer[i] - srcbuffer[i-skip]);
+		}
+		
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+	template<typename T>
+	void wavelet_level<T>::SynthesisFilterSubsampHaar (T * srcLo, T * srcHi, T * dst, int pitch, int dstlen) {
+		
+		/* Basic convolution code
+		 * Applies a Haar filter 
+		 * Input was subsampled by two
+		 */
+		
+		
+		// calculate coefficients
+				
+		//TODO: this code is buggy...
+		for (int n=0; n<skip; n++) {
+			for (size_t i = m_pad; i < (dstlen+m_pad-2*skip); i+=2*skip) {
+				dst[pitch*(i-m_pad+n)] = srcLo[i/2+n]+srcHi[i/2+n];
+				dst[pitch*(i-m_pad+skip+n)] = srcLo[i/2+n]-srcHi[i/2+n];
+			}
+		}
+		
+		if ((dstlen+m_pad-2*skip)<dstlen-1) {
+			for (int n=0; n<skip; n++) {
+				for (size_t i=(dstlen+m_pad-2*skip); i<dstlen+m_pad-1; i++) {
+					dst[pitch*(dstlen-m_pad+n)] = srcLo[i/2+n]+srcHi[i/2+n];
+				}
+			}
+		}
+		
+		
+	}
+	
+
+
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	
+	
+	template<typename T> template<typename E>
+	void wavelet_level<T>::decompose_level(E *src, float *filterV, float *filterH, int taps, int offset) { 
+		
+		T *tmpLo = new T[m_w*m_h2];
+		T *tmpHi = new T[m_w*m_h2];
+		
+		T *buffer = new T[MAX(m_w,m_h)+2*m_pad+skip];
+		
+		/* filter along columns */
+//OpenMP here		
+		for (int j=0; j<m_w; j++) {
+			loadbuffer(src+j, buffer, m_w/*pitch*/, m_h/*srclen*/);//pad a column of data and load it to buffer
+			if (subsamp_out) {
+				AnalysisFilterSubsamp (buffer, tmpLo+j, tmpHi+j, filterV, filterV+taps, taps, offset, m_w/*output_pitch*/, m_h/*srclen*/);
+				//AnalysisFilterSubsampHaar (buffer, tmpLo+j, tmpHi+j, m_w, m_h);
+			} else {
+				//AnalysisFilter (buffer, tmpLo+j, tmpHi+j, filterV, filterV+taps, taps, offset, m_w/*output_pitch*/, m_h/*srclen*/);
+				AnalysisFilterHaar (buffer, tmpLo+j, tmpHi+j, m_w, m_h);
+			}
+		}
+		
+		/* filter along rows */
+//OpenMP here		
+		for (int i=0; i<m_h2; i++) {
+			loadbuffer(tmpLo+i*m_w, buffer, 1/*pitch*/, m_w/*srclen*/);//pad a row of data and load it to buffer
+			if (subsamp_out) {
+				AnalysisFilterSubsamp (buffer, wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
+				//AnalysisFilterSubsampHaar (buffer, wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, 1, m_w);
+			} else {
+				//AnalysisFilter (buffer, wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
+				AnalysisFilterHaar (buffer, wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, 1, m_w);
+			}
+			loadbuffer(tmpHi+i*m_w, buffer, 1/*pitch*/, m_w/*srclen*/);
+			if (subsamp_out) {
+				AnalysisFilterSubsamp (buffer, wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
+				//AnalysisFilterSubsampHaar (buffer, wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, 1, m_w);
+			} else {
+				//AnalysisFilter (buffer, wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, filterH, filterH+taps, taps, offset, 1/*output_pitch*/, m_w/*srclen*/);
+				AnalysisFilterHaar (buffer, wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, 1, m_w);
+			}
+		}
+				
+		delete[] tmpLo;
+		delete[] tmpHi;
+		delete[] buffer;
+	}
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+	
+	template<typename T> template<typename E>
+	void wavelet_level<T>::reconstruct_level(E *dst, float *filterV, float *filterH, int taps, int offset) { 
+		
+		T *tmpLo = new T[m_w*m_h2];
+		T *tmpHi = new T[m_w*m_h2];
+		
+		int buflen = MAX(m_w2,m_h2);
+		float *bufferLo = new float[buflen];
+		float *bufferHi = new float[buflen];
+		
+		/* filter along rows */
+//OpenMP here		
+		for (int i=0; i<m_h2; i++) {
+			
+			if (subsamp_out) {
+				SynthesisFilterSubsamp (wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, tmpLo+i*m_w, bufferLo, bufferHi,  \
+										filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
+				SynthesisFilterSubsamp (wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, tmpHi+i*m_w, bufferLo, bufferHi,  \
+										filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
+				//SynthesisFilterSubsampHaar (wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, tmpLo+i*m_w, 1, m_w);//TODO: this is buggy
+				//SynthesisFilterSubsampHaar (wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, tmpHi+i*m_w, 1, m_w);
+			} else {
+				//SynthesisFilter (wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, tmpLo+i*m_w, bufferLo, bufferHi, \
+								 filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
+				//SynthesisFilter (wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, tmpHi+i*m_w, bufferLo, bufferHi, \
+								 filterH, filterH+taps, taps, offset, 1/*pitch*/, m_w/*dstlen*/);
+				SynthesisFilterHaar (wavcoeffs[0]+i*m_w2, wavcoeffs[1]+i*m_w2, tmpLo+i*m_w, bufferLo, bufferHi, 1, m_w);
+				SynthesisFilterHaar (wavcoeffs[2]+i*m_w2, wavcoeffs[3]+i*m_w2, tmpHi+i*m_w, bufferLo, bufferHi, 1, m_w);
+			}
+		}
+		
+		/* filter along columns */
+//OpenMP here		
+		for (int j=0; j<m_w; j++) {
+			if (subsamp_out) {
+				SynthesisFilterSubsamp (tmpLo+j, tmpHi+j, dst+j, bufferLo, bufferHi, \
+										filterV, filterV+taps, taps, offset, m_w/*pitch*/, m_h/*dstlen*/);
+				//SynthesisFilterSubsampHaar (tmpLo+j, tmpHi+j, dst+j, m_w, m_h);
+			} else {
+				//SynthesisFilter (tmpLo+j, tmpHi+j, dst+j, bufferLo, bufferHi, \
+								 filterV, filterV+taps, taps, offset, m_w/*pitch*/, m_h/*dstlen*/);
+				SynthesisFilterHaar (tmpLo+j, tmpHi+j, dst+j, bufferLo, bufferHi, m_w, m_h);
+			}
+		}
+		
+		
+		delete[] tmpLo;
+		delete[] tmpHi;
+		delete[] bufferLo;
+		delete[] bufferHi;
+		
+	}
+
+	
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	
+};
+
+#endif
diff --git a/rtengine/dcrop.cc b/rtengine/dcrop.cc
index ebfe0a5fe..4221b35d4 100644
--- a/rtengine/dcrop.cc
+++ b/rtengine/dcrop.cc
@@ -86,8 +86,8 @@ void Crop::update (int todo) {
 
     bool needstransform  = parent->ipf.needsTransform();
 
-    if (todo & M_INIT) {
-        Glib::Mutex::Lock lock(parent->minit);  // Also used in inproccoord
+    if (todo & (M_INIT|M_LINDENOISE)) {
+        Glib::Mutex::Lock lock(parent->minit);  // Also used in improccoord
 
         int tr = TR_NONE;
         if (params.coarse.rotate==90)  tr |= TR_R90;
@@ -101,7 +101,16 @@ void Crop::update (int todo) {
             setCropSizes (rqcropx, rqcropy, rqcropw, rqcroph, skip, true);
         PreviewProps pp (trafx, trafy, trafw*skip, trafh*skip, skip);
         parent->imgsrc->getImage (parent->currWB, tr, origCrop, pp, params.hlrecovery, params.icm, params.raw );
-    }
+
+		//parent->imgsrc->convertColorSpace(origCrop, params.icm);
+
+        if (todo & M_LINDENOISE) {
+			if (skip==1 && params.dirpyrDenoise.enabled) {
+				parent->ipf.RGB_denoise(origCrop, origCrop, /*Roffset,*/ params.dirpyrDenoise, params.defringe);
+			}
+        }
+        parent->imgsrc->convertColorSpace(origCrop, params.icm, params.raw);
+}
 
     // transform
     if ((!needstransform && transCrop) || (transCrop && (transCrop->width!=cropw || transCrop->height!=croph))) {
@@ -182,7 +191,6 @@ void Crop::update (int todo) {
 		if (skip==1) {
 			parent->ipf.impulsedenoise (labnCrop);
 			parent->ipf.defringe (labnCrop);
-			parent->ipf.dirpyrdenoise (labnCrop);
 			parent->ipf.MLsharpen (labnCrop);
 			parent->ipf.MLmicrocontrast (labnCrop);
 			//parent->ipf.MLmicrocontrast (labnCrop);
diff --git a/rtengine/demosaic_algos.cc b/rtengine/demosaic_algos.cc
index 47a34c8d5..202da7b11 100644
--- a/rtengine/demosaic_algos.cc
+++ b/rtengine/demosaic_algos.cc
@@ -57,18 +57,18 @@ void RawImageSource::eahd_demosaic () {
   }
 
   // prepare cache and constants for cielab conversion
-	//TODO: revisit after conversion to D50 illuminant
-	lc00 = (0.412453 * rgb_cam[0][0] + 0.357580 * rgb_cam[0][1] + 0.180423 * rgb_cam[0][2]) ;// / 0.950456;
-  lc01 = (0.412453 * rgb_cam[1][0] + 0.357580 * rgb_cam[1][1] + 0.180423 * rgb_cam[1][2]) ;// / 0.950456;
-  lc02 = (0.412453 * rgb_cam[2][0] + 0.357580 * rgb_cam[2][1] + 0.180423 * rgb_cam[2][2]) ;// / 0.950456;
+  //TODO: revisit after conversion to D50 illuminant
+  lc00 = (0.412453 * imatrices.rgb_cam[0][0] + 0.357580 * imatrices.rgb_cam[0][1] + 0.180423 * imatrices.rgb_cam[0][2]) ;// / 0.950456;
+  lc01 = (0.412453 * imatrices.rgb_cam[1][0] + 0.357580 * imatrices.rgb_cam[1][1] + 0.180423 * imatrices.rgb_cam[1][2]) ;// / 0.950456;
+  lc02 = (0.412453 * imatrices.rgb_cam[2][0] + 0.357580 * imatrices.rgb_cam[2][1] + 0.180423 * imatrices.rgb_cam[2][2]) ;// / 0.950456;
 
-  lc10 = 0.212671 * rgb_cam[0][0] + 0.715160 * rgb_cam[0][1] + 0.072169 * rgb_cam[0][2];
-  lc11 = 0.212671 * rgb_cam[1][0] + 0.715160 * rgb_cam[1][1] + 0.072169 * rgb_cam[1][2];
-  lc12 = 0.212671 * rgb_cam[2][0] + 0.715160 * rgb_cam[2][1] + 0.072169 * rgb_cam[2][2];
+  lc10 = 0.212671 * imatrices.rgb_cam[0][0] + 0.715160 * imatrices.rgb_cam[0][1] + 0.072169 * imatrices.rgb_cam[0][2];
+  lc11 = 0.212671 * imatrices.rgb_cam[1][0] + 0.715160 * imatrices.rgb_cam[1][1] + 0.072169 * imatrices.rgb_cam[1][2];
+  lc12 = 0.212671 * imatrices.rgb_cam[2][0] + 0.715160 * imatrices.rgb_cam[2][1] + 0.072169 * imatrices.rgb_cam[2][2];
 
-  lc20 = (0.019334 * rgb_cam[0][0] + 0.119193 * rgb_cam[0][1] + 0.950227 * rgb_cam[0][2]) ;// / 1.088754;
-  lc21 = (0.019334 * rgb_cam[1][0] + 0.119193 * rgb_cam[1][1] + 0.950227 * rgb_cam[1][2]) ;// / 1.088754;
-  lc22 = (0.019334 * rgb_cam[2][0] + 0.119193 * rgb_cam[2][1] + 0.950227 * rgb_cam[2][2]) ;// / 1.088754;
+  lc20 = (0.019334 * imatrices.rgb_cam[0][0] + 0.119193 * imatrices.rgb_cam[0][1] + 0.950227 * imatrices.rgb_cam[0][2]) ;// / 1.088754;
+  lc21 = (0.019334 * imatrices.rgb_cam[1][0] + 0.119193 * imatrices.rgb_cam[1][1] + 0.950227 * imatrices.rgb_cam[1][2]) ;// / 1.088754;
+  lc22 = (0.019334 * imatrices.rgb_cam[2][0] + 0.119193 * imatrices.rgb_cam[2][1] + 0.950227 * imatrices.rgb_cam[2][2]) ;// / 1.088754;
 
   int maxindex = 2*65536;
   cache = new double[maxindex];
@@ -530,7 +530,7 @@ void RawImageSource::hphd_demosaic () {
 #endif
 
   hphd_green (hpmap);
-  freeArray<float>(hpmap, H);
+  freeArray<float>(hpmap, H);//TODO: seems to cause sigabrt ???  why???
 
   if (plistener)
     plistener->setProgress (0.66);
@@ -820,22 +820,19 @@ void RawImageSource::ppg_demosaic()
     if(plistener) plistener->setProgress(0.67 + 0.33*row/(height-1));
   }
 
-  red = new float*[H];
-  for (int i=0; i<H; i++) {
-    red[i] = new float[W];
-    for (int j=0; j<W; j++)
+  red(W,H);
+  for (int i=0; i<H; i++) 
+	  for (int j=0; j<W; j++) {
         red[i][j] = image[i*W+j][0];
   }
-  green = new float*[H];
-  for (int i=0; i<H; i++) {
-    green[i] = new float[W];
-    for (int j=0; j<W; j++)
+  green(W,H);
+  for (int i=0; i<H; i++) 
+	  for (int j=0; j<W; j++) {
         green[i][j] = image[i*W+j][1];
   }
-  blue = new float*[H];
-  for (int i=0; i<H; i++) {
-    blue[i] = new float[W];
-    for (int j=0; j<W; j++)
+  blue(W,H);
+  for (int i=0; i<H; i++) 
+	  for (int j=0; j<W; j++) {
         blue[i][j] = image[i*W+j][2];
   }
   free (image);
@@ -919,7 +916,7 @@ void RawImageSource::ahd_demosaic(int winx, int winy, int winw, int winh)
     for (i=0; i < 3; i++)
         for (j=0; j < colors; j++)
             for (xyz_cam[i][j] = k=0; k < 3; k++)
-	            xyz_cam[i][j] += xyz_rgb[i][k] * rgb_cam[k][j] / d65_white[i];
+	            xyz_cam[i][j] += xyz_rgb[i][k] * imatrices.rgb_cam[k][j] / d65_white[i];
 
     border_interpolate(5, image);
     buffer = (char *) malloc (13*TS*TS*sizeof(float));		/* 1664 kB */
@@ -1055,21 +1052,15 @@ void RawImageSource::ahd_demosaic(int winx, int winy, int winw, int winh)
 
 void RawImageSource::nodemosaic()
 {
-    red = new float*[H];
-    green = new float*[H];
-    blue = new float*[H];
+    red(W,H);
+    green(W,H);
+    blue(W,H);
     for (int i=0; i<H; i++) {
-        red[i] = new float[W];
-        green[i] = new float[W];
-        blue[i] = new float[W];
         for (int j=0; j<W; j++){
         	switch( FC(i,j)){
 				case 0: red[i][j] = rawData[i][j]; green[i][j]=blue[i][j]=0; break;
 				case 1: green[i][j] = rawData[i][j]; red[i][j]=blue[i][j]=0; break;
 				case 2: blue[i][j] = rawData[i][j]; red[i][j]=green[i][j]=0; break;
-					//red[i][j] = rawData[i][j];
-					//green[i][j] = rawData[i][j];
-					//blue[i][j] = rawData[i][j];
         	}
         }
     }
diff --git a/rtengine/gamutbdy.cc b/rtengine/gamutbdy.cc
new file mode 100644
index 000000000..2e3ba23d9
--- /dev/null
+++ b/rtengine/gamutbdy.cc
@@ -0,0 +1 @@
+// File will be deleted on commit
\ No newline at end of file
diff --git a/rtengine/gauss.h b/rtengine/gauss.h
index a98cf67eb..bf9cc271d 100644
--- a/rtengine/gauss.h
+++ b/rtengine/gauss.h
@@ -217,7 +217,193 @@ template<class T> void gaussVertical (T** src, T** dst, AlignedBufferMP<double>
         double temp2Hm1 = src[H-1][i] + M[0][0]*(temp2[H-1] - src[H-1][i]) + M[0][1]*(temp2[H-2] - src[H-1][i]) + M[0][2]*(temp2[H-3] - src[H-1][i]);
         double temp2H   = src[H-1][i] + M[1][0]*(temp2[H-1] - src[H-1][i]) + M[1][1]*(temp2[H-2] - src[H-1][i]) + M[1][2]*(temp2[H-3] - src[H-1][i]);
         double temp2Hp1 = src[H-1][i] + M[2][0]*(temp2[H-1] - src[H-1][i]) + M[2][1]*(temp2[H-2] - src[H-1][i]) + M[2][2]*(temp2[H-3] - src[H-1][i]);
+		
+        temp2[H-1] = temp2Hm1;
+        temp2[H-2] = B * temp2[H-2] + b1*temp2[H-1] + b2*temp2H + b3*temp2Hp1;
+        temp2[H-3] = B * temp2[H-3] + b1*temp2[H-2] + b2*temp2[H-1] + b3*temp2H;
+        
+        for (int j=H-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<H; j++)
+            dst[j][i] = (T)temp2[j];
+    }
+}   
 
+
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+template<class T> void gaussDerivH (T** src, T** dst, AlignedBufferMP<double> &buffer, int W, int H, double sigma) {
+
+	
+    if (sigma<0.6) {
+        // apply symmetric derivative
+#ifdef _OPENMP
+#pragma omp for
+#endif
+		for (int i=0; i<H; i++) {
+			AlignedBuffer<double>* pBuf = buffer.acquire();
+			T* temp = (T*)pBuf->data;
+			// double* temp = buffer->data;// replaced by 2 lines above
+			for (int j=1; j<W-1; j++)
+				temp[j] = (0.5 * (src[i][j+1] - src[i][j-1]) );
+			dst[i][0] = (src[i][1]-src[i][0]);
+			//memcpy (dst[i]+1, temp+1, (W-2)*sizeof(T));
+			for (int j=1; j<W-1; j++)
+				dst[i][j] = temp[j];
+			
+			buffer.release(pBuf);
+			dst[i][W-1] = (src[i][W-1]-src[i][W-2]);
+		}
+        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);
+
+#pragma omp for
+    for (int i=0; i<H; i++) {
+		AlignedBuffer<double>* pBuf = buffer.acquire();
+		T* temp2 = (T*)pBuf->data;
+		// double* temp2 = buffer->data;// replaced by 2 lines above
+		
+		double src0 = (src[i][1]-src[i][0]);
+		
+        temp2[0] = B * src0 + b1*src0 + b2*src0 + b3*src0;
+        temp2[1] = B * 0.5*(src[i][2]-src[i][0]) + b1*temp2[0]  + b2*src0 + b3*src0;
+        temp2[2] = B * 0.5*(src[i][3]-src[i][1]) + b1*temp2[1]  + b2*temp2[0]  + b3*src0;
+		
+        for (int j=3; j<W-1; j++)
+            temp2[j] = B * 0.5*(src[i][j+1]-src[i][j-1]) + b1*temp2[j-1] + b2*temp2[j-2] + b3*temp2[j-3];
+		
+		double srcWm1 = (src[i][W-1]-src[i][W-2]);
+
+		temp2[W-1] = B * srcWm1 + b1*temp2[W-2] + b2*temp2[W-3] + b3*temp2[W-4];
+		
+        double temp2Wm1 = srcWm1 + M[0][0]*(temp2[W-1] - srcWm1) + M[0][1]*(temp2[W-2] - srcWm1) + M[0][2]*(temp2[W-3] - srcWm1);
+        double temp2W   = srcWm1 + M[1][0]*(temp2[W-1] - srcWm1) + M[1][1]*(temp2[W-2] - srcWm1) + M[1][2]*(temp2[W-3] - srcWm1);
+        double temp2Wp1 = srcWm1 + M[2][0]*(temp2[W-1] - srcWm1) + M[2][1]*(temp2[W-2] - srcWm1) + M[2][2]*(temp2[W-3] - srcWm1);
+		
+        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];
+            
+        buffer.release(pBuf);
+    }
+}
+
+//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+template<class T> void gaussDerivV (T** src, T** dst, AlignedBufferMP<double> &buffer, int W, int H, double sigma) {
+
+	if (sigma<0.6) {
+        // apply symmetric derivative
+#ifdef _OPENMP
+#pragma omp for
+#endif
+		for (int j=0; j<W; j++) {
+			AlignedBuffer<double>* pBuf = buffer.acquire();
+			T* temp = (T*)pBuf->data;
+			// double* temp = buffer->data;// replaced by 2 lines above
+			for (int i = 1; i<H-1; i++) 
+				temp[i] = (0.5 * (src[i+1][j] - src[i-1][j]) );
+			dst[0][j] = (src[1][j]-src[0][j]);
+			for (int i=1; i<H-1; i++)
+				dst[i][j] = temp[i];
+			
+			buffer.release(pBuf);
+			
+			dst[H-1][j] = (src[H-1][j]-src[H-2][j]);
+		}
+        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);
+#ifdef _OPENMP
+#pragma omp for
+#endif
+    for (int i=0; i<W; i++) {
+		AlignedBuffer<double>* pBuf = buffer.acquire();
+		T* temp2 = (T*)pBuf->data;
+		// double* temp2 = buffer->data;// replaced by 2 lines above
+		
+		double src0 = 0.5*(src[1][i]-src[0][i]);
+		
+    	temp2[0] = B * src0 + b1*src0 + b2*src0 + b3*src0;
+        temp2[1] = B * 0.5*(src[2][i]-src[0][i]) + b1*temp2[0]  + b2*src0 + b3*src0;
+        temp2[2] = B * 0.5*(src[3][i]-src[1][i]) + b1*temp2[1]  + b2*temp2[0]  + b3*src0;
+		
+        for (int j=3; j<H-1; j++)
+            temp2[j] = B * 0.5*(src[j+1][i]-src[j-1][i]) + b1*temp2[j-1] + b2*temp2[j-2] + b3*temp2[j-3];
+		
+		double srcHm1 = 0.5*(src[H-1][i]-src[H-2][i]);
+		
+		temp2[H-1] = B * srcHm1 + b1*temp2[H-2] + b2*temp2[H-3] + b3*temp2[H-4];
+		
+        double temp2Hm1 = srcHm1 + M[0][0]*(temp2[H-1] - srcHm1) + M[0][1]*(temp2[H-2] - srcHm1) + M[0][2]*(temp2[H-3] - srcHm1);
+        double temp2H   = srcHm1 + M[1][0]*(temp2[H-1] - srcHm1) + M[1][1]*(temp2[H-2] - srcHm1) + M[1][2]*(temp2[H-3] - srcHm1);
+        double temp2Hp1 = srcHm1 + M[2][0]*(temp2[H-1] - srcHm1) + M[2][1]*(temp2[H-2] - srcHm1) + M[2][2]*(temp2[H-3] - srcHm1);
+		
         temp2[H-1] = temp2Hm1;
         temp2[H-2] = B * temp2[H-2] + b1*temp2[H-1] + b2*temp2H + b3*temp2Hp1;
         temp2[H-3] = B * temp2[H-3] + b1*temp2[H-2] + b2*temp2[H-1] + b3*temp2H;
diff --git a/rtengine/green_equil_RT.cc b/rtengine/green_equil_RT.cc
index 66ff90ffe..da3255c23 100644
--- a/rtengine/green_equil_RT.cc
+++ b/rtengine/green_equil_RT.cc
@@ -40,7 +40,8 @@ void RawImageSource::green_equilibrate(float thresh)
 	int height=H, width=W; 
 
 	// local variables
-	array2D<float> cfa (width,height,rawData);
+	float** rawptr = rawData;
+	array2D<float> cfa (width,height,rawptr);
 	//array2D<int> checker (width,height,ARRAY2D_CLEAR_DATA);
 	
 	
diff --git a/rtengine/imagesource.h b/rtengine/imagesource.h
index 528bee047..83c9f2eff 100644
--- a/rtengine/imagesource.h
+++ b/rtengine/imagesource.h
@@ -53,6 +53,15 @@ class PreviewProps {
 
 };
 
+class ImageMatrices {
+
+public:
+    double rgb_cam[3][3];
+    double cam_rgb[3][3];
+    double xyz_cam[3][3];
+    double cam_xyz[3][3];
+};
+
 class ImageSource : public InitialImage {
 
     private:
@@ -62,6 +71,7 @@ class ImageSource : public InitialImage {
         cmsHPROFILE embProfile;
         Glib::ustring fileName;
         ImageData* idata;
+        ImageMatrices imatrices;
 
     public:
                             ImageSource () : references (1), embProfile(NULL), idata(NULL) {}
@@ -77,10 +87,12 @@ class ImageSource : public InitialImage {
 
         virtual bool        IsrgbSourceModified() =0; // tracks whether cached rgb output of demosaic has been modified
 
-        // use the right after demosaicing image, add coarse transformation and put the result in the provided Imagefloat*
+        // use right after demosaicing image, add coarse transformation and put the result in the provided Imagefloat*
         virtual void        getImage    (ColorTemp ctemp, int tran, Imagefloat* image, PreviewProps pp, HRecParams hlp, ColorManagementParams cmp, RAWParams raw) {}
         // true is ready to provide the AutoWB, i.e. when the image has been demosaiced for RawImageSource
         virtual bool        isWBProviderReady () =0;
+
+		virtual void        convertColorSpace(Imagefloat* image, ColorManagementParams cmp, RAWParams raw) =0;// DIRTY HACK: this method is derived in rawimagesource and strimagesource, but (...,RAWParams raw) will be used ONLY for raw images
         virtual ColorTemp   getWB       () =0;
         virtual ColorTemp   getAutoWB   () =0;
         virtual ColorTemp   getSpotWB   (std::vector<Coord2D> red, std::vector<Coord2D> green, std::vector<Coord2D>& blue, int tran) =0;
@@ -93,7 +105,9 @@ class ImageSource : public InitialImage {
         virtual void        getSize     (int tran, PreviewProps pp, int& w, int& h) {}
         virtual int         getRotateDegree() const { return 0; }
 
-        virtual ImageData*  getImageData () =0;
+        virtual ImageData*     getImageData () =0;
+        virtual ImageMatrices* getImageMatrices () =0;
+
         virtual void        setProgressListener (ProgressListener* pl) {}
 
                 void        increaseRef () { references++; }
diff --git a/rtengine/improccoordinator.cc b/rtengine/improccoordinator.cc
index 25bda07ca..9804b60b7 100644
--- a/rtengine/improccoordinator.cc
+++ b/rtengine/improccoordinator.cc
@@ -150,7 +150,23 @@ void ImProcCoordinator::updatePreviewImage (int todo, Crop* cropCall) {
     {
 
         if (settings->verbose) printf("Demosaic %s\n",rp.dmethod.c_str());
-        imgsrc->demosaic( rp );
+		
+		//TODO - denoise branch - is this code for WB params still necessary?
+		currWB = ColorTemp (params.wb.temperature, params.wb.green, params.wb.method);
+        if (params.wb.method=="Camera")
+            currWB = imgsrc->getWB ();
+        else if (params.wb.method=="Auto") {
+            if (!awbComputed) {
+                autoWB = imgsrc->getAutoWB ();
+                awbComputed = true;
+            }
+            currWB = autoWB;
+        }
+        params.wb.temperature = currWB.getTemp ();
+        params.wb.green = currWB.getGreen ();
+        
+    	imgsrc->demosaic( rp );
+    	
         if (highDetailNeeded) {
             highDetailRawComputed = true;
             if (params.hlrecovery.enabled && params.hlrecovery.method=="Color") {
@@ -159,10 +175,18 @@ void ImProcCoordinator::updatePreviewImage (int todo, Crop* cropCall) {
         }
         else
             highDetailRawComputed = false;
+		
+		
+		LUTu aehist; int aehistcompr;
+		double clip;
+		int brightness, contrast, black, hlcompr, hlcomprthresh;
+		
+		imgsrc->getAutoExpHistogram (aehist, aehistcompr);
+		ipf.getAutoExp (aehist, aehistcompr, imgsrc->getDefGain(), clip, params.dirpyrDenoise.expcomp, brightness, contrast, black, hlcompr, hlcomprthresh);	
     }
 
 
-    if (todo & M_INIT) {
+    if (todo & (M_INIT|M_LINDENOISE)) {
         Glib::Mutex::Lock lock(minit);  // Also used in crop window
 
         imgsrc->HLRecovery_Global( params.hlrecovery ); // this handles Color HLRecovery
@@ -193,6 +217,18 @@ void ImProcCoordinator::updatePreviewImage (int todo, Crop* cropCall) {
         PreviewProps pp (0, 0, fw, fh, scale);
         setScale (scale);
         imgsrc->getImage (currWB, tr, orig_prev, pp, params.hlrecovery, params.icm, params.raw);
+
+		//imgsrc->convertColorSpace(orig_prev, params.icm, params.raw);
+
+        if (todo & M_LINDENOISE) {
+        	//printf("denoising!\n");
+			if (scale==1 && params.dirpyrDenoise.enabled) {
+				ipf.RGB_denoise(orig_prev, orig_prev, params.dirpyrDenoise, params.defringe);
+			}
+			ImageMatrices* imatrices = imgsrc->getImageMatrices ();
+        }
+        imgsrc->convertColorSpace(orig_prev, params.icm, params.raw);
+
         ipf.firstAnalysis (orig_prev, &params, vhist16, imgsrc->getGamma());
     }
     readyphase++;
@@ -306,9 +342,6 @@ void ImProcCoordinator::updatePreviewImage (int todo, Crop* cropCall) {
             readyphase++;
 			progress ("Defringing...",100*readyphase/numofphases);
             ipf.defringe (nprevl);
-            readyphase++;
-            progress ("Denoising luma/chroma...",100*readyphase/numofphases);
-            ipf.dirpyrdenoise (nprevl);
             readyphase++;
 			if (params.sharpenEdge.enabled) {
                 progress ("Edge sharpening...",100*readyphase/numofphases);
@@ -620,6 +653,7 @@ void ImProcCoordinator::saveInputICCReference (const Glib::ustring& fname) {
 	params.wb.temperature = currWB.getTemp ();
 	params.wb.green = currWB.getGreen ();
 	imgsrc->getImage (currWB, 0, im, pp, ppar.hlrecovery, ppar.icm, ppar.raw);
+	imgsrc->convertColorSpace(im, ppar.icm, params.raw);
 	im16 = im->to16();
 	im16->saveTIFF (fname,16,true);
 	//im->saveJPEG (fname, 85);
diff --git a/rtengine/improcfun.cc b/rtengine/improcfun.cc
index 41cca39c7..0ce21cd38 100644
--- a/rtengine/improcfun.cc
+++ b/rtengine/improcfun.cc
@@ -22,6 +22,7 @@
 
 #include "rtengine.h"
 #include "improcfun.h"
+#include "curves.h"
 #include "colorclip.h"
 #include "gauss.h"
 #include "bilateral2.h"
@@ -34,6 +35,8 @@
 #include "iccmatrices.h"
 #include "color.h"
 #include "calc_distort.h"
+#include "cplx_wavelet_dec.h"
+#include "boxblur.h"
 #include "rt_math.h"
 
 #ifdef _OPENMP
@@ -907,12 +910,6 @@ void ImProcFunctions::colorCurve (LabImage* lold, LabImage* lnew) {
 			PF_correct_RT(lab, lab, params->defringe.radius, params->defringe.threshold);
 	}
 	
-	void ImProcFunctions::dirpyrdenoise (LabImage* lab) {
-		
-		if (params->dirpyrDenoise.enabled && lab->W>=8 && lab->H>=8)
-			
-			dirpyrLab_denoise(lab, lab, params->dirpyrDenoise );
-	}
 	
 	void ImProcFunctions::dirpyrequalizer (LabImage* lab) {
 		
diff --git a/rtengine/improcfun.h b/rtengine/improcfun.h
index ec659e6b3..f50b4b26b 100644
--- a/rtengine/improcfun.h
+++ b/rtengine/improcfun.h
@@ -30,6 +30,8 @@
 #include "LUT.h"
 #include "lcp.h"
 #include "curves.h"
+#include <fftw3.h>
+#include "cplx_wavelet_dec.h"
 
 namespace rtengine {
 
@@ -103,22 +105,40 @@ class ImProcFunctions {
 
 		void impulsedenoise   (LabImage* lab);//Emil's impulse denoise
 		void impulse_nr (LabImage* lab, double thresh);
-		void dirpyrdenoise    (LabImage* lab);//Emil's pyramid denoise
+		void dirpyrdenoise    (LabImage* src);//Emil's pyramid denoise
 		void dirpyrequalizer  (LabImage* lab);//Emil's equalizer
-
+		
 		void EPDToneMap(LabImage *lab, unsigned int Iterates = 0, int skip = 1);
 		
+		// pyramid denoise
 		procparams::DirPyrDenoiseParams dnparams;
 		void dirpyrLab_denoise(LabImage * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams );//Emil's directional pyramid denoise
 		void dirpyr           (LabImage* data_fine, LabImage* data_coarse, int level, LUTf &rangefn_L, LUTf &rangefn_ab,
 							   int pitch, int scale, const int luma, int chroma );
 		void idirpyr          (LabImage* data_coarse, LabImage* data_fine, int level, LUTf &rangefn_L, LUTf & nrwt_l, LUTf & nrwt_ab,
 							   int pitch, int scale, const int luma, const int chroma/*, LUTf & Lcurve, LUTf & abcurve*/ );
-
-		void dirpyrLab_equalizer (LabImage * src, LabImage * dst, const double * mult );//Emil's directional pyramid equalizer
-		void dirpyr_eq           (LabImage* data_coarse, LabImage* data_fine, LUTf & rangefn, int level, int pitch, int scale, const double * mult );
-		void idirpyr_eq          (LabImage* data_coarse, LabImage* data_fine, int *** buffer, int level, int pitch, int scale, const double * mult );
-
+		
+		// FT denoise
+		//void RGB_InputTransf(Imagefloat * src, LabImage * dst, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe);
+		//void RGB_OutputTransf(LabImage * src, Imagefloat * dst, const procparams::DirPyrDenoiseParams & dnparams);
+		//void output_tile_row (float *Lbloxrow, float ** Lhipassdn, float ** tilemask, int height, int width, int top, int blkrad );
+		void RGB_denoise(Imagefloat * src, Imagefloat * dst, const procparams::DirPyrDenoiseParams & dnparams, const procparams::DefringeParams & defringe);
+		void RGBtile_denoise (float * fLblox, int vblproc, int hblproc, int numblox_H, int numblox_W, float noisevar_L );	//for DCT
+		void RGBoutput_tile_row (float *Lbloxrow, float ** Ldetail, float ** tilemask_out, int height, int width, int top );
+		//void WaveletDenoise(cplx_wavelet_decomposition &DualTreeCoeffs, float noisevar );
+		//void WaveletDenoise(wavelet_decomposition &WaveletCoeffs, float noisevar );
+		void WaveletDenoiseAll(wavelet_decomposition &WaveletCoeffs_L, wavelet_decomposition &WaveletCoeffs_a, 
+							   wavelet_decomposition &WaveletCoeffs_b, float noisevar_L, float noisevar_ab );
+		void WaveletDenoiseAll_BiShrink(wavelet_decomposition &WaveletCoeffs_L, wavelet_decomposition &WaveletCoeffs_a, 
+										wavelet_decomposition &WaveletCoeffs_b, float noisevar_L, float noisevar_ab );
+		//void BiShrink(float * ReCoeffs, float * ImCoeffs, float * ReParents, float * ImParents, 
+		//			  int W, int H, int level, int padding, float noisevar);
+		//void Shrink(float ** WavCoeffs, int W, int H, int level, float noisevar);
+		void ShrinkAll(float ** WavCoeffs_L, float ** WavCoeffs_a, float ** WavCoeffs_b, int level, 
+					   int W_L, int H_L, int W_ab, int H_ab, int skip_L, int skip_ab, float noisevar_L, float noisevar_ab);
+		float MadMax(float * HH_Coeffs, int &max, int datalen);
+		
+		// pyramid equalizer
 		void dirpyr_equalizer    (float ** src, float ** dst, int srcwidth, int srcheight, const double * mult );//Emil's directional pyramid equalizer
 		void dirpyr_channel      (float ** data_fine, float ** data_coarse, int width, int height, LUTf & rangefn, int level, int scale, const double * mult  );
 		void idirpyr_eq_channel  (float ** data_coarse, float ** data_fine, float ** buffer, int width, int height, int level, const double * mult );
diff --git a/rtengine/labimage.h b/rtengine/labimage.h
index cc600a22a..7389e0fbc 100644
--- a/rtengine/labimage.h
+++ b/rtengine/labimage.h
@@ -26,10 +26,10 @@ namespace rtengine {
 class LabImage {
 private:
 	bool fromImage;
-	float * data;
 
 public:
 	int W, H;
+	float * data;
 	float** L;
 	float** a;
 	float** b;
diff --git a/rtengine/procevents.h b/rtengine/procevents.h
index de45f62bb..6e9fb1311 100644
--- a/rtengine/procevents.h
+++ b/rtengine/procevents.h
@@ -194,8 +194,9 @@ enum ProcEvent {
 	EvVibranceSkinTonesCurve=169,
 	EvLLCCurve=170,
 	EvLLCredsk=171,
+	EvDPDNLdetail=172,
 	
-	NUMOFEVENTS=172
+	NUMOFEVENTS=173
 };
 }    
 #endif    
diff --git a/rtengine/procparams.cc b/rtengine/procparams.cc
index 89caff871..49cc8255a 100644
--- a/rtengine/procparams.cc
+++ b/rtengine/procparams.cc
@@ -232,9 +232,11 @@ void ProcParams::setDefaults () {
     defringe.threshold          = 25;
 
     dirpyrDenoise.enabled       = false;
-    dirpyrDenoise.luma          = 10;
-    dirpyrDenoise.chroma        = 10;
-    dirpyrDenoise.gamma         = 2.0;
+    dirpyrDenoise.luma          = 30;
+    dirpyrDenoise.Ldetail       = 50;
+    dirpyrDenoise.chroma        = 30;
+    dirpyrDenoise.gamma         = 1.7;
+	dirpyrDenoise.expcomp       = 0.0;
 
     edgePreservingDecompositionUI.enabled = false;
     edgePreservingDecompositionUI.Strength = 0.25;
@@ -537,8 +539,9 @@ int ProcParams::save (Glib::ustring fname, Glib::ustring fname2, ParamsEdited* p
 
     // save dirpyrDenoise
     if (!pedited || pedited->dirpyrDenoise.enabled) keyFile.set_boolean ("Directional Pyramid Denoising", "Enabled", dirpyrDenoise.enabled);
-    if (!pedited || pedited->dirpyrDenoise.luma)    keyFile.set_integer ("Directional Pyramid Denoising", "Luma",    dirpyrDenoise.luma);
-    if (!pedited || pedited->dirpyrDenoise.chroma)  keyFile.set_integer ("Directional Pyramid Denoising", "Chroma",  dirpyrDenoise.chroma);
+    if (!pedited || pedited->dirpyrDenoise.luma)    keyFile.set_double ("Directional Pyramid Denoising", "Luma",    dirpyrDenoise.luma);
+    if (!pedited || pedited->dirpyrDenoise.Ldetail) keyFile.set_double ("Directional Pyramid Denoising", "Ldetail", dirpyrDenoise.Ldetail);
+    if (!pedited || pedited->dirpyrDenoise.chroma)  keyFile.set_double ("Directional Pyramid Denoising", "Chroma",  dirpyrDenoise.chroma);
     if (!pedited || pedited->dirpyrDenoise.gamma)   keyFile.set_double  ("Directional Pyramid Denoising", "Gamma",   dirpyrDenoise.gamma);
 
     //Save edgePreservingDecompositionUI.
@@ -961,11 +964,12 @@ if (keyFile.has_group ("Impulse Denoising")) {
 }
 		
 	// load dirpyrDenoise
-if (keyFile.has_group ("Directional Pyramid Denoising")) {
+if (keyFile.has_group ("Directional Pyramid Denoising")) {//TODO: No longer an accurate description for FT denoise
 	if (keyFile.has_key ("Directional Pyramid Denoising", "Enabled"))    { dirpyrDenoise.enabled = keyFile.get_boolean ("Directional Pyramid Denoising", "Enabled"); if (pedited) pedited->dirpyrDenoise.enabled = true; }
-	if (keyFile.has_key ("Directional Pyramid Denoising", "Luma"))       { dirpyrDenoise.luma    = keyFile.get_integer ("Directional Pyramid Denoising", "Luma"); if (pedited) pedited->dirpyrDenoise.luma = true; }
-	if (keyFile.has_key ("Directional Pyramid Denoising", "Chroma"))     { dirpyrDenoise.chroma  = keyFile.get_integer ("Directional Pyramid Denoising", "Chroma"); if (pedited) pedited->dirpyrDenoise.chroma = true; }
-	if (keyFile.has_key ("Directional Pyramid Denoising", "Gamma"))      { dirpyrDenoise.gamma  = keyFile.get_double ("Directional Pyramid Denoising", "Gamma"); if (pedited) pedited->dirpyrDenoise.gamma = true; }
+	if (keyFile.has_key ("Directional Pyramid Denoising", "Luma"))       { dirpyrDenoise.luma    = keyFile.get_double ("Directional Pyramid Denoising", "Luma"); if (pedited) pedited->dirpyrDenoise.luma = true; }
+	if (keyFile.has_key ("Directional Pyramid Denoising", "Ldetail"))    { dirpyrDenoise.Ldetail = keyFile.get_double ("Directional Pyramid Denoising", "Ldetail"); if (pedited) pedited->dirpyrDenoise.Ldetail = true; }
+	if (keyFile.has_key ("Directional Pyramid Denoising", "Chroma"))     { dirpyrDenoise.chroma  = keyFile.get_double ("Directional Pyramid Denoising", "Chroma"); if (pedited) pedited->dirpyrDenoise.chroma = true; }
+	if (keyFile.has_key ("Directional Pyramid Denoising", "Gamma"))      { dirpyrDenoise.gamma   = keyFile.get_double ("Directional Pyramid Denoising", "Gamma"); if (pedited) pedited->dirpyrDenoise.gamma = true; }
 }
 
 //Load EPD.
@@ -1315,6 +1319,7 @@ bool ProcParams::operator== (const ProcParams& other) {
 		&& impulseDenoise.thresh == other.impulseDenoise.thresh
 		&& dirpyrDenoise.enabled == other.dirpyrDenoise.enabled
 		&& dirpyrDenoise.luma == other.dirpyrDenoise.luma
+		&& dirpyrDenoise.Ldetail == other.dirpyrDenoise.Ldetail
 		&& dirpyrDenoise.chroma == other.dirpyrDenoise.chroma
 		&& dirpyrDenoise.gamma == other.dirpyrDenoise.gamma
 		&& edgePreservingDecompositionUI.enabled == other.edgePreservingDecompositionUI.enabled
diff --git a/rtengine/procparams.h b/rtengine/procparams.h
index 728e92af0..ba7462d1e 100644
--- a/rtengine/procparams.h
+++ b/rtengine/procparams.h
@@ -405,9 +405,11 @@ class WBParams {
 		
     public:
         bool    enabled;
-        int		luma;
-        int     chroma;
-		float	gamma;
+        double	luma;
+        double	Ldetail;
+        double  chroma;
+		double	gamma;
+		double	expcomp;
 	};
 
 //EPD related parameters.
diff --git a/rtengine/rawimagesource.cc b/rtengine/rawimagesource.cc
index 4a44a5750..1ad945f84 100644
--- a/rtengine/rawimagesource.cc
+++ b/rtengine/rawimagesource.cc
@@ -75,10 +75,10 @@ RawImageSource::RawImageSource ()
 ,border(4)
 ,ri(NULL)
 ,cache(NULL)
-,rawData(NULL)
-,green(NULL)
-,red(NULL)
-,blue(NULL)
+,rawData(0,0)
+,green(0,0)
+,red(0,0)
+,blue(0,0)
 {
     hrmap[0] = NULL;
     hrmap[1] = NULL;
@@ -99,14 +99,9 @@ RawImageSource::~RawImageSource () {
         delete ri;
     }
 
-    if (green)
-        freeArray<float>(green, H);
-    if (red)
-        freeArray<float>(red, H);
-    if (blue)
-        freeArray<float>(blue, H);
-    if(rawData)
-    	freeArray<float>(rawData, H);
+    flushRGB();
+    flushRawData();
+
     if( cache )
         delete [] cache;
     if (hrmap[0]!=NULL) {
@@ -218,9 +213,9 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
     double r, g, b;
     float rm, gm, bm;
     ctemp.getMultipliers (r, g, b);
-    rm = cam_rgb[0][0]*r + cam_rgb[0][1]*g + cam_rgb[0][2]*b;
-    gm = cam_rgb[1][0]*r + cam_rgb[1][1]*g + cam_rgb[1][2]*b;
-    bm = cam_rgb[2][0]*r + cam_rgb[2][1]*g + cam_rgb[2][2]*b;
+    rm = imatrices.cam_rgb[0][0]*r + imatrices.cam_rgb[0][1]*g + imatrices.cam_rgb[0][2]*b;
+    gm = imatrices.cam_rgb[1][0]*r + imatrices.cam_rgb[1][1]*g + imatrices.cam_rgb[1][2]*b;
+    bm = imatrices.cam_rgb[2][0]*r + imatrices.cam_rgb[2][1]*g + imatrices.cam_rgb[2][2]*b;
     rm = camwb_red / rm;
     gm = camwb_green / gm;
     bm = camwb_blue / bm;
@@ -427,7 +422,12 @@ void RawImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
     // Color correction (only when running on full resolution)
     if (ri->isBayer() && pp.skip==1)
         processFalseColorCorrection (image, raw.ccSteps);
-    colorSpaceConversion (image, cmp, raw, embProfile, camProfile, xyz_cam, (static_cast<const ImageData*>(getMetaData()))->getCamera());
+    // *** colorSpaceConversion was here ***
+    //colorSpaceConversion (image, cmp, raw, embProfile, camProfile, xyz_cam, (static_cast<const ImageData*>(getMetaData()))->getCamera());
+}
+
+void RawImageSource::convertColorSpace(Imagefloat* image, ColorManagementParams cmp, RAWParams raw) {
+    colorSpaceConversion (image, cmp, raw, embProfile, camProfile, imatrices.xyz_cam, (static_cast<const ImageData*>(getMetaData()))->getCamera());
 }
 	
 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -892,10 +892,10 @@ int RawImageSource::load (Glib::ustring fname, bool batch) {
     fuji = ri->get_FujiWidth()!=0;
     for (int i=0; i<3; i++)
         for (int j=0; j<3; j++)
-            rgb_cam[i][j] = ri->get_rgb_cam(i,j);
+        	imatrices.rgb_cam[i][j] = ri->get_rgb_cam(i,j);
     // compute inverse of the color transformation matrix
 	// first arg is matrix, second arg is inverse
-    inverse33 (rgb_cam, cam_rgb);
+    inverse33 (imatrices.rgb_cam, imatrices.cam_rgb);
 
     d1x  = ! ri->get_model().compare("D1X");
     if (d1x)
@@ -904,13 +904,13 @@ int RawImageSource::load (Glib::ustring fname, bool batch) {
         embProfile = cmsOpenProfileFromMem (ri->get_profile(), ri->get_profileLen());
 
     // create profile
-    memset (xyz_cam, 0, sizeof(xyz_cam));
+    memset (imatrices.xyz_cam, 0, sizeof(imatrices.xyz_cam));
     for (int i=0; i<3; i++)
         for (int j=0; j<3; j++)
             for (int k=0; k<3; k++)
-                xyz_cam[i][j] += xyz_sRGB[i][k] * rgb_cam[k][j]; 
-    camProfile = iccStore->createFromMatrix (xyz_cam, false, "Camera");
-    inverse33 (xyz_cam, cam_xyz);
+            	imatrices.xyz_cam[i][j] += xyz_sRGB[i][k] * imatrices.rgb_cam[k][j];
+    camProfile = iccStore->createFromMatrix (imatrices.xyz_cam, false, "Camera");
+    inverse33 (imatrices.xyz_cam, imatrices.cam_xyz);
 
 	float pre_mul[4];
 	ri->get_colorsCoeff( pre_mul, scale_mul, c_black);//modify  for black level
@@ -920,9 +920,9 @@ int RawImageSource::load (Glib::ustring fname, bool batch) {
 	camwb_blue = ri->get_pre_mul(2) / pre_mul[2];
 	initialGain = 1.0 / min(pre_mul[0], pre_mul[1], pre_mul[2]);
 
-    double cam_r = rgb_cam[0][0]*camwb_red + rgb_cam[0][1]*camwb_green + rgb_cam[0][2]*camwb_blue;
-    double cam_g = rgb_cam[1][0]*camwb_red + rgb_cam[1][1]*camwb_green + rgb_cam[1][2]*camwb_blue;
-    double cam_b = rgb_cam[2][0]*camwb_red + rgb_cam[2][1]*camwb_green + rgb_cam[2][2]*camwb_blue;
+    double cam_r = imatrices.rgb_cam[0][0]*camwb_red + imatrices.rgb_cam[0][1]*camwb_green + imatrices.rgb_cam[0][2]*camwb_blue;
+    double cam_g = imatrices.rgb_cam[1][0]*camwb_red + imatrices.rgb_cam[1][1]*camwb_green + imatrices.rgb_cam[1][2]*camwb_blue;
+    double cam_b = imatrices.rgb_cam[2][0]*camwb_red + imatrices.rgb_cam[2][1]*camwb_green + imatrices.rgb_cam[2][2]*camwb_blue;
 
     wb = ColorTemp (cam_r, cam_g, cam_b);
 
@@ -935,9 +935,9 @@ int RawImageSource::load (Glib::ustring fname, bool batch) {
     rml.ciffLength = ri->get_ciffLen();
     idata = new ImageData (fname, &rml);
 
-    green = allocArray<float>(W,H);
-    red   = allocArray<float>(W,H);
-    blue  = allocArray<float>(W,H);
+    green(W,H);
+    red(W,H);
+    blue(W,H);
     //hpmap = allocArray<char>(W, H);
 
     if (plistener) {
@@ -1156,23 +1156,19 @@ void RawImageSource::flushRawData() {
         cache = 0;
     }
     if (rawData) {
-        freeArray<float>(rawData, H);
-        rawData = 0;
+        rawData(0,0);
     }
 }
 
 void RawImageSource::flushRGB() {
     if (green) {
-        freeArray<float>(green, H);
-        green = 0;
+        green(0,0);
     }
     if (red) {
-        freeArray<float>(red, H);
-        red = 0;
+        red(0,0);
     }
     if (blue) {
-        freeArray<float>(blue, H);
-        blue = 0;
+        blue(0,0);
     }
 }
 
@@ -1201,7 +1197,7 @@ void RawImageSource::copyOriginalPixels(const RAWParams &raw, RawImage *src, Raw
 
 	if (ri->isBayer()) {
 		if (!rawData)
-			rawData = allocArray<float>(W,H);
+			rawData(W,H);
 		if (riDark && W == riDark->get_width() && H == riDark->get_height()) {
 			for (int row = 0; row < H; row++) {
 				for (int col = 0; col < W; col++) {
@@ -1286,7 +1282,7 @@ void RawImageSource::copyOriginalPixels(const RAWParams &raw, RawImage *src, Raw
 	} else {
         // No bayer pattern
         // TODO: Is there a flat field correction possible?
-		if (!rawData) rawData = allocArray<float>(3*W,H);
+		if (!rawData) rawData(3*W,H);
 
 		if (riDark && W == riDark->get_width() && H == riDark->get_height()) {
 			for (int row = 0; row < H; row++) {
@@ -1496,7 +1492,7 @@ void RawImageSource::processFalseColorCorrectionThread  (Imagefloat* im, int row
   float middle_Q[6];
   float* tmp;
 
-  int ppx, px=(row_from-1)%3, cx=row_from%3, nx=0;
+  int ppx=0, px=(row_from-1)%3, cx=row_from%3, nx=0;
   
     convert_row_to_YIQ (im->r[row_from-1], im->g[row_from-1], im->b[row_from-1], rbconv_Y[px], rbconv_I[px], rbconv_Q[px], W);
     convert_row_to_YIQ (im->r[row_from], im->g[row_from], im->b[row_from], rbconv_Y[cx], rbconv_I[cx], rbconv_Q[cx], W);
@@ -1658,7 +1654,7 @@ void RawImageSource::colorSpaceConversion (Imagefloat* im, ColorManagementParams
         for (int i=0; i<3; i++)
             for (int j=0; j<3; j++) 
                 for (int k=0; k<3; k++) 
-                    mat[i][j] += work[i][k] * camMatrix[k][j]; // rgb_xyz * xyz_cam
+                    mat[i][j] += work[i][k] * camMatrix[k][j]; // rgb_xyz * imatrices.xyz_cam
 
 
     if (in==NULL) {
@@ -1687,9 +1683,9 @@ void RawImageSource::colorSpaceConversion (Imagefloat* im, ColorManagementParams
         #pragma omp parallel for
 		for ( int h = 0; h < im->height; ++h )
 			for ( int w = 0; w < im->width; ++w ) {
-				im->r[h][w] /= 65535.0;
-				im->g[h][w] /= 65535.0;
-				im->b[h][w] /= 65535.0;
+				im->r[h][w] /= 65535.0f ;
+				im->g[h][w] /= 65535.0f ;
+				im->b[h][w] /= 65535.0f ;
 			}
 
 
@@ -1884,7 +1880,7 @@ TMatrix work = iccStore->workingSpaceInverseMatrix (cmp.working);
 		for (int i=0; i<3; i++)
 			for (int j=0; j<3; j++) 
 				for (int k=0; k<3; k++) 
-					mat[i][j] += work[i][k] * camMatrix[k][j]; // rgb_xyz * xyz_cam
+					mat[i][j] += work[i][k] * camMatrix[k][j]; // rgb_xyz * imatrices.xyz_cam
 		
 #pragma omp parallel for
 		for (int i=0; i<im->height; i++)
@@ -2192,7 +2188,7 @@ void RawImageSource::hlRecovery (std::string method, float* red, float* green, f
     if (method=="Luminance")
         HLRecovery_Luminance (red, green, blue, red, green, blue, width, 65535.0);
     else if (method=="CIELab blending")
-        HLRecovery_CIELab (red, green, blue, red, green, blue, width, 65535.0, xyz_cam, cam_xyz);
+        HLRecovery_CIELab (red, green, blue, red, green, blue, width, 65535.0, imatrices.xyz_cam, imatrices.cam_xyz);
     /*else if (method=="Color")
         HLRecovery_ColorPropagation (red, green, blue, i, sx1, width, skip);*/
 	else if (method=="Blend")	// derived from Dcraw
@@ -2236,6 +2232,9 @@ void RawImageSource::getRAWHistogram (LUTu & histRedRaw, LUTu & histGreenRaw, LU
 
     histRedRaw.clear(); histGreenRaw.clear(); histBlueRaw.clear();
 	float mult = 65535.0 / ri->get_white();
+
+	// WARNING: This parallelization is not thread-safe because it R/W in histRedRaw, histGreenRaw, histBlueRaw
+	// which are defined before the parallel section and must survive after it
     #pragma omp parallel for
     for (int i=border; i<H-border; i++) {
         int start, end, idx;
@@ -2385,9 +2384,9 @@ void RawImageSource::getRowStartEnd (int x, int &start, int &end) {
 		double greens = avg_g/gn * camwb_green;
 		double blues  = avg_b/bn * camwb_blue;
 		
-		double rm = rgb_cam[0][0]*reds + rgb_cam[0][1]*greens + rgb_cam[0][2]*blues;
-		double gm = rgb_cam[1][0]*reds + rgb_cam[1][1]*greens + rgb_cam[1][2]*blues;
-		double bm = rgb_cam[2][0]*reds + rgb_cam[2][1]*greens + rgb_cam[2][2]*blues;
+		double rm = imatrices.rgb_cam[0][0]*reds + imatrices.rgb_cam[0][1]*greens + imatrices.rgb_cam[0][2]*blues;
+		double gm = imatrices.rgb_cam[1][0]*reds + imatrices.rgb_cam[1][1]*greens + imatrices.rgb_cam[1][2]*blues;
+		double bm = imatrices.rgb_cam[2][0]*reds + imatrices.rgb_cam[2][1]*greens + imatrices.rgb_cam[2][2]*blues;
 		
 		return ColorTemp (rm, gm, bm);
 	}
@@ -2466,9 +2465,9 @@ void RawImageSource::getRowStartEnd (int x, int &start, int &end) {
 			greens = greens/rn * camwb_green;
 			blues = blues/rn * camwb_blue;
 			
-			double rm = rgb_cam[0][0]*reds + rgb_cam[0][1]*greens + rgb_cam[0][2]*blues;
-			double gm = rgb_cam[1][0]*reds + rgb_cam[1][1]*greens + rgb_cam[1][2]*blues;
-			double bm = rgb_cam[2][0]*reds + rgb_cam[2][1]*greens + rgb_cam[2][2]*blues;
+			double rm = imatrices.rgb_cam[0][0]*reds + imatrices.rgb_cam[0][1]*greens + imatrices.rgb_cam[0][2]*blues;
+			double gm = imatrices.rgb_cam[1][0]*reds + imatrices.rgb_cam[1][1]*greens + imatrices.rgb_cam[1][2]*blues;
+			double bm = imatrices.rgb_cam[2][0]*reds + imatrices.rgb_cam[2][1]*greens + imatrices.rgb_cam[2][2]*blues;
 			
 			return ColorTemp (rm, gm, bm);
 		}
diff --git a/rtengine/rawimagesource.h b/rtengine/rawimagesource.h
index 68b3ac8aa..2c29ad396 100644
--- a/rtengine/rawimagesource.h
+++ b/rtengine/rawimagesource.h
@@ -24,6 +24,8 @@
 #include "dcp.h"
 #include "array2D.h"
 #include "curves.h"
+#include <fftw3.h>
+#include "color.h"
 #include "../rtgui/cacheimagedata.h"
 
 #define HR_SCALE 2
@@ -101,14 +103,14 @@ class RawImageSource : public ImageSource {
         double* cache;
         int threshold;
 
-        float** rawData;  // holds preprocessed pixel values, rowData[i][j] corresponds to the ith row and jth column
+        array2D<float> rawData;  // holds preprocessed pixel values, rowData[i][j] corresponds to the ith row and jth column
 
         // the interpolated green plane:
-        float** green; 
+        array2D<float> green; 
         // the interpolated red plane:
-        float** red;
+        array2D<float> red;
         // the interpolated blue plane:
-        float** blue;
+        array2D<float> blue;
 
 
         void hphd_vertical       (float** hpmap, int col_from, int col_to);
@@ -150,21 +152,23 @@ class RawImageSource : public ImageSource {
         ColorTemp   getWB       () { return wb; }
         ColorTemp   getAutoWB   ();
         ColorTemp   getSpotWB   (std::vector<Coord2D> red, std::vector<Coord2D> green, std::vector<Coord2D>& blue, int tran);
-        bool        isWBProviderReady () { return rawData != (float**)NULL; };
+        bool        isWBProviderReady () { return rawData; }
 
         double      getDefGain  () { return defGain; }
 
-        double      getGamma    () { return CurveFactory::sRGBGamma; }
+        double      getGamma    () { return Color::sRGBGamma; }
 
         void        getFullSize (int& w, int& h, int tr = TR_NONE);
         void        getSize     (int tran, PreviewProps pp, int& w, int& h);
         int         getRotateDegree() const { return ri->get_rotateDegree(); }
 
         ImageData*  getImageData () { return idata; }
+        ImageMatrices* getImageMatrices () { return &imatrices; }
         void        setProgressListener (ProgressListener* pl) { plistener = pl; }
         void        getAutoExpHistogram (LUTu & histogram, int& histcompr);
         void        getRAWHistogram (LUTu & histRedRaw, LUTu & histGreenRaw, LUTu & histBlueRaw);
 
+        void convertColorSpace(Imagefloat* image, ColorManagementParams cmp, RAWParams raw);
         static void colorSpaceConversion16 (Image16* im, ColorManagementParams cmp, cmsHPROFILE embedded, cmsHPROFILE camprofile, double cam[3][3], std::string camName);
         static void colorSpaceConversion (Imagefloat* im, ColorManagementParams cmp, RAWParams raw, cmsHPROFILE embedded, cmsHPROFILE camprofile, double cam[3][3], std::string camName);
         static void inverse33 (const double (*coeff)[3], double (*icoeff)[3]);
@@ -201,6 +205,8 @@ class RawImageSource : public ImageSource {
         int  findHotDeadPixel( PixelsMap &bpMap, float thresh);
 
         void cfa_linedn (float linenoiselevel);//Emil's line denoise
+        void cfa_tile_denoise (fftwf_complex * fcfablox, int vblk, int hblk, int numblox_H, int numblox_W, float noisevar, float * rolloff );
+        void cfa_output_tile_row (float *cfabloxrow, float ** cfahipassdn, float ** tilemask_out, int height, int width, int top, int blkrad);
 
         void green_equilibrate (float greenthresh);//Emil's green equilibration
 
diff --git a/rtengine/refreshmap.cc b/rtengine/refreshmap.cc
index c7dfbf519..656c3f1a2 100644
--- a/rtengine/refreshmap.cc
+++ b/rtengine/refreshmap.cc
@@ -108,10 +108,10 @@ DARKFRAME,        // EvLCPFile
 0,                // EvEqlEnabled:obsolete
 IMPULSEDENOISE,   // EvIDNEnabled,
 IMPULSEDENOISE,   // EvIDNThresh,
-DIRPYRDENOISE,    // EvDPDNEnabled,
-DIRPYRDENOISE,    // EvDPDNLuma,
-DIRPYRDENOISE,    // EvDPDNChroma,
-DIRPYRDENOISE,    // EvDPDNGamma,
+ALLNORAW,         // EvDPDNEnabled,
+ALLNORAW,         // EvDPDNLuma,
+ALLNORAW,         // EvDPDNChroma,
+ALLNORAW,         // EvDPDNGamma,
 DIRPYREQUALIZER,  // EvDirPyrEqualizer,
 DIRPYREQUALIZER,  // EvDirPyrEqlEnabled,
 LUMINANCECURVE,   // EvLSaturation,
@@ -150,8 +150,8 @@ FLATFIELD,        // EvFlatFieldAutoSelect,
 FLATFIELD,        // EvFlatFieldBlurRadius,
 FLATFIELD,        // EvFlatFieldBlurType,
 TRANSFORM,        // EvAutoDIST,
-DIRPYRDENOISE,    // EvDPDNLumCurve,
-DIRPYRDENOISE,    // EvDPDNChromCurve,
+ALLNORAW,         // EvDPDNLumCurve,
+ALLNORAW,         // EvDPDNChromCurve,
 GAMMA,            // EvGAMMA
 GAMMA,            // EvGAMPOS
 GAMMA,            // EvGAMFREE
@@ -191,6 +191,7 @@ LUMINANCECURVE,   // EvLCCurve
 LUMINANCECURVE,   // EvLCHCurve
 RGBCURVE,         // EvVibranceSkinTonesCurve
 LUMINANCECURVE,   // EvLLCCurve
-LUMINANCECURVE    // EvLLCredsk
+LUMINANCECURVE,   // EvLLCredsk
+ALLNORAW          // EvDPDNLdetail
 };
 
diff --git a/rtengine/refreshmap.h b/rtengine/refreshmap.h
index 8594d2ac5..00bdabada 100644
--- a/rtengine/refreshmap.h
+++ b/rtengine/refreshmap.h
@@ -30,21 +30,22 @@
 
 // Elementary functions that can be done to
 // the preview image when an event occurs
-#define M_PREPROC    (1<<9)
-#define M_RAW        (1<<8)
-#define M_INIT       (1<<7)
-#define M_TRANSFORM  (1<<6)
-#define M_BLURMAP    (1<<5)
-#define M_AUTOEXP    (1<<4)
-#define M_RGBCURVE   (1<<3)
-#define M_LUMACURVE  (1<<2)
-#define M_LUMINANCE  (1<<1)
-#define M_COLOR      (1<<0)
+#define M_PREPROC     (1<<10)
+#define M_RAW         (1<<9)
+#define M_INIT        (1<<8)
+#define M_LINDENOISE  (1<<7)
+#define M_TRANSFORM   (1<<6)
+#define M_BLURMAP     (1<<5)
+#define M_AUTOEXP     (1<<4)
+#define M_RGBCURVE    (1<<3)
+#define M_LUMACURVE   (1<<2)
+#define M_LUMINANCE   (1<<1)
+#define M_COLOR       (1<<0)
 
 // Bitfield of functions to do to the preview image when an event occurs
 // Use those or create new ones for your new events
-#define FIRST           (M_PREPROC|M_RAW|M_INIT|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)  // without HIGHQUAL
-#define ALL             (M_PREPROC|M_RAW|M_INIT|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)  // without HIGHQUAL
+#define FIRST           (M_PREPROC|M_RAW|M_INIT|M_LINDENOISE|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)  // without HIGHQUAL
+#define ALL             (M_PREPROC|M_RAW|M_INIT|M_LINDENOISE|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)  // without HIGHQUAL
 #define TRANSFORM       (M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
 #define RETINEX         (M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
 #define AUTOEXP         (M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
@@ -53,10 +54,10 @@
 #define SHARPENING       M_LUMINANCE
 #define IMPULSEDENOISE   M_LUMINANCE
 #define DEFRINGE         M_LUMINANCE
-#define WHITEBALANCE    (M_INIT|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
-#define DEMOSAIC        (M_RAW|M_INIT|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
-#define DARKFRAME       (M_PREPROC|M_RAW|M_INIT|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
-#define FLATFIELD       (M_PREPROC|M_RAW|M_INIT|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
+#define WHITEBALANCE    (M_INIT|M_LINDENOISE|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
+#define DEMOSAIC        (M_RAW|M_INIT|M_LINDENOISE|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
+#define DARKFRAME       (M_PREPROC|M_RAW|M_INIT|M_LINDENOISE|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
+#define FLATFIELD       (M_PREPROC|M_RAW|M_INIT|M_LINDENOISE|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
 #define DIRPYRDENOISE   (M_COLOR|M_LUMINANCE)
 #define CROP             M_MINUPDATE
 #define RESIZE           M_VOID
@@ -66,7 +67,7 @@
 #define OUTPUTPROFIL    (M_COLOR|M_LUMINANCE)
 #define GAMMA           (M_COLOR|M_LUMINANCE)
 #define NONE             0
-#define ALLNORAW        (M_INIT|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
+#define ALLNORAW        (M_INIT|M_LINDENOISE|M_TRANSFORM|M_BLURMAP|M_AUTOEXP|M_RGBCURVE|M_LUMACURVE|M_LUMINANCE|M_COLOR)
 
 extern int refreshmap[];
 #endif    
diff --git a/rtengine/rtthumbnail.cc b/rtengine/rtthumbnail.cc
index 67e93381c..538f50f26 100644
--- a/rtengine/rtthumbnail.cc
+++ b/rtengine/rtthumbnail.cc
@@ -678,6 +678,9 @@ IImage8* Thumbnail::processImage (const procparams::ProcParams& params, int rhei
         }
     }
 */
+
+    // if luma denoise has to be done for thumbnails, it should be right here
+
     // perform color space transformation
     if (isRaw)
         RawImageSource::colorSpaceConversion16 (resImg, params.icm, embProfile, camProfile, cam2xyz, camName );
@@ -981,7 +984,7 @@ unsigned char* Thumbnail::getGrayscaleHistEQ (int trim_width) {
         unsigned long cutoff = thumbImg->height * thumbImg->height * 4 * BurnOffPct;
 
         int max_;
-	unsigned long sum=0;
+        unsigned long sum=0;
         for (max_=65535; max_>16384 && sum<cutoff; max_--) sum+=hist16[max_];
 
         delete[] hist16;
@@ -992,8 +995,8 @@ unsigned char* Thumbnail::getGrayscaleHistEQ (int trim_width) {
         for (int i=0; i<thumbImg->height; i++)
             for (int j=(thumbImg->width-trim_width)/2; j<trim_width+(thumbImg->width-trim_width)/2; j++) {
                 int r= gammatab[min(thumbImg->r[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
-		int g= gammatab[min(thumbImg->g[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
-		int b= gammatab[min(thumbImg->b[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
+                int g= gammatab[min(thumbImg->g[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
+                int b= gammatab[min(thumbImg->b[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
                 tmpdata[ix++] = (r*19595+g*38469+b*7472) >> 16;
             }
     }
@@ -1087,8 +1090,8 @@ bool Thumbnail::writeImage (const Glib::ustring& fname, int format) {
             for (int i=0; i<thumbImg->height; i++)
                 for (int j=0; j<thumbImg->width; j++) {
                     tmpdata[ix++] = gammatab[min(thumbImg->r[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
-		    tmpdata[ix++] = gammatab[min(thumbImg->g[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
-		    tmpdata[ix++] = gammatab[min(thumbImg->b[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
+                    tmpdata[ix++] = gammatab[min(thumbImg->g[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
+                    tmpdata[ix++] = gammatab[min(thumbImg->b[i][j],static_cast<unsigned short>(max_)) * scaleForSave >> 13];
                 }
         }
         else {
diff --git a/rtengine/simpleprocess.cc b/rtengine/simpleprocess.cc
index 49eddd282..5c437f4d7 100644
--- a/rtengine/simpleprocess.cc
+++ b/rtengine/simpleprocess.cc
@@ -107,6 +107,12 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p
     imgsrc->getImage (currWB, tr, baseImg, pp, params.hlrecovery, params.icm, params.raw);
     if (pl) pl->setProgress (0.45);
 
+    // perform luma/chroma denoise
+    LabImage* labView = new LabImage (fw,fh);
+    if (params.dirpyrDenoise.enabled) {
+		ipf.RGB_denoise(baseImg, baseImg, params.dirpyrDenoise, params.defringe);
+    }
+    imgsrc->convertColorSpace(baseImg, params.icm, params.raw);
 
     // perform first analysis
     LUTu hist16 (65536);
@@ -146,7 +152,7 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p
     }
 
     // at this stage, we can flush the raw data to free up quite an important amount of memory
-    // commented out because it make the application crash when batch processing...
+    // commented out because it makes the application crash when batch processing...
     // TODO: find a better place to flush rawData and rawRGB
     //imgsrc->flushRawData();
     //imgsrc->flushRGB();
@@ -172,8 +178,6 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p
 	CurveFactory::RGBCurve (params.rgbCurves.gcurve, gCurve, 1);
 	CurveFactory::RGBCurve (params.rgbCurves.bcurve, bCurve, 1);
 
-	LabImage* labView = new LabImage (fw,fh);
-
     ipf.rgbProc (baseImg, labView, curve1, curve2, curve, shmap, params.toneCurve.saturation, rCurve, gCurve, bCurve, nonStandardCurve, expcomp, hlcompr, hlcomprthresh);
 
     // Freeing baseImg because not used anymore
@@ -187,6 +191,9 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p
     if (pl) 
         pl->setProgress (0.5);
 
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	// start tile processing...???
 
     // luminance histogram update
     hist16.clear();
@@ -212,7 +219,6 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p
 
 	ipf.impulsedenoise (labView);
 	ipf.defringe (labView);
-	ipf.dirpyrdenoise (labView);
 	if (params.sharpenEdge.enabled) {
 		 ipf.MLsharpen(labView);
 	}
@@ -234,6 +240,9 @@ IImage16* processImage (ProcessingJob* pjob, int& errorCode, ProgressListener* p
 	// directional pyramid equalizer
     ipf.dirpyrequalizer (labView);//TODO: this is the luminance tonecurve, not the RGB one
 
+	// end tile processing...???
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
     if (pl) pl->setProgress (0.60);
 
diff --git a/rtengine/stdimagesource.cc b/rtengine/stdimagesource.cc
index 138a553da..6ce81ba22 100644
--- a/rtengine/stdimagesource.cc
+++ b/rtengine/stdimagesource.cc
@@ -290,7 +290,8 @@ void StdImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
 	//Image16* tmpim = new Image16 (image->width,image->height);
     getImage_ (ctemp, tran, image, pp, true, hrp);
 
-	colorSpaceConversion (image, cmp, embProfile);
+    // *** colorSpaceConversion was there ***
+/*    	colorSpaceConversion (image, cmp, embProfile);
 	
 	for ( int h = 0; h < image->height; ++h )
 		for ( int w = 0; w < image->width; ++w ) {
@@ -299,17 +300,20 @@ void StdImageSource::getImage (ColorTemp ctemp, int tran, Imagefloat* image, Pre
 			image->b[h][w] *= 65535.0f;
 			//if (h==100 && w==100) printf("stdimsrc after R= %f  G= %f  B= %f  \n",image->r[h][w],image->g[h][w],image->b[h][w]);
 		}
-	
+*/
     // Flip if needed
     if (tran & TR_HFLIP)
-	 hflip (image);
-	 if (tran & TR_VFLIP)
-	 vflip (image);
-	
-	
+        hflip (image);
+    if (tran & TR_VFLIP)
+        vflip (image);
+
     t2.set ();
 }
 	
+void StdImageSource::convertColorSpace(Imagefloat* image, ColorManagementParams cmp, RAWParams raw) {
+    colorSpaceConversion (image, cmp, embProfile);
+}
+
 void StdImageSource::colorSpaceConversion (Imagefloat* im, ColorManagementParams cmp, cmsHPROFILE embedded) {
 	
 	cmsHPROFILE in;
@@ -339,6 +343,15 @@ void StdImageSource::colorSpaceConversion (Imagefloat* im, ColorManagementParams
 		
         cmsDeleteTransform(hTransform);
 	}
+
+	// Code moved from getImage to change the range of the float value from [0.;1.] to [0.;65535.]
+	for ( int h = 0; h < im->height; ++h )
+		for ( int w = 0; w < im->width; ++w ) {
+			im->r[h][w] *= 65535.0f;
+			im->g[h][w] *= 65535.0f;
+			im->b[h][w] *= 65535.0f;
+			//if (h==100 && w==100) printf("stdimsrc after R= %f  G= %f  B= %f  \n",im->r[h][w],im->g[h][w],im->b[h][w]);
+		}
 }
 	
 
@@ -370,6 +383,8 @@ void StdImageSource::colorSpaceConversion16 (Image16* im, ColorManagementParams
         
         cmsDeleteTransform(hTransform);
     }
+
+    // WARNING: A range update may be missing here (see colorSpaceConversion above)
 }
 
 void StdImageSource::getFullSize (int& w, int& h, int tr) {
diff --git a/rtengine/stdimagesource.h b/rtengine/stdimagesource.h
index 22e810d09..220fa8ff1 100644
--- a/rtengine/stdimagesource.h
+++ b/rtengine/stdimagesource.h
@@ -57,7 +57,10 @@ class StdImageSource : public ImageSource {
         void        getSize     (int tran, PreviewProps pp, int& w, int& h);
 
         ImageData*  getImageData () { return idata; }
+        ImageMatrices* getImageMatrices () { return (ImageMatrices*)NULL; }
         void        setProgressListener (ProgressListener* pl) { plistener = pl; }
+
+        void convertColorSpace(Imagefloat* image, ColorManagementParams cmp, RAWParams raw);// RAWParams raw will not be used for non-raw files (see imagesource.h)
         static void colorSpaceConversion (Imagefloat* im, ColorManagementParams cmp, cmsHPROFILE embedded);
 		static void colorSpaceConversion16 (Image16* im, ColorManagementParams cmp, cmsHPROFILE embedded);
 
diff --git a/rtgui/CMakeLists.txt b/rtgui/CMakeLists.txt
index 9fb7d7e34..0f63bceab 100644
--- a/rtgui/CMakeLists.txt
+++ b/rtgui/CMakeLists.txt
@@ -42,9 +42,9 @@ if (WIN32)
     #set_target_properties (rth PROPERTIES LINK_FLAGS "-mwindows")
 else (WIN32)
     include_directories (${EXTRA_INCDIR} ${GLIB2_INCLUDE_DIRS} ${GLIBMM_INCLUDE_DIRS}
-        ${GTK_INCLUDE_DIRS} ${GTKMM_INCLUDE_DIRS} ${GIO_INCLUDE_DIRS} ${GIOMM_INCLUDE_DIRS} ${IPTCDATA_INCLUDE_DIRS} ${LCMS_INCLUDE_DIRS} ${EXPAT_INCLUDE_DIRS} ${GTHREAD_INCLUDE_DIRS} ${GOBJECT_INCLUDE_DIRS} )
+        ${GTK_INCLUDE_DIRS} ${GTKMM_INCLUDE_DIRS} ${GIO_INCLUDE_DIRS} ${GIOMM_INCLUDE_DIRS} ${IPTCDATA_INCLUDE_DIRS} ${LCMS_INCLUDE_DIRS} ${EXPAT_INCLUDE_DIRS} ${FFTW3F_LIBRARY_DIRS} ${GTHREAD_INCLUDE_DIRS} ${GOBJECT_INCLUDE_DIRS} )
     link_directories (${EXTRA_LIBDIR} ${GLIB2_LIBRARY_DIRS} ${GLIBMM_LIBRARY_DIRS}
-        ${GTK_LIBRARY_DIRS} ${GTKMM_LIBRARY_DIRS} ${GIO_LIBRARY_DIRS} ${GIOMM_LIBRARY_DIRS} ${IPTCDATA_LIBRARY_DIRS} ${LCMS_LIBRARY_DIRS} ${EXPAT_LIBRARY_DIRS} ${GTHREAD_LIBRARY_DIRS} ${GOBJECT_LIBRARY_DIRS})
+        ${GTK_LIBRARY_DIRS} ${GTKMM_LIBRARY_DIRS} ${GIO_LIBRARY_DIRS} ${GIOMM_LIBRARY_DIRS} ${IPTCDATA_LIBRARY_DIRS} ${LCMS_LIBRARY_DIRS} ${EXPAT_LIBRARY_DIRS} ${FFTW3F_LIBRARY_DIRS} ${GTHREAD_LIBRARY_DIRS} ${GOBJECT_LIBRARY_DIRS})
 endif (WIN32)
 # create config.h which defines where data are stored
 configure_file ("${CMAKE_CURRENT_SOURCE_DIR}/config.h.in" "${CMAKE_CURRENT_BINARY_DIR}/config.h")
@@ -55,6 +55,6 @@ set_target_properties (rth PROPERTIES COMPILE_FLAGS "${CMAKE_CXX_FLAGS}" OUTPUT_
 #target_link_libraries (rth rtengine ${JPEG_LIBRARIES} ${PNG_LIBRARIES} ${ZLIB_LIBRARIES} ${TIFF_LIBRARIES} ${EXTRA_LIB} ${GOBJECT_LIBRARIES} ${GTHREAD_LIBRARIES}
 #   ${GLIB2_LIBRARIES} ${GLIBMM_LIBRARIES} ${GTK_LIBRARIES} ${GTKMM_LIBRARIES} ${GIO_LIBRARIES} ${GIOMM_LIBRARIES} ${LCMS_LIBRARIES} ${IPTCDATA_LIBRARIES})
 target_link_libraries (rth rtengine ${JPEG_LIBRARIES} ${PNG_LIBRARIES} ${ZLIB_LIBRARIES} ${TIFF_LIBRARIES} ${GOBJECT_LIBRARIES} ${GTHREAD_LIBRARIES}
-   ${GLIB2_LIBRARIES} ${GLIBMM_LIBRARIES} ${GTK_LIBRARIES} ${GTKMM_LIBRARIES} ${GIO_LIBRARIES} ${GIOMM_LIBRARIES} ${LCMS_LIBRARIES} ${EXPAT_LIBRARIES} ${IPTCDATA_LIBRARIES} ${EXTRA_LIB_RTGUI})
+   ${GLIB2_LIBRARIES} ${GLIBMM_LIBRARIES} ${GTK_LIBRARIES} ${GTKMM_LIBRARIES} ${GIO_LIBRARIES} ${GIOMM_LIBRARIES} ${LCMS_LIBRARIES} ${EXPAT_LIBRARIES} ${FFTW3F_LIBRARIES} ${IPTCDATA_LIBRARIES} ${EXTRA_LIB_RTGUI})
 install (TARGETS rth DESTINATION ${BINDIR})
 
diff --git a/rtgui/batchtoolpanelcoord.cc b/rtgui/batchtoolpanelcoord.cc
index 71608ba3e..4b4457054 100644
--- a/rtgui/batchtoolpanelcoord.cc
+++ b/rtgui/batchtoolpanelcoord.cc
@@ -216,7 +216,7 @@ void BatchToolPanelCoordinator::initSession () {
 			if (options.baBehav[ADDSET_VIGN_CENTER])  pparams.vignetting.centerY = 0;
 
 			if (options.baBehav[ADDSET_DIRPYREQ]) for (int i=0; i<5; i++) pparams.dirpyrequalizer.mult[i] = 0;
-			if (options.baBehav[ADDSET_DIRPYRDN_CHLUM])  pparams.dirpyrDenoise.luma = pparams.dirpyrDenoise.chroma = 0;
+			if (options.baBehav[ADDSET_DIRPYRDN_CHLUM])  pparams.dirpyrDenoise.Ldetail = pparams.dirpyrDenoise.luma = pparams.dirpyrDenoise.chroma = 0;
 			if (options.baBehav[ADDSET_DIRPYRDN_GAMMA])  pparams.dirpyrDenoise.gamma = 0;
 
 			if (options.baBehav[ADDSET_PREPROCESS_GREENEQUIL])  pparams.raw.greenthresh = 0;
diff --git a/rtgui/dirpyrdenoise.cc b/rtgui/dirpyrdenoise.cc
index ccf8db930..0da2f8a11 100644
--- a/rtgui/dirpyrdenoise.cc
+++ b/rtgui/dirpyrdenoise.cc
@@ -36,19 +36,25 @@ DirPyrDenoise::DirPyrDenoise () : Gtk::VBox(), FoldableToolPanel(this)  {
 	
 	enaConn = enabled->signal_toggled().connect( sigc::mem_fun(*this, &DirPyrDenoise::enabledChanged) );
 	
-	luma  = Gtk::manage (new Adjuster (M("TP_DIRPYRDENOISE_LUMA"), 0, 100, 1, 5));
-	chroma    = Gtk::manage (new Adjuster (M("TP_DIRPYRDENOISE_CHROMA"), 0, 100, 1, 5));
-	gamma	= Gtk::manage (new Adjuster (M("TP_DIRPYRDENOISE_GAMMA"), 1.0, 3.0, 0.01, 2.0));
+
+	luma  = Gtk::manage (new Adjuster (M("TP_DIRPYRDENOISE_LUMA"), 0, 100, 0.01, 0));
+	Ldetail  = Gtk::manage (new Adjuster (M("TP_DIRPYRDENOISE_LDETAIL"), 0, 100, 0.01, 80));
+	chroma    = Gtk::manage (new Adjuster (M("TP_DIRPYRDENOISE_CHROMA"), 0, 100, 0.01, 15));
+	gamma	= Gtk::manage (new Adjuster (M("TP_DIRPYRDENOISE_GAMMA"), 1.0, 3.0, 0.01, 1.7));
+
 	
 	luma->setAdjusterListener (this);
+	Ldetail->setAdjusterListener (this);
 	chroma->setAdjusterListener (this); 
     gamma->setAdjusterListener (this); 
 	
-	luma->show();
-	chroma->show();
+    luma->show();
+    Ldetail->show();
+    chroma->show();
 	gamma->show();
 	
 	pack_start (*luma);
+	pack_start (*Ldetail);
 	pack_start (*chroma);
 	pack_start (*gamma);
 
@@ -59,9 +65,10 @@ void DirPyrDenoise::read (const ProcParams* pp, const ParamsEdited* pedited) {
     disableListener ();
 	
     if (pedited) {
-        luma->setEditedState    (pedited->dirpyrDenoise.luma ? Edited : UnEdited);
-        chroma->setEditedState      (pedited->dirpyrDenoise.chroma ? Edited : UnEdited);
-        gamma->setEditedState      (pedited->dirpyrDenoise.gamma ? Edited : UnEdited);
+    	luma->setEditedState      (pedited->dirpyrDenoise.luma ? Edited : UnEdited);
+    	Ldetail->setEditedState   (pedited->dirpyrDenoise.Ldetail ? Edited : UnEdited);
+    	chroma->setEditedState    (pedited->dirpyrDenoise.chroma ? Edited : UnEdited);
+        gamma->setEditedState     (pedited->dirpyrDenoise.gamma ? Edited : UnEdited);
         enabled->set_inconsistent (!pedited->dirpyrDenoise.enabled);
     }
 	
@@ -72,8 +79,9 @@ void DirPyrDenoise::read (const ProcParams* pp, const ParamsEdited* pedited) {
     lastEnabled = pp->dirpyrDenoise.enabled;
 	
     luma->setValue    (pp->dirpyrDenoise.luma);
-    chroma->setValue      (pp->dirpyrDenoise.chroma);
-	gamma->setValue      (pp->dirpyrDenoise.gamma);
+    Ldetail->setValue (pp->dirpyrDenoise.Ldetail);
+    chroma->setValue  (pp->dirpyrDenoise.chroma);
+    gamma->setValue   (pp->dirpyrDenoise.gamma);
 
     enableListener ();
 }
@@ -81,46 +89,56 @@ void DirPyrDenoise::read (const ProcParams* pp, const ParamsEdited* pedited) {
 void DirPyrDenoise::write (ProcParams* pp, ParamsEdited* pedited) {
 	
 	pp->dirpyrDenoise.luma      = luma->getValue ();
+	pp->dirpyrDenoise.Ldetail   = Ldetail->getValue ();
 	pp->dirpyrDenoise.chroma	= chroma->getValue ();
 	pp->dirpyrDenoise.gamma		= gamma->getValue ();
 	pp->dirpyrDenoise.enabled   = enabled->get_active();
 	
     if (pedited) {
         pedited->dirpyrDenoise.luma     = luma->getEditedState ();
-        pedited->dirpyrDenoise.chroma	= chroma->getEditedState ();
-		pedited->dirpyrDenoise.gamma	= gamma->getEditedState ();
-		pedited->dirpyrDenoise.enabled  = !enabled->get_inconsistent();
+        pedited->dirpyrDenoise.Ldetail  = Ldetail->getEditedState ();
+        pedited->dirpyrDenoise.chroma   = chroma->getEditedState ();
+        pedited->dirpyrDenoise.gamma    = gamma->getEditedState ();
+        pedited->dirpyrDenoise.enabled  = !enabled->get_inconsistent();
     }
 }
 
 void DirPyrDenoise::setDefaults (const ProcParams* defParams, const ParamsEdited* pedited) {
 	
-    luma->setDefault (defParams->dirpyrDenoise.luma);
-    chroma->setDefault (defParams->dirpyrDenoise.chroma);
-    gamma->setDefault (defParams->dirpyrDenoise.gamma);
+	luma->setDefault    (defParams->dirpyrDenoise.luma);
+	Ldetail->setDefault (defParams->dirpyrDenoise.Ldetail);
+	chroma->setDefault  (defParams->dirpyrDenoise.chroma);
+	gamma->setDefault   (defParams->dirpyrDenoise.gamma);
 
     if (pedited) {
-        luma->setDefaultEditedState		(pedited->dirpyrDenoise.luma ? Edited : UnEdited);
-        chroma->setDefaultEditedState   (pedited->dirpyrDenoise.chroma ? Edited : UnEdited);
-		gamma->setDefaultEditedState   (pedited->dirpyrDenoise.gamma ? Edited : UnEdited);
+    	luma->setDefaultEditedState     (pedited->dirpyrDenoise.luma ? Edited : UnEdited);
+    	Ldetail->setDefaultEditedState	(pedited->dirpyrDenoise.Ldetail ? Edited : UnEdited);
+    	chroma->setDefaultEditedState   (pedited->dirpyrDenoise.chroma ? Edited : UnEdited);
+		gamma->setDefaultEditedState    (pedited->dirpyrDenoise.gamma ? Edited : UnEdited);
    }
     else {
-        luma->setDefaultEditedState (Irrelevant);
-        chroma->setDefaultEditedState   (Irrelevant);
-        gamma->setDefaultEditedState   (Irrelevant);
+    	luma->setDefaultEditedState     (Irrelevant);
+    	Ldetail->setDefaultEditedState  (Irrelevant);
+    	chroma->setDefaultEditedState   (Irrelevant);
+        gamma->setDefaultEditedState    (Irrelevant);
     }
 }
 
 void DirPyrDenoise::adjusterChanged (Adjuster* a, double newval) {
 	
+	Glib::ustring costr;
+    costr = Glib::ustring::format (std::setw(3), std::fixed, std::setprecision(2), a->getValue());
+        
     if (listener && enabled->get_active()) {
-        if (a==luma)
-            listener->panelChanged (EvDPDNLuma, Glib::ustring::format ((int)a->getValue()));
+        if (a==Ldetail)
+            listener->panelChanged (EvDPDNLdetail, costr);
+        else if (a==luma)
+            listener->panelChanged (EvDPDNLuma, costr);
         else if (a==chroma)
-             listener->panelChanged (EvDPDNChroma, Glib::ustring::format ((int)a->getValue()));
-        else if (a==gamma) 
-             listener->panelChanged (EvDPDNGamma, Glib::ustring::format (std::setw(2), std::fixed, std::setprecision(2), a->getValue()));
-    }
+            listener->panelChanged (EvDPDNChroma, costr);
+        else if (a==gamma)
+			listener->panelChanged (EvDPDNGamma, costr);
+	}
 }
 
 void DirPyrDenoise::enabledChanged () {
@@ -150,6 +168,7 @@ void DirPyrDenoise::setBatchMode (bool batchMode) {
 	
     ToolPanel::setBatchMode (batchMode);
     luma->showEditedCB ();
+    Ldetail->showEditedCB ();
     chroma->showEditedCB ();
     gamma->showEditedCB ();
 }
@@ -157,6 +176,7 @@ void DirPyrDenoise::setBatchMode (bool batchMode) {
 void DirPyrDenoise::setAdjusterBehavior (bool chrolumaadd, bool gammaadd) {
 
 	luma->setAddMode(chrolumaadd);
+	Ldetail->setAddMode(chrolumaadd);
 	chroma->setAddMode(chrolumaadd);
 	gamma->setAddMode(gammaadd);
 }
@@ -164,6 +184,7 @@ void DirPyrDenoise::setAdjusterBehavior (bool chrolumaadd, bool gammaadd) {
 void DirPyrDenoise::trimValues (rtengine::procparams::ProcParams* pp) {
 
 	luma->trimValue(pp->dirpyrDenoise.luma);
+	Ldetail->trimValue(pp->dirpyrDenoise.Ldetail);
 	chroma->trimValue(pp->dirpyrDenoise.chroma);
 	gamma->trimValue(pp->dirpyrDenoise.gamma);
 }
diff --git a/rtgui/dirpyrdenoise.h b/rtgui/dirpyrdenoise.h
index 8ffbf294f..0aad76286 100644
--- a/rtgui/dirpyrdenoise.h
+++ b/rtgui/dirpyrdenoise.h
@@ -26,8 +26,9 @@
 class DirPyrDenoise : public Gtk::VBox, public AdjusterListener, public FoldableToolPanel {
 
   protected:
-    Adjuster* luma;
-    Adjuster* chroma;
+	Adjuster* luma;
+	Adjuster* Ldetail;
+	Adjuster* chroma;
 	Adjuster* gamma;
 
     Gtk::CheckButton* enabled;
diff --git a/rtgui/dirpyrequalizer.cc b/rtgui/dirpyrequalizer.cc
index 1ff0540b6..7b8bc7e61 100644
--- a/rtgui/dirpyrequalizer.cc
+++ b/rtgui/dirpyrequalizer.cc
@@ -58,7 +58,7 @@ DirPyrEqualizer::DirPyrEqualizer () : Gtk::VBox(), FoldableToolPanel(this) {
         ss = Glib::ustring::format(i);
         if(i == 0) {
             ss += Glib::ustring::compose(" (%1)", M("TP_DIRPYREQUALIZER_LUMAFINEST"));
-            multiplier[i] = Gtk::manage ( new Adjuster (ss, 0, 4, 0.01, 3.0) );
+            multiplier[i] = Gtk::manage ( new Adjuster (ss, 0, 4, 0.01, 1.0) );
         } else {
 	    if(i == 3)
                 ss += Glib::ustring::compose(" (%1)", M("TP_DIRPYREQUALIZER_LUMACOARSEST"));    
diff --git a/rtgui/paramsedited.cc b/rtgui/paramsedited.cc
index 288f32ddc..6e7390054 100644
--- a/rtgui/paramsedited.cc
+++ b/rtgui/paramsedited.cc
@@ -111,6 +111,7 @@ void ParamsEdited::set (bool v) {
 	impulseDenoise.thresh      = v;
 	dirpyrDenoise.enabled      = v;
 	dirpyrDenoise.luma         = v;
+	dirpyrDenoise.Ldetail      = v;
 	dirpyrDenoise.chroma       = v;
 	dirpyrDenoise.gamma        = v;
 	edgePreservingDecompositionUI.enabled = v;
@@ -312,6 +313,7 @@ void ParamsEdited::initFrom (const std::vector<rtengine::procparams::ProcParams>
 
         dirpyrDenoise.enabled = dirpyrDenoise.enabled && p.dirpyrDenoise.enabled == other.dirpyrDenoise.enabled;
         dirpyrDenoise.luma = dirpyrDenoise.luma && p.dirpyrDenoise.luma == other.dirpyrDenoise.luma;
+        dirpyrDenoise.Ldetail = dirpyrDenoise.Ldetail && p.dirpyrDenoise.Ldetail == other.dirpyrDenoise.Ldetail;
         dirpyrDenoise.chroma = dirpyrDenoise.chroma && p.dirpyrDenoise.chroma == other.dirpyrDenoise.chroma;
         dirpyrDenoise.gamma = dirpyrDenoise.gamma && p.dirpyrDenoise.gamma == other.dirpyrDenoise.gamma;
 
@@ -522,6 +524,7 @@ void ParamsEdited::combine (rtengine::procparams::ProcParams& toEdit, const rten
 
 	if (dirpyrDenoise.enabled)				toEdit.dirpyrDenoise.enabled 	= mods.dirpyrDenoise.enabled;
 	if (dirpyrDenoise.luma)					toEdit.dirpyrDenoise.luma		= dontforceSet && options.baBehav[ADDSET_DIRPYRDN_CHLUM] ? toEdit.dirpyrDenoise.luma + mods.dirpyrDenoise.luma : mods.dirpyrDenoise.luma;
+	if (dirpyrDenoise.Ldetail)				toEdit.dirpyrDenoise.Ldetail	= dontforceSet && options.baBehav[ADDSET_DIRPYRDN_CHLUM] ? toEdit.dirpyrDenoise.Ldetail + mods.dirpyrDenoise.Ldetail : mods.dirpyrDenoise.Ldetail;
 	if (dirpyrDenoise.chroma)				toEdit.dirpyrDenoise.chroma		= dontforceSet && options.baBehav[ADDSET_DIRPYRDN_CHLUM] ? toEdit.dirpyrDenoise.chroma + mods.dirpyrDenoise.chroma : mods.dirpyrDenoise.chroma;
 	if (dirpyrDenoise.gamma)				toEdit.dirpyrDenoise.gamma		= dontforceSet && options.baBehav[ADDSET_DIRPYRDN_GAMMA] ? toEdit.dirpyrDenoise.gamma + mods.dirpyrDenoise.gamma : mods.dirpyrDenoise.gamma;
 
diff --git a/rtgui/paramsedited.h b/rtgui/paramsedited.h
index 08ad94d87..9ad8fb9dd 100644
--- a/rtgui/paramsedited.h
+++ b/rtgui/paramsedited.h
@@ -184,6 +184,7 @@ class DirPyrDenoiseParamsEdited {
 	
 public:
 	bool enabled;
+	bool Ldetail;
 	bool luma;
 	bool chroma;
 	bool gamma;