2D wavelet transform with symmetric border extension

[yoavnaveh]

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);

}

[yoavnaveh]

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

[Chao_Y_Intel]

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

[yoavnaveh]

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

[Chao_Y_Intel]

Yoav,

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

thanks,
Chao

[Le_Guelvouit__Gaëtan]

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!