2D wavelet transform with symmetric border extension

yoavnaveh — Tue, 09 Nov 2010 15:52:57 GMT

Hello,

I'm using 2D wavelet transforms for image compression. for some reason, IPP does not have symmetric border extension, but only wrap-around border extension. this decrease the effectivness of the compression and creates a lot of non-zero coefficients. I've implemented a symmetric extension but it doesn't work - I can't reconstruct my image (forward and then inverse wavelet transforms differs from the original image).
below I've added both my symmetric extension function as well as a test code. in the test code, the wrap-around border extension calls are commented (if you un-comment them, you get the reconstruction).

please help me, as this is crucial to my project.

thanks,

Yoav

----------------------------------------------------------------------------------------------------------------------------
symmetric extension function:
----------------------------------------------------------------------------------------------------------------------------

bool CopySymmBorder_32f_C1R(const float* pSrc, int srcStep, IppiSize srcRoiSize, float* pDst, int dstStep, IppiSize dstRoiSize,
int topBorderHeight, int leftBorderWidth) {
// parameters consistency checks:
if (!pSrc || !pDst) return false;
int sizeofFloat = sizeof(float);
if (srcStep <= 0 || dstStep <= 0 || srcStep % sizeofFloat != 0 || dstStep % sizeofFloat != 0)
return false;
if (srcRoiSize.width <= 0 || srcRoiSize.height <= 0 || dstRoiSize.width <= 0 || dstRoiSize.height <= 0)
return false;
if (topBorderHeight < 0 || leftBorderWidth < 0)
return false;
if (dstRoiSize.width < srcRoiSize.width + leftBorderWidth)
return false;
if (dstRoiSize.height < srcRoiSize.height + topBorderHeight)
return false;

int rightBorderWidth = dstRoiSize.width - srcRoiSize.width - leftBorderWidth;
int botBorderHeight = dstRoiSize.height - srcRoiSize.height - topBorderHeight;

// copy into the appropriate position in :
float* pCurrDst = pDst + topBorderHeight*dstStep/sizeofFloat + leftBorderWidth;
if (ippStsNoErr != ippiCopy_32f_C1R (pSrc, srcStep, pCurrDst, dstStep, srcRoiSize))
return false;

bool needFurtherExtension = false;
const float* ps = NULL;
float* pd = NULL;
int N, M;
int i, j, ind;
int currTop, currBot, currLeft, currRight;

// extending the top border:
if (topBorderHeight > srcRoiSize.height) {
needFurtherExtension = true;
currTop = topBorderHeight > srcRoiSize.height;
}
else
currTop = 0;
N = std::min(topBorderHeight, srcRoiSize.height);

for (i = 0; i < N; ++i) {
pd = pDst + (N-i-1)*dstStep/sizeofFloat + leftBorderWidth;
ps = pSrc + i*srcStep/sizeofFloat;
memcpy (pd, ps, srcRoiSize.width*sizeofFloat);
}

// extending the bottom border:
if (botBorderHeight > srcRoiSize.height) {
needFurtherExtension = true;
currBot = dstRoiSize.height - (botBorderHeight > srcRoiSize.height);
}
else
currBot = dstRoiSize.height;
N = std::min(botBorderHeight, srcRoiSize.height);

for (i = 0; i < N; ++i) {
pd = pDst + (i+topBorderHeight+srcRoiSize.height)*dstStep/sizeofFloat + leftBorderWidth;
ps = pSrc + (srcRoiSize.height-i-1)*srcStep/sizeofFloat;
memcpy (pd, ps, srcRoiSize.width*sizeofFloat);
}

// extending left border:
if (leftBorderWidth > srcRoiSize.width) {
needFurtherExtension = true;
currLeft = leftBorderWidth - srcRoiSize.width;
}
else
currLeft = 0;
M = std::min(leftBorderWidth, srcRoiSize.width);
N = srcRoiSize.height;

if (M > 0) {
ind = topBorderHeight*dstStep/sizeofFloat + leftBorderWidth;
for (i = 0; i < N; ++i) {
ps = pSrc + i*srcStep/sizeofFloat;
for (j = 0; j < M; ++j) {
pDst[ind-j-1] = ps;
}
ind += dstStep/sizeofFloat;
}

// extending the (top-left, bottom-left) diagonal regions: flipping the left border
N = std::min(topBorderHeight, srcRoiSize.height);
for (i = 0; i < N; ++i) {
pd = pDst + (topBorderHeight-i-1)*dstStep/sizeofFloat;
ps = pDst + (topBorderHeight+i)*dstStep/sizeofFloat;
memcpy (pd, ps, M*sizeofFloat);
}

N = std::min(botBorderHeight, srcRoiSize.height);
for (i = 0; i < N; ++i) {
pd = pDst + (topBorderHeight+srcRoiSize.height+i)*dstStep/sizeofFloat;
ps = pDst + (topBorderHeight+srcRoiSize.height-i-1)*dstStep/sizeofFloat;
memcpy (pd, ps, M*sizeofFloat);
}
}

// extending right border:
if (rightBorderWidth > srcRoiSize.width) {
needFurtherExtension = true;
currRight = dstRoiSize.width - (rightBorderWidth - srcRoiSize.width);
}
else
currRight = dstRoiSize.width;
M = std::min(rightBorderWidth, srcRoiSize.width);
N = srcRoiSize.height;

if (M > 0) {
ind = topBorderHeight*dstStep/sizeofFloat + leftBorderWidth + srcRoiSize.width;
for (i = 0; i < N; ++i) {
ps = pSrc + i*srcStep/sizeofFloat;
for (j = 0; j < M; ++j) {
pDst[ind+j] = ps[srcRoiSize.width-j-1];
}
ind += dstStep/sizeofFloat;
}

// extending the (top-right, bottom-right) diagonal regions: flipping the right border
N = std::min(topBorderHeight, srcRoiSize.height);
for (i = 0; i < N; ++i) {
pd = pDst + (topBorderHeight-i-1)*dstStep/sizeofFloat + leftBorderWidth + srcRoiSize.width;
ps = pDst + (topBorderHeight+i)*dstStep/sizeofFloat + leftBorderWidth + srcRoiSize.width;
memcpy (pd, ps, M*sizeofFloat);
}

N = std::min(botBorderHeight, srcRoiSize.height);
for (i = 0; i < N; ++i) {
pd = pDst + (topBorderHeight+srcRoiSize.height+i)*dstStep/sizeofFloat + leftBorderWidth + srcRoiSize.width;
ps = pDst + (topBorderHeight+srcRoiSize.height-i-1)*dstStep/sizeofFloat + leftBorderWidth + srcRoiSize.width;
memcpy (pd, ps, M*sizeofFloat);
}
}

if (needFurtherExtension) {
//IppiSize middleRoi = {currRight - currLeft, currBot - currTop};
//if (ippStsNoErr != ippiCopyReplicateBorder_32f_C1IR (pDst, dstStep, middleRoi, dstRoiSize, currTop, currLeft))
// return false;
return false;
}

return true;
}

----------------------------------------------------------------------------------------------------------------------------
test code:
----------------------------------------------------------------------------------------------------------------------------

void main() {
IppiWTFwdSpec_32f_C1R* pSpec;
IppiWTInvSpec_32f_C1R* pSpecInv;

// Biorthogonal 4.4 (2) decomposition filters:
Ipp32f pTapsLowF[9] = { 0.02674875741080976, -0.01686411844287495, -0.07822326652898785, 0.2668641184428723, 0.6029490182363579, 0.2668641184428723, -0.07822326652898785, -0.01686411844287495, 0.02674875741080976};
int lenLowF = 9;
int anchorLowF = 4;
Ipp32f pTapsHighF[7] = {0.09127176311424948, -0.05754352622849957, -0.5912717631142470, 1.115087052456994, -0.5912717631142470, -0.05754352622849957, 0.09127176311424948};
int lenHighF = 7;
int anchorHighF = 4;

Ipp32f pSrc[8*8] = { 0.0, 0.1, 0.0, 11.0, 11.9, 0.0, 0.0, 0.0,
0.0, 0.0, 0.222, 11.0, 11.9, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 11.0, 11.0, 0.0, 0.0, 0.0,
11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0,
11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0, 11.0,
0.0, 0.0, 7.8, 11.0, 11.0, 0.0, 0.0, 0.0,
0.0, 2.5, 0.0, 11.0, 11.0, 0.0, 0.0, 0.0,
7.8, 0.0, 0.0, 11.0, 11.0, 0.0, 0.0, 0.0};

// borders manipulations:
int leftBorderLow = lenLowF - 1 - anchorLowF; // 4
int leftBorderHigh = lenHighF - 1 - anchorHighF; // 2
int rightBorderLow = lenLowF - 2 - leftBorderLow; // 3
int rightBorderHigh = lenHighF - 2 - leftBorderHigh; // 3
int leftTopBorder = IPP_MAX(leftBorderLow, leftBorderHigh);
int rightBottomBorder = IPP_MAX(rightBorderLow, rightBorderHigh);

// allocation the buffer with borders:

int srcWidthWithBorders = 8 + leftTopBorder + rightBottomBorder;
int srcHeightWithBorders = 8 + leftTopBorder + rightBottomBorder;
int srcStepWBorder;
Ipp32f* pSrcB = ippiMalloc_32f_C1 (srcWidthWithBorders, srcHeightWithBorders, &srcStepWBorder); //ippsMalloc_32f(bufSize);

//Ipp32f pSrcB[9*9];
int srcStep = 8*sizeof(Ipp32f);
IppiSize roiSize = {8, 8};
//int srcStepB = srcWidthWithBorders*sizeof(Ipp32f);
IppiSize roiSizeB = {srcWidthWithBorders, srcHeightWithBorders};

Ipp32f pDetailXDst[4*4];
Ipp32f pDetailYDst[4*4];
Ipp32f pDetailXYDst[4*4];
Ipp32f pApproxDst[4*4];
IppiSize dstRoiSize = {4, 4};
Ipp8u* pBuffer;

int approxStep, detailXStep, detailYStep, detailXYStep;
approxStep = detailXStep = detailYStep = detailXYStep = 4*sizeof(Ipp32f);

// adds border to the source image
//ippiCopyWrapBorder_32f_C1R(pSrc, srcStep, roiSize, pSrcB, srcStepWBorder, roiSizeB, leftTopBorder, leftTopBorder);
CopySymmBorder_32f_C1R (pSrc, srcStep, roiSize, pSrcB, srcStepWBorder, roiSizeB, leftTopBorder, leftTopBorder);
//----- performs forward wavelet transform
// (1) allocate the transform context & memory
ippiWTFwdInitAlloc_32f_C1R ( &pSpec, pTapsLowF, lenLowF, anchorLowF, pTapsHighF, lenHighF, anchorHighF);

// (2) allocate the auxiliary buffer
int bufSize;
ippiWTFwdGetBufSize_C1R(pSpec, &bufSize);
pBuffer = ippsMalloc_8u(bufSize);

// (3) perform the forward transform
int srcROIOffset = leftTopBorder * srcStepWBorder/4 + leftTopBorder;
if (ippiWTFwd_32f_C1R (pSrcB + srcROIOffset, srcStepWBorder, pApproxDst, approxStep, pDetailXDst,
detailXStep, pDetailYDst, detailYStep, pDetailXYDst, detailXYStep,
dstRoiSize, pSpec, pBuffer) != ippStsNoErr) {
printf("something's wrong in ippiWTFwd_32f_C1R");
}

//------------------------------------

Ipp8u* pBufferInv;
Ipp32f pDstInv[8*8];
IppiSize roiInvSize = {4, 4};
int stepDstInv = 8*sizeof(Ipp32f);
int bufSizeInv;

// Biorthogonal 4.4 (2) reconstruction filters:
Ipp32f pTapsLowI[7] = {-0.09127176311424948, -0.05754352622849957, 0.5912717631142470, 1.115087052456994, 0.5912717631142470, -0.05754352622849957, -0.09127176311424948};
int lenLowI = 7;
int anchorLowI = 3;
Ipp32f pTapsHighI[9] = {0.02674875741080976, 0.01686411844287495, -0.07822326652898785, -0.2668641184428723, 0.6029490182363579, -0.2668641184428723, -0.07822326652898785, 0.01686411844287495, 0.02674875741080976};
int lenHighI = 9;
int anchorHighI = 3;

//----- performs inverse wavelet transform
// (1) allocate the transform context & memory
ippiWTInvInitAlloc_32f_C1R (&pSpecInv, pTapsLowI, lenLowI, anchorLowI, pTapsHighI, lenHighI, anchorHighI);

// (2) allocate the auxiliary buffer
ippiWTInvGetBufSize_C1R(pSpecInv, &bufSizeInv);
pBufferInv = ippsMalloc_8u(bufSizeInv);

// borders computation
leftBorderLow = (lenLowI - 1 - anchorLowI) / 2;
leftBorderHigh = (lenHighI - 1 - anchorHighI) / 2;
rightBorderLow = (anchorLowI + 1) / 2;
rightBorderHigh = (anchorHighI + 1) / 2;
int apprLeftBorder = leftBorderLow; // 1
int apprRightBorder = rightBorderLow;
int apprTopBorder = leftBorderLow; // 1
int apprBottomBorder = rightBorderLow;
int detxLeftBorder = leftBorderLow;
int detxRightBorder = rightBorderLow;
int detxTopBorder = leftBorderHigh; // 2
int detxBottomBorder = rightBorderHigh;
int detyLeftBorder = leftBorderHigh;
int detyRightBorder = rightBorderHigh;
int detyTopBorder = leftBorderLow; // 1
int detyBottomBorder = rightBorderLow;
int detxyLeftBorder = leftBorderHigh;
int detxyRightBorder = rightBorderHigh;
int detxyTopBorder = leftBorderHigh; // 2
int detxyBottomBorder = rightBorderHigh;

// allocate the approx image with borders
int appBStep;
IppiSize roiInvSizeBApp = {4 + apprLeftBorder + apprRightBorder, 4 + apprTopBorder + apprBottomBorder};
Ipp32f* pAppB = ippiMalloc_32f_C1 (4 + apprLeftBorder + apprRightBorder, 4 + apprTopBorder + apprBottomBorder, &appBStep);

//ippiCopyWrapBorder_32f_C1R (pApproxDst, approxStep, dstRoiSize, pAppB, appBStep, roiInvSizeBApp, apprTopBorder, apprLeftBorder);
CopySymmBorder_32f_C1R (pApproxDst, approxStep, dstRoiSize, pAppB, appBStep, roiInvSizeBApp, apprTopBorder, apprLeftBorder);
// allocate the X-detail image with borders
int detXBStep;
IppiSize roiInvSizeBDetX = {4 + detxLeftBorder + detxRightBorder, 4 + detxTopBorder + detxBottomBorder};
Ipp32f* pXB = ippiMalloc_32f_C1 (4 + detxLeftBorder + detxRightBorder, 4 + detxTopBorder + detxBottomBorder, &detXBStep);

//ippiCopyWrapBorder_32f_C1R(pDetailXDst, detailXStep, dstRoiSize, pXB, detXBStep, roiInvSizeBDetX, detxTopBorder, detxLeftBorder);
CopySymmBorder_32f_C1R(pDetailXDst, detailXStep, dstRoiSize, pXB, detXBStep, roiInvSizeBDetX, detxTopBorder, detxLeftBorder);

// allocate the Y-detail image with borders
int detYBStep;
IppiSize roiInvSizeBDetY = {4 + detyLeftBorder + detyRightBorder, 4 + detyTopBorder + detyBottomBorder};
Ipp32f* pYB = ippiMalloc_32f_C1 (4 + detyLeftBorder + detyRightBorder, 4 + detyTopBorder + detyBottomBorder, &detYBStep);

// ippiCopyWrapBorder_32f_C1R(pDetailYDst, detailYStep, dstRoiSize, pYB, detYBStep, roiInvSizeBDetY, detyTopBorder, detyLeftBorder);
CopySymmBorder_32f_C1R(pDetailYDst, detailYStep, dstRoiSize, pYB, detYBStep, roiInvSizeBDetY, detyTopBorder, detyLeftBorder);

// allocate the XY-detail image with borders
int detXYBStep;
IppiSize roiInvSizeBDetXY = {4 + detxyLeftBorder + detxyRightBorder, 4 + detxyTopBorder + detxyBottomBorder};
Ipp32f* pXYB = ippiMalloc_32f_C1 (4 + detxyLeftBorder + detxyRightBorder, 4 + detxyTopBorder + detxyBottomBorder, &detXYBStep);

// ippiCopyWrapBorder_32f_C1R(pDetailXYDst, detailXYStep, dstRoiSize, pXYB, detXYBStep, roiInvSizeBDetXY, detxyTopBorder, detxyLeftBorder);
CopySymmBorder_32f_C1R(pDetailXYDst, detailXYStep, dstRoiSize, pXYB, detXYBStep, roiInvSizeBDetXY, detxyTopBorder, detxyLeftBorder);

//performs inverse wavelet transform
if (ippiWTInv_32f_C1R(pAppB+apprTopBorder*(appBStep/4)+apprLeftBorder, appBStep,
pXB+detxTopBorder*(detXBStep/4)+detxLeftBorder, detXBStep,
pYB+detyTopBorder*(detYBStep/4)+detyLeftBorder, detYBStep,
pXYB+detxyTopBorder*(detXYBStep/4)+detxyLeftBorder, detXYBStep,
roiInvSize, pDstInv, stepDstInv, pSpecInv, pBufferInv) != ippStsNoErr) {
printf("something's wrong in ippiWTInv_32f_C1R");
}

// free context memory
ippiWTInvFree_32f_C1R (pSpecInv);
ippiWTFwdFree_32f_C1R (pSpec);

}

2D wavelet transform with symmetric border extension

yoavnaveh — Thu, 11 Nov 2010 08:17:01 GMT

Hello again,

I tried another test with symmetric extensions: actually, there are 2 ways to extend, lets say, the 1D signal {1,2,3,4} in a symmetric fashion. (1) {...3,2,1,1,2,3,4,4,3,2...}. (2) {...3,2,1,2,3,4,3,2...}.
I tries both ways with the 2D wavelet transform and could not reconstruct the image.

hope you can help me, as this problem, and the compression performance that I get as a result, put my whole project with IPP wavelets in danger.

thanks,

Yoav

2D wavelet transform with symmetric border extension

Chao_Y_Intel — Fri, 12 Nov 2010 07:28:05 GMT

Yoav

I read the paper here:

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.99.3753&rep=rep1&type=pdf

2.3.3 part on Symmetric extension.

For Symmetric extension, it does not just simply use 1) , or 2) method. For low frequency part, and high frequency part, it need to combine 1) and 2) extension in different cases.

Based on that paper, I created an example code. It looks to be working here. See attached.

Thanks,
Chao

2D wavelet transform with symmetric border extension

yoavnaveh — Sun, 14 Nov 2010 15:45:02 GMT

Hello Chao,

thank you very much for taking the time and effort to takle (and solve) my problem. I read the paper you refered me to, looked at your implementation, implemented it in my code and it worked on the first run!!

As you can see, this was not a trivial problem, and I would have had to work very hard in order to solve it by myself. Your help was invaluable and I thank you again for doing a great job.

Yoav

2D wavelet transform with symmetric border extension

Chao_Y_Intel — Mon, 15 Nov 2010 02:08:45 GMT

Yoav,

You are welcome. Feel free to share your questions or suggestion on using IPP.

thanks,
Chao

Hi Yoav and Chao,

Le_Guelvouit__Gaëtan — Fri, 26 Apr 2019 11:51:57 GMT

Hi Yoav and Chao,

I have exactly the same problem. Unfortunately, the example source code Chao provided is not available any more. Would it be possible to repost it?

Thanks!

topic 2D wavelet transform with symmetric border extension in Intel® Integrated Performance Primitives

2D wavelet transform with symmetric border extension

2D wavelet transform with symmetric border extension

2D wavelet transform with symmetric border extension

2D wavelet transform with symmetric border extension

2D wavelet transform with symmetric border extension

Hi Yoav and Chao,