Welcome to PySABER’s documentation!

Introduction

PySABER is a python package for characterizing the X-ray source and detector blurs in cone-beam X-ray imaging systems. SABER is an abbreviation for systems approach to blur estimation and reduction. Note that even parallel beam X-rays in synchrotrons are in fact cone beams albeit with a large source to object distance (SOD). X-ray images, also called radiographs, are simultaneously blurred by both the X-ray source spot blur and detector blur. This python package uses a numerical optimization algorithm to disentangle and estimate both forms of blur simultaneously. The point spread function (PSF) of X-ray source blur is modeled using a density function with two parameters. The first parameter is the full width half maximum (FWHM) of the PSF along the x-axis (row-wise) and second is the FWHM along the y-axis (column-axis). The PSF of detector blur is modeled as the sum of two density functions, each with its own FWHM parameter, that are mixed together by a mixture (or weight) parameter. All these parameters are then estimated using numerical optimization from normalized radiographs of a sharp edge such as a thick Tungsten plate rollbar. To simultaneously estimate the PSFs of both source and detector blurs, radiographs must be acquired at two different values for the ratio of the source to object distance (SOD) and object to detector distance (ODD). If each radiograph has a single straight edge, then the measurements must be repeated for two different, preferably perpendicular, orientations of the edge. If the radiograph consists of two intersecting perpendicular edges, then a single radiograph at each specified SOD/ODD is sufficient.

Once the parameters of both source and detector blurs are estimated, this package is also useful to reduce blur in radiographs using deblurring algorithms. Currently, Wiener filtering and regularized least squares deconvolution are two deblurring algorithms that are supported for deblurring. Both these techniques use the estimated blur parameters to deblur radiographs. For more detailed explanation on the experimental methodology and theory of PySABER, please read the paper in References.

References

If you use pysaber Package, we request you to cite the following article -

  • K. Aditya Mohan, Robert M. Panas, and Jefferson A. Cuadra. “SABER: A Systems Approach to Blur Estimation and Reduction in X-ray Imaging.” arXiv preprint arXiv:1905.03935 (2019) [pdf]

License

LLNL-CODE-766837. This project is licensed under the MIT License.

MIT License

Copyright (c) 2018, Lawrence Livermore National Security, LLC

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

Installation

pysaber is installed using the python package manager pip. To install pysaber, run the following command in a terminal:

pip install pysaber

Alternatively, to install using the source code in the github repository pysaber, first download the repository using the download link in the top right corner of the webpage. Or, you can also git clone the repository directly from github. In a terminal, change the current directory to the outermost folder of this downloaded repository, which contains the README file, and run the following command:

pip install .

It is recommended to install the pysaber Package within a python virtual environment.

Tutorial

The steps involved in estimating blur PSFs and deblurring radiographs are outlined in the links below.

Input Sanity Check

  • The first step to estimating the blur model involves computation of the transmission function, which is the ideal radiograph image that is formed in the absence of X-ray source and detector blurs. This computation of transmission function is performed internally in the function pysaber.estimate_blur(), which is used to estimate the blur model by computing the parameters of X-ray source and detector blurs. However, this computation of transmission function is not fail-proof and may result in inaccurate edge localization if certain assumptions made when computing the transmission function are not satisfied.

  • Before using pysaber.estimate_blur() to estimate blur model, it is recommended to check for accurate edge localization in the transmission function. The transmission function can be computed using the function pysaber.get_trans_masks().

  • The function pysaber.get_trans_masks() also returns the mask arrays for the transmission function and radiograph, which are used to include or exclude certain pixels from blur estimation. By default, the radiograph mask only excludes a small number of pixels along the boundary of the radiograph from blur estimation. Additional pixels can be excluded from blur estimation by appropriately setting the input arguments of the functions pysaber.get_trans_masks() and pysaber.estimate_blur(). The mask for transmission function should also exclude the padded pixels in addition to those pixels excluded by the radiograph mask. Hence, pysaber.get_trans_masks() is also useful to check if user expectations for the mask arrays are satisfied.

  • Example python scripts that demonstrate the above procedure are shown below. To obtain the data that is required to run this script, download and unzip the zip file at the link data. To run the script as is within the current working directory, the files in the zip file must be placed within a folder called data.

Verify transmission function and masks for horizontal edge radiograph.
import numpy as np
from PIL import Image #To read images in TIFF format
from pysaber import get_trans_masks #To compute transmission function and masks
import matplotlib.pyplot as plt #To display images
    
pix_wid = 0.675 #Width of each pixel in micrometers

#Read a horizontal edge radiograph for which the transmission function must be computed
rad = Image.open('data/horz_edge_25mm.tif') #Read radiograph
rad = np.asarray(rad) #Convert to numpy array
bright = Image.open('data/horz_bright.tif') #Read bright field
bright = np.asarray(bright) #Convert to numpy array
dark = Image.open('data/horz_dark.tif') #Read dark field
dark = np.asarray(dark) #Convert to numpy array
nrad = (rad-dark)/(bright-dark) #Normalize radiograph

#Get the transmission function and masks 
trans,trans_mask,rad_mask = get_trans_masks(nrad,edge='straight',pad=[3,3])
#Use pad = [1,1] if you do not want padding

#Display the array trans
#Visually check for inaccurate localization of the sharp-edge
plt.imshow(trans,cmap='gray')
plt.colorbar()
plt.show()

#Display and inspect the mask for transmission function
plt.imshow(trans_mask,cmap='gray')
plt.show()

#Display and inspect the mask for radiograph
plt.imshow(rad_mask,cmap='gray')
plt.show()

#Show a line plot comparing the measured radiograph and transmission function
sz = nrad.shape
coords = np.arange(-(sz[0]//2),sz[0]//2,1)*pix_wid
mid = (trans.shape[0]//2,trans.shape[1]//2)

plt.plot(coords,nrad[:,sz[1]//2])
#Due to padding, trans is three times the size of nrad in each dimension
#For proper alignment in the presence of padding, both nrad and trans are center aligned
#Center alignment is used since an equal amount of padding is applied at both ends of each axis
plt.plot(coords,trans[mid[0]-(sz[0]//2):mid[0]+(sz[0]//2),mid[1]])
plt.xlabel('micrometers')
plt.legend(['Measured','Transmission'])
plt.show()
Verify transmission function and masks for vertical edge radiograph.
import numpy as np
from PIL import Image #To read images in TIFF format
from pysaber import get_trans_masks #To compute transmission function and masks
import matplotlib.pyplot as plt #To display images
    
pix_wid = 0.675 #Width of each pixel in micrometers

#Read a vertical edge radiograph for which the transmission function must be computed
rad = Image.open('data/vert_edge_25mm.tif') #Read radiograph
rad = np.asarray(rad) #Convert to numpy array 
bright = Image.open('data/vert_bright.tif') #Read bright field
bright = np.asarray(bright) #Convert to numpy array
dark = Image.open('data/vert_dark.tif') #Read dark field
dark = np.asarray(dark) #Convert to numpy array
nrad = (rad-dark)/(bright-dark) #Normalize radiograph

#Get the transmission function and masks 
trans,trans_mask,rad_mask = get_trans_masks(nrad,edge='straight',pad=[3,3])
#Use pad = [1,1] if you do not want padding

#Display the array trans
#Visually check for inaccurate localization of the sharp-edge
plt.imshow(trans,cmap='gray')
plt.colorbar()
plt.show()

#Display and inspect the mask for transmission function
plt.imshow(trans_mask,cmap='gray')
plt.show()

#Display and inspect the mask for radiograph
plt.imshow(rad_mask,cmap='gray')
plt.show()

#Show a line plot comparing the measured radiograph and transmission function
sz = nrad.shape
coords = np.arange(-(sz[1]//2),sz[1]//2,1)*pix_wid
mid = (trans.shape[0]//2,trans.shape[1]//2)

plt.plot(coords,nrad[sz[0]//2,:])
#Due to padding, trans is three times the size of nrad in each dimension
#For proper alignment in the presence of padding, both nrad and trans are center aligned
#Center alignment is used since an equal amount of padding is applied at both ends of each axis
plt.plot(coords,trans[mid[0],mid[1]-(sz[1]//2):mid[1]+(sz[1]//2)])
plt.xlabel('micrometers')
plt.legend(['Measured','Transmission'])
plt.show()

Estimate Blur Model

  • To estimate the blur model, we must acquire radiographs of a sharp edge such as a Tungsten rollbar. Each sharp edge radiograph can either contain a single straight edge or two mutually perpendicular intersecting edges. If imaging a single straight edge, then radiographs must be acquired at two different perpendicular orientations of the straight edge. Also, radiographs must be acquired at two different values of SOD/ODD, where SOD is the source to object distance and ODD is the object to detector distance.

  • Next, the radiographs must be appropriately normalized. For each radiograph, acquire a bright field image (measurements with X-rays but no sample) and a dark field image (measurements without X-rays). Then, compute the normalized radiograph by dividing the difference between the radiograph and the dark field image with the difference between the bright field and the dark field image.

  • Using the normalized radiographs, estimate parameters of X-ray source blur and detector blur using the function pysaber.estimate_blur().

  • An example python script that demonstrates blur estimation using radiographs of a single straight edge at various orientations and SOD/ODD values is shown below. To obtain the data that is required to run this script, download and unzip the zip file at the link data. To run the script as is within the current working directory, the files in the zip file must be placed within a folder called data.

import numpy as np #For mathematics on vectors
from PIL import Image #To read images in TIFF format
from pysaber import estimate_blur #To estimate blur PSF parameters

pix_wid = 0.675 #Width of each pixel in micrometers

#Horizontal and vertical edge radiographs
edge_files = ['data/horz_edge_25mm.tif','data/horz_edge_50mm.tif',
                'data/vert_edge_25mm.tif','data/vert_edge_50mm.tif'] 
#Filenames of bright field images for normalization
bright_files = ['data/horz_bright.tif','data/horz_bright.tif',
                    'data/vert_bright.tif','data/vert_bright.tif']
#Filenames of dark field images for normalization
dark_files = ['data/horz_dark.tif','data/horz_dark.tif',
                'data/vert_dark.tif','data/vert_dark.tif']
#Source to object (SOD) distances in micrometers for each radiograph in edge_files
sod = [24751.89,50251.79,24753.05,50253.35] 
#Source to detector (SDD) distances in micrometers for each radiograph in edge_files
sdd = [71003.08,71003.08,71010.86,71010.86] 

rads = [] #List that will contain normalized radiographs
odd = [] #Object to detector distance (ODD) for each radiograph in rads
for i in range(len(edge_files)): #Loop through all the radiograph files
    rad = Image.open(edge_files[i]) #Read radiograph
    rad = np.asarray(rad) #Convert to numpy array
    bright = Image.open(bright_files[i]) #Read bright field
    bright = np.asarray(bright) #Convert to numpy array
    dark = Image.open(dark_files[i]) #Read dark field
    dark = np.asarray(dark) #Convert to numpy array
    nrad = (rad-dark)/(bright-dark) #Normalize radiograph
    rads.append(nrad) #Add normalized radiograph to the list rads
    odd.append(sdd[i]-sod[i]) #Add corresponding ODD to the list odd

#Estimate X-ray source blur, detector blur, and every radiograph's transmission function
#To reduce run time and quickly produce a result, the value for argument thresh can be reduced.
#However, reducing thresh may produce an inaccurate blur model that does not fit the measured data  
src_params,det_params,trans_params = estimate_blur(rads,sod,odd,pix_wid,
                edge='straight',thresh=1e-6,pad=[3,3],power=1.0,save_dir='./')
#src_params is a python dictionary of parameters that quantify X-ray source blur
#det_params is a python dictionary of parameters that quantify blur from the detector panel
#Both src_params and det_params characterize the effective blur and are used during deblurring 
#trans_params is a list of lists, each of which contains the low/high values of transmission function
#trans_params is useful to check accuracy of fit and not useful for deblurring. 

#help(estimate_blur)
#Uncomment above line to get help on using the function estimate_blur

print("-------------------------------------------------------")
print("Source blur model parameters are {}".format(src_params)) 
#Print parameters of source blur
print("Detector blur model parameters are {}".format(det_params)) 
#Print parameters of detector blur
print("Transmission function parameters are {}".format(trans_params)) 
#Print parameters of transmission functions

Validate Blur Model

  • We must ensure that the estimated parameters are indeed a good fit for the measured data. This is done by comparing line profiles across the sharp edge between the measured radiograph and the predicted radiograph from the blur model. The output of the blur model given parameters of source blur, detector blur, and transmission function is computed using the function pysaber.get_trans_fit(). The fit must be evaluated for every sharp-edge radiograph that is input to pysaber.estimate_blur().

  • Verify the agreement between the measured radiograph and the blur model prediction. Carefully zoom into the region containing the sharp edge and verify if the predicted blur matches with the blur in the measured radiograph. Also, verify the agreement between the measured radiograph values and blur model prediction in regions further away from the sharp edge. The predicted radiograph that is output by the blur model contains additional padding. Hence, it is necessary to account for this padding when comparing with the measured radiograph.

  • If the fit is not tight, consider reducing the value of the input argument thresh of the function pysaber.estimate_blur() to obtain a better fit. Reducing the convergence threshold, thresh, can improve the agreement between the measured radiograph and the blur model prediction, but will inevitably result in longer run times. A good fit indicates that the blur model is able to accurately model the X-ray source and detector blurs.

  • Example python scripts for line profile comparisons between the blur model prediction and measured radiograph are shown below. To obtain the data that is required to run this script, download and unzip the zip file at the link data. To run the script as is within the current working directory, the files in the zip file must be placed within a folder called data.

Line profile comparisons across a horizontal edge radiograph.
import numpy as np
from PIL import Image #To read images in TIFF format
from pysaber import get_trans_fit #To get blurred radiograph as predicted by the blur model 
import matplotlib.pyplot as plt #To display images
    
pix_wid = 0.675 #Width of each pixel in micrometers

#Read a horizontal edge radiograph for which accuracy of fit must be analyzed
rad = Image.open('data/horz_edge_25mm.tif') #Read radiograph
rad = np.asarray(rad) #Convert to numpy array
bright = Image.open('data/horz_bright.tif') #Read bright field
bright = np.asarray(bright) #Convert to numpy array
dark = Image.open('data/horz_dark.tif') #Read dark field
dark = np.asarray(dark) #Convert to numpy array
nrad = (rad-dark)/(bright-dark) #Normalize radiograph
sod = 24751.89 #Source to object distance (SOD) of radiograph
sdd = 71003.08 #Source to detector distance (SDD) of radiograph

#Parameters of X-ray source blur
src_params = {'source_FWHM_x_axis':2.69,
                'source_FWHM_y_axis':3.01, 
                'norm_power':1.0,
                'cutoff_FWHM_multiplier':10}
#Parameters of detector blur
det_params = {'detector_FWHM_1':1.85, 
                'detector_FWHM_2':126.5, 
                'detector_weight_1':0.916, 
                'norm_power':1.0, 
                'cutoff_FWHM_1_multiplier':10, 
                'cutoff_FWHM_2_multiplier':10}
#Transmission function parameters
trans_params = [0.015,0.98]

#Get the blurred radiograph as predicted by the blur model
pred_nrad,_ = get_trans_fit(nrad,sod,sdd-sod,pix_wid,src_params,det_params,trans_params,edge='straight',pad=[3,3])

#Show a line plot comparing the measured radiograph and the predicted blurred radiograph
sz = nrad.shape
coords = np.arange(-(sz[0]//2),sz[0]//2,1)*pix_wid
mid = (pred_nrad.shape[0]//2,pred_nrad.shape[1]//2)

plt.plot(coords,nrad[:,sz[1]//2])
#Due to padding, pred_nrad is three times the size of nrad in each dimension
#For proper alignment in the presence of padding, both nrad and pred_nrad are center aligned
#Center alignment is used since an equal amount of padding is applied at both ends of each axis
plt.plot(coords,pred_nrad[mid[0]-(sz[0]//2):mid[0]+(sz[0]//2),mid[1]])
plt.xlabel('micrometers')
plt.legend(['Measured','Prediction'])
plt.show()
Line profile comparisons across a vertical edge radiograph.
import numpy as np
from PIL import Image #To read images in TIFF format
from pysaber import get_trans_fit #To get blurred radiograph as predicted by the blur model 
import matplotlib.pyplot as plt #To display images
    
pix_wid = 0.675 #Width of each pixel in micrometers

#Read a vertical edge radiograph for which accuracy of fit must be analyzed
rad = Image.open('data/vert_edge_25mm.tif') #Read radiograph
rad = np.asarray(rad) #Convert to numpy array 
bright = Image.open('data/vert_bright.tif') #Read bright field
bright = np.asarray(bright) #Convert to numpy array
dark = Image.open('data/vert_dark.tif') #Read dark field
dark = np.asarray(dark) #Convert to numpy array
nrad = (rad-dark)/(bright-dark) #Normalize radiograph
sod = 24753.05 #Source to object distance (SOD) of radiograph
sdd = 71010.86 #Source to detector distance (SDD) of radiograph

#Parameters of X-ray source blur
src_params = {'source_FWHM_x_axis':2.69,
                'source_FWHM_y_axis':3.01, 
                'norm_power':1.0,
                'cutoff_FWHM_multiplier':10}
#Parameters of detector blur
det_params = {'detector_FWHM_1':1.85, 
                'detector_FWHM_2':126.5, 
                'detector_weight_1':0.916, 
                'norm_power':1.0, 
                'cutoff_FWHM_1_multiplier':10, 
                'cutoff_FWHM_2_multiplier':10}
#Transmission function parameters
trans_params = [0.015,0.98]

#Get the blurred radiograph as predicted by the blur model
pred_nrad,_ = get_trans_fit(nrad,sod,sdd-sod,pix_wid,src_params,det_params,trans_params,edge='straight',pad=[3,3])

#Show a line plot comparing the measured radiograph and the predicted blurred radiograph
sz = nrad.shape
coords = np.arange(-(sz[1]//2),sz[1]//2,1)*pix_wid
mid = (pred_nrad.shape[0]//2,pred_nrad.shape[1]//2)

plt.plot(coords,nrad[sz[0]//2,:])
#Due to padding, pred_nrad is three times the size of nrad in each dimension
#For proper alignment in the presence of padding, both nrad and pred_nrad are center aligned
#Center alignment is used since an equal amount of padding is applied at both ends of each axis
plt.plot(coords,pred_nrad[mid[0],mid[1]-(sz[1]//2):mid[1]+(sz[1]//2)])
plt.xlabel('micrometers')
plt.legend(['Measured','Prediction'])
plt.show()

Visualize Blur PSF

  • For further analysis and visualization, we can also compute the point spread functions (PSF) of source blur and detector blur. The PSF of source blur is computed using the function pysaber.get_source_psf() and PSF of detector blur is computed using pysaber.get_detector_psf().

  • The function pysaber.get_source_psf() is useful to compute PSF either in the plane of the X-ray source or the plane of the detector. Since source blur PSF on the detector plane is a function of the object’s source to object distance (SOD) and object to detector distance (ODD), SOD and ODD must be specified when computing source PSF in the plane of the detector. To compute source blur PSF in the source plane, it is sufficient to use the default values for SOD and ODD in pysaber.get_source_psf().

  • The detector blur PSF obtained using pysaber.get_detector_psf() models blur due to the scintillator and detector panel. Hence, it is independent of SOD and ODD.

  • Example python scripts that demonstrate visualization of source and detector PSFs are shown below.

Plot X-ray source blur PSF
import numpy as np #For mathematics on vectors
import matplotlib.pyplot as plt #For plotting and showing images
from pysaber import get_source_psf #To compute PSF of source blur 

pix_wid = 0.675 #Width of each pixel in micrometers
#Parameters of X-ray source blur
src_params = {'source_FWHM_x_axis':2.69,
                'source_FWHM_y_axis':3.01, 
                'norm_power':1.0,
                'cutoff_FWHM_multiplier':10}

#Get point spread function (PSF) of source blur in the plane of the X-ray source.
#Do not supply SOD and SDD if you need PSF in the source plane.
source_psf = get_source_psf(pix_wid,src_params)

#Display the source blur PSF on the source plane as an image
sz = source_psf.shape
x = np.arange(-(sz[1]//2),(sz[1]//2)+1,1)*pix_wid
y = np.arange(-(sz[0]//2),(sz[0]//2)+1,1)*pix_wid
plt.pcolormesh(x,y,source_psf,cmap='gray')
plt.xlabel('micrometers')
plt.ylabel('micrometers')
plt.title('X-ray source PSF at source plane')
plt.colorbar()
plt.show()

#To get source blur PSF on the detector plane, supply SOD and ODD = SDD - SOD to the function get_source_psf
sod = 25000 #in micrometers
sdd = 71000 #in micrometers
source_psf = get_source_psf(pix_wid,src_params,sod,sdd-sod)

#help(get_source_psf)
#Uncomment the above line to get help on using the function get_source_psf

#Display the source blur PSF on the detector plane as an image
sz = source_psf.shape
x = np.arange(-(sz[1]//2),(sz[1]//2)+1,1)*pix_wid
y = np.arange(-(sz[0]//2),(sz[0]//2)+1,1)*pix_wid
plt.pcolormesh(x,y,source_psf,cmap='gray')
plt.xlabel('micrometers')
plt.ylabel('micrometers')
plt.title('X-ray source PSF at detector plane')
plt.colorbar()
plt.show()
Plot X-ray detector blur PSF
import numpy as np #For mathematics on vectors
import matplotlib.pyplot as plt #For displaying images
from matplotlib.colors import LogNorm #To display image values in logarithm scale 
from pysaber import get_detector_psf #To compute PSF of detector blur

pix_wid = 0.675 #Width of each pixel in micrometers
#Parameters of detector blur
det_params = {'detector_FWHM_1':1.85, 
                'detector_FWHM_2':126.5, 
                'detector_weight_1':0.916, 
                'norm_power':1.0, 
                'cutoff_FWHM_1_multiplier':10, 
                'cutoff_FWHM_2_multiplier':10}

#Get point spread function (PSF) of detector blur as a 2D numpy array.
detector_psf = get_detector_psf(pix_wid,det_params)

#help(get_detector_psf)
#Uncomment the above line to get help in using the function get_detector_psf

#Display the PSF of detector blur
sz = detector_psf.shape
x = np.arange(-(sz[1]//2),(sz[1]//2)+1,1)*pix_wid
y = np.arange(-(sz[0]//2),(sz[0]//2)+1,1)*pix_wid
plt.pcolormesh(x,y,detector_psf,cmap='gray',norm=LogNorm())
plt.xlabel('micrometers')
plt.ylabel('micrometers')
plt.title('Detector PSF')
plt.colorbar()
plt.show()

Deblur Radiographs

  • Once the parameters of source and detector PSFs are estimated, radiographs of any arbitrary sample acquired at any source to object distance (SOD) and source to detector distance (SDD) can be deblurred using various techniques.

  • To deblur a radiograph using Wiener filter, the function pysaber.wiener_deblur() is used. To deblur using regularized least squares deconvolution (RLSD), use the function pysaber.least_squares_deblur().

  • Deblurring increases sharpness and resolution. However, it also introduces ringing artifacts and increases noise. To reduce noise and ringing artifacts, the regularization parameter can be increased. Ringing artifacts also increase with increasing inaccuracy of the blur model. Thus, it is essential to obtain a good fit between the measured radiograph and the blur model prediction as explained in the section Validate Blur Model.

  • The python scripts shown below demonstrate deblurring of radiographs using Wiener filter and RLSD. To obtain the data that is required to run this script, download and unzip the zip file at the link data. To run the script as is within the current working directory, the files in the zip file must be placed within a folder called data.

Deblurring using Wiener filter
import numpy as np
from pysaber import wiener_deblur #To deblur using Wiener filtering
import matplotlib.pyplot as plt #To display images
from PIL import Image #To read images in TIFF format

rad_file = 'data/horz_edge_25mm.tif' #Filename of radiograph
bright_file = 'data/horz_bright.tif' #Bright field
dark_file = 'data/horz_dark.tif' #Dark field

sdd = 71003.08 #Source to detector distance (SDD) in micrometers
sod = 24751.89 #Source to object distance (SOD) in micrometers
pix_wid = 0.675 #Pixel width in micrometers
reg_param = 0.1 #Regularization parameter

rad = np.asarray(Image.open(rad_file)) #Read radiograph and convert to numpy array
bright = np.asarray(Image.open(bright_file)) #Read bright field image and convert to numpy array
dark = np.asarray(Image.open(dark_file)) #Read dark field image and convert to numpy array
norm_rad = (rad-dark)/(bright-dark) #Normalize radiograph

#X-ray source blur parameters
src_params = {'source_FWHM_x_axis':2.69,
                'source_FWHM_y_axis':3.01, 
                'norm_power':1.0,
                'cutoff_FWHM_multiplier':10}

#Detector blur parameters
det_params = {'detector_FWHM_1':1.85, 
                'detector_FWHM_2':126.5, 
                'detector_weight_1':0.916, 
                'norm_power':1.0, 
                'cutoff_FWHM_1_multiplier':10, 
                'cutoff_FWHM_2_multiplier':10}

#Deblur the radiograph using Wiener filter
wiener_rad = wiener_deblur(norm_rad,sod,sdd-sod,pix_wid,src_params,det_params,reg_param)

#Display deblurred radiograph
sz = wiener_rad.shape
x = np.arange(-(sz[1]//2),(sz[1]//2)+1,1)*pix_wid
y = np.arange(-(sz[0]//2),(sz[0]//2)+1,1)*pix_wid
plt.pcolormesh(x,y,wiener_rad,cmap='gray')
plt.xlabel('micrometers')
plt.ylabel('micrometers')
plt.title('Wiener deblur')
plt.colorbar()
plt.show()
Deblurring using Regularized Least Squares Deconvolution
import numpy as np
from pysaber import least_squares_deblur #To deblur using regularized least squares deconvolution (RLSD)
import matplotlib.pyplot as plt #To display images
from PIL import Image #To read images in TIFF format

rad_file = 'data/horz_edge_25mm.tif' #Filename of radiograph
bright_file = 'data/horz_bright.tif' #Bright field image
dark_file = 'data/horz_dark.tif' #Dark field image

sdd = 71003.08 #Source to detector distance (SDD) in micrometers
sod = 24751.89 #Source to object distance (SOD) in micrometers
pix_wid = 0.675 #Pixel width in micrometers
reg_param = 0.001 #Regularization parameter

rad = np.asarray(Image.open(rad_file)) #Read radiograph and convert to numpy array
bright = np.asarray(Image.open(bright_file)) #Read bright field image and convert to numpy array
dark = np.asarray(Image.open(dark_file)) #Read dark field image and convert to numpy array
norm_rad = (rad-dark)/(bright-dark) #Normalize radiograph

#X-ray source blur parameters
src_params = {'source_FWHM_x_axis':2.69,
                'source_FWHM_y_axis':3.01, 
                'norm_power':1.0,
                'cutoff_FWHM_multiplier':10}

#Detector blur parameters
det_params = {'detector_FWHM_1':1.85, 
                'detector_FWHM_2':126.5, 
                'detector_weight_1':0.916, 
                'norm_power':1.0, 
                'cutoff_FWHM_1_multiplier':10, 
                'cutoff_FWHM_2_multiplier':10}

#Deblur the radiograph using regularized least squares deconvolution (RLSD)
rlsd_rad = least_squares_deblur(norm_rad,sod,sdd-sod,pix_wid,src_params,det_params,reg_param,thresh=2e-4)

#Display deblurred radiograph of artifact
sz = rlsd_rad.shape
x = np.arange(-(sz[1]//2),(sz[1]//2)+1,1)*pix_wid
y = np.arange(-(sz[0]//2),(sz[0]//2)+1,1)*pix_wid
plt.pcolormesh(x,y,rlsd_rad,cmap='gray')
plt.xlabel('micrometers')
plt.ylabel('micrometers')
plt.title('RLSD deblur')
plt.colorbar()
plt.show()

pysaber Package

Functions

apply_blur_psfs(rad, sod, odd, pix_wid, ...)

Function to blur the input radiograph with point spread functions (PSF) of X-ray source and detector blurs.

estimate_blur(rads, sod, odd, pix_wid, edge)

Estimate parameters of point spread functions (PSF) that model X-ray source blur and/or detector blur from normalized radiographs of a straight sharp edge or mutually perpendicular intersecting pair of sharp edges.

get_detector_psf(pix_wid, det_pars)

Function to compute point spread function (PSF) of detector blur.

get_effective_psf(pix_wid, src_pars, det_pars)

Function to compute the effective point spread function (PSF), which is the convolution of X-ray source and detector PSFs.

get_source_psf(pix_wid, src_pars[, sod, odd])

Function to compute the point spread function (PSF) of X-ray source blur in the plane of the X-ray source or the detector.

get_trans_fit(rad, sod, odd, pix_wid, ...[, pad])

Function to compute the blur model prediction and ideal transmission function for a radiograph with a single straight edge or two mutually perpendicular edges.

get_trans_masks(rad, edge[, tran_pars, pad, ...])

Function to compute transmission function and masks for a radiograph with a single straight edge or two mutually perpendicular edges.

least_squares_deblur(rad, sod, odd, pix_wid, ...)

Function to reduce blur (deblur) in radiographs using a regularized least squares iterative algorithm.

wiener_deblur(rad, sod, odd, pix_wid, ...)

Function to reduce blur (deblur) in a radiograph using Wiener filtering.

Feedback

This software is under development and may have bugs. If you run into any problems, please raise a issue on github at the link issues. There is a lot of scope to improve the performance and functionality of this python package. Furthermore, since this package solves a non-convex optimization problem, there is a remote possibility that the final solution may be a local optima that does not properly fit the data. If there is sufficient interest, we will invest time to significantly reduce the run time, improve convergence and usability, and add additional features and functionalities.

Indices and tables