Return dft in polar (or log-polar) units, the angle step
    (and the log base)

    Parameters
    ----------
    im: 2d array
        The image
    nangle: number, optional
        The number of angles in the polar representation
    radiimax: number, optional
        The number of radius in the polar representation
    isshiftdft: boolean, default False
        True if the image is pre processed (DFT + fftshift)
    truesize: 2 numbers, required if isshiftdft is True
        The true size of the image
    logpolar: boolean, default False
        True if want the log polar representation instead of polar
    logoutput: boolean, default False
        True if want the log of the output
    interpolation: string, default 'bilinear'
        ('bicubic', 'bilinear', 'nearest') The interpolation
        technique. (For now, avoid bicubic)

    Returns
    -------
    im: 2d array
        The (log) polar representation of the input image
    log_base: number, only if logpolar is True
        the log base if this is log polar representation

    Notes
    -----
    radiimax is the maximal radius (log of radius if logpolar is true).
    if not provided, it is deduced from the image size

    To get log-polar, set logpolar to True
    log_base is the base of the log. It is deduced from radiimax.
    Two images that will be compared should therefore have the same radiimax.

def polar_fft(im, nangle=None, radiimax=None, *, isshiftdft=False,
              truesize=None, logpolar=False, logoutput=False,
              interpolation='bilinear'):
    """Return dft in polar (or log-polar) units, the angle step
    (and the log base)

    Parameters
    ----------
    im: 2d array
        The image
    nangle: number, optional
        The number of angles in the polar representation
    radiimax: number, optional
        The number of radius in the polar representation
    isshiftdft: boolean, default False
        True if the image is pre processed (DFT + fftshift)
    truesize: 2 numbers, required if isshiftdft is True
        The true size of the image
    logpolar: boolean, default False
        True if want the log polar representation instead of polar
    logoutput: boolean, default False
        True if want the log of the output
    interpolation: string, default 'bilinear'
        ('bicubic', 'bilinear', 'nearest') The interpolation
        technique. (For now, avoid bicubic)

    Returns
    -------
    im: 2d array
        The (log) polar representation of the input image
    log_base: number, only if logpolar is True
        the log base if this is log polar representation

    Notes
    -----
    radiimax is the maximal radius (log of radius if logpolar is true).
    if not provided, it is deduced from the image size

    To get log-polar, set logpolar to True
    log_base is the base of the log. It is deduced from radiimax.
    Two images that will be compared should therefore have the same radiimax.
    """
    im = np.asarray(im, dtype=np.float32)
    # get dft if not already done
    if not isshiftdft:
        truesize = im.shape
        # substract mean to avoid having large central value
        im = im - im.mean()
        im = centered_mag_sq_ccs(dft_optsize(im))
    # We need truesize! otherwise border effects.
    assert(truesize is not None)

    # the center is shifted from 0,0 to the ~ center
    #(eg. if 4x4 image, the center is at [2,2], as if 5x5)
    qshape = np.asarray([im.shape[0] // 2, im.shape[1]])
    center = np.asarray([qshape[0], 0])

    # if the angle Step is not given, take the number of pixel
    # on the perimeter as the target #=range/step

    if nangle is None:
        # TODO: understand why nangle need to be exactly truesize
        nangle = np.min(truesize)  # range is pi, nbangle = 2r =~pi r
#        nangle-=2

    # get the theta range
    theta = np.linspace(-np.pi / 2, np.pi / 2, nangle,
                        endpoint=False, dtype=np.float32)

    # For the radii, the units are comparable if the log_base and radiimax are
    # the same. Therefore, log_base is deduced from radiimax
    # The step is assumed to be 1
    if radiimax is None:
        radiimax = qshape.min()

    # also as the circle is an ellipse in the image,
    # we want the radius to be from 0 to 1
    if logpolar:
        # The log base solves log_radii_max=log_{log_base}(linear_radii_max)
        # where we decided arbitrarely that linear_radii_max=log_radii_max
        log_base = np.exp(np.log(radiimax) / radiimax)
        radius = ((log_base ** np.arange(0, radiimax, dtype=np.float32))
                  / radiimax)
    else:
        radius = np.linspace(0, 1, radiimax, endpoint=False,
                             dtype=np.float32)

    # get x y coordinates matrix (The units are 0 to 1, therefore a circle is
    # represented as an ellipse)
    y = cv2.gemm(np.sin(theta), radius, qshape[0], 0, 0,
                 flags=cv2.GEMM_2_T) + center[0]
    x = cv2.gemm(np.cos(theta), radius, qshape[1], 0, 0,
                 flags=cv2.GEMM_2_T) + center[1]
    interp = cv2.INTER_LINEAR
    if interpolation == 'bicubic':
        interp = cv2.INTER_CUBIC
    if interpolation == 'nearest':
        interp = cv2.INTER_NEAREST

    # get output
    output = cv2.remap(im, x, y, interp)  # LINEAR, CUBIC,LANCZOS4
    # apply log
    if logoutput:
        output = cv2.log(output)

    if logpolar:
        return output, log_base
    else:
        return output