def __log_gabor_3d(image, num_scales=4, num_phi_orientations=6,
                   num_theta_orientations=4, min_wavelength=3,
                   scaling_constant=2, center_sigma=0.65, d_theta_sigma=1.5,
                   d_phi_sigma=1.5):
    # Pre-compute sigma values
    theta_sigma = np.pi / num_theta_orientations / d_theta_sigma
    phi_sigma = (2 * np.pi) / num_phi_orientations / d_phi_sigma

    # Allocate space for return structures
    bandpass = np.empty([num_scales, image.shape[0], image.shape[1],
                         image.shape[2]], dtype=np.complex)
    log_gabor = np.empty([num_scales, image.shape[0], image.shape[1],
                          image.shape[2]])
    S = np.zeros(image.shape)
    complex_conv = np.empty([num_scales, num_theta_orientations,
                             num_phi_orientations, image.shape[0],
                             image.shape[1], image.shape[2]], dtype=np.complex)
    tmp_complex_conv = np.empty([num_scales, image.shape[0], image.shape[1],
                                 image.shape[2]], dtype=np.complex)

    # Pre-compute fourier values
    image_fft = np.fft.fftn(image)

    axis0, axis1, axis2 = __adjusted_meshgrid(image.shape)

    radius = np.sqrt(axis0 ** 2 + axis1 ** 2 + axis2 ** 2)
    theta = np.arctan2(axis0, axis1)
    # TODO: Is adding the mean REALLY a good idea?
    m_ab = np.abs(np.mean(radius))
    phi = np.arccos(axis2 / (radius + m_ab))

    radius = np.fft.ifftshift(radius)
    radius[0, 0, 0] = 1.0
    theta = np.fft.ifftshift(theta)
    phi = np.fft.ifftshift(phi)

    sin_theta = np.sin(theta)
    cos_theta = np.cos(theta)
    sin_phi = np.sin(phi)
    cos_phi = np.cos(phi)

    # Compute the lowpass filter
    butterworth_filter = __frequency_butterworth_filter(image.shape, 0.45, 15)

    # Compute radial component of filter
    for s in range(num_scales):
        wavelength = min_wavelength * scaling_constant ** s
        fo = 1.0 / wavelength

        l = np.exp((-np.log(radius / fo) ** 2) /
                   (2.0 * np.log(center_sigma) ** 2))
        l = l * butterworth_filter
        l[0, 0, 0] = 0.0

        log_gabor[s, :, :, :] = l
        bandpass[s, :, :, :] = np.fft.ifft2(image_fft * l)

    # Computer angular component of filter
    for e in range(num_theta_orientations):
        # Pre-compute filter data specific to this orientation
        elevation_angle = e * np.pi / num_theta_orientations

        d_theta_sin = (sin_theta * np.cos(elevation_angle) -
                       cos_theta * np.sin(elevation_angle))
        d_theta_cos = (cos_theta * np.cos(elevation_angle) +
                       sin_theta * np.sin(elevation_angle))
        d_theta = np.abs(np.arctan2(d_theta_sin, d_theta_cos))

        for a in range(num_phi_orientations):
            azimuth_angle = a * 2 * np.pi / num_phi_orientations
            d_phi_sin = (sin_phi * np.cos(azimuth_angle) -
                         cos_phi * np.sin(azimuth_angle))
            d_phi_cos = (cos_phi * np.cos(azimuth_angle) +
                         sin_phi * np.sin(azimuth_angle))
            d_phi = np.abs(np.arctan2(d_phi_sin, d_phi_cos))

            phi_spread = (-d_phi ** 2) / (2 * phi_sigma ** 2)
            theta_spread = (-d_theta ** 2) / (2 * theta_sigma ** 2)
            spread = np.exp(phi_spread + theta_spread)

            # For each scale, multiply by the angular spread
            for s in range(0, num_scales):
                filter_bank = log_gabor[s] * spread

                shifted_filter = np.fft.fftshift(filter_bank)
                S += shifted_filter * np.conjugate(shifted_filter)

                tmp_complex_conv[s, :, :] = np.fft.ifft2(image_fft *
                                                         filter_bank)

            complex_conv[:, e, a, :, :] = tmp_complex_conv[None, None, ...]

    # TODO: Do we need to flip S as in the 2D version?
    return complex_conv, bandpass, S
dataset/ETHPy150Open menpo/menpo/menpo/math/convolution.py/__log_gabor_3d/original