text: import numpy as np
  from tqdm import trange
  """
  Created by Mikhaylov Andrey. 
  Last update 04.02.2022.
  
  andrey.mikhaylov@kit.edu
  
  
  
  """
  
  __version__ = '1.01_04.02.2022'
  
  
  
  
  
  
  
  
  
  
  def custom_sino_filter(sinogram,ramp = False):
      """
      Hilber or Hilbert+Ramp filtering of DPC sinogram
  
      Parameters
      ----------
      sinogram : 2D numpy array
          Your sinogram. Index 0 - angles, index 1 - detector size, horizontal.
      ramp : bool, optional
          If you want ramp filtering after Hilbert. The default is False.
  
      Returns
      -------
      2D numpy array
          Filtered sinogram.
  
      """
      new_sino = np.zeros_like(sinogram)    
      n_angles,ncolums = sinogram.shape
  
      'Filter creation'
      custom_filter_space = np.arange(-ncolums/2,ncolums/2,1)
      custom_filter = np.sign(custom_filter_space)/(2*np.pi*1j)
      custom_filter[0] = 0
      custom_filter[-1] = 0
      
      'Filterng'
      for i in range(n_angles):
          new_sino[i] = np.real(np.fft.ifft(np.fft.fft(sinogram[i]) * custom_filter))    
      'Ramp filtering'
      if ramp:
          ramp_filter = np.abs(custom_filter_space)
          ramp_filter[0] = 0
          ramp_filter[-1] = 0
          new_new_sino = np.zeros_like(sinogram)         
          for i in range(n_angles):
              new_new_sino[i] = np.real(np.fft.ifft(np.fft.fft(new_sino[i]) * ramp_filter))
          
          return new_new_sino
          
      else:
          return new_sino
  
  def ful_rec_phase_hor(stack,ramp = False,pad = 100,angle = 180):
      """
      Full tomography reconstruction from the DPC data
  
      Parameters
      ----------
      stack : 3D numpy array
          Projection stack. Index 0 - angles, index 1 - nr. of slices/vertical detector size, index 2 - detector size, horizontal.
      ramp : bool, optional
          f you want ramp filtering after Hilbert. The default is False.
      pad : int, optional
          Padding size. Will add this amount of pixels during rotation. The default is 100.
      angle : int, optional
          biggest angle value in deg. The default is 180.
  
      Returns
      -------
      stack_final : 3D numpy array
          reconstructed volume.
  
      """
      import skimage.transform as tra
      n_angles,ncolums,nrow = stack.shape
      ang_step = 180/n_angles
     
      stack_filtered = np.zeros_like(stack)
      
      for i in range(ncolums):
          single_slice = stack[:,i,:]
          stack_filtered[:,i,:] = custom_sino_filter(single_slice,ramp)
          
      stack_final = np.zeros((nrow,nrow,ncolums))        
      
      for horisontal_slice in trange(ncolums):
          
          final_sino = stack_filtered[:,horisontal_slice,:]
          
          final_slice = np.zeros((nrow+2*pad,nrow+2*pad))
          
          for angle in range(n_angles):
              angle_n = ang_step*angle
  
              a = np.pad(np.tile(final_sino[angle],(nrow,1)),(pad,pad))
              b = tra.rotate(a,angle_n)
              final_slice += b
          
          stack_final[:,:,horisontal_slice] = final_slice[pad:nrow+pad,pad:nrow+pad]
      
      return stack_final
  
  
  
  "Example"
  
  
  if __name__ == "__main__":
      from skimage.transform import radon,rescale
      from skimage.data import shepp_logan_phantom
      import matplotlib.pyplot as plt
      
      image = shepp_logan_phantom()
      image = rescale(image, scale=0.3, mode='reflect')
      theta = np.linspace(0., 180., 360, endpoint=False)
      sinogram = radon(image, theta=theta)
      
      
      stack = np.tile(sinogram,(2,1,1))
      stack = np.swapaxes(stack, 0, 2)
      stack = np.swapaxes(stack, 2, 1)
      
      plt.imshow(stack[:,0,:])
      plt.title('initial sinogram')
      plt.colorbar()
      plt.show()        
      
      phase_stack = np.zeros_like(stack)
      
             
      for i in trange(360):
          phase_stack[i] = np.gradient(stack[i],axis=1)
      
      plt.imshow(phase_stack[:,1,:])
      plt.title('DPC sinogram')
      plt.colorbar()
      plt.show()        
      
      
      your_reconstruction = ful_rec_phase_hor(phase_stack)
      
      plt.imshow(your_reconstruction[:,:,0])
      plt.colorbar()
      plt.show()  
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      