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_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 check_iradon_center(size, theta, circle):
    debug = False
    # Create a test sinogram corresponding to a single projection
    # with a single non-zero pixel at the rotation center
    if circle:
        sinogram = np.zeros((size, 1), dtype=np.float)
        sinogram[size // 2, 0] = 1.0
    else:
        diagonal = int(np.ceil(np.sqrt(2) * size))
        sinogram = np.zeros((diagonal, 1), dtype=np.float)
        sinogram[sinogram.shape[0] // 2, 0] = 1.0
    maxpoint = np.unravel_index(np.argmax(sinogram), sinogram.shape)
    print("shape of generated sinogram", sinogram.shape)
    print("maximum in generated sinogram", maxpoint)
    # Compare reconstructions for theta=angle and theta=angle + 180;
    # these should be exactly equal
    reconstruction = iradon(sinogram, theta=[theta], circle=circle)
    reconstruction_opposite = iradon(sinogram, theta=[theta + 180], circle=circle)
    print("rms deviance:", np.sqrt(np.mean((reconstruction_opposite - reconstruction) ** 2)))
    if debug:
        import matplotlib.pyplot as plt

        imkwargs = dict(cmap="gray", interpolation="nearest")
        plt.figure()
        plt.subplot(221)
        plt.imshow(sinogram, **imkwargs)
        plt.subplot(222)
        plt.imshow(reconstruction_opposite - reconstruction, **imkwargs)
        plt.subplot(223)
        plt.imshow(reconstruction, **imkwargs)
        plt.subplot(224)
        plt.imshow(reconstruction_opposite, **imkwargs)
        plt.show()

    assert np.allclose(reconstruction, reconstruction_opposite)

def tiltaxisalign(im_series,tilt_angles,shift_and_tilt=('hold','hold')):
    series_shape = np.shape(im_series)
    
    new_series = im_series.copy()
    final_series = im_series.copy()
    
    #deg0_int = input('Which image is the 0 degree image? ')
    midy = int(series_shape[1]/2)
    
    axis_shift = shift_and_tilt[0]
    axis_tilt = shift_and_tilt[1]
    
    if axis_shift == 'hold':
        shift_continue = 1
    
        while shift_continue == 1:
            plt.imshow(iradon(np.rot90(new_series[:,midy,:]),  # rot
            theta = tilt_angles,output_size = series_shape[2]))# anti-clokwise
            plt.show()
            
            axis_shift = float(input('By how many pixels from the original mid-point should the tilt axis be shifted? '))
            
            for i in range(series_shape[0]):
                new_series[i,:,:] = interpolation.shift(im_series.copy()[i,:,:],(0,axis_shift)) # shift along np x-axis
                
            shift_continue = int(input('Would you like to apply further image shifts (1 for yes, 0 for no)? '))
                
    for i in range(series_shape[0]):
        final_series[i,:,:] = interpolation.shift(final_series[i,:,:],(0,axis_shift))
        
    topy = int(3*series_shape[1]/8)
    bottomy = int(5*series_shape[1]/8)

    if axis_tilt == 'hold':
        tilt_series = final_series.copy()
        tilt_continue = 1
        while tilt_continue == 1:
            plt.imshow(iradon(np.rot90(new_series[:,topy,:]), theta = tilt_angles,output_size = series_shape[2]))
            plt.show()
            plt.imshow(iradon(np.rot90(new_series[:,bottomy,:]), theta = tilt_angles,output_size = series_shape[2]))
            plt.show()
            
            axis_tilt = float(input('By what angle from the original y axis (in degrees) should the tilt axis be rotated? '))
            
            for i in range(series_shape[0]):
                new_series[i,:,:] = interpolation.rotate(tilt_series.copy()[i,:,:],axis_tilt,reshape=False)
                    
            tilt_continue = int(input('Would you like to try another tilt angle (1 for yes, 0 for no)? '))
                    
    for i in range(series_shape[0]):
        final_series[i,:,:] = interpolation.rotate(final_series[i,:,:],axis_tilt,reshape=False)
            
    shift_and_tilt = (axis_shift,axis_tilt)
            
    return(final_series, shift_and_tilt)

def tomo_reconstruct_layer(rad_stack,cross_sectional_dim,layer_row=1024,start_tomo_ang=0., end_tomo_ang=360.,tomo_num_imgs=360, center=0.,pixel_size=0.00148):
    sinogram=np.squeeze(rad_stack[:,layer_row,:])

    rotation_axis_pos=-int(np.round(center/pixel_size))
    #rotation_axis_pos=13
    
    theta = np.linspace(start_tomo_ang, end_tomo_ang, tomo_num_imgs, endpoint=False)
    
    max_rad=int(cross_sectional_dim/pixel_size/2.*np.sqrt(2.))
    
    if rotation_axis_pos>=0:
        sinogram_cut=sinogram[:,2*rotation_axis_pos:]
    else:
        sinogram_cut=sinogram[:,:(2*rotation_axis_pos)]
    
    dist_from_edge=np.round(sinogram_cut.shape[1]/2.).astype(int)-max_rad                          
    
    sinogram_cut=sinogram_cut[:,dist_from_edge:-dist_from_edge]
    
    print('Inverting Sinogram....')
    reconstruction_fbp = iradon(sinogram_cut.T, theta=theta, circle=True)
    
    reconstruction_fbp=np.rot90(reconstruction_fbp,3)#Rotation to get the result consistent with hexrd, needs to be checked
    
    return reconstruction_fbp

    def reconstruct(self, sinogram, centre_of_rotations, angles, vol_shape):
        in_pData = self.get_plugin_in_datasets()[0]
        sinogram = np.swapaxes(sinogram, 0, 1)
        sinogram = self._shift(sinogram, centre_of_rotations)
        sino = sinogram.astype(np.float64)
        theta = np.linspace(0, 180, sinogram.shape[1])
        result = transform.iradon(
            sino,
            theta=theta,
            output_size=(in_pData.get_shape()[1]),
            # self.parameters['output_size'],
            filter="ramp",  # self.parameters['filter'],
            interpolation="linear",
            # self.parameters['linear'],
            circle=False,
        )  # self.parameters[False])

        for i in range(self.parameters["iterations"]):
            print "Iteration %i" % i
            result = transform.iradon_sart(
                sino,
                theta=theta,
                image=result,
                # self.parameters['result'],
                projection_shifts=None,
                # self.parameters['None'],
                clip=None,
                # self.parameters[None],
                relaxation=0.15
                # self.parameters[0.15])
            )
        return result

 def reconstruct(self, sinogram, centre_of_rotation,
                 angles, shape, center):
     print sinogram.shape
     sinogram = np.swapaxes(sinogram, 0, 1)
     sinogram = self._shift(sinogram, centre_of_rotation)
     sino = np.nan_to_num(sinogram)
     theta = np.linspace(0, 180, sinogram.shape[1])
     result = \
         transform.iradon(sino, theta=theta,
                          output_size=(sinogram.shape[0]),
                          # self.parameters['output_size'],
                          filter='ramp',  # self.parameters['filter'],
                          interpolation='linear',
                          # self.parameters['linear'],
                          circle=False)  # self.parameters[False])
     for i in range(self.parameters["iterations"]):
         print "Iteration %i" % i
         result = transform.iradon_sart(sino, theta=theta, image=result,
                                        # self.parameters['result'],
                                        projection_shifts=None,
                                        # self.parameters['None'],
                                        clip=None,
                                        # self.parameters[None],
                                        relaxation=0.15
                                        # self.parameters[0.15])
                                        )
     return result

def centreshift(data,angles,tiltrange,increment):
    series_shape = np.shape(data)
    new_series = data.copy()
    midy = int(series_shape[1]/2)
    for i in range(series_shape[0]):
        new_series[i,:,:] = interpolation.shift(data.copy()[i,:,:],(0,axis_shift))
        
    
    iradon(np.rot90(data[:,midy,:]), theta = angles,output_size = series_shape[2])
        
    
#if __name__ == "__main__":
    #x = hspy.load('C:/Users/Tom/Documents/TEM data/20150521_SS316needle/EDX Tomo/Cr/Signed_zero/cr_sub_hdr0_2.ali')
    '''x = h5py.File('C:/Users/Tom/Documents/TEM data/20150521_SS316needle/EDX Tomo/Cr/Signed_zero/cr_sub_hdr0_2_bin8.h5', "r")
    tilt_angles = np.linspace(-90.,90.,37)
    x_data = x['/0']
    reconstruction = reconstruct(x_data,tilt_angles)'''
    '''

def PD_denoise(sinogram,rec_FBP,thetas,outerloops,innerloops,
                           CGloops,image_res,convtestnorm,lambdapen,slicenum):

   """
   Solves the problem via the following formulation.  We minimize F(Ku)+G(u),
   where K:= (grad,R)^T and F(x):= ||x_1||_1+lambda/2||x_2-g||^2, and G(u):=0.
   Then, F^*(y) = \delta_p(y_1)+1/(2lambda)||y_2||^2+<g,y_2>.
   We then perform the primal dual algorithm of Chambolle-Pock '11.
   We assume it is better to have repeated G evaluations and fewer thing stored.
   Do not use; needs significant tuning/possible debugging.
   """
   uk = (TV_Split_utilities.generateSheppLogan(Phant_size))
   radius = min(uk.shape) // 2
   c0, c1 = np.ogrid[0:uk.shape[0], 0:uk.shape[1]]
   Circle_Mask = ((c0 - uk.shape[0] // 2) ** 2+ (c1 - uk.shape[1] // 2) ** 2)
   Circle_Mask = Circle_Mask <= 0.95*(radius ** 2)
   Circle_Mask = np.matrix(Circle_Mask)
   uk[0==Circle_Mask]=0.
   y_onek = np.gradient(0.*uk)
   y_twok = 0.*sinogram
   sigma = 1.e-2
   tau = 1.e-2
   theta = 1.
   ubark = uk
   Err_diff = np.zeros((int(.5*rec_FBP.size)))
   for i in range(int(.5*rec_FBP.size)):
        y_onek[0],y_onek[1], y_twok = resolvent_Fstar(y_onek[0]+sigma*
                  grad_one(ubark),y_onek[1]+sigma*grad_two(ubark),
                  y_twok+sigma*radon(ubark,thetas,circle=True),sigma,lambdapen,
                  sinogram)
        ubark = (1+theta)*resolvent_G(uk+tau*(-1.*Divu(y_onek)+iradon(y_twok,
                 thetas,output_size=image_res,circle=True,filter = 'ramp' )))-theta*uk
        ubark[0==Circle_Mask]=0.
        uk = resolvent_G(uk-tau*(-1.*Divu(y_onek)+iradon(y_twok,
                 thetas,output_size=image_res,circle=True,filter = 'ramp' )))
        uk[0==Circle_Mask]=0.
        Err_diff[i] = np.linalg.norm((ubark-uk)/theta)
        plt.imshow(uk);plt.pause(.0001)
        print('Slice %d]' %(slicenum))
        print('Outer Iterate = ' ,i)
        print('Update = ' ,Err_diff[i])
        if Err_diff[i]<1.e-5:
            break
   return uk,Err_diff

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 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 inverAbel(x,y):
	x = numpy.array(x)
	y = list(y)
	dx = abs(x[1]-x[0])
	length = len(y)
	#~ print len(numpy.array(y*100))
	sino = numpy.array(y*100).reshape(100,length).transpose()
	reconstruction = iradon(sino)
	iry = reconstruction[:,math.ceil(len(reconstruction)*0.5)]/dx
	cropxs = int((length - len(iry))*0.5)-1
	cropxe = cropxs + len(iry)
	irx = x[cropxs:cropxe]
	#~ print len(irx),len(iry)
	return  irx,iry

def reconstruct(tilt_data,tilt_angles,algorithm='ASTRA_SIRT',iterations=20,geometry='parallel3d'):
    '''Function to reconstruct a tilt series.
    
    Data should be in the format (Angles,X,Y), where the tilt axis is situated about the mid column of the data.'''
    
    t0 = time.time()
    data_shape = np.shape(tilt_data)
    y_size = data_shape[1]
    
    recon = np.zeros((data_shape[2],y_size,data_shape[2]))
    #sino = np.zeros(data_shape)
    
    #Reconstruction using Filtered/Weighted Backprojection from skimage (check how filtering is done)
    if algorithm == 'SKI_FBP' or 'SKI_WBP':
        for y in range(0,y_size-1):
            recon[:,y,:] = iradon(np.rot90(tilt_data[:,y,:]), theta = tilt_angles,output_size = data_shape[2])
            
            #This is supposed to reorder axis to orientation for projection but not sure the x and y are correct
            for y in range(0,y_size-1):
                recon[:,y,:] = np.rot90(recon[:,y,:])
                recon[:,y,:] = np.rot90(recon[:,y,:])
                #recon[:,y,:] = np.rot90(recon[:,y,:])
            
    if algorithm == 'SKI_SART':
        for y in range(0,y_size-1):
            recon[:,y,:] = iradon_sart(np.rot90(tilt_data[:,y,:]), theta = tilt_angles,clip=(0,np.max(tilt_data)))
            
        if iterations > 1:
            for it in range(iterations-1):
                print("Iteration number "+str(it+2)+" in progress.")
                for y in range(0,y_size-1):
                    recon[:,y,:] = iradon_sart(np.rot90(tilt_data[:,y,:]), theta = tilt_angles,image=recon[:,y,:],clip=(0,np.max(tilt_data)))
                    
        #This is supposed to reorder axis to orientation for projection but not sure the x and y are correct
        for y in range(0,y_size-1):
            recon[:,y,:] = np.rot90(recon[:,y,:])
            recon[:,y,:] = np.rot90(recon[:,y,:])
            #recon[:,y,:] = np.rot90(recon[:,y,:])
                    
    if algorithm == 'ASTRA_SIRT':
        recon = astrarecon(tilt_data,tilt_angles,iterations,geometry)
      
        '''for z in xrange(0,data_shape[2]):
            recon[:,:,z] = np.rot90(recon[:,:,z])
            recon[:,:,z] = np.rot90(recon[:,:,z])
            recon[:,:,z] = np.rot90(recon[:,:,z])'''
        
    print("Reconstruction completed in {} seconds.".format(time.time() - t0))
            
    return(recon)

 def reconstruct(self, sino, centre_of_rotations, angles, vol_shape, init):
     in_pData = self.get_plugin_in_datasets()[0]
     in_meta_data = self.get_in_meta_data()[0]
     sinogram = np.swapaxes(sino, 0, 1)
     sinogram = self._shift(sinogram, centre_of_rotations)
     theta = in_meta_data.get_meta_data('rotation_angle')
     result = \
         transform.iradon(sinogram, theta=theta,
                          output_size=(in_pData.get_shape()[1]),
                          # self.parameters['output_size'],
                          filter='ramp',  # self.parameters['filter'],
                          interpolation='linear',
                          # self.parameters['linear'],
                          circle=False)  # self.parameters[False])
     return result

def iridgelet(ridge):
    lvl, nr,nt = np.shape(ridge)
 #   for l in np.linspace(0,lvl-1,lvl):
#        x = ski.radon(ridge[l,:,:],theta = np.linspace(0,180-1,nt))
#        if l ==0:
#            nr,nt = np.shape(x)
#            rad_ridgelet = np.zeros((lvl,nr,nt))
#        rad_ridgelet[l,:,:] = x
    rad_ridgelet = ridge
    rad = np.zeros((nr,nt))
        
    for i in np.linspace(0,nt-1,nt):
        rad[:,i] = iuwt_1D(rad_ridgelet[:,:,i])
    img = ski.iradon(rad)
    return img

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 reconstruct(self, sinogram, centre_of_rotation,
                 angles, shape, center):
     print sinogram.shape
     sinogram = np.swapaxes(sinogram, 0, 1)
     sinogram = self._shift(sinogram, centre_of_rotation)
     sino = np.nan_to_num(sinogram)
     theta = np.linspace(0, 180, sinogram.shape[1])
     result = \
         transform.iradon(sino, theta=theta,
                          output_size=(sinogram.shape[0]),
                          # self.parameters['output_size'],
                          filter='ramp',  # self.parameters['filter'],
                          interpolation='linear',
                          # self.parameters['linear'],
                          circle=False)  # self.parameters[False])
     return result

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

def test_iradon_bias_circular_phantom():
    """
    test that a uniform circular phantom has a small reconstruction bias
    """
    pixels = 128
    xy = np.arange(-pixels / 2, pixels / 2) + 0.5
    x, y = np.meshgrid(xy, xy)
    image = x**2 + y**2 <= (pixels/4)**2

    theta = np.linspace(0., 180., max(image.shape), endpoint=False)
    sinogram = radon(image, theta=theta)

    reconstruction_fbp = iradon(sinogram, theta=theta)
    error = reconstruction_fbp - image
    
    tol = 5e-5
    roi_err = np.abs(np.mean(error))
    assert( roi_err < tol )