def check_radon_iradon_circle(interpolation, shape, output_size):
    # Forward and inverse radon on synthetic data
    image = _random_circle(shape)
    radius = min(shape) // 2
    sinogram_rectangle = radon(image, circle=False)
    reconstruction_rectangle = iradon(sinogram_rectangle,
                                      output_size=output_size,
                                      interpolation=interpolation,
                                      circle=False)
    sinogram_circle = radon(image, circle=True)
    reconstruction_circle = iradon(sinogram_circle,
                                   output_size=output_size,
                                   interpolation=interpolation,
                                   circle=True)
    # Crop rectangular reconstruction to match circle=True reconstruction
    width = reconstruction_circle.shape[0]
    excess = int(np.ceil((reconstruction_rectangle.shape[0] - width) / 2))
    s = np.s_[excess:width + excess, excess:width + excess]
    reconstruction_rectangle = reconstruction_rectangle[s]
    # Find the reconstruction circle, set reconstruction to zero outside
    c0, c1 = np.ogrid[0:width, 0:width]
    r = np.sqrt((c0 - width // 2)**2 + (c1 - width // 2)**2)
    reconstruction_rectangle[r > radius] = 0.
    print(reconstruction_circle.shape)
    print(reconstruction_rectangle.shape)
    np.allclose(reconstruction_rectangle, reconstruction_circle)

def test_iradon_angles():
    """
    Test with different number of projections
    """
    size = 100
    # Synthetic data
    image = np.tri(size) + np.tri(size)[::-1]
    # Large number of projections: a good quality is expected
    nb_angles = 200
    radon_image_200 = radon(image, theta=np.linspace(0, 180, nb_angles,
                                                     endpoint=False))
    reconstructed = iradon(radon_image_200)
    delta_200 = np.mean(abs(_rescale_intensity(image) - _rescale_intensity(reconstructed)))
    assert delta_200 < 0.03
    # Lower number of projections
    nb_angles = 80
    radon_image_80 = radon(image, theta=np.linspace(0, 180, nb_angles,
                                                    endpoint=False))
    # Test whether the sum of all projections is approximately the same
    s = radon_image_80.sum(axis=0)
    assert np.allclose(s, s[0], rtol=0.01)
    reconstructed = iradon(radon_image_80)
    delta_80 = np.mean(abs(image / np.max(image) -
                           reconstructed / np.max(reconstructed)))
    # Loss of quality when the number of projections is reduced
    assert delta_80 > delta_200

def check_sinogram_circle_to_square(size):
    from skimage.transform.radon_transform import _sinogram_circle_to_square

    image = _random_circle((size, size))
    theta = np.linspace(0.0, 180.0, size, False)
    sinogram_circle = radon(image, theta, circle=True)
    argmax_shape = lambda a: np.unravel_index(np.argmax(a), a.shape)
    print("\n\targmax of circle:", argmax_shape(sinogram_circle))
    sinogram_square = radon(image, theta, circle=False)
    print("\targmax of square:", argmax_shape(sinogram_square))
    sinogram_circle_to_square = _sinogram_circle_to_square(sinogram_circle)
    print("\targmax of circle to square:", argmax_shape(sinogram_circle_to_square))
    error = abs(sinogram_square - sinogram_circle_to_square)
    print(np.mean(error), np.max(error))
    assert argmax_shape(sinogram_square) == argmax_shape(sinogram_circle_to_square)

def test_iradon_sart():
    debug = False

    image = rescale(PHANTOM, 0.8)
    theta_ordered = np.linspace(0., 180., image.shape[0], endpoint=False)
    theta_missing_wedge = np.linspace(0., 150., image.shape[0], endpoint=True)
    for theta, error_factor in ((theta_ordered, 1.),
                                (theta_missing_wedge, 2.)):
        sinogram = radon(image, theta, circle=True)
        reconstructed = iradon_sart(sinogram, theta)

        if debug:
            from matplotlib import pyplot as plt
            plt.figure()
            plt.subplot(221)
            plt.imshow(image, interpolation='nearest')
            plt.subplot(222)
            plt.imshow(sinogram, interpolation='nearest')
            plt.subplot(223)
            plt.imshow(reconstructed, interpolation='nearest')
            plt.subplot(224)
            plt.imshow(reconstructed - image, interpolation='nearest')
            plt.show()

        delta = np.mean(np.abs(reconstructed - image))
        print('delta (1 iteration) =', delta)
        assert delta < 0.02 * error_factor
        reconstructed = iradon_sart(sinogram, theta, reconstructed)
        delta = np.mean(np.abs(reconstructed - image))
        print('delta (2 iterations) =', delta)
        assert delta < 0.014 * error_factor
        reconstructed = iradon_sart(sinogram, theta, clip=(0, 1))
        delta = np.mean(np.abs(reconstructed - image))
        print('delta (1 iteration, clip) =', delta)
        assert delta < 0.018 * error_factor

        np.random.seed(1239867)
        shifts = np.random.uniform(-3, 3, sinogram.shape[1])
        x = np.arange(sinogram.shape[0])
        sinogram_shifted = np.vstack(np.interp(x + shifts[i], x,
                                               sinogram[:, i])
                                     for i in range(sinogram.shape[1])).T
        reconstructed = iradon_sart(sinogram_shifted, theta,
                                    projection_shifts=shifts)
        if debug:
            from matplotlib import pyplot as plt
            plt.figure()
            plt.subplot(221)
            plt.imshow(image, interpolation='nearest')
            plt.subplot(222)
            plt.imshow(sinogram_shifted, interpolation='nearest')
            plt.subplot(223)
            plt.imshow(reconstructed, interpolation='nearest')
            plt.subplot(224)
            plt.imshow(reconstructed - image, interpolation='nearest')
            plt.show()

        delta = np.mean(np.abs(reconstructed - image))
        print('delta (1 iteration, shifted sinogram) =', delta)
        assert delta < 0.022 * error_factor

def wu_hash(fxy):
    # compute radon hash, use 180 deg w/sampl. intervall 1
    fxy_rad = radon(fxy)
    # divide into 40x10 blocks
    bl = blocks(fxy_rad,40,10)
    # compute mean values of the blocks
    ms = []
    for x in xrange(0,len(bl)):
        els = []
        for y in xrange(0,len(bl[x])):
            els.append(np.mean(bl[x][y]))
        ms.append(els)
    # wavelet decomposition with haar wavelet
    # for each column resulting in (approx, detail)
    # approx is thrown away, resulting in a
    # list of 40 lists with each 5 higher order elements
    dec = []
    for x in xrange(0,len(ms)):
        dec.append(pywt.dwt(ms[x],"haar")[1])
    # apply fft to each component and throw imaginary
    # components away
    ffts = map(np.fft.fft,dec)
    reals = []
    for x in xrange(0,len(ffts)):
        reals_of_x = []
        for c in ffts[x]:
            reals_of_x.append(c.real)
        reals.append(reals_of_x)    
    return reals

 def extract(self, image):
     """ Applies the radon transform to the image.
     """
     gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
     try:
         sinogram = radon(gray, theta=self.theta, circle=False)
     except:
         print(traceback.format_exc())
     return sinogram

def _radon_costf(frame, cent, radint, coords):
    """ Radon cost function used in frame_center_radon().
    """
    frame_shifted = frame_shift(frame, coords[0], coords[1])
    frame_shifted_ann = get_annulus(frame_shifted, radint, cent-radint)
    theta = np.linspace(start=0., stop=360., num=frame_shifted_ann.shape[0],    
                    endpoint=False)
    sinogram = radon(frame_shifted_ann, theta=theta, circle=True)
    costf = np.sum(np.abs(sinogram[cent,:]))
    return costf

def low_L2_objfun(u,sinogram, lambdapen,dkx,dky,bkx,bky,thetas,image_res,
                  Circle_Mask,TV,TVT,TVTTV):
    u = np.reshape(u,(image_res,image_res))
    u[0==Circle_Mask] = 0.
    dtheta = (thetas[1]-thetas[0])/180.     #Assume regular spaced angles
    dx = 1./image_res #Assume square image
    J = .5*np.sum((radon(u,thetas,circle=True)-sinogram)**2)*dtheta*dx
    J += .5*lambdapen*np.sum( (dkx-bkx-u*TV)**2 )*dx*dx
    J += .5*lambdapen*np.sum( (dky-bky-TVT*u)**2 )*dx*dx
    return J

def check_radon_center(shape, circle):
    # Create a test image with only a single non-zero pixel at the origin
    image = np.zeros(shape, dtype=np.float)
    image[(shape[0] // 2, shape[1] // 2)] = 1.
    # Calculate the sinogram
    theta = np.linspace(0., 180., max(shape), endpoint=False)
    sinogram = radon(image, theta=theta, circle=circle)
    # The sinogram should be a straight, horizontal line
    sinogram_max = np.argmax(sinogram, axis=0)
    print(sinogram_max)
    assert np.std(sinogram_max) < 1e-6

    def sinogram_radon(self, img, detectors=0, alpha=0, scans=0):
        self.__setup(img, detectors, alpha, scans)

        self._original_image = img
        scale = 1.0*self._detectors/len(img)
        img = rescale(img, scale=scale)

        theta = np.linspace(0., float(self._alpha), self._scans, endpoint=False)
        # print theta
        self._sinogram = radon(img, theta=theta, circle=True)
        return self._sinogram

def tiff_sinos_pad(sinogram,flag,thetas):
    """
    Pads Octopus sinograms (after shape has been extracted elsewhere) so that
    applying radon to other data won't need a mask and will still be the right
    shape.  NOTE: Very Very hacked together.
    """
    from skimage.transform import radon, iradon
    if flag=='FBP':#apply Radon transform to FBP of sinogram to use as data
        imres=sinogram.shape[0]
        sinogram = radon(iradon(sinogram,thetas,output_size=imres,circle=True),
                         thetas,circle=False)
    elif flag=='pad':#Insert 0's into sinogram on either side
        imres=sinogram.shape[0]
        temp = radon(iradon(0.*sinogram,thetas,output_size=imres,circle=True),
                         thetas,circle=False)
        sizediff = abs(sinogram.shape[0]-temp.shape[0])
        if sizediff>0: #padding is needed
            sinogram = np.concatenate((temp[0:np.ceil(sizediff/2.),:],sinogram,
                                           temp[0:np.floor(sizediff/2.),:]))
    return sinogram

    def compute_skew(cls, image):
        image = image - np.mean(image)  # Demean; make the brightness extend above and below zero

        # Do the radon transform and display the result
        sinogram = radon(image)

        # Find the RMS value of each row and find "busiest" rotation,
        # where the transform is lined up perfectly with the alternating dark
        # text and white lines
        r = np.array([rms_flat(line) for line in sinogram.transpose()])
        rotation = np.argmax(r)
        return (90 - rotation) / 100

def low_L2_objgrad(u,sinogram, lambdapen,dkx,dky,bkx,bky,thetas,image_res,
                  Circle_Mask,TV,TVT,TVTTV):
    u = np.reshape(u,(image_res,image_res))
    u[0==Circle_Mask] = 0.
    dtheta = (thetas[1]-thetas[0])/180.
    dx = 1./image_res
    grad = (iradon(radon(u,thetas,circle = True)-sinogram,thetas,
                        output_size = image_res, circle = True,filter = None
                        ))*dtheta*dx*(2.*len(thetas))/np.pi
    grad -= dx*dx*((lambdapen*TVTTV*u+lambdapen*u*TVTTV) +lambdapen*(TVT*(dkx-
                bkx)+(dky-bky)*TV))
    return np.ravel(grad)

def check_radon_iradon_minimal(shape, slices):
    debug = False
    theta = np.arange(180)
    image = np.zeros(shape, dtype=np.float)
    image[slices] = 1.
    sinogram = radon(image, theta)
    reconstructed = iradon(sinogram, theta)
    print('\n\tMaximum deviation:', np.max(np.abs(image - reconstructed)))
    if debug:
        _debug_plot(image, reconstructed, sinogram)
    if image.sum() == 1:
        assert (np.unravel_index(np.argmax(reconstructed), image.shape)
                == np.unravel_index(np.argmax(image), image.shape))

def _radon_costf2(frame, cent, radint, coords):
    """ Radon cost function used in frame_center_radon().
    """
    frame_shifted = frame_shift(frame, coords[0], coords[1])
    frame_shifted_ann = get_annulus(frame_shifted, radint, cent-radint)
    samples = 10
    theta = np.hstack((np.linspace(start=40, stop=50, num=samples, endpoint=False), 
                   np.linspace(start=130, stop=140, num=samples, endpoint=False),
                   np.linspace(start=220, stop=230, num=samples, endpoint=False),
                   np.linspace(start=310, stop=320, num=samples, endpoint=False)))
    
    sinogram = radon(frame_shifted_ann, theta=theta, circle=True)
    costf = np.sum(np.abs(sinogram[cent,:]))
    return costf

def test_radon_circle():
    a = np.ones((10, 10))
    assert_raises(ValueError, radon, a, circle=True)

    # Synthetic data, circular symmetry
    shape = (61, 79)
    c0, c1 = np.ogrid[0:shape[0], 0:shape[1]]
    r = np.sqrt((c0 - shape[0] // 2)**2 + (c1 - shape[1] // 2)**2)
    radius = min(shape) // 2
    image = np.clip(radius - r, 0, np.inf)
    image = _rescale_intensity(image)
    angles = np.linspace(0, 180, min(shape), endpoint=False)
    sinogram = radon(image, theta=angles, circle=True)
    assert np.all(sinogram.std(axis=1) < 1e-2)

    # Synthetic data, random
    image = _random_circle(shape)
    sinogram = radon(image, theta=angles, circle=True)
    mass = sinogram.sum(axis=0)
    average_mass = mass.mean()
    relative_error = np.abs(mass - average_mass) / average_mass
    print(relative_error.max(), relative_error.mean())
    assert np.all(relative_error < 3.2e-3)

def ridgelet(img):
    n1,n2 = np.shape(img)
    lvl = np.int(np.log2(n1))
    rad = ski.radon(img, theta = np.linspace(0,180,n1))
    nr,nt = np.shape(rad)
    rad_ridgelet = np.zeros((lvl,nr,nt))
    for i in np.linspace(0,nt-1,nt):
        rad_ridgelet[:,:,i] = np.reshape(wavelet_1D(rad[:,i],lvl),(lvl,nr))
#    ridgelet = np.zeros((lvl,n1,n2))
    ##Not right
   # for l in np.linspace(0,lvl-1,lvl):
 #       ridgelet[l,:,:] = ski.iradon(rad_ridgelet[l,:,:])
    ridgelet = rad_ridgelet
    return ridgelet

def data_load(dataflag,Phant_size,TIFF_dir,TIFF_slice,res_freq,theta_freq,
                     theta_count,padflag,FITS_dir,FITS_slice,corruptflag,
                     shift_amount,dp_perc):
        if dataflag == 1:
            image = TV_Split_utilities.generateSheppLogan(Phant_size)
            image_res = image.shape[0]
            thetas = np.linspace(0.,180.,theta_count, endpoint = False)
            image = np.expand_dims(image,axis = 2)
            sinogram = np.empty((image.shape[0],theta_count,image.shape[2]))
            for i in range(image.shape[2]):
                sinogram[:,:,i] = radon(image[:,:,i], thetas,circle = True)
                if corruptflag>0:
                    sinogram[:,:,i] = TV_Split_utilities.add_deadpix(
                                        sinogram[:,:,i],dp_perc)

        if dataflag == 2:
            sinogram = TV_Split_utilities.FITS_to_sinos(FITS_dir,FITS_slice,
                        theta_freq)
            sinores = sinogram.shape
            #Downsample
            ind = np.linspace(0,sinores[0]-1,num = (sinores[0]-1)/res_freq
                    ).astype('int')
            sinogram = sinogram[ind,:]
            ind = np.linspace(0,sinores[1]-1,num =
                    (sinores[1]-1)/theta_freq).astype('int')
            sinogram = sinogram[:,ind]
            sinogram = sinogram.astype(float)
            image_res = sinogram.shape[0]
            theta_count = sinogram.shape[1]
            thetas = np.linspace(0.,180.,theta_count, endpoint = False)
        if dataflag == 4:
            sinogram = TV_Split_utilities.tiff_sino_to_image_slice(TIFF_dir,
                                                                      TIFF_slice)
            sinores = sinogram.shape
            #Downsample
            ind = np.linspace(0,sinores[0]-1,num = (sinores[0]-1)/res_freq
                    ).astype('int')
            sinogram = sinogram[ind,:]
            ind = np.linspace(0,sinores[1]-1,num =
                    (sinores[1]-1)/theta_freq).astype('int')
            sinogram = sinogram[:,ind]
            sinogram = sinogram.astype(float)
            image_res = sinogram.shape[0]
            theta_count = sinogram.shape[1]
            thetas = np.linspace(0.,180.,theta_count, endpoint = False)

        #Normalize for all sinograms
        sinogram = sinogram/np.max(np.abs(sinogram))
        return sinogram,thetas,image_res

def check_radon_iradon(interpolation_type, filter_type):
    debug = False
    image = PHANTOM
    reconstructed = iradon(radon(image), filter=filter_type, interpolation=interpolation_type)
    delta = np.mean(np.abs(image - reconstructed))
    print("\n\tmean error:", delta)
    if debug:
        _debug_plot(image, reconstructed)
    if filter_type in ("ramp", "shepp-logan"):
        if interpolation_type == "nearest":
            allowed_delta = 0.03
        else:
            allowed_delta = 0.025
    else:
        allowed_delta = 0.05
    assert delta < allowed_delta

def check_radon_iradon(interpolation_type, filter_type):
    debug = False
    image = PHANTOM
    reconstructed = iradon(radon(image, circle=False), filter=filter_type,
                           interpolation=interpolation_type)
    delta = np.mean(np.abs(image - reconstructed))
    print('\n\tmean error:', delta)
    if debug:
        _debug_plot(image, reconstructed)
    if filter_type in ('ramp', 'shepp-logan'):
        if interpolation_type == 'nearest':
            allowed_delta = 0.03
        else:
            allowed_delta = 0.025
    else:
        allowed_delta = 0.05
    assert delta < allowed_delta