ipprWarp (3D affine transformation) failures (test code attached)

[Ken_Thornton]

(See similar post about ipprResize)

I am running 2020.0.2 IPP in 32-bit mode on 64-bit Windows 10 Enterprise (10.0.19042).

I have been using ipprWarp to resize a (dimIn x dimIn x dimIn) input volume to size (dimOut x dimOut x dimOut). I have been testing all combinations of sizes 2 <= dimIn, dimOut <= 50 where dimIn != dimOut.

Out of 2352 tests, I get 899 failures. A failure means that an Intel interpolated voxel value differs by more than 2 from what I think it should be. 

I assume that ipprWarp works as follows:

1. For each output volume voxel, get the 3D coordinates p of the center of that voxel
2. Determine corresponding 3D coordinates q in the input volume
3. Do a trilinear interpolation at q in the input volume and assign that value to the output voxel

I have attached files ToIntelWarp.h / ToIntelWarp.cpp to replicate my test and ipprWarp_failures.txt which give my results.

[ShanmukhS_Intel]

Hi, 

Thanks for reaching out to us. 

We would like to recommend you to try executing the code in latest version of oneAPI(2021.3) and get back to us if issue persists. 

Best Regards,

Shanmukh.SS

[Ken_Thornton]

I just installed 2021.3.0 and re-ran the test and got the same results.

Below is the library version information.

Please look at my test code and see if you can replicate the problems.

Thanks,

Ken Thornton

[Gennady_F_Intel]

Yes, we confirmed the problem you reported and will investigate the case. The thread will be updated.

[Andrey_B_Intel]

Hi Ken.

I am attaching simple and minimized code to demonstrate formulas of affine linear transform in IPP for 2-D. It is easy to extended for 3d.

void reference_affine_linear_16u(const Ipp16u *pSrc, int srcStep, IppiSize srcSize,
                            Ipp16u *pDst, int dstStep, IppiSize dstSize,
    double mat[2][3])
{
    double inv[2][3];
    int i, j;
    srcStep >>= 1;
    dstStep >>= 1;
    
    //inverse matrix
    double det =  mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1];
    
    inv[0][0]  =  mat[1][1] / det;
    inv[0][1]  = -mat[0][1] / det;
    inv[1][0]  = -mat[1][0] / det;
    inv[1][1]  =  mat[0][0] / det;

    inv[0][2] = -mat[0][2];
    inv[1][2] = -mat[1][2];
    
    //check if in source image
    #define IN(X,Y) ((X)>= 0 && (X) < srcSize.width && (Y) >= 0 && (Y) < srcSize.height)
    for (i = 0; i < dstSize.height; i++) {
        for (j = 0; j < dstSize.width; j++) {
            double x, y;
            x = (j)* inv[0][0] + (i)* inv[0][1] + inv[0][2]; 
            y = (j)* inv[1][0] + (i)* inv[1][1] + inv[1][2];
            int x0 = (int)x;
            int x1 = (int)x+1;
            int y0 = (int)y;
            int y1 = (int)y+1;
            //if all 4 pixels inside source
            if (IN(x0, y0) && IN(x1, y0) && IN(x0, y1) && IN(x1, y1)){
                double kx = x - floor(x);
                double ky = y - floor(y);
                int s00 = pSrc[srcStep*y0 + x0];
                int s01 = pSrc[srcStep*y0 + x1];
                int s10 = pSrc[srcStep*y1 + x0];
                int s11 = pSrc[srcStep*y1 + x1];
                //interpolate linear by x
                double u = kx * (s01 - s00) + s00;
                double b = kx * (s11 - s10) + s10;
                //by y
                double v = ky * (b - u) + u;
                //cvt from double and write to dst
                pDst[dstStep*i + j] = (Ipp16u)v; /* just to simplify*/
                                                 /* actually IPP rounds to nearest*/
              }
        }
    }
}

 I hope this code provides more details how center of pixels(voxels) is calculated in IPP.

Thanks. 

[Ken_Thornton]

Andrew,

The 2D code above doesn't really help as my issues are with 3D resizing/warping. 

Have you run my test code?

I showed that out of 2,352 resizing tests (i.e. resizing a (dimIn,dimIn,dimIn) volume to (dimOut,dimOut,dimOut) for various combinations of dimIn and dimOut), 186 tests failed.

For example, in your 2D code above, (0,0) is the center of the first pixel. In the Intel 3D resizing code, (0.5,0.5,0.5) is the center of the first voxel.

Also, I don't think the 3D code requires all 8 voxels to be inside the src volume (the 3D equivalent to your 2D test on line 33).

Regards,

Ken Thornton