Accuracy 2D wavelet transform and reconstruction

[Bouchan__Thomas]

Hi, 

I have been working on the implementation of the wavelet transform in 2D and I am quite surprised to see that the reconstruction does not give the same result as the original image when I do a decomposition followed by a reconstruction, using Haar or DB4 wavelet. In Matlab the reconstruction is almost perfect so I was expecting the same here. 

The main issue I have is that my reconstructed image looks shifted by 1 pixel in both dimensions. I have checked my pointers and they look correct.

I wondered if it was due to my implementation or if it is an expected behavior. 

I can provide my code if needed. 

Thanks, 

Thomas

[Gennady_F_Intel]

Do you mean then you use Matlab, the output result are bit-to-bit the same as input one?

[Bouchan__Thomas]

Thanks for your help,

When I use Matlab, I have a very similar image after the reconstruction (in terms of alignment). When using IPP library, I don't manage to get the same results. I guess that is it due to my understanding of the anchor parameter. To me, the result should be the same, no matter its value, but with my implementation, it's not true. See examples in attached document.

I also provided my code as it is probably easier to spot the error. There are two functions (1 for the forward wavelet transform, one for the reconstruction). I used to proof check my, the files written in the debug sections (at the end of each functions)

I would really appreciate some advice on that.

Regards

Thomas

void ComputeWaveletTransform(Ipp32f* src, IppiSize srcRoi, int srcStep, Ipp32f* dst, IppiSize* dstRoi, int* dstStep) 
{
	/*
	INPUT src, srcRoi, dst*
	OUTPUT dst, dstRoi, dstStep 
	Are provided only pointer to the outputs
	*/
	IppStatus status;
	
	// Initialisation 
	int lenLow = 8, lenHigh = 8, anchorLow = 0, anchorHigh = 0;
	// coefficients for DB4 forward
	Ipp32f pTapsLow[8] = { -0.010597401784997278, 0.032883011666982945, 0.030841381835986965, -0.18703481171888114, -0.02798376941698385, 0.6308807679295904, 0.7148465705525415, 0.23037781330885523 };
	Ipp32f pTapsHigh[8] = { -0.23037781330885523, 0.7148465705525415, -0.6308807679295904, -0.02798376941698385, 0.18703481171888114, 0.030841381835986965, -0.032883011666982945, -0.010597401784997278 };

	//int lenLow = 2, lenHigh = 2, anchorLow = 0, anchorHigh = 0;
	//// coefficients for Haar forward
	//Ipp32f pTapsLow[2] = { 0.7071067811865476, 0.7071067811865476 };
	//Ipp32f pTapsHigh[2] = { -0.7071067811865476, 0.7071067811865476 };


	int specSize, bufSize;
	
	status = ippiWTFwdGetSize_32f(1, lenLow, anchorLow, lenHigh, anchorHigh, &specSize, &bufSize);
	if (status != 0) { printf("ippiWTFwdGetSize_32f result : %s\n", ippGetStatusString(status)); }

	IppiWTFwdSpec_32f_C1R* pSpec = (IppiWTFwdSpec_32f_C1R*)ippMalloc(specSize);
	Ipp8u* buffer = (Ipp8u*)ippMalloc(bufSize);

	status = ippiWTFwdInit_32f_C1R(pSpec, pTapsLow, lenLow, anchorLow, pTapsHigh, lenHigh, anchorHigh);
	if (status != 0) { printf("ippiWTFwdInit_32f_C1R result : %s\n", ippGetStatusString(status)); }


	// Image preprocessing (add borders etc)
	int leftBorderLow = lenLow - 1 - anchorLow;
	int leftBorderHigh = lenHigh - 1 - anchorHigh;
	int rightBorderLow = anchorLow;
	int rightBorderHigh = anchorHigh;
	int leftTopBorder = IPP_MAX(leftBorderLow, leftBorderHigh);
	int rightBottomBorder = IPP_MAX(rightBorderLow, rightBorderHigh);

	IppiSize roiBorderSize = { srcRoi.width + leftTopBorder + rightBottomBorder, srcRoi.height + leftTopBorder + rightBottomBorder };
	int sourceBorderStep = roiBorderSize.width * sizeof(Ipp32f);
	size_t sizeSourceBorder = roiBorderSize.height * roiBorderSize.width * sizeof(Ipp32f);
	Ipp32f* bSourceImage = (Ipp32f*)ippMalloc(sizeSourceBorder);

	//printf("\nPixel addition: \nTop/Left : %d \nBottom/Right : %d\n\n", leftTopBorder, rightBottomBorder);
	status = ippiCopyMirrorBorder_32f_C1R(src, srcStep, srcRoi, bSourceImage, sourceBorderStep, roiBorderSize, leftTopBorder, leftTopBorder);
	if (status != 0) { printf("ippiCopyMirrorBorder_32f_C1R result : %s\n", ippGetStatusString(status)); }

	IppiSize detailsSize = { srcRoi.width / 2, srcRoi.height / 2 };
	size_t detailsBuf = detailsSize.height * detailsSize.width * sizeof(Ipp32f);
	int detailNPix = detailsSize.height * detailsSize.width;
	int detailStep = detailsSize.width * sizeof(Ipp32f);

	Ipp32f* approx = (Ipp32f*)ippMalloc(detailsBuf);
	Ipp32f* vertical = (Ipp32f*)ippMalloc(detailsBuf);
	Ipp32f* horizontal = (Ipp32f*)ippMalloc(detailsBuf);
	Ipp32f* diagonal = (Ipp32f*)ippMalloc(detailsBuf);

	Ipp32f* pSource = bSourceImage + (leftTopBorder) * (roiBorderSize.width + 1);

	status = ippiWTFwd_32f_C1R(pSource, sourceBorderStep, approx, detailStep, horizontal, detailStep, vertical, detailStep, diagonal, detailStep, detailsSize, pSpec, buffer);
	if (status != 0) { printf("ippiWTFwd_32f_C1R result : %s\n", ippGetStatusString(status)); }

	// Concatenating 4 components
	
	dstRoi->width = srcRoi.width;
	dstRoi->height = srcRoi.height;
	size_t dstNElem = dstRoi->height * dstRoi->width;
	*dstStep = dstRoi->width * sizeof(Ipp32f);
	
	status = ippiCopy_32f_C1R(approx, detailStep, dst, *dstStep, detailsSize);
	status = ippiCopy_32f_C1R(horizontal, detailStep, dst + detailsSize.width, *dstStep, detailsSize);
	status = ippiCopy_32f_C1R(vertical, detailStep, dst + dstRoi->width * detailsSize.height, *dstStep, detailsSize);
	status = ippiCopy_32f_C1R(diagonal, detailStep, dst + dstRoi->width * detailsSize.height + detailsSize.width, *dstStep, detailsSize);
	
	// DEBUG 
	FILE *fileToWrite;
	fileToWrite = fopen("outputLenaDebugWT.bin", "wb");
	if (fileToWrite == NULL)
	{
		printf("File has not been opened correctly to write \n");
	}
	else {
		fwrite(src, 4, 512*512, fileToWrite);
		fwrite(bSourceImage, 4, roiBorderSize.height * roiBorderSize.width, fileToWrite);

		fwrite(approx, 4, detailNPix, fileToWrite);
		fwrite(horizontal, 4, detailNPix, fileToWrite);
		fwrite(vertical, 4, detailNPix, fileToWrite);
		fwrite(diagonal, 4, detailNPix, fileToWrite);
		fwrite(dst, 4, dstRoi->height * dstRoi->width, fileToWrite);
		fclose(fileToWrite);
	}

	ippFree(pSpec); ippFree(buffer); ippFree(bSourceImage); ippFree(approx); ippFree(diagonal); ippFree(horizontal); ippFree(vertical);
}


void ComputeInverseWT(const Ipp32f* wt, int wtStep, IppiSize wtRoi, Ipp32f* dst, int dstStep)
{
	// ATTENTION :
	// dst has to be initialized with ippMalloc before calling this function
	IppStatus s1, s2, s3, s4;

	int lenLow = 8, lenHigh = 8, anchorLow = 0, anchorHigh = 0;
	Ipp32f pTapsLow[8] = { 0.23037781330885523, 0.7148465705525415, 0.6308807679295904, -0.02798376941698385, -0.18703481171888114, 0.030841381835986965, 0.032883011666982945, -0.010597401784997278 };
	Ipp32f pTapsHigh[8] = { -0.010597401784997278, -0.032883011666982945, 0.030841381835986965, 0.18703481171888114, -0.02798376941698385, -0.6308807679295904, 0.7148465705525415, -0.23037781330885523 };

	//int lenLow = 2, lenHigh = 2, anchorLow = 0, anchorHigh = 0;
	//Ipp32f pTapsLow[2] = { 0.7071067811865476, 0.7071067811865476 };
	//Ipp32f pTapsHigh[2] = { 0.7071067811865476, -0.7071067811865476 };


	int pSpecSize, pBufSize;

	IppiSize detailRoi = { wtRoi.width / 2, wtRoi.height / 2 };

	int leftBorderLow = (lenLow - 1 - anchorLow) / 2;
	int leftBorderHigh = (lenHigh - 1 - anchorHigh) / 2;
	int rightBorderLow = (anchorLow + 1) / 2;
	int rightBorderHigh = (anchorHigh + 1) / 2;
	int apprLeftBorder = leftBorderLow;
	int apprRightBorder = rightBorderLow;
	int apprTopBorder = leftBorderLow;
	int apprBottomBorder = rightBorderLow;
	int detxLeftBorder = leftBorderLow;
	int detxRightBorder = rightBorderLow;
	int detxTopBorder = leftBorderHigh;
	int detxBottomBorder = rightBorderHigh;
	int detyLeftBorder = leftBorderHigh;
	int detyRightBorder = rightBorderHigh;
	int detyTopBorder = leftBorderLow;
	int detyBottomBorder = rightBorderLow;
	int detxyLeftBorder = leftBorderHigh;
	int detxyRightBorder = rightBorderHigh;
	int detxyTopBorder = leftBorderHigh;
	int detxyBottomBorder = rightBorderHigh;

	IppiSize apprSizeB = { detailRoi.width + apprLeftBorder + apprRightBorder, detailRoi.height + apprTopBorder + apprBottomBorder };
	IppiSize hzSizeB = { detailRoi.width + detxLeftBorder + detxRightBorder, detailRoi.height + detxTopBorder + detxBottomBorder };
	IppiSize vertSizeB = { detailRoi.width + detyLeftBorder + detyRightBorder , detailRoi.height + detyTopBorder + detyBottomBorder };
	IppiSize diagSizeB = { detailRoi.width + detxyLeftBorder + detxyRightBorder, detailRoi.height + detxyTopBorder + detxyBottomBorder };
	
	int approxStep = apprSizeB.width * sizeof(Ipp32f);
	int hzStep = hzSizeB.width * sizeof(Ipp32f);
	int vertStep = vertSizeB.width * sizeof(Ipp32f);
	int diagStep = diagSizeB.width * sizeof(Ipp32f);

	int detailStep = detailRoi.width * sizeof(Ipp32f);

	Ipp32f* approxB = ippMalloc(apprSizeB.height * apprSizeB.width * sizeof(Ipp32f));
	Ipp32f* verticalB = ippMalloc(hzSizeB.height * hzSizeB.width * sizeof(Ipp32f));
	Ipp32f* horizontalB = ippMalloc(vertSizeB.height * vertSizeB.width * sizeof(Ipp32f));
	Ipp32f* diagonalB = ippMalloc(diagSizeB.height * diagSizeB.width * sizeof(Ipp32f));


	// Define pointers to the first element of the image (without border)
	Ipp32f* pApprox = approxB + apprSizeB.width * apprTopBorder + apprLeftBorder;
	Ipp32f* pHztal = horizontalB + hzSizeB.width * detxTopBorder + detxLeftBorder;
	Ipp32f* pVert = verticalB + vertSizeB.width * detyTopBorder + detyLeftBorder;
	Ipp32f* pDiag = diagonalB + diagSizeB.width * detxyTopBorder + detxyLeftBorder;

	s1 = ippiCopyMirrorBorder_32f_C1R(wt, wtStep, detailRoi, approxB, approxStep, apprSizeB, apprTopBorder, apprLeftBorder);
	s2 = ippiCopyMirrorBorder_32f_C1R(wt + detailRoi.width, wtStep, detailRoi, horizontalB, hzStep, hzSizeB, detxTopBorder, detxLeftBorder);
	s3 = ippiCopyMirrorBorder_32f_C1R(wt + (wtStep / sizeof(Ipp32f))*detailRoi.height, wtStep, detailRoi, verticalB, vertStep, vertSizeB, detyTopBorder, detyLeftBorder);
	s4 = ippiCopyMirrorBorder_32f_C1R(wt + (wtStep / sizeof(Ipp32f))*detailRoi.height + detailRoi.width, wtStep, detailRoi, diagonalB, diagStep, diagSizeB, detxyTopBorder, detxyLeftBorder);
	if (s1 + s2 + s3 + s4 != 0) // Detect error and display if needed
	{
		printf("Extending borders \nApprox : %s\nHorizontal : %s\nVertical : %s\nDiagonal : %s\n", ippGetStatusString(s1), ippGetStatusString(s2), ippGetStatusString(s3), ippGetStatusString(s4));
	}

		
	// Initialisation of the Inv WT
	ippiWTInvGetSize_32f(1, lenLow, anchorLow, lenHigh, anchorHigh, &pSpecSize, &pBufSize);
	Ipp8u* pBuffer = (Ipp8u*)ippMalloc(pBufSize);
	IppiWTInvSpec_32f_C1R* pSpec = (IppiWTInvSpec_32f_C1R*)ippMalloc(pSpecSize);
	
	s1 = ippiWTInvInit_32f_C1R(pSpec, pTapsLow, lenLow, anchorLow, pTapsHigh, lenHigh, anchorHigh);

	if (s1 != 0)
	{
		printf("ippiWTInvInit_32f_C1R result : %s\n", ippGetStatusString(s1));
	}
	
	s1 = ippiWTInv_32f_C1R(pApprox, approxStep, pHztal, hzStep, pVert, vertStep, pDiag, diagStep, detailRoi, dst, dstStep, pSpec, pBuffer);
	if (s1 != 0)
	{
		printf("ippiWTInv_32f_C1R result : %s\n", ippGetStatusString(s1));
	}

	//// DEBUG 
	FILE *fileToWrite;
	fileToWrite = fopen("outputLenaDebug.bin", "wb");
	if (fileToWrite == NULL)
	{
		printf("File has not been opened correctly to write \n");
	}
	else {
		//fwrite(wt, 4, wtRoi.height * wtRoi.width, fileToWrite);
		fwrite(wt, 4, 512*512, fileToWrite);

		fwrite(approxB, 4, apprSizeB.height * apprSizeB.width, fileToWrite);
		fwrite(horizontalB, 4, hzSizeB.height * hzSizeB.width, fileToWrite);
		fwrite(verticalB, 4, vertSizeB.height * vertSizeB.width, fileToWrite);
		fwrite(diagonalB, 4, diagSizeB.height * diagSizeB.width, fileToWrite);
		fwrite(dst, 4, 256*256, fileToWrite);
		fclose(fileToWrite);
	}

	//// END DEBUG
	ippFree(approxB); 
	ippFree(verticalB); 
	ippFree(horizontalB); 
	ippFree(diagonalB); 

	ippFree(pBuffer); 
	ippFree(pSpec);
}

[Bouchan__Thomas]

Hi, 

Does anyone has an idea of my mistake?

Regards

Thomas