Functions operating on packed format, real to complex conversion, absent ippiInv/ippiDivCRev, etc.

[kavermeer]

Hi,

We're doing some simple image registration, where we calculate:

r = F^-1 ( ( F(A) * conj(F(B)) ) / magn( F(A) * conj(F(B)) ) )

where A is one image, B is the other image, F() is the Fourier transform, F^-1 is the inverse Fourier transform, conj() is the conjugate, magn() is the magnitude and all multiplications and divisions are per-element.

The DFT and FFT functions in IPP return packed data. The library also provides some function to work directly on packed data. So, we can calculate F(A) * conj(F(B)) by ippiMulPackConj (which also returns a packed image - I'll call this t1) and the magnitude of t1 by ippiMagnitudePack (which returns a normal image - I'll call this t2).

At this stage, I have no other option than to start doing ugly things. Two alternatives:
1. Calculate 1/t2, convert that to a complex image, convert that to a packed image and multiply it by t1 (ippiMulPack), then do the inverse Fourier.
2. Convert t2 to a complex image, convert t1 to a complex image, do the division, convert the result to a packed image, then do the inverse Fourier.

In the first case, I'm looking for a way to calculate 1/t2 (problem 1) and to do the conversion from a real image to a complex image (problem 2). In the second case, the conversion from a real image to a complex image (problem 1) is also needed.

For problem 1, I can create an image filled with ones, and do ippiDiv, which seems inefficient. Alternatively, I can use ippsInv or ippsDivCRev, but these are not ippi functions, so I have to handle the step size at the end of the line myself.

For problem 2, I again have to revert to ipps and use ippsRealToCplx, as there seems to be no ippi function to do this.

Given the problems that I encounter, I'm wondering whether I'm using the right approach. Any comments?

Best,
Koen

[Vladimir_Dudnik]

Hello,

there is comment from our experts:

[cpp]For (1) and (2) we have supporting functionality (ippSP functions are not suitable for that as Pack2D format doesnt correspond to Pack 1D):

/* /////////////////////////////////////////////////////////////////////////////
//  Name:           ippiPackToCplxExtend
//
//  Purpose:        Converts an image in RCPack2D format to a complex data image.
//
//  Returns:
//      ippStsNoErr            No errors
//      ippStsNullPtrErr       pSrc == NULL, or pDst == NULL
//      ippStsStepErr          One of the step values is less zero or negative
//      ippStsSizeErr          The srcSize has a field with zero or negative value
//
//  Parameters:
//    pSrc        Pointer to the source image data (point to pixel (0,0))
//    srcSize     Size of the source image
//    srcStep     Step through  the source image
//    pDst        Pointer to the destination image
//    dstStep     Step through the destination image
//  Notes:
*/

IPPAPI (IppStatus, ippiPackToCplxExtend_32f32fc_C1R, (const Ipp32f* pSrc,
        IppiSize srcSize, int srcStep,
        Ipp32fc* pDst, int dstStep ))


/* /////////////////////////////////////////////////////////////////////////////
//  Name:           ippiCplxExtendToPack
//
//  Purpose:        Converts an image in complex data format to RCPack2D image.
//
//  Returns:
//      ippStsNoErr            No errors
//      ippStsNullPtrErr       pSrc == NULL, or pDst == NULL
//      ippStsStepErr          One of the step values is less zero or negative
//      ippStsSizeErr          The srcSize has a field with zero or negative value
//
//  Parameters:
//    pSrc        Pointer to the source image data (point to pixel (0,0))
//    srcSize     Size of the source image
//    srcStep     Step through  the source image
//    pDst        Pointer to the destination image
//    dstStep     Step through the destination image
//  Notes:
*/
IPPAPI(IppStatus,ippiCplxExtendToPack_32fc32f_C1R,(const Ipp32fc* pSrc, int srcStep, IppiSize srcSize,
                                                                           Ipp32f* pDst, int dstStep ))

After numerator/denominator conversion to cplx format ippiDiv_32fc can be used and then  conversion back to Pack format.
[/cpp]

Regards,
Vladimir

[kavermeer]

Hi Vladimir,

Thanks for your answer.

Quoting - Vladimir Dudnik (Intel)

Hello,

there is comment from our experts:

[cpp]For (1) and (2) we have supporting functionality (ippSP functions are not suitable for that as Pack2D format doesnt correspond to Pack 1D):

/* /////////////////////////////////////////////////////////////////////////////
//  Name:           ippiPackToCplxExtend
//
//  Purpose:        Converts an image in RCPack2D format to a complex data image.

/* /////////////////////////////////////////////////////////////////////////////
//  Name:           ippiCplxExtendToPack
//
//  Purpose:        Converts an image in complex data format to RCPack2D image.

After numerator/denominator conversion to cplx format ippiDiv_32fc can be used and then  conversion back to Pack format.
[/cpp]

Regards,
Vladimir

So, if I'm not mistaken, you and/or your experts suggest the following algorithm:

FA = Fourier transform of A (FFTFwd), resulting in a packed image
FB = Fourier transform of B (FFTFwd), resulting in a packed image
t1 = multiplication of FA with the conjugate of FB (MulPackConj), resulting in a packed image
t2 = magnitude of t1 (MagnitudePack), resulting in a real image
t2c = conversion of t2 to complex (how?), resulting in a complex image
t1c = conversion of t1 to complex (PackToCplxExtend), resulting in a complex image
t3c = division of t1 by t2 (IppiDiv), resulting in a complex image
t3 = conversion of t3 to packed (CplxExtendToPack), resulting in a packed image
r = inverse Fourier transform of t3 (FFTInv), resulting in a real image

There are two things I don't like: The conversions and the complex division.

First, regarding the conversions: I understand the conversions between packed and complex, as indicated by your experts. What I was looking for is the 'right' way to convert a real image to a complex image (to create image t2c). I cannot find an ippi function to build a complex image from a real image. I only found a ipps function to do so, as indicated in the original message. What is the best way to convert a real image to a complex image?

Second, regarding the complex division: I only need to calculate a complex image divided by a real image. It seems to me that replacing that by a complex image divided by another complex image is a rather inefficient way to do the calculation. (a+bi) / c = (a/c) + (b/c)i vs. (a+bi) / (c+di ) = (ac+bd) / (c^2+d^2) + (bc-ad) / (c^2+d^2)i.

Best,
Koen

[Naveen_G_Intel]

Hi Koen,

Sorry for silence on this thread for a long time, I would like to know what value you like to have for zero/zero 1, 0 or NaN?

Regards,

Naveen Gv

[igorastakhov]

Hi all,

as this is not the first request for such kind functionality, use the next reference code for division in Pack format (this function performs division of Packed FFT output by its Magnitude - the corresponding IPP function will be available in the next IPP release):

#include

///////////////////////////////////////////////////////////////////////////// */
/* Name: NormalizePack_32f_C1R
// Purpose: Divides pixel values of an image in RCPack2D format
// by its magnitude values and store the results
// in a destination image in RCPack2D format
// Parameters:
// pSrc Pointer to the source image in RCPack2D format
// srcStep Step through the source image
// pDst Pointer to the destination image in RCPack2D format
// dstStep Step through the destination image
// roiSize Size of the source and destination ROI
// Returns:
// ippStsNoErr No errors
// ippStsNullPtrErr One of the pointers is NULL
// ippStsStepErr One of the step values is zero or negative
// ippStsSizeErr The roiSize has a field with negative or zero value
// Notes:
//double precision is used for magnitude calculation
//if magnitude==0, then normalized value is set to 0
//thereare no any special checks for NaNor Inf
*/

#define MAGNREC(x) x = (!x) ? 0 : (1./sqrt(x))

static void ownCalcMagnRec_OneCol (
Ipp32f *src0, Ipp32f *src1, int sline,
Ipp32f *dst0, Ipp32f *dst1, int dline,
int length, int flag)
{
int i;
float re, im;
double val;

for (i = 1; i < length; i += 2) {
re = *src0; im = *src1;
val = re * re + im * im;
MAGNREC(val);
*dst0 = (float)(re * val);
*dst1 = (float)(im * val);
src0 = (Ipp32f*)((Ipp8u*)src0 + sline);
src1 = (Ipp32f*)((Ipp8u*)src1 + sline);
dst0 = (Ipp32f*)((Ipp8u*)dst0 + dline);
dst1 = (Ipp32f*)((Ipp8u*)dst1 + dline);
}
if (!flag) dst0[0] = (src0[0] >= 0) ? 1.f : -1.f;
}

static void ownCalcMagnRec_32f_C1 (
Ipp32fc *src, Ipp32fc *dst, int length)
{
int i;
float re, im;
double val;

for (i = 0; i < length; i ++) {
re = src.re;
im = src.im;
val = re * re + im * im;
MAGNREC(val);
dst.re = (float)(re * val);
dst.im = (float)(im * val);
}
}

IppStatus NormalizePack_32f_C1R (
const Ipp32f* pSrc, int srcStep, Ipp32f* pDst, int dstStep,
IppiSize roiSize)
{
int dw, dh, w, h, j;
Ipp32f *src0, *src1, *dst0, *dst1;
Ipp32fc *sfc, *dfc;

dw = roiSize.width;
dh = roiSize.height;
w = (dw & 1) ? 1 : 0;
h = (dh & 1) ? 1 : 0;

/* 1st column */
pDst[0] = (pSrc[0] >= 0) ? 1.f : -1.f;
src0 = (Ipp32f*)((Ipp8u*)pSrc + srcStep);
dst0 = (Ipp32f*)((Ipp8u*)pDst + dstStep);
src1 = (Ipp32f*)((Ipp8u*)src0 + srcStep);
dst1 = (Ipp32f*)((Ipp8u*)dst0 + dstStep);
ownCalcMagnRec_OneCol (
src0, src1, srcStep+srcStep,
dst0, dst1, dstStep+dstStep,
dh-1, h);

/* last column */
if (!w) {
pDst[dw-1] = (pSrc[dw-1] >= 0) ? 1.f : -1.f;
src0 = (Ipp32f*)((Ipp8u*)pSrc + srcStep);
dst0 = (Ipp32f*)((Ipp8u*)pDst + dstStep);
src1 = (Ipp32f*)((Ipp8u*)src0 + srcStep);
dst1 = (Ipp32f*)((Ipp8u*)dst0 + dstStep);
ownCalcMagnRec_OneCol (
src0+dw-1, src1+dw-1, srcStep+srcStep,
dst0+dw-1, dst1+dw-1, dstStep+dstStep,
dh-1, h);
}

/* region of complex numbers */
dw = (dw - 1) >> 1;
sfc = (Ipp32fc*)(pSrc + 1);
dfc = (Ipp32fc*)(pDst + 1);
for (j = 0; j < dh; j ++) {
ownCalcMagnRec_32f_C1(sfc, dfc, dw);
sfc = (Ipp32fc*)((Ipp8u*)sfc + srcStep);
dfc = (Ipp32fc*)((Ipp8u*)dfc + dstStep);
}

return ippStsNoErr;
}

regards,
Igor