diff --git a/rtengine/FTblockDN.cc b/rtengine/FTblockDN.cc
index f9d809dfb..44077d996 100644
--- a/rtengine/FTblockDN.cc
+++ b/rtengine/FTblockDN.cc
@@ -35,7 +35,7 @@
 #include "rt_math.h"
 #include "mytime.h"
 #include "sleef.c"
-#include "opthelper.h"
+#include "opthelper.h"
 
 #ifdef _OPENMP
 #include <omp.h>
@@ -244,7 +244,7 @@ void ImProcFunctions::Tile_calc (int tilesize, int overlap, int kall, int imwidt
 int denoiseNestedLevels = 1;
 
 
-void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagefloat * dst,Imagefloat * calclum, float * ch_M, float *max_r, float *max_b, bool isRAW, const procparams::DirPyrDenoiseParams & dnparams, const double expcomp, const NoiseCurve & noiseLCurve, const NoiseCurve & noiseCCurve, float &chaut, float &redaut, float &blueaut, float &maxredaut, float &maxblueaut, float &nresi, float &highresi)
+SSEFUNCTION void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagefloat * dst,Imagefloat * calclum, float * ch_M, float *max_r, float *max_b, bool isRAW, const procparams::DirPyrDenoiseParams & dnparams, const double expcomp, const NoiseCurve & noiseLCurve, const NoiseCurve & noiseCCurve, float &chaut, float &redaut, float &blueaut, float &maxredaut, float &maxblueaut, float &nresi, float &highresi)
 {
 //#ifdef _DEBUG
 	MyTime t1e,t2e;
@@ -278,20 +278,20 @@ void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagefloat * dst,I
 			{wprofi[1][0],wprofi[1][1],wprofi[1][2]},
 			{wprofi[2][0],wprofi[2][1],wprofi[2][2]}
 			};
-		lumcalc = new float*[hei];
-			for (int i=0; i<hei; i++)
-				lumcalc[i] = new float[wid];
-		acalc = new float*[hei];
-			for (int i=0; i<hei; i++)
-				acalc[i] = new float[wid];
-		bcalc = new float*[hei];
-			for (int i=0; i<hei; i++)
-				bcalc[i] = new float[wid];
+		lumcalc = new float*[(hei+1)/2];
+			for (int i=0; i<(hei+1)/2; i++)
+				lumcalc[i] = new float[(wid+1)/2];
+		acalc = new float*[(hei+1)/2];
+			for (int i=0; i<(hei+1)/2; i++)
+				acalc[i] = new float[(wid+1)/2];
+		bcalc = new float*[(hei+1)/2];
+			for (int i=0; i<(hei+1)/2; i++)
+				bcalc[i] = new float[(wid+1)/2];
 #ifdef _OPENMP
 #pragma omp parallel for
 #endif		
-		for(int ii=0;ii<hei;ii++){
-			for(int jj=0;jj<wid;jj++){
+		for(int ii=0;ii<hei;ii+=2){
+			for(int jj=0;jj<wid;jj+=2){
 				float LLum,AAum,BBum;
 			
 				float RL = calclum->r(ii,jj);
@@ -301,9 +301,9 @@ void ImProcFunctions::RGB_denoise(int kall, Imagefloat * src, Imagefloat * dst,I
 				float XL,YL,ZL;
 				Color::rgbxyz(RL,GL,BL,XL,YL,ZL,wpi);
 				Color::XYZ2Lab(XL, YL, ZL, LLum, AAum, BBum);							
-				lumcalc[ii][jj]=LLum;
-				acalc[ii][jj]=AAum;
-				bcalc[ii][jj]=BBum;
+				lumcalc[ii>>1][jj>>1]=LLum;
+				acalc[ii>>1][jj>>1]=AAum;
+				bcalc[ii>>1][jj>>1]=BBum;
 			}
 		}
 	delete calclum;
@@ -416,7 +416,7 @@ if(settings->leveldnti ==1) {
 
 	int numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip;
 	
-	Tile_calc (tilesize, overlap, 2, imwidth, imheight, numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip);
+	Tile_calc (tilesize, overlap, 2, imwidth, imheight, numtiles_W, numtiles_H, tilewidth, tileheight, tileWskip, tileHskip);
 	//now we have tile dimensions, overlaps
 	//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -472,14 +472,12 @@ if(settings->leveldnti ==1) {
 		printf("RGB_denoise uses %d main thread(s) and up to %d nested thread(s) for each main thread\n",numthreads,denoiseNestedLevels);
 #endif
 
-	float *nbrwtArray[denoiseNestedLevels*numthreads];
 	float *LbloxArray[denoiseNestedLevels*numthreads];
 	float *fLbloxArray[denoiseNestedLevels*numthreads];
 
 	for(int i=0;i<denoiseNestedLevels*numthreads;i++) {
 		LbloxArray[i]  = (float*) fftwf_malloc(max_numblox_W*TS*TS*sizeof(float));
 		fLbloxArray[i] = (float*) fftwf_malloc(max_numblox_W*TS*TS*sizeof(float));
-		nbrwtArray[i] = new float[TS*TS];
 	}
 
 	TMatrix wiprof = iccStore->workingSpaceInverseMatrix (params->icm.working);
@@ -511,6 +509,12 @@ if(settings->leveldnti ==1) {
 	float residblue=0.f;
 	
 	int pos;
+	float** noisevarlum = new float*[(tileheight+1)/2];
+	for (int i=0; i<(tileheight+1)/2; i++)
+		noisevarlum[i] = new float[(tilewidth+1)/2];
+	float** noisevarchrom = new float*[(tileheight+1)/2];
+	for (int i=0; i<(tileheight+1)/2; i++)
+		noisevarchrom[i] = new float[(tilewidth+1)/2];
 
 
 #ifdef _OPENMP
@@ -562,17 +566,9 @@ if(settings->leveldnti ==1) {
 			//wavelet denoised image
 			LabImage * labdn = new LabImage(width,height);
 			float* mad_LL = new float [height*width];
-			float** noisevarlum;
-					noisevarlum = new float*[height];
-						for (int i=0; i<height; i++)
-						noisevarlum[i] = new float[width];
 
 			float* mad_aa = new float [height*width];
 			float* mad_bb = new float [height*width];
-			float** noisevarchrom;
-					noisevarchrom = new float*[height];
-						for (int i=0; i<height; i++)
-						noisevarchrom[i] = new float[width];
 
 			//residual between input and denoised L channel
 			array2D<float> Ldetail(width,height,ARRAY2D_CLEAR_DATA);
@@ -626,6 +622,7 @@ if(settings->leveldnti ==1) {
 						R_ = R_<65535.0f ? gamcurve[R_] : (Color::gammanf(R_/65535.f, gam)*32768.0f);
 						G_ = G_<65535.0f ? gamcurve[G_] : (Color::gammanf(G_/65535.f, gam)*32768.0f);
 						B_ = B_<65535.0f ? gamcurve[B_] : (Color::gammanf(B_/65535.f, gam)*32768.0f);
+						
 						//true conversion xyz=>Lab
 						float L,a,b;
 						float X,Y,Z;
@@ -633,35 +630,38 @@ if(settings->leveldnti ==1) {
 
 						//convert to Lab
 						Color::XYZ2Lab(X, Y, Z, L, a, b);
-						
-						if(noiseLCurve) {
-							float kN = lumcalc[i][j]; //with no gamma and take into account working profile
-							float epsi = 0.01f;
-							if(kN<2.f)
-								kN = 2.f;//avoid divided by zero
-							if(kN>32768.f)
-								kN = 32768.f;	// not strictly necessary							
-							float kinterm = epsi + noiseLCurve[xdivf(kN,15)*500.f];
-							float ki = kinterm*100.f;
-							ki += noiseluma;
-							noisevarlum[i1][j1] = SQR((ki/125.f)*(1.f+ki/25.f));
-						}
-						if(noiseCCurve) {
-							float aN = acalc[i][j];
-							float bN = bcalc[i][j];
-							float cN = sqrtf(SQR(aN)+SQR(bN));
-							if(cN < 100.f)
-								cN = 100.f;//avoid divided by zero
-							float Cinterm = 1.f + ponderCC*4.f*noiseCCurve[cN/60.f];//C=f(C)
-							noisevarchrom[i1][j1] = max(noisevarab_b,noisevarab_r)*SQR(Cinterm);
-						}
-						//end chroma
 
 						labdn->L[i1][j1] = L;
 						labdn->a[i1][j1] = a;
 						labdn->b[i1][j1] = b;
 					
-						Lin[i1][j1] = L;
+						Lin[i1][j1] = L;
+
+						if(((i1|j1)&1) == 0) {
+							if(noiseLCurve) {
+								float kN = lumcalc[i>>1][j>>1]; //with no gamma and take into account working profile
+								float epsi = 0.01f;
+								if(kN<2.f)
+									kN = 2.f;//avoid divided by zero
+								if(kN>32768.f)
+									kN = 32768.f;	// not strictly necessary							
+								float kinterm = epsi + noiseLCurve[xdivf(kN,15)*500.f];
+								float ki = kinterm*100.f;
+								ki += noiseluma;
+								noisevarlum[i1>>1][j1>>1] = SQR((ki/125.f)*(1.f+ki/25.f));
+							}
+							if(noiseCCurve) {
+								float aN = acalc[i>>1][j>>1];
+								float bN = bcalc[i>>1][j>>1];
+								float cN = sqrtf(SQR(aN)+SQR(bN));
+								if(cN < 100.f)
+									cN = 100.f;//avoid divided by zero
+								float Cinterm = 1.f + ponderCC*4.f*noiseCCurve[cN/60.f];//C=f(C)
+								noisevarchrom[i1>>1][j1>>1] = max(noisevarab_b,noisevarab_r)*SQR(Cinterm);
+							}
+						}
+						//end chroma
+
 					}
 				}
 				}
@@ -680,29 +680,30 @@ if(settings->leveldnti ==1) {
 						//conversion colorspace to determine luminance with no gamma
 						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);
-						if(noiseLCurve) {
-						//	float noiseluma=(float) dnparams.luma;	
-							float kN = lumcalc[i][j];
-							float epsi = 0.01f;
-							if (kN<2.f) 
-								kN = 2.f;
-							if (kN>32768.f)
-								kN = 32768.f;
-							float kinterm = epsi + noiseLCurve[xdivf(kN,15)*500.f];
-							float ki = kinterm*100.f;
-							ki += noiseluma;
-							noisevarlum[i1][j1] = SQR((ki/125.f)*(1.f+ki/25.f));
-						}
-						if(noiseCCurve) {
-							float aN = acalc[i][j];
-							float bN = bcalc[i][j];
-							
-							float cN = sqrtf(SQR(aN)+SQR(bN));
-							if(cN < 100.f)
-								cN = 100.f;//avoid divided by zero
-							float Cinterm = 1.f + ponderCC*4.f*noiseCCurve[cN/60.f];
-							noisevarchrom[i1][j1] = max(noisevarab_b,noisevarab_r)*SQR(Cinterm);
+						Z = Z<65535.0f ? gamcurve[Z] : (Color::gamma((double)Z/65535.0, gam, gamthresh, gamslope, 1.0, 0.0)*32768.0f);
+						if(((i1|j1)&1) == 0) {
+							if(noiseLCurve) {
+							//	float noiseluma=(float) dnparams.luma;	
+								float kN = lumcalc[i>>1][j>>1];
+								float epsi = 0.01f;
+								if (kN<2.f) 
+									kN = 2.f;
+								if (kN>32768.f)
+									kN = 32768.f;
+								float kinterm = epsi + noiseLCurve[xdivf(kN,15)*500.f];
+								float ki = kinterm*100.f;
+								ki += noiseluma;
+								noisevarlum[i1>>1][j1>>1] = SQR((ki/125.f)*(1.f+ki/25.f));
+							}
+							if(noiseCCurve) {
+								float aN = acalc[i>>1][j>>1];
+								float bN = bcalc[i>>1][j>>1];
+								float cN = sqrtf(SQR(aN)+SQR(bN));
+								if(cN < 100.f)
+									cN = 100.f;//avoid divided by zero
+								float Cinterm = 1.f + ponderCC*4.f*noiseCCurve[cN/60.f];
+								noisevarchrom[i1>>1][j1>>1] = max(noisevarab_b,noisevarab_r)*SQR(Cinterm);
+							}
 						}
 						//end chroma					
 						labdn->L[i1][j1] = Y;
@@ -745,35 +746,35 @@ if(settings->leveldnti ==1) {
 
 						//convert Lab
 						Color::XYZ2Lab(X, Y, Z, L, a, b);
-						float Llum,alum,blum;
-						if(noiseLCurve || noiseCCurve) {
-							float XL,YL,ZL;
-							Color::rgbxyz(rLum,gLum,bLum,XL,YL,ZL,wp);
-							Color::XYZ2Lab(XL, YL, ZL, Llum, alum, blum);
+						float Llum,alum,blum;
+						if(((i1|j1)&1) == 0) {
+							if(noiseLCurve || noiseCCurve) {
+								float XL,YL,ZL;
+								Color::rgbxyz(rLum,gLum,bLum,XL,YL,ZL,wp);
+								Color::XYZ2Lab(XL, YL, ZL, Llum, alum, blum);
+							}
+					
+							if(noiseLCurve) {
+								float kN = Llum;
+								float epsi=0.01f;
+								if(kN<2.f) kN=2.f;
+								if(kN>32768.f) kN=32768.f;								
+								float kinterm=epsi + noiseLCurve[xdivf(kN,15)*500.f];
+								float ki=kinterm*100.f;
+								ki+=noiseluma;
+								noiseluma += 1.f;
+								noisevarlum[i1>>1][j1>>1]=SQR((ki/125.f)*(1.f+ki/25.f));
+							}	
+							if(noiseCCurve) {
+								float aN=alum;
+								float bN=blum;
+								float cN=sqrtf(SQR(aN)+SQR(bN));
+								if(cN < 100.f)
+									cN=100.f;//avoid divided by zero ???
+								float Cinterm=1.f + ponderCC*4.f*noiseCCurve[cN/60.f];
+								noisevarchrom[i1>>1][j1>>1]=max(noisevarab_b,noisevarab_r)*SQR(Cinterm);
+							}
 						}
-				
-						if(noiseLCurve) {
-							float kN = Llum;
-							float epsi=0.01f;
-
-							if(kN<2.f) kN=2.f;
-							if(kN>32768.f) kN=32768.f;								
-							float kinterm=epsi + noiseLCurve[xdivf(kN,15)*500.f];
-							float ki=kinterm*100.f;
-							ki+=noiseluma;
-							noiseluma += 1.f;
-							noisevarlum[i1][j1]=SQR((ki/125.f)*(1.f+ki/25.f));
-						}	
-						if(noiseCCurve) {
-							float aN=alum;
-							float bN=blum;
-							float cN=sqrtf(SQR(aN)+SQR(bN));
-							if(cN < 100.f)
-								cN=100.f;//avoid divided by zero ???
-							float Cinterm=1.f + ponderCC*4.f*noiseCCurve[cN/60.f];
-							noisevarchrom[i1][j1]=max(noisevarab_b,noisevarab_r)*SQR(Cinterm);
-						}
-						
 						labdn->L[i1][j1] = L;
 						labdn->a[i1][j1] = a;
 						labdn->b[i1][j1] = b;
@@ -838,13 +839,18 @@ if(settings->leveldnti ==1) {
 #pragma omp section
 #endif
 {
-			Ldecomp = new wavelet_decomposition (labdn->data, labdn->W, labdn->H, levwav/*maxlevels*/, 0/*subsampling*/ );
+			Ldecomp = new wavelet_decomposition (labdn->data, labdn->W, labdn->H, levwav/*maxlevels*/, 1/*subsampling*/ );
 }
 #ifdef _OPENMP
 #pragma omp section
 #endif
-{			// only one section for both, because they are much faster then Ldecomp
+{
 			adecomp = new wavelet_decomposition (labdn->data+datalen, labdn->W, labdn->H,levwav, 1 );
+}
+#ifdef _OPENMP
+#pragma omp section
+#endif
+{
 			bdecomp = new wavelet_decomposition (labdn->data+2*datalen, labdn->W, labdn->H, levwav, 1 );
 }
 }
@@ -1231,10 +1237,6 @@ if(settings->leveldnti ==1) {
 		}
 
 			//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
@@ -1269,21 +1271,19 @@ if(settings->leveldnti ==1) {
 #else
 			int subThread = 0;
 #endif
-			AlignedBuffer<float> pBuf(width + TS + 2*blkrad*offset);
-			float * nbrwt = nbrwtArray[subThread];
+			float blurbuffer[TS*TS] ALIGNED64;
 			float * Lblox = LbloxArray[subThread];
-			float * fLblox = fLbloxArray[subThread];;
-			float * blurbuffer = new float[TS*TS];
+			float * fLblox = fLbloxArray[subThread];
+			float pBuf[width + TS + 2*blkrad*offset] ALIGNED16;
+			float nbrwt[TS*TS] ALIGNED64;
 #ifdef _OPENMP
 #pragma omp for
 #endif
 			for (int vblk=0; vblk<numblox_H; vblk++) {
-				//printf("vblock=%d",vblk);
 
 				int top = (vblk-blkrad)*offset;
-				float * datarow = (float*)pBuf.data +blkrad*offset;
+				float * datarow = pBuf +blkrad*offset;
 
-				//TODO: implement using AlignedBufferMP
 				for (int i=0/*, row=top*/; i<TS; i++/*, row++*/) {
 					int row = top + i;
 					int rr = row;
@@ -1294,7 +1294,7 @@ if(settings->leveldnti ==1) {
 					}
 
 					for (int j=0; j<labdn->W; j++) {
-						datarow[j] = (Lin[rr][j]-Lwavdn[rr][j]);
+						datarow[j] = (Lin[rr][j]-labdn->L[rr][j]);
 					}
 
 					for (int j=-blkrad*offset; j<0; j++) {
@@ -1305,13 +1305,12 @@ if(settings->leveldnti ==1) {
 					}//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
+							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];
 							}
@@ -1350,7 +1349,6 @@ if(settings->leveldnti ==1) {
 
 			}//end of vertical block loop
 			//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-			delete [] blurbuffer;
 
 }
 			}
@@ -1360,7 +1358,7 @@ if(settings->leveldnti ==1) {
 					//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;
+					labdn->L[i][j] += hpdn;
 
 				}
 			}
@@ -1369,8 +1367,8 @@ if(settings->leveldnti ==1) {
 			// transform denoised "Lab" to output RGB
 
 			//calculate mask for feathering output tile overlaps
-			float * Vmask = new float [height+1];
-			float * Hmask = new float [width+1];
+			float Vmask[height+1] ALIGNED16;
+			float Hmask[width+1] ALIGNED16;
 
 			for (int i=0; i<height; i++) {
 				Vmask[i] = 1;
@@ -1383,7 +1381,7 @@ if(settings->leveldnti ==1) {
 			}
 				
 			for (int i=0; i<overlap; i++) {
-				float mask = SQR(sin((M_PI*i)/(2*overlap)));
+				float mask = SQR(xsinf((M_PI*i)/(2*overlap)));
 				if (tiletop>0) Vmask[i] = mask;
 				if (tilebottom<imheight) Vmask[height-i] = mask;
 				if (tileleft>0) Hmask[i] = mask/newGain;
@@ -1514,28 +1512,23 @@ if(settings->leveldnti ==1) {
 			//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 			delete labdn;
-			for (int i=0; i<height; i++)
-				delete [] noisevarlum[i];
-			delete [] noisevarlum;
 			delete [] mad_LL;
-			for (int i=0; i<height; i++)
-				delete [] noisevarchrom[i];
-			delete [] noisevarchrom;
 			delete [] mad_aa;
 			delete [] mad_bb;
 
-			delete[] Vmask;
-			delete[] Hmask;
-
-
-
 		}//end of tile row
 	}//end of tile loop
+	for (int i=0; i<(tileheight+1)/2; i++)
+		delete [] noisevarlum[i];
+	delete [] noisevarlum;
+	for (int i=0; i<(tileheight+1)/2; i++)
+		delete [] noisevarchrom[i];
+	delete [] noisevarchrom;
+
 	}
 	for(int i=0;i<denoiseNestedLevels*numthreads;i++) {
 		fftwf_free(LbloxArray[i]);
 		fftwf_free(fLbloxArray[i]);
-		delete [] nbrwtArray[i];
 	}
 
 #ifdef _OPENMP
@@ -1733,13 +1726,13 @@ for(int iteration=1;iteration<=dnparams.passes;iteration++){
 }		
 	//end median
 	if(noiseLCurve || noiseCCurve) {
-			for (int i=0; i<hei; i++)
+			for (int i=0; i<(hei+1)/2; i++)
 				delete [] lumcalc[i];
 			delete [] lumcalc;
-			for (int i=0; i<hei; i++)
+			for (int i=0; i<(hei+1)/2; i++)
 				delete [] acalc[i];
 			delete [] acalc;
-			for (int i=0; i<hei; i++)
+			for (int i=0; i<(hei+1)/2; i++)
 				delete [] bcalc[i];
 			delete [] bcalc;
 			
@@ -1771,7 +1764,7 @@ SSEFUNCTION	void ImProcFunctions::RGBtile_denoise (float * fLblox, int hblproc,
 		__m128	noisevar_Ldetailv = _mm_set1_ps( noisevar_Ldetail );
 		__m128	onev = _mm_set1_ps( 1.0f );
 		for (int n=0; n<TS*TS; n+=4) {		//for DCT
-			tempv  = onev - xexpf( -SQRV( LVFU(nbrwt[n]))/noisevar_Ldetailv);
+			tempv  = onev - xexpf( -SQRV( LVF(nbrwt[n]))/noisevar_Ldetailv);
 			_mm_storeu_ps( &fLblox[blkstart+n], LVFU(fLblox[blkstart+n]) * tempv );
 		}//output neighbor averaged result
 #else
@@ -1960,29 +1953,24 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 
 		float madL[8][3], mada[8][3], madb[8][3];
 		
-		int maxWL = 0, maxHL = 0, maxWab = 0, maxHab = 0;
+		int maxWL = 0, maxHL = 0;
 		for (int lvl=0; lvl<maxlvl; lvl++) {
 			if(WaveletCoeffs_L.level_W(lvl) > maxWL)
 				maxWL = WaveletCoeffs_L.level_W(lvl);
 			if(WaveletCoeffs_L.level_H(lvl) > maxHL)
 				maxHL = WaveletCoeffs_L.level_H(lvl);
-			if(WaveletCoeffs_a.level_W(lvl) > maxWab)
-				maxWab = WaveletCoeffs_a.level_W(lvl);
-			if(WaveletCoeffs_a.level_H(lvl) > maxHab)
-				maxHab = WaveletCoeffs_a.level_H(lvl);
 		}
 
 #ifdef _OPENMP
 #pragma omp parallel num_threads(denoiseNestedLevels) if(denoiseNestedLevels>1)
 #endif
 {
-		float *buffer[6];
+		float *buffer[5];
 		buffer[0] = new float[maxWL*maxHL+32];
 		buffer[1] = new float[maxWL*maxHL+64];
 		buffer[2] = new float[maxWL*maxHL+96];
-		buffer[3] = new float[maxWab*maxHab+128];
+		buffer[3] = new float[maxWL*maxHL+128];
 		buffer[4] = new float[maxWL*maxHL+160];
-		buffer[5] = new float[maxWL*maxHL+192];
 		
 
 #ifdef _OPENMP
@@ -1994,9 +1982,6 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 			int Wlvl_L = WaveletCoeffs_L.level_W(lvl);
 			int Hlvl_L = WaveletCoeffs_L.level_H(lvl);
 			
-			if(Wlvl_L!=width) Wlvl_L=width;
-			if(Hlvl_L!=height) Hlvl_L=height;
-			
 			int Wlvl_ab = WaveletCoeffs_a.level_W(lvl);
 			int Hlvl_ab = WaveletCoeffs_a.level_H(lvl);
 
@@ -2024,9 +2009,6 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 			int Wlvl_L = WaveletCoeffs_L.level_W(lvl);
 			int Hlvl_L = WaveletCoeffs_L.level_H(lvl);
 			
-			if(Wlvl_L!=width) Wlvl_L=width;
-			if(Hlvl_L!=height) Hlvl_L=height;
-			
 			int Wlvl_ab = WaveletCoeffs_a.level_W(lvl);
 			int Hlvl_ab = WaveletCoeffs_a.level_H(lvl);
 
@@ -2044,8 +2026,7 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 			} else {
 				//simple wavelet shrinkage
 				float * sfave = buffer[0]+32;
-				float * sfaved = buffer[4]+160;
-				float * WavCoeffsLtemp = buffer[3]+128;
+				float * sfaved = buffer[3]+128;
 				float * blurBuffer = buffer[1]+64;
 				for (int dir=1; dir<4; dir++) {
 					float mad_Lr = madL[lvl][dir-1];
@@ -2070,48 +2051,42 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 						}
 					}
 					if (noisevar_abr>0.001f  || noisevar_abb>0.001f) {
-						for(int i=0;i<Hlvl_ab;i++)
-							for(int j=0;j<Wlvl_ab;j++)
-								WavCoeffsLtemp[i*Wlvl_ab+j] = WavCoeffs_L[dir][((i*skip_L)/skip_ab)*Wlvl_L + ((j*skip_L)/skip_ab)];
 					
 #ifdef __SSE2__
-						int j;
 						__m128 onev = _mm_set1_ps(1.f);
 						__m128 rmad_Lm9v = onev / _mm_set1_ps(mad_Lr * 9.f);
 						__m128 mad_av;
 						__m128 mad_bv;
 						__m128 mag_Lv, mag_av, mag_bv;
-						__m128 tempav, tempbv;
-						for (int i=0; i<Hlvl_ab; i++) {
-							int coeffloc_ab = i*Wlvl_ab;
-							for (j=0; j<Wlvl_ab-3; j+=4, coeffloc_ab+=4) {
-								mad_av = LVFU(mad_aa[coeffloc_ab]);
-								mad_bv = LVFU(mad_bb[coeffloc_ab]);
-							
-								tempav = LVFU(WavCoeffs_a[dir][coeffloc_ab]);
-								tempbv = LVFU(WavCoeffs_b[dir][coeffloc_ab]);
-								mag_Lv = LVFU(WavCoeffsLtemp[coeffloc_ab]);
-								mag_av = SQRV(tempav);
-								mag_bv = SQRV(tempbv);
-								mag_Lv = SQRV(mag_Lv) * rmad_Lm9v;
-								_mm_storeu_ps(&WavCoeffs_a[dir][coeffloc_ab], tempav * SQRV((onev-xexpf(-(mag_av/mad_av)-(mag_Lv)))));
-								_mm_storeu_ps(&WavCoeffs_b[dir][coeffloc_ab], tempbv * SQRV((onev-xexpf(-(mag_bv/mad_bv)-(mag_Lv)))));
-							}
-							for (; j<Wlvl_ab; j++,coeffloc_ab++) {
-								float mag_L = SQR(WavCoeffsLtemp[coeffloc_ab ]);
-								float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab]);
-								float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab]);
-								WavCoeffs_a[dir][coeffloc_ab] *= SQR(1.f-xexpf(-(mag_a/mad_aa[coeffloc_ab])-(mag_L/(9.f*mad_Lr)))/*satfactor_a*/);
-								WavCoeffs_b[dir][coeffloc_ab] *= SQR(1.f-xexpf(-(mag_b/mad_bb[coeffloc_ab])-(mag_L/(9.f*mad_Lr)))/*satfactor_b*/);
-							}
+						__m128 tempav, tempbv;
+						int coeffloc_ab;
+						for (coeffloc_ab=0; coeffloc_ab<Hlvl_ab*Wlvl_ab-3; coeffloc_ab+=4) {
+							mad_av = LVFU(mad_aa[coeffloc_ab]);
+							mad_bv = LVFU(mad_bb[coeffloc_ab]);
+						
+							tempav = LVFU(WavCoeffs_a[dir][coeffloc_ab]);
+							tempbv = LVFU(WavCoeffs_b[dir][coeffloc_ab]);
+							mag_Lv = LVFU(WavCoeffs_L[dir][coeffloc_ab]);
+							mag_av = SQRV(tempav);
+							mag_bv = SQRV(tempbv);
+							mag_Lv = SQRV(mag_Lv) * rmad_Lm9v;
+							_mm_storeu_ps(&WavCoeffs_a[dir][coeffloc_ab], tempav * SQRV((onev-xexpf(-(mag_av/mad_av)-(mag_Lv)))));
+							_mm_storeu_ps(&WavCoeffs_b[dir][coeffloc_ab], tempbv * SQRV((onev-xexpf(-(mag_bv/mad_bv)-(mag_Lv)))));
+						}
+						// few remaining pixels
+						for (; coeffloc_ab<Hlvl_ab*Wlvl_ab; coeffloc_ab++) {
+							float mag_L = SQR(WavCoeffs_L[dir][coeffloc_ab ]);
+							float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab]);
+							float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab]);
+							WavCoeffs_a[dir][coeffloc_ab] *= SQR(1.f-xexpf(-(mag_a/mad_aa[coeffloc_ab])-(mag_L/(9.f*mad_Lr)))/*satfactor_a*/);
+							WavCoeffs_b[dir][coeffloc_ab] *= SQR(1.f-xexpf(-(mag_b/mad_bb[coeffloc_ab])-(mag_L/(9.f*mad_Lr)))/*satfactor_b*/);
 						}//now chrominance coefficients are denoised
 #else
 						for (int i=0; i<Hlvl_ab; i++) {
 							for (int j=0; j<Wlvl_ab; j++) {
-
 								int coeffloc_ab = i*Wlvl_ab+j;
 
-								float mag_L = SQR(WavCoeffsLtemp[coeffloc_ab ])+eps;
+								float mag_L = SQR(WavCoeffs_L[dir][coeffloc_ab ])+eps;
 								float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab])+eps;
 								float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab])+eps;
 
@@ -2142,23 +2117,19 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 							}
 						}
 #ifdef __SSE2__
-						int j;
 						__m128 mad_Lv;
 						__m128 ninev = _mm_set1_ps( 9.0f );
 						__m128 epsv = _mm_set1_ps(eps);
-						__m128 mag_Lv;
-						for (int i=0; i<Hlvl_L; i++) {
-							int coeffloc_L = i*Wlvl_L;
-							for (j=0; j<Wlvl_L-3; j+=4,coeffloc_L+=4) {
-								mad_Lv = LVFU(mad_LL[coeffloc_L]);
-								mag_Lv = SQRV(LVFU(WavCoeffs_L[dir][coeffloc_L]));
-								_mm_storeu_ps(&sfave[coeffloc_L], mag_Lv / ( mag_Lv + mad_Lv * xexpf(-mag_Lv/(mad_Lv*ninev) )+ epsv));
-								
-							}	
-							for (; j<Wlvl_L; j++, coeffloc_L++) {
-								float mag_L = SQR(WavCoeffs_L[dir][coeffloc_L]);
-								sfave[coeffloc_L] = mag_L/(mag_L+mad_LL[coeffloc_L]*xexpf(-mag_L/(9.f*mad_LL[coeffloc_L]))+eps);
-							}
+						__m128 mag_Lv;
+						int coeffloc_L;
+						for (coeffloc_L=0; coeffloc_L<Hlvl_L*Wlvl_L-3; coeffloc_L+=4) {
+							mad_Lv = LVFU(mad_LL[coeffloc_L]);
+							mag_Lv = SQRV(LVFU(WavCoeffs_L[dir][coeffloc_L]));
+							_mm_storeu_ps(&sfave[coeffloc_L], mag_Lv / ( mag_Lv + mad_Lv * xexpf(-mag_Lv/(mad_Lv*ninev) )+ epsv));
+						}	
+						for (; coeffloc_L<Hlvl_L*Wlvl_L; coeffloc_L++) {
+							float mag_L = SQR(WavCoeffs_L[dir][coeffloc_L]);
+							sfave[coeffloc_L] = mag_L/(mag_L+mad_LL[coeffloc_L]*xexpf(-mag_L/(9.f*mad_LL[coeffloc_L]))+eps);
 						}
 #else
 						for (int i=0; i<Hlvl_L; i++)
@@ -2173,20 +2144,18 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 #ifdef __SSE2__
 						__m128 sfavev;
 						__m128 sf_Lv;
-						for (int i=0; i<Hlvl_L; i++) {
-							int coeffloc_L = i*Wlvl_L;
-							for (j=0; j<Wlvl_L-3; j+=4,coeffloc_L+=4) {
-								sfavev = LVFU(sfaved[coeffloc_L]);
-								sf_Lv = LVFU(sfave[coeffloc_L]);
-								_mm_storeu_ps(&WavCoeffs_L[dir][coeffloc_L], LVFU(WavCoeffs_L[dir][coeffloc_L]) * (SQRV(sfavev)+SQRV(sf_Lv))/(sfavev+sf_Lv+epsv));
-								//use smoothed shrinkage unless local shrinkage is much less
-							}//now luminance coeffs are denoised
-							for (; j<Wlvl_L; j++,coeffloc_L++) {
-								float sf_L = sfave[coeffloc_L];
-								//use smoothed shrinkage unless local shrinkage is much less
-								WavCoeffs_L[dir][coeffloc_L] *= (SQR(sfaved[coeffloc_L])+SQR(sf_L))/(sfaved[coeffloc_L]+sf_L+eps);
-							}//now luminance coeffs are denoised
-						}
+						for (coeffloc_L=0; coeffloc_L<Hlvl_L*Wlvl_L-3; coeffloc_L+=4) {
+							sfavev = LVFU(sfaved[coeffloc_L]);
+							sf_Lv = LVFU(sfave[coeffloc_L]);
+							_mm_storeu_ps(&WavCoeffs_L[dir][coeffloc_L], LVFU(WavCoeffs_L[dir][coeffloc_L]) * (SQRV(sfavev)+SQRV(sf_Lv))/(sfavev+sf_Lv+epsv));
+							//use smoothed shrinkage unless local shrinkage is much less
+						}
+						// few remaining pixels
+						for (; coeffloc_L<Hlvl_L*Wlvl_L; coeffloc_L++) {
+							float sf_L = sfave[coeffloc_L];
+							//use smoothed shrinkage unless local shrinkage is much less
+							WavCoeffs_L[dir][coeffloc_L] *= (SQR(sfaved[coeffloc_L])+SQR(sf_L))/(sfaved[coeffloc_L]+sf_L+eps);
+						}//now luminance coeffs are denoised
 #else
 						for (int i=0; i<Hlvl_L; i++)
 							for (int j=0; j<Wlvl_L; j++) {
@@ -2194,14 +2163,13 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 								float sf_L = sfave[coeffloc_L];
 								//use smoothed shrinkage unless local shrinkage is much less
 								WavCoeffs_L[dir][coeffloc_L] *= (SQR(sfaved[coeffloc_L])+SQR(sf_L))/(sfaved[coeffloc_L]+sf_L+eps);
-
 							}//now luminance coeffs are denoised
 #endif
 					}
 				}
 			}
 		}
-		for(int i=5;i>=0;i--)
+		for(int i=4;i>=0;i--)
 			delete [] buffer[i];
 
 }
@@ -2214,16 +2182,12 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 	{
 
 		int maxlvl = WaveletCoeffs_L.maxlevel();
-		int maxWL = 0, maxHL = 0, maxWab = 0, maxHab = 0;
+		int maxWL = 0, maxHL = 0;
 		for (int lvl=0; lvl<maxlvl; lvl++) {
 			if(WaveletCoeffs_L.level_W(lvl) > maxWL)
 				maxWL = WaveletCoeffs_L.level_W(lvl);
 			if(WaveletCoeffs_L.level_H(lvl) > maxHL)
 				maxHL = WaveletCoeffs_L.level_H(lvl);
-			if(WaveletCoeffs_a.level_W(lvl) > maxWab)
-				maxWab = WaveletCoeffs_a.level_W(lvl);
-			if(WaveletCoeffs_a.level_H(lvl) > maxHab)
-				maxHab = WaveletCoeffs_a.level_H(lvl);
 		}
 
 
@@ -2231,13 +2195,12 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 #pragma omp parallel num_threads(denoiseNestedLevels) if(denoiseNestedLevels>1)
 #endif
 {
-		float *buffer[6];
+		float *buffer[5];
 		buffer[0] = new float[maxWL*maxHL+32];
 		buffer[1] = new float[maxWL*maxHL+64];
 		buffer[2] = new float[maxWL*maxHL+96];
-		buffer[3] = new float[maxWab*maxHab+128];
+		buffer[3] = new float[maxWL*maxHL+128];
 		buffer[4] = new float[maxWL*maxHL+160];
-		buffer[5] = new float[maxWL*maxHL+192];
 		
 		
 #ifdef _OPENMP
@@ -2268,7 +2231,7 @@ SSEFUNCTION	void ImProcFunctions::WaveletDenoiseAll_BiShrink(wavelet_decompositi
 					  skip_L, skip_ab, noisevar_L, noisevarlum,  noisevarchrom, width, height, mad_LL, mad_aa, mad_bb, noisevar_abr, noisevar_abb, noi, noiseLCurve, noiseCCurve, chaut, redaut, blueaut, maxredaut, maxblueaut, callby, autoch, perf);
 
 		}
-		for(int i=5;i>=0;i--)
+		for(int i=4;i>=0;i--)
 			delete [] buffer[i];
 }
 	}
@@ -2282,18 +2245,17 @@ SSEFUNCTION	void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
 									{
 		//simple wavelet shrinkage
 		const float eps = 0.01f;
-		W_L = width;
-		H_L = height;
-		if(autoch && noisevar_abr <=0.001f) noisevar_abr=0.02f;
-		if(autoch && noisevar_abb <=0.001f) noisevar_abb=0.02f;
+		if(autoch && noisevar_abr <=0.001f)
+			noisevar_abr=0.02f;
+		if(autoch && noisevar_abb <=0.001f)
+			noisevar_abb=0.02f;
 		
 		float * sfavea = buffer[0]+32;
 		float * sfaveb = buffer[1]+64;
-		float * sfavead = buffer[4]+160;
-		float * sfavebd = buffer[5]+192;
+		float * sfavead = buffer[3]+128;
+		float * sfavebd = buffer[4]+160;
 		float * sfave = sfavea; // we can safely reuse sfavea here, because they are not used together
 		float * sfaved = sfavead; // we can safely reuse sfavead here, because they are not used together
-		float * WavCoeffsLtemp = buffer[3]+128;
 		float * blurBuffer = buffer[2]+96;
 
 		for (int dir=1; dir<4; dir++) {
@@ -2332,39 +2294,31 @@ SSEFUNCTION	void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
 			}
 			
 			if (noisevar_abr>0.001f  ||  noisevar_abb>0.001f ) {
-
-				for(int i=0;i<H_ab;i++)
-					for(int j=0;j<W_ab;j++)
-						WavCoeffsLtemp[i*W_ab+j] = WavCoeffs_L[dir][((i*skip_L)/skip_ab)*W_L + ((j*skip_L)/skip_ab)];
-
 #ifdef __SSE2__
-				int j;
 				__m128 onev = _mm_set1_ps(1.f);
 				__m128 rmadLm9v = onev / _mm_set1_ps(madL * 9.f);
 				__m128 mad_av ;
 				__m128 mad_bv;
-				__m128 mag_Lv, mag_av, mag_bv;
-				for (int i=0; i<H_ab; i++) {
-					int coeffloc_ab = i*W_ab;
-					for (j=0; j<W_ab-3; j+=4, coeffloc_ab+=4) {
-						mad_av = LVFU(mad_aa[coeffloc_ab]);
-						mad_bv = LVFU(mad_bb[coeffloc_ab]);
+				__m128 mag_Lv, mag_av, mag_bv;
+				int coeffloc_ab;
+				for (coeffloc_ab=0; coeffloc_ab<H_ab*W_ab-3; coeffloc_ab+=4) {
+					mad_av = LVFU(mad_aa[coeffloc_ab]);
+					mad_bv = LVFU(mad_bb[coeffloc_ab]);
 
-						mag_Lv = LVFU(WavCoeffsLtemp[coeffloc_ab]);
-						mag_av = SQRV(LVFU(WavCoeffs_a[dir][coeffloc_ab]));
-						mag_bv = SQRV(LVFU(WavCoeffs_b[dir][coeffloc_ab]));
-						mag_Lv = (SQRV(mag_Lv)) * rmadLm9v;
-						_mm_storeu_ps(&sfavea[coeffloc_ab], (onev-xexpf(-(mag_av/mad_av)-(mag_Lv))));
-						_mm_storeu_ps(&sfaveb[coeffloc_ab], (onev-xexpf(-(mag_bv/mad_bv)-(mag_Lv))));
-
-					}
-					for (; j<W_ab; j++,coeffloc_ab++) {
-						float mag_L = SQR(WavCoeffsLtemp[coeffloc_ab]);
-						float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab]);
-						float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab]);
-						sfavea[coeffloc_ab] = (1.f-xexpf(-(mag_a/mad_aa[coeffloc_ab])-(mag_L/(9.f*madL))));
-						sfaveb[coeffloc_ab] = (1.f-xexpf(-(mag_b/mad_bb[coeffloc_ab])-(mag_L/(9.f*madL))));
-					}
+					mag_Lv = LVFU(WavCoeffs_L[dir][coeffloc_ab]);
+					mag_av = SQRV(LVFU(WavCoeffs_a[dir][coeffloc_ab]));
+					mag_bv = SQRV(LVFU(WavCoeffs_b[dir][coeffloc_ab]));
+					mag_Lv = (SQRV(mag_Lv)) * rmadLm9v;
+					_mm_storeu_ps(&sfavea[coeffloc_ab], (onev-xexpf(-(mag_av/mad_av)-(mag_Lv))));
+					_mm_storeu_ps(&sfaveb[coeffloc_ab], (onev-xexpf(-(mag_bv/mad_bv)-(mag_Lv))));
+				}
+				// few remaining pixels
+				for (; coeffloc_ab<H_ab*W_ab; coeffloc_ab++) {
+					float mag_L = SQR(WavCoeffs_L[dir][coeffloc_ab]);
+					float mag_a = SQR(WavCoeffs_a[dir][coeffloc_ab]);
+					float mag_b = SQR(WavCoeffs_b[dir][coeffloc_ab]);
+					sfavea[coeffloc_ab] = (1.f-xexpf(-(mag_a/mad_aa[coeffloc_ab])-(mag_L/(9.f*madL))));
+					sfaveb[coeffloc_ab] = (1.f-xexpf(-(mag_b/mad_bb[coeffloc_ab])-(mag_L/(9.f*madL))));
 				}//now chrominance coefficients are denoised
 
 #else
@@ -2388,30 +2342,26 @@ SSEFUNCTION	void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
 				__m128 epsv = _mm_set1_ps(eps);
 				__m128 sfav, sfbv;
 				__m128 sfaveav, sfavebv;
-				for (int i=0; i<H_ab; i++) {
-					int coeffloc_ab = i*W_ab;
-					for (j=0; j<W_ab-3; j+=4,coeffloc_ab+=4) {
-						sfav = LVFU(sfavea[coeffloc_ab]);
-						sfbv = LVFU(sfaveb[coeffloc_ab]);
-						sfaveav = LVFU(sfavead[coeffloc_ab]);
-						sfavebv = LVFU(sfavebd[coeffloc_ab]);
+				for (coeffloc_ab=0; coeffloc_ab<H_ab*W_ab-3; coeffloc_ab+=4) {
+					sfav = LVFU(sfavea[coeffloc_ab]);
+					sfbv = LVFU(sfaveb[coeffloc_ab]);
+					sfaveav = LVFU(sfavead[coeffloc_ab]);
+					sfavebv = LVFU(sfavebd[coeffloc_ab]);
 
-						//use smoothed shrinkage unless local shrinkage is much less
-						_mm_storeu_ps( &WavCoeffs_a[dir][coeffloc_ab], LVFU(WavCoeffs_a[dir][coeffloc_ab]) * (SQRV(sfaveav)+SQRV(sfav))/(sfaveav+sfav+epsv));
-						_mm_storeu_ps( &WavCoeffs_b[dir][coeffloc_ab], LVFU(WavCoeffs_b[dir][coeffloc_ab]) * (SQRV(sfavebv)+SQRV(sfbv))/(sfavebv+sfbv+epsv));
-
-					}//now chrominance coefficients are denoised
-					for (; j<W_ab; j++,coeffloc_ab++) {
-						//modification Jacques feb 2013
-						float sfa = sfavea[coeffloc_ab];
-						float sfb = sfaveb[coeffloc_ab];
-
-						//use smoothed shrinkage unless local shrinkage is much less
-						WavCoeffs_a[dir][coeffloc_ab] *= (SQR(sfavead[coeffloc_ab])+SQR(sfa))/(sfavead[coeffloc_ab]+sfa+eps);
-						WavCoeffs_b[dir][coeffloc_ab] *= (SQR(sfavebd[coeffloc_ab])+SQR(sfb))/(sfavebd[coeffloc_ab]+sfb+eps);
-
-					}//now chrominance coefficients are denoised
+					//use smoothed shrinkage unless local shrinkage is much less
+					_mm_storeu_ps( &WavCoeffs_a[dir][coeffloc_ab], LVFU(WavCoeffs_a[dir][coeffloc_ab]) * (SQRV(sfaveav)+SQRV(sfav))/(sfaveav+sfav+epsv));
+					_mm_storeu_ps( &WavCoeffs_b[dir][coeffloc_ab], LVFU(WavCoeffs_b[dir][coeffloc_ab]) * (SQRV(sfavebv)+SQRV(sfbv))/(sfavebv+sfbv+epsv));
 				}
+				// few remaining pixels
+				for (; coeffloc_ab<H_ab*W_ab; coeffloc_ab++) {
+					//modification Jacques feb 2013
+					float sfa = sfavea[coeffloc_ab];
+					float sfb = sfaveb[coeffloc_ab];
+
+					//use smoothed shrinkage unless local shrinkage is much less
+					WavCoeffs_a[dir][coeffloc_ab] *= (SQR(sfavead[coeffloc_ab])+SQR(sfa))/(sfavead[coeffloc_ab]+sfa+eps);
+					WavCoeffs_b[dir][coeffloc_ab] *= (SQR(sfavebd[coeffloc_ab])+SQR(sfb))/(sfavebd[coeffloc_ab]+sfb+eps);
+				}//now chrominance coefficients are denoised
 #else
 				for (int i=0; i<H_ab; i++) {
 					for (int j=0; j<W_ab; j++) {
@@ -2422,7 +2372,6 @@ SSEFUNCTION	void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
 						//use smoothed shrinkage unless local shrinkage is much less
 						WavCoeffs_a[dir][coeffloc_ab] *= (SQR(sfavead[coeffloc_ab])+SQR(sfa))/(sfavead[coeffloc_ab]+sfa+eps);
 						WavCoeffs_b[dir][coeffloc_ab] *= (SQR(sfavebd[coeffloc_ab])+SQR(sfb))/(sfavebd[coeffloc_ab]+sfb+eps);
-
 					}//now chrominance coefficients are denoised
 				}
 #endif
@@ -2455,6 +2404,7 @@ SSEFUNCTION	void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
 					magv = SQRV(LVFU(WavCoeffs_L[dir][i]));
 					_mm_storeu_ps( &sfave[i], magv / (magv + mad_Lv*xexpf(-magv/(ninev * mad_Lv)) + epsv));
 				}
+				// few remaining pixels
 				for (; i<W_L*H_L; i++) {
 					float mag = SQR(WavCoeffs_L[dir][i]);
 					sfave[i] = mag/(mag+mad_LL[i]*xexpf(-mag/(9*mad_LL[i]))+eps);
@@ -2471,12 +2421,13 @@ SSEFUNCTION	void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
 
 #ifdef __SSE2__ 
 				__m128	sfv;
-				for (int i=0; i<W_L*H_L-3; i+=4) {
+				for (i=0; i<W_L*H_L-3; i+=4) {
 					sfv = LVFU(sfave[i]);
 					//use smoothed shrinkage unless local shrinkage is much less
 					_mm_storeu_ps( &WavCoeffs_L[dir][i], _mm_loadu_ps( &WavCoeffs_L[dir][i]) * (SQRV( LVFU(sfaved[i] )) + SQRV(sfv)) / (LVFU(sfaved[i])+sfv+epsv));
 				}
-				for (int i=(W_L*H_L)-((W_L*H_L)%4); i<W_L*H_L; i++) {
+				// few remaining pixels
+				for (; i<W_L*H_L; i++) {
 					float sf = sfave[i];
 
 					//use smoothed shrinkage unless local shrinkage is much less
@@ -2493,10 +2444,7 @@ SSEFUNCTION	void ImProcFunctions::ShrinkAll(float ** WavCoeffs_L, float ** WavCo
 				}//now luminance coefficients are denoised
 #endif
 			}
-
-
 		}
-
 	}
 
 
@@ -2721,7 +2669,7 @@ SSEFUNCTION void ImProcFunctions::RGB_denoise_info(Imagefloat * src, Imagefloat
 		//nothing to do
 		return;
 	}
-
+
 	int hei,wid;
 	float** lumcalc;
 	float** acalc;
diff --git a/rtengine/boxblur.h b/rtengine/boxblur.h
index 9871d43cc..dcd35fdfe 100644
--- a/rtengine/boxblur.h
+++ b/rtengine/boxblur.h
@@ -708,20 +708,25 @@ template<class T, class A> SSEFUNCTION void boxabsblur (T* src, A* dst, int radx
 		//horizontal blur
 		for (int row = 0; row < H; row++) {
 			int len = radx + 1;
-			temp[row*W+0] = fabsf((float)src[row*W+0]);
+			float tempval = fabsf((float)src[row*W+0]);
 			for (int j=1; j<=radx; j++) {
-				temp[row*W+0] += fabsf((float)src[row*W+j]);
+				tempval += fabsf((float)src[row*W+j]);
 			}
-			temp[row*W+0] = temp[row*W+0] / len;
+			tempval /= len;
+			temp[row*W+0] = tempval;
 			for (int col=1; col<=radx; col++) {
-				temp[row*W+col] = (temp[row*W+col-1]*len + fabsf(src[row*W+col+radx]))/(len+1);
+				tempval = (tempval*len + fabsf(src[row*W+col+radx]))/(len+1);
+				temp[row*W+col] = tempval;
 				len ++;
 			}
+			float rlen = 1.f/(float)len;
 			for (int col = radx+1; col < W-radx; col++) {
-				temp[row*W+col] = temp[row*W+col-1] + ((float)(fabsf(src[row*W+col+radx]) - fabsf(src[row*W+col-radx-1])))/len;
+				tempval = tempval + ((float)(fabsf(src[row*W+col+radx]) - fabsf(src[row*W+col-radx-1])))*rlen;
+				temp[row*W+col] = tempval;
 			}
 			for (int col=W-radx; col<W; col++) {
-				temp[row*W+col] = (temp[row*W+col-1]*len - fabsf(src[row*W+col-radx-1]))/(len-1);
+				tempval = (tempval*len - fabsf(src[row*W+col-radx-1]))/(len-1);
+				temp[row*W+col] = tempval;
 				len --;
 			}
 		}
@@ -737,29 +742,30 @@ template<class T, class A> SSEFUNCTION void boxabsblur (T* src, A* dst, int radx
 #ifdef __SSE2__
 		__m128	leninitv = _mm_set1_ps( (float)(rady+1));
 		__m128 	onev = _mm_set1_ps( 1.0f );
-		__m128	tempv,lenv,lenp1v,lenm1v;
+		__m128	tempv,lenv,lenp1v,lenm1v,rlenv;
 		for (int col = 0; col < W-3; col+=4) {
 			lenv = leninitv;
-			tempv = LVFU(temp[0*W+col]);
+			tempv = LVF(temp[0*W+col]);
 			for (int i=1; i<=rady; i++) {
-				tempv = tempv + LVFU(temp[i*W+col]);
+				tempv = tempv + LVF(temp[i*W+col]);
 			}
 			tempv = tempv / lenv;
-			_mm_storeu_ps( &dst[0*W+col], tempv );
+			_mm_store_ps( &dst[0*W+col], tempv );
 			for (int row=1; row<=rady; row++) {
 				lenp1v = lenv + onev;
-				tempv = (tempv*lenv + LVFU(temp[(row+rady)*W+col]))/lenp1v;
-				_mm_storeu_ps( &dst[row*W+col],tempv);
+				tempv = (tempv*lenv + LVF(temp[(row+rady)*W+col]))/lenp1v;
+				_mm_store_ps( &dst[row*W+col],tempv);
 				lenv = lenp1v;
 			}
+			rlenv = onev / lenv;
 			for (int row = rady+1; row < H-rady; row++) {
-				tempv = tempv + (LVFU(temp[(row+rady)*W+col])- LVFU(temp[(row-rady-1)*W+col]))/lenv;
-				_mm_storeu_ps( &dst[row*W+col], tempv);
+				tempv = tempv + (LVF(temp[(row+rady)*W+col])- LVF(temp[(row-rady-1)*W+col]))*rlenv;
+				_mm_store_ps( &dst[row*W+col], tempv);
 			}
 			for (int row=H-rady; row<H; row++) {
 				lenm1v = lenv - onev;
-				tempv = (tempv*lenv - LVFU(temp[(row-rady-1)*W+col]))/lenm1v;
-				_mm_storeu_ps( &dst[row*W+col], tempv);
+				tempv = (tempv*lenv - LVF(temp[(row-rady-1)*W+col]))/lenm1v;
+				_mm_store_ps( &dst[row*W+col], tempv);
 				lenv = lenm1v;
 			}
 	}
diff --git a/rtengine/cplx_wavelet_filter_coeffs.h b/rtengine/cplx_wavelet_filter_coeffs.h
index cfe3619fe..8f2bcbc49 100644
--- a/rtengine/cplx_wavelet_filter_coeffs.h
+++ b/rtengine/cplx_wavelet_filter_coeffs.h
@@ -22,7 +22,7 @@ namespace rtengine {
 
 const int Daub4_len=6;
 const int Daub4_offset=2;
-const float Daub4_anal[2][6] = {//analysis filter
+const float Daub4_anal[2][6] ALIGNED16 = {//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
index a5454a0ec..58603e804 100644
--- a/rtengine/cplx_wavelet_level.h
+++ b/rtengine/cplx_wavelet_level.h
@@ -25,7 +25,6 @@
 #include <cstddef>
 #include "rt_math.h"
 #include "opthelper.h"
-
 namespace rtengine {
 
 	template<typename T>
@@ -50,15 +49,15 @@ namespace rtengine {
 
 		// load a row/column of input data, possibly with padding
 
-		void AnalysisFilterHaarVertical (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen);
-		void AnalysisFilterHaarHorizontal (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen);
+		void AnalysisFilterHaarVertical (T * srcbuffer, T * dstLo, T * dstHi, int pitch, int srclen, int row);
+		void AnalysisFilterHaarHorizontal (T * srcbuffer, T * dstLo, T * dstHi, int srclen, int row);
 		void SynthesisFilterHaarHorizontal (T * srcLo, T * srcHi, T * dst, int dstlen);
 		void SynthesisFilterHaarVertical (T * srcLo, T * srcHi, T * dst, int pitch, int dstlen);
 
 		void AnalysisFilterSubsampHorizontal (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, 
-							 int taps, int offset, int pitch, int srclen, int m_w2);
+							 int taps, int offset, int pitch, int srclen, int m_w2, int row);
 		void AnalysisFilterSubsampVertical (T * srcbuffer, T * dstLo, T * dstHi, float *filterLo, float *filterHi, 
-							 int taps, int offset, int pitch, int srclen);
+							 int taps, int offset, int pitch, int srclen, int row);
 		void SynthesisFilterSubsampHorizontal (T * srcLo, T * srcHi, T * dst,
 							  float *filterLo, float *filterHi, int taps, int offset, int dstlen);
 		void SynthesisFilterSubsampVertical (T * srcLo, T * srcHi, T * dst, float *filterLo, float *filterHi, int taps, int offset, int pitch, int dstlen);
@@ -66,7 +65,7 @@ namespace rtengine {
 	public:
 
 		T ** wavcoeffs;
-				// full size
+		// full size
 		size_t m_w, m_h;
 
 		// size of low frequency part
@@ -152,38 +151,33 @@ namespace rtengine {
 	}
 
 	template<typename T>
-	void wavelet_level<T>::AnalysisFilterHaarHorizontal (T * RESTRICT srcbuffer, T * RESTRICT dstLo, T * RESTRICT dstHi, int pitch, int srclen) {
-
+	void wavelet_level<T>::AnalysisFilterHaarHorizontal (T * RESTRICT srcbuffer, T * RESTRICT dstLo, T * RESTRICT dstHi, int srclen, int row) {
 		/* Basic convolution code
 		 * Applies a Haar filter 
-		*/										
-		for(int j=0;j<pitch;j++) {
+		*/							
 			for(int i = 0; i < (srclen - skip); i++) {
-				dstLo[j*srclen+i] = 0.5f*(srcbuffer[j*srclen+i] + srcbuffer[j*srclen+i+skip]);
-				dstHi[j*srclen+i] = 0.5f*(srcbuffer[j*srclen+i] - srcbuffer[j*srclen+i+skip]);
+				dstLo[row*srclen+i] = (srcbuffer[i] + srcbuffer[i+skip]);
+				dstHi[row*srclen+i] = (srcbuffer[i] - srcbuffer[i+skip]);
 			}
 			for(size_t i = max(srclen-skip,skip); i < (srclen); i++) {
-				dstLo[j*srclen+i] = 0.5f*(srcbuffer[j*srclen+i] + srcbuffer[j*srclen+i-skip]);
-				dstHi[j*srclen+i] = 0.5f*(srcbuffer[j*srclen+i] - srcbuffer[j*srclen+i-skip]);
+				dstLo[row*srclen+i] = (srcbuffer[i] + srcbuffer[i-skip]);
+				dstHi[row*srclen+i] = (srcbuffer[i] - srcbuffer[i-skip]);
 			}
-		}
 	}
 
-	template<typename T> void wavelet_level<T>::AnalysisFilterHaarVertical (T * RESTRICT srcbuffer, T * RESTRICT dstLo, T * RESTRICT dstHi, int pitch, int srclen) {
-
+	template<typename T> void wavelet_level<T>::AnalysisFilterHaarVertical (T * RESTRICT srcbuffer, T * RESTRICT dstLo, T * RESTRICT dstHi, int pitch, int srclen, int row) {
 	/* Basic convolution code
 	 * Applies a Haar filter 
 	*/
-		for(int i = 0; i < (srclen - skip); i++) {
+		if(row < (srclen - skip)) {
 			for(int j=0;j<pitch;j++) {
-				dstLo[(pitch*(i))+j] = 0.5f*(srcbuffer[i*pitch+j] + srcbuffer[(i+skip)*pitch+j]);
-				dstHi[(pitch*(i))+j] = 0.5f*(srcbuffer[i*pitch+j] - srcbuffer[(i+skip)*pitch+j]);
+				dstLo[j] = 0.25f*(srcbuffer[row*pitch+j] + srcbuffer[(row+skip)*pitch+j]);
+				dstHi[j] = 0.25f*(srcbuffer[row*pitch+j] - srcbuffer[(row+skip)*pitch+j]);
 			}
-		}
-		for(size_t i = max(srclen-skip,skip); i < (srclen); i++) {
+		} else if(row>=max(srclen-skip,skip)) {
 			for(int j=0;j<pitch;j++) {
-				dstLo[(pitch*(i))+j] = 0.5f*(srcbuffer[i*pitch+j] + srcbuffer[(i-skip)*pitch+j]);
-				dstHi[(pitch*(i))+j] = 0.5f*(srcbuffer[i*pitch+j] - srcbuffer[(i-skip)*pitch+j]);
+				dstLo[j] = 0.25f*(srcbuffer[row*pitch+j] + srcbuffer[(row-skip)*pitch+j]);
+				dstHi[j] = 0.25f*(srcbuffer[row*pitch+j] - srcbuffer[(row-skip)*pitch+j]);
 			}
 		}
 	}
@@ -227,40 +221,37 @@ namespace rtengine {
 
 	template<typename T>
 	void wavelet_level<T>::AnalysisFilterSubsampHorizontal (T * RESTRICT srcbuffer, T * RESTRICT dstLo, T * RESTRICT dstHi, float * RESTRICT filterLo, float *filterHi, 
-												int taps, int offset, int pitch, int srclen, int m_w2) {
-
+												int taps, int offset, int pitch, int srclen, int m_w2, int row) {
 		/* 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 k=0;k<pitch;k++) {
+		// calculate coefficients
 			for(int i = 0; i < srclen; i+=2) {
 				float lo = 0.f, hi = 0.f;
 				if (LIKELY(i>skip*taps && i<srclen-skip*taps)) {//bulk
 					for (int j=0, l=-skip*offset; j<taps; j++, l+=skip) {
-						float src = srcbuffer[k*srclen+i-l];
+						float src = srcbuffer[i-l];
 						lo += filterLo[j] * src;//lopass channel
 						hi += filterHi[j] * src;//hipass channel
 					}
 				} else {
 					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[k*srclen+arg];//lopass channel
-						hi += filterHi[j] * srcbuffer[k*srclen+arg];//hipass channel
+						lo += filterLo[j] * srcbuffer[arg];//lopass channel
+						hi += filterHi[j] * srcbuffer[arg];//hipass channel
 					}
 				}
-				dstLo[k*m_w2+((i/2))] = lo;
-				dstHi[k*m_w2+((i/2))] = hi;
+				dstLo[row*m_w2+((i/2))] = lo;
+				dstHi[row*m_w2+((i/2))] = hi;
 			}
-		}
 	}
 
 	template<typename T> void wavelet_level<T>::AnalysisFilterSubsampVertical (T * RESTRICT srcbuffer, T * RESTRICT dstLo, T * RESTRICT dstHi, float * RESTRICT filterLo, float * RESTRICT filterHi, 
-													int taps, int offset, int pitch, int srclen) {
+													int taps, int offset, int pitch, int srclen, int row) {
+
 			/* Basic convolution code
 			 * Applies an FIR filter 'filter' with filter length 'taps', 
 			 * aligning the 'offset' element of the filter with
@@ -269,31 +260,28 @@ namespace rtengine {
 			 */
 
 			// calculate coefficients
-
-		for(int i = 0; i < srclen; i+=2) {
-			if (LIKELY(i>skip*taps && i<srclen-skip*taps)) {//bulk
+			if (LIKELY(row>skip*taps && row<srclen-skip*taps)) {//bulk
 				for (int k=0; k<pitch; k++) {
 					float lo = 0.f, hi = 0.f;
 					for (int j=0, l=-skip*offset; j<taps; j++, l+=skip) {
-						lo += filterLo[j] * srcbuffer[(i-l)*pitch+k];//lopass channel
-						hi += filterHi[j] * srcbuffer[(i-l)*pitch+k];//hipass channel
+						lo += filterLo[j] * srcbuffer[(row-l)*pitch+k];//lopass channel
+						hi += filterHi[j] * srcbuffer[(row-l)*pitch+k];//hipass channel
 					}
-					dstLo[(pitch*(i/2))+k] = lo;
-					dstHi[(pitch*(i/2))+k] = hi;
+					dstLo[k] = lo;
+					dstHi[k] = hi;
 				}
 			} else {//boundary
 				for (int k=0; k<pitch; k++) {
 					float lo = 0.f, hi = 0.f;
 					for (int j=0; j<taps; j++) {
-						int arg = max(0,min(i+skip*(offset-j),srclen-1))*pitch+k;//clamped BC's
+						int arg = max(0,min(row+skip*(offset-j),srclen-1))*pitch+k;//clamped BC's
 						lo += filterLo[j] * srcbuffer[arg];//lopass channel
 						hi += filterHi[j] * srcbuffer[arg];//hipass channel
 					}
-					dstLo[(pitch*(i/2))+k] = lo;
-					dstHi[(pitch*(i/2))+k] = hi;
+					dstLo[k] = lo;
+					dstHi[k] = hi;
 				}
 			}
-		}
 	}
 
 	template<typename T> void wavelet_level<T>::SynthesisFilterSubsampHorizontal (T * RESTRICT srcLo, T * RESTRICT srcHi, T * RESTRICT dst, float * RESTRICT filterLo, float * RESTRICT filterHi, int taps, int offset, int dstlen) {
@@ -325,7 +313,7 @@ namespace rtengine {
 						tot += ((filterLo[j] * srcLo[k*m_w2+arg] + filterHi[j] * srcHi[k*m_w2+arg]));
 					}
 				}
-				dst[k*m_w+(i-m_pad)] = 2.f * tot;
+				dst[k*m_w+(i-m_pad)] = tot;
 			}
 		}
 	}
@@ -353,7 +341,7 @@ namespace rtengine {
 					for (int j=begin, l=0; j<taps; j+=2, l+=skip) {
 						tot += ((filterLo[j] * srcLo[(i_src-l)*pitch+k] + filterHi[j] * srcHi[(i_src-l)*pitch+k]));
 					}
-					dst[pitch*(i-m_pad)+k] = 2.f * tot;
+					dst[pitch*(i-m_pad)+k] = 4.f * tot;
 				}
 			} else {//boundary
 				for (int k=0; k<pitch; k++) {
@@ -362,7 +350,7 @@ namespace rtengine {
 						int arg = max(0,min((i_src-l),srclen-1))*pitch+k;//clamped BC's
 						tot += ((filterLo[j] * srcLo[arg] + filterHi[j] * srcHi[arg]));
 					}
-					dst[pitch*(i-m_pad)+k] = 2.f * tot;
+					dst[pitch*(i-m_pad)+k] = 4.f * tot;
 				}
 			}
 		}
@@ -370,22 +358,22 @@ namespace rtengine {
 
 	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 tmpLo[m_w] ALIGNED64;
+		T tmpHi[m_w] ALIGNED64;
 		/* filter along rows and columns */
 		if(subsamp_out) {
-			AnalysisFilterSubsampVertical (src, tmpLo, tmpHi, filterV, filterV+taps, taps, offset, m_w/*output_pitch*/, m_h/*srclen*/);
-			AnalysisFilterSubsampHorizontal (tmpLo, wavcoeffs[0], wavcoeffs[1], filterH, filterH+taps, taps, offset, m_h2/*output_pitch*/, m_w/*srclen*/, m_w2);
-			AnalysisFilterSubsampHorizontal (tmpHi, wavcoeffs[2], wavcoeffs[3], filterH, filterH+taps, taps, offset, m_h2/*output_pitch*/, m_w/*srclen*/, m_w2);
+			for(int row=0;row<m_h;row+=2) {
+				AnalysisFilterSubsampVertical (src, tmpLo, tmpHi, filterV, filterV+taps, taps, offset, m_w/*output_pitch*/, m_h/*srclen*/, row);
+				AnalysisFilterSubsampHorizontal (tmpLo, wavcoeffs[0], wavcoeffs[1], filterH, filterH+taps, taps, offset, m_h2/*output_pitch*/, m_w/*srclen*/, m_w2, row/2);
+				AnalysisFilterSubsampHorizontal (tmpHi, wavcoeffs[2], wavcoeffs[3], filterH, filterH+taps, taps, offset, m_h2/*output_pitch*/, m_w/*srclen*/, m_w2, row/2);
+			}
 		} else {
-			AnalysisFilterHaarVertical (src, tmpLo, tmpHi, m_w, m_h);
-			AnalysisFilterHaarHorizontal (tmpLo, wavcoeffs[0], wavcoeffs[1], m_h, m_w);
-			AnalysisFilterHaarHorizontal (tmpHi, wavcoeffs[2], wavcoeffs[3], m_h, m_w);
+			for(int row=0;row<m_h;row++) {
+				AnalysisFilterHaarVertical (src, tmpLo, tmpHi, m_w, m_h, row);
+				AnalysisFilterHaarHorizontal (tmpLo, wavcoeffs[0], wavcoeffs[1], m_w, row);
+				AnalysisFilterHaarHorizontal (tmpHi, wavcoeffs[2], wavcoeffs[3], m_w, row);
+			}
 		}
-
-		delete[] tmpLo;
-		delete[] tmpHi;
 	}
 
 	template<typename T> template<typename E> void wavelet_level<T>::reconstruct_level(E* tmpLo, E* tmpHi, E *dst, float *filterV, float *filterH, int taps, int offset) {