Comparing 2D Opencv DFT transform and FINUFFT: Non uniform grid for forward FFT and uniform grid for INVERSE-FFT

[astroteo]

Dear all,

Premise: I am not an expert in fourier 2D transform and I might need some help in understanding some results

What I am trying to achieve is the following:

1. create an "Image": a set of 2D spatial values on a grid where the value of each pixel represent the intesity of a field. Such grid/Image is supposed to be widely composed by zero values. I create the grid/image by mean of the following function:

import finufft
import numpy as np
import matplotlib.pyplot as plt
import cv2
def create_img():
    BOX_H = 100
    pt_intensity = 1
    pc_width = 160
    pc_height = 100
    pc_res = 2 # [m/px]
    box_w = 20 # [m] width of the parallelipiped
    box_h = BOX_H # [m] height of the parallelepiped
    box_d = 50 # [m] depth of the parallelipiped
    box_origin = [20,30] # x,z coordinates of the box origin (corner)

    # point data
    pc_grid = np.zeros((pc_height, pc_width), dtype=np.float32)

    #Fill box-1
    i_max = int(min((box_origin[1]*pc_res+box_d)/pc_res, pc_height))
    j_max = int(min((box_origin[0]*pc_res+box_w)/pc_res, pc_width))
    print(str(i_max) + ' ' + str(j_max))
    for i in range(box_origin[1], i_max):
        for j in range(box_origin[0], j_max):
            pc_grid[i,j] = box_h

    #Fill box-2
    box_w = 50 # [m] width of the parallelipiped
    box_h = BOX_H -0.5 * BOX_H # [m] height of the parallelepiped
    box_d = 20 # [m] depth of the parallelipiped
    box_origin = [80,50] # x,z coordinates of the box origin (corner)
    i_max = int(min((box_origin[1]*pc_res+box_d)/pc_res, pc_height))
    j_max = int(min((box_origin[0]*pc_res+box_w)/pc_res, pc_width))

    print(str(i_max) + ' ' + str(j_max))
    for i in range(box_origin[1], i_max):
        for j in range(box_origin[0], j_max):
            pc_grid[i,j] = box_h
    return pc_grid

to visualize the img i use pyplot as follow:

raster_img = create_img()
plt.figure()
plt.title('raster-image') 
plt.imshow(raster_img)

2. compute the NU-FFT transform via FINUFFT only on non zero points of the image/grid:

therefore to select only NON-ZERO points from the original image I apply the following loop:

fx_s = []
fy_s = []
field =  []

cnt = 0

for i in range(0,raster_img.shape[0]):
    for j in range(0,raster_img.shape[1]):

        if raster_img[i][j] > 1e-6:
            fx_s.append( i/raster_img.shape[1] * 2 * np.pi)# (i - RASTER_HEIGHT/2)/RASTER_HEIGHT * np.pi)
            fy_s.append( j/raster_img.shape[0]  * 2 * np.pi)# (j - RASTER_WIDTH/2)/RASTER_WIDTH  * np.pi)

            field.append(raster_img[i][j]) #* (RASTER_HEIGHT * RASTER_WIDTH) )# raster_img[i][j] )

fx_s = np.array(fx_s).astype(np.float64)
fy_s = np.array(fy_s).astype(np.float64)  

field = np.array(field).astype(np.complex128)

Then I compute the NUFFT by using the 'type-1' function, forcing FINUFFT to return the Fourier coefficients on a grid of the same size of the same size of the original image ( and visualize the spectra):

result = finufft.nufft2d1(
                              fx_s, 
                              fy_s,
                              field,
                              (raster_img.shape[0],raster_img.shape[1]),
                              isign=+1)
plt.figure()
plt.title('fft-image [NUFFT]') 
plt.imshow(np.abs(result))

3. Finally I compute the INVERSE-FFT using both FINUFFT (2d type-1 again) and OPENCV and comparing the result:

First I create the spatial coordinate (coincident with the original grid, scaled in between [0,2pi]):

sx = []
sy = []
spec =  []

for i in range(0,raster_img.shape[0]):
    spec_row = []
    for j in range(0,raster_img.shape[1]):
        sx.append( i/raster_img.shape[0] * 2 *np.pi)
        sy.append( j/raster_img.shape[1]* 2 * np.pi) 
        spec_row.append(result[i][j])
    spec.append(spec_row)

sx = np.array(sx)
sy = np.array(sy)
spec = np.array(spec)

Then I compute the inverse transform via FINUFFT:

i_result = finufft.nufft2d1(
                              sx.flatten(), 
                              sy.flatten(), 
                              spec.flatten(),
                              (raster_img.shape[0] ,raster_img.shape[1]),
                              isign= -1)
i_result = np.fft.fftshift(i_result)

plt.figure()
plt.title('inverse-fft-image [via FINUFFT]:') 
plt.imshow(np.abs(i_result))

Second I compute the inverse transform via DFT algorithm from OPENCV on the same spectra spec / result (which is uniform):

cv_img = cv2.merge([np.real(result),np.imag(result)])
cv_img_back = np.fft.ifftshift(nufft_img)
img_back_cv = cv2.idft(cv_img_back)
img_back_cv_mag = cv2.magnitude(img_back_cv[:,:,0], img_back_cv[:,:,1])

plt.figure()
plt.title('inverse-fft-image [via OPEN-CV]:') 
plt.imshow(img_back_cv_mag)

If I compute the difference between the 2 inverse transform I obtain almost the same result:

cv_img_cplx = img_back_cv[:,:,0] + 1.j *img_back_cv[:,:,1]

max_error_finufft_opencv = np.max(np.abs(cv_img_cplx - i_result))
avg_error_finufft_opencv = np.mean(np.abs(cv_img_cplx - i_result))

plt.figure()
plt.title('difference between OPENCV-INVERSE and FINUFFT-INVERSE,\n max_error = ' + str(max_error_finufft_opencv) +
                                          '\n average-error = ' +str(max_error_finufft_opencv)) 
plt.imshow(np.abs(np.abs(cv_img_cplx - i_result)))

Except from a scaling the resulting images are equal.

4. It also possible to observe that the spectra computed via OPENCV via DFT on the uniform (full) grid and the one I showed before are almost equal ( again except from the scaling) :

   
   cv_fft_img = cv2.dft(raster_img,flags = cv2.DFT_COMPLEX_OUTPUT)
   cv_fft_img_shift = np.fft.fftshift(cv_fft_img)
   cv_fft_img_mag = cv2.magnitude(cv_fft_img_shift[:,:,0], cv_fft_img_shift[:,:,1])

plt.figure() plt.title('fft-image[OPENCV]') plt.imshow(cv_fft_img_mag)

plt.figure() plt.title('fft-image [NUFFT]') plt.imshow(np.abs(result))

cv_fft_cplx = cv_fft_img_shift[:,:,0] + 1.j *cv_fft_img_shift[:,:,1]

max_error_finufft_opencv_direct = np.max(np.abs(cv_fft_cplx - result)) avg_error_finufft_opencv_direct = np.mean(np.abs(cv_fft_cplx -result))

plt.figure() plt.title('difference between OPENCV-DIRECT and FINUFFT-DIRECT,\n max_error = ' + str(max_error_finufft_opencv_direct) + '\n average-error = ' +str(max_error_finufft_opencv_direct)) plt.imshow(np.abs(np.abs(cv_fft_cplx - result)))

![fftnufft](https://user-images.githubusercontent.com/19655904/195860644-acde003a-b942-4989-911b-6620de5c8705.png)
![fftcv](https://user-images.githubusercontent.com/19655904/195860650-ef981781-247b-4d0e-9260-4fc11f2bebe8.png)

I am wondering how it is possible that the transform results almost equal even if via FINUFFT i am computing the fourier transform only on the non-zero points of the original field.

If I am missing something please let me know, thanks in advance for yout support! 

[ahbarnett]

Hi Matteo, Please see other Issue response. This code has too much stuff not relevant to the task. Several points:

• The code create_img seems to fill a rectangle with one value and s subrectangle with another; the array is not sparse.
• Unless the sparsity is 1% or less, t's probably not worht using FINUFFT; just use an FFT.
• The inverse of a type 1 is not a type 2. Inverses are not as with FFTs. I hope that helps to get you started. The math of what the NUFFT does is very clear in our docs - please study them. Best, Alex

On Fri, Oct 14, 2022 at 8:08 AM matteo baiguera @.***> wrote:

Dear all,

Premise: I am not an expert in fourier 2D transform and I might need some help in understanding some results

What I am trying to achieve is the following:

1.

create an "Image": a set of 2D spatial values on a grid where the value of each pixel represent the intesity of a fiels. Such grid/Image is supposed to be widely composed by zero values. I create the grid/image by mean of the following function:

import finufftimport numpy as npimport matplotlib.pyplot as pltimport cv2def create_img(): BOX_H = 100 pt_intensity = 1 pc_width = 160 pc_height = 100 pc_res = 2 # [m/px] box_w = 20 # [m] width of the parallelipiped box_h = BOX_H # [m] height of the parallelepiped box_d = 50 # [m] depth of the parallelipiped box_origin = [20,30] # x,z coordinates of the box origin (corner)

# point data
pc_grid = np.zeros((pc_height, pc_width), dtype=np.float32)

#Fill box-1
i_max = int(min((box_origin[1]*pc_res+box_d)/pc_res, pc_height))
j_max = int(min((box_origin[0]*pc_res+box_w)/pc_res, pc_width))
print(str(i_max) + ' ' + str(j_max))
for i in range(box_origin[1], i_max):
    for j in range(box_origin[0], j_max):
        pc_grid[i,j] = box_h

#Fill box-2
box_w = 50 # [m] width of the parallelipiped
box_h = BOX_H -0.5 * BOX_H # [m] height of the parallelepiped
box_d = 20 # [m] depth of the parallelipiped
box_origin = [80,50] # x,z coordinates of the box origin (corner)
i_max = int(min((box_origin[1]*pc_res+box_d)/pc_res, pc_height))
j_max = int(min((box_origin[0]*pc_res+box_w)/pc_res, pc_width))

print(str(i_max) + ' ' + str(j_max))
for i in range(box_origin[1], i_max):
    for j in range(box_origin[0], j_max):
        pc_grid[i,j] = box_h
return pc_grid

to visualize the img i use pyplot as follow:

raster_img = create_img()plt.figure()plt.title('raster-image') plt.imshow(raster_img)

2.> compute the NU-FFT transform via FINUFFT only on non zero points of the image/grid:

therefore to select only NON-ZERO points from the original image I apply the following loop:

fx_s = []fy_s = []field = [] cnt = 0 for i in range(0,raster_img.shape[0]): for j in range(0,raster_img.shape[1]):

    if raster_img[i][j] > 1e-6:
        fx_s.append( i/raster_img.shape[1] * 2 * np.pi)# (i - RASTER_HEIGHT/2)/RASTER_HEIGHT * np.pi)
        fy_s.append( j/raster_img.shape[0]  * 2 * np.pi)# (j - RASTER_WIDTH/2)/RASTER_WIDTH  * np.pi)

        field.append(raster_img[i][j]) #* (RASTER_HEIGHT * RASTER_WIDTH) )# raster_img[i][j] )

        fx_s = np.array(fx_s).astype(np.float64)fy_s = np.array(fy_s).astype(np.float64)

field = np.array(field).astype(np.complex128)

Then I compute the NUFFT by using the 'type-1' function, forcing FINUFFT to return the Fourier coefficients on a grid of the same size of the same size of the original image ( and visualize the spectra):

result = finufft.nufft2d1( fx_s, fy_s, field, (raster_img.shape[0],raster_img.shape[1]), isign=+1)plt.figure()plt.title('fft-image [NUFFT]') plt.imshow(np.abs(result))

3.> Finally I compute the INVERSE-FFT using both FINUFFT (2d type-1 again) and OPENCV and comparing the result:

First I create the spatial coordinate (coincident with the original grid, scaled in between [0,2pi]):

sx = []sy = []spec = []

for i in range(0,raster_img.shape[0]): spec_row = [] for j in range(0,raster_img.shape[1]): sx.append( i/raster_img.shape[0] 2 np.pi) sy.append( j/raster_img.shape[1] 2 np.pi) spec_row.append(result[i][j]) spec.append(spec_row)

  sx = np.array(sx)sy = np.array(sy)spec = np.array(spec)

Then I compute the inverse transform via FINUFFT:

i_result = finufft.nufft2d1( sx.flatten(), sy.flatten(), spec.flatten(), (raster_img.shape[0] ,raster_img.shape[1]), isign= -1)i_result = np.fft.fftshift(i_result) plt.figure()plt.title('inverse-fft-image [via FINUFFT]:') plt.imshow(np.abs(i_result))

Second I comoute the inverse transform via DFT algorithm from OPENCV on the same spectra 'spec/result' (which is uniform):

cv_img = cv2.merge([np.real(result),np.imag(result)])cv_img_back = np.fft.ifftshift(nufft_img)img_back_cv = cv2.idft(cv_img_back)img_back_cv_mag = cv2.magnitude(img_back_cv[:,:,0], img_back_cv[:,:,1]) plt.figure()plt.title('inverse-fft-image [via OPEN-CV]:') plt.imshow(img_back_cv_mag)

If I compute the difference between the 2 inverse transform I obtain almost the same result:

cv_img_cplx = img_back_cv[:,:,0] + 1.j *img_back_cv[:,:,1] max_error_finufft_opencv = np.max(np.abs(cv_img_cplx - i_result))avg_error_finufft_opencv = np.mean(np.abs(cv_img_cplx - i_result))

plt.figure()plt.title('difference between OPENCV-INVERSE and FINUFFT-INVERSE,\n max_error = ' + str(max_error_finufft_opencv) + '\n average-error = ' +str(max_error_finufft_opencv)) plt.imshow(np.abs(np.abs(cv_img_cplx - i_result)))

Except from a scaling the resulting images are equal.

4.>It also possible to observe that the spectra computed via OPENCV via DFT on the uniform (full) grid and the one I showed before are almost equal ( again except from the scaling) :

cv_fft_img = cv2.dft(raster_img,flags = cv2.DFT_COMPLEX_OUTPUT)cv_fft_img_shift = np.fft.fftshift(cv_fft_img)cv_fft_img_mag = cv2.magnitude(cv_fft_img_shift[:,:,0], cv_fft_img_shift[:,:,1]) plt.figure()plt.title('fft-image[OPENCV]') plt.imshow(cv_fft_img_mag) plt.figure()plt.title('fft-image [NUFFT]') plt.imshow(np.abs(result)) cv_fft_cplx = cv_fft_img_shift[:,:,0] + 1.j *cv_fft_img_shift[:,:,1] max_error_finufft_opencv_direct = np.max(np.abs(cv_fft_cplx - result))avg_error_finufft_opencv_direct = np.mean(np.abs(cv_fft_cplx -result))

plt.figure()plt.title('difference between OPENCV-DIRECT and FINUFFT-DIRECT,\n max_error = ' + str(max_error_finufft_opencv_direct) + '\n average-error = ' +str(max_error_finufft_opencv_direct)) plt.imshow(np.abs(np.abs(cv_fft_cplx - result)))

I am wondering how it is possible that the transform results almost equal even if via FINUFFT i am computing the fourier transform only on the non-zero points of the original field.

If I am missing something please let me know, thanks in advance for yout support!

— Reply to this email directly, view it on GitHub https://github.com/flatironinstitute/finufft/issues/242, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACNZRSUHH7TJKAJBDQM4NM3WDFENFANCNFSM6AAAAAARFFXPRM . You are receiving this because you are subscribed to this thread.Message ID: @.***>

-- *-------------------------------------------------------------------~^`^~._.~' |\ Alex Barnett Center for Computational Mathematics, Flatiron Institute | \ http://users.flatironinstitute.org/~ahb 646-876-5942

[astroteo]

Dear Alex,

As you might have already seen, I just closed the other issue as now I can correctly transform and inverse-transform via finufft2d3 .

I am trying to improve some of the results metioned in this reference: FFT-based Terrain Segmentation for Underwater Mapping where in section 3A. is stated that:

Interpolation in the data sets processed here is not likely to modify the frequency content of the original elevation signal since empty cells represent at most 10% of the elevation grid. However, this is not true in more sparsescans, such as side-looking Velodyne ..... . Alternatively, non-uniform DFTs or Polar FFTs could be used.

In short what we want to achieve is to:

1. fit a 3d point-cloud from a laser scan   [we first squash the 3d-points into a plane obtaining 2D raster image where the intensity of each pixel represnt the elevetion of each 3D point],
2. transform from spatial domain to frequency domain
3. filter-out high frequencies content and transform back into spatial domain.

What the authors states is that if the 10% of the raster image is "empty" [meaning that no point is present there] a non uniform transform is recommended. Since only a small % of the beams from the laser are reflected ( because of emptyness in the scene) the raster we obtain can be empty for wide regions.

Maybe the word "sparse" i previously used is not 100% correct from a mathematical perspective: we are actually dealing with a 2D-spatial field empty for the 10%--90% of a 2D-grid, which actually creates a non-uniform set of points.

Moreover at my "level of sparsity" I can see some differences in the spectra between FINUFFT on non-unifrom and DFT (from opencv) including the full grid.

raster_img = create_img()

plt.figure()
plt.title('raster img \n [INPUT FIELD TO FINUFFT & OPENCV] ') 
plt.imshow(raster_img)
plt.figure()

fx_s = []
fy_s = []

sx_s = []
sy_s = []

field = []

for i in range(0,raster_img.shape[0]):
    for j in range(0,raster_img.shape[1]):
        if raster_img[i][j] > 1e-6:

            fx_s.append( i * (2*np.pi / raster_img.shape[0]))
            fy_s.append( j * (2*np.pi / raster_img.shape[1]))

            field.append(raster_img[i][j])

        sx_s.append( i - raster_img.shape[0]/2 )
        sy_s.append( j - raster_img.shape[1]/2 )

fx_s = np.array(fx_s).astype(np.float64)
fy_s = np.array(fy_s).astype(np.float64)

sx_s = np.array(sx_s).astype(np.float64)
sy_s = np.array(sy_s).astype(np.float64)  

field = np.array(field).astype(np.complex128)

# Spatial-domain-> Freq-domain : FINUFFT
fft_img = finufft.nufft2d3(fx_s, 
                           fy_s, 
                           field.flatten(), 
                           sx_s,
                           sy_s,
                           isign=1
                           )
# Spatial-domain-> Freq-domain : DFT/OPENCV
opencv_fft_img =  cv2.dft(raster_img,flags = cv2.DFT_COMPLEX_OUTPUT)
opencv_fft_img = np.fft.fftshift(opencv_fft_img)
opencv_fft_img = cv2.magnitude(opencv_fft_img[:,:,0], opencv_fft_img[:,:,1])

plt.figure()
plt.title('DFT [OPENCV] ') 
plt.imshow(opencv_fft_img)
plt.figure()

plt.figure()
plt.title('fft-image \n [NUFFT TYPE-3] ') 
plt.imshow(np.abs(fft_img.reshape(raster_img.shape[0],raster_img.shape[1])))
plt.figure()

#normalize the 2 spectra
fft_img_s = fft_img.reshape(raster_img.shape[0],raster_img.shape[1])/np.max(np.abs(fft_img)) 
opencv_fft_img_s = opencv_fft_img/np.max(opencv_fft_img)

#subtract the 2 spectra
diff_fft_imgs = np.abs(np.abs(fft_img_s) - opencv_fft_img_s)

plt.figure()
plt.title(' [NUFFT] - DFT [OPENCV]    ') 
plt.imshow(np.abs(diff_fft_imgs ))
plt.figure()

Thanks in advance if you might be willling to support us on such comparison between DFT and NU-FFT and many thanks again for having solved the issue regarding the scaling.

Best regards, Matteo Baiguera.