Simulation

phasetorch.mono.sim.simdata(delta_projs, beta_projs, pix_wid, energy, prop_dists, processes=1, device='cpu', dtype='float', mult_FWHM=5, ret_sqroot=False, pad=True)

Simulate normalized radiographs that contain both absorption and phase contrast from projections of absorption index and refractive index decrement.

Simulates normalized radiographs (pixel values that range from 0 to 1 approx.) for X-ray phase contrast CT. As input, it takes the projections of the refractive index decrement (delta), \({\delta}\), and absorption index (beta), \({\beta}\), at various views and object to detector propagation distances. To simulate phase and absorption contrast in the radiographs, the function also needs the pixel width, X-ray energy, and propagation distances as additional inputs. The returned array is the normalized radiographs at all the views and propagation distances. The same length units should be used for all input parameters.

Parameters

• delta_projs (numpy.ndarray) – Refractive index decrement projections.

• beta_projs (numpy.ndarray) – Absorption index projections.

• pix_wid (float) – Pixel width.

• energy (float) – X-ray energy (in keV).

• prop_dists (list of float) – Propagation distances from object to detector (or scintillator).

• processes (int, default=1) – Number of parallel processes. If running on GPUs, this parameter is the total number of available GPUs.

• device ({'cpu', 'cuda'}, default='cpu') – Target device for running this function. String ‘cuda’ specifies GPU execution.

• dtype ({'float', 'double'}, default='float') – Specifies the floating point precision during computation. ‘float’ uses 32-bit single precision floats and ‘double’ uses 64-bit double precision floats.

• mult_FWHM (float, default=5) – Specifies the dominant width of the Fresnel impulse response function in the spatial domain as a multiple of its full width half maximum (FWHM). Used to determine padding.

• ret_sqroot (bool, detault=False) – If True, return the square root of the X-ray intensity measurements i.e., square root of the measured radiographs.

• pad (bool, default=True) – If True, use appropriate amount of padding for the projection arrays during simulation.

Returns

rads – Simulated normalized radiographs.

Return type

numpy.ndarray

Notes

The array shape for delta_projs and beta_projs is \((N_{views}, N_y, N_x)\), where \(N_{views}\) is the total number of views, \(N_y\) is the number of rows, and \(N_x\) is the number of columns. Optionally, the shape can also be \((N_y, N_x)\), which is equivalent to \(N_{views}=1\). The size of prop_dists is \(N_{dists}\), which is the total number of propagation distances. The returned normalized radiographs rads is of shape \((N_{views}, N_{dists}, N_y, N_x)\).

Examples

if __name__ == '__main__': #Required for parallel compute
    import numpy as np
    import matplotlib.pyplot as plt
    from phasetorch.sim import simdata

    x, y = np.meshgrid(np.arange(0, 128.0)*0.001, np.arange(0, 128.0)*0.001) #Grid shape 128 x 128. Pixel 0.001mm.
    projs = np.zeros_like(x)[np.newaxis] #Ground-truth projections. Unitary dimension indicates one view.
    mask = (0.032**2 - (x-0.064)**2) > 0 #Mask for non-zero projections of 0.032mm cylinder.
    projs[0][mask] = 2*np.sqrt((0.032**2 - (x-0.064)**2)[mask]) #Generate path lengths.
    delta_projs, beta_projs = projs*1e-7, projs*1e-8 #delta = 1e-7, beta=1e-8

    rads = simdata(delta_projs, beta_projs, 0.001, 20, [100, 200], processes=1, device='cuda') #Radiographs at distances of 100mm and 200mm for 20keV X-rays

    plt.imshow(rads[0, 0]) #Radiograph at 100mm
    plt.show()

    plt.imshow(rads[0, 1]) #Radiograph at 200mm
    plt.show()

Phase Retrieval

phasetorch.mono.pr.paganin(rads, pix_wid, energy, prop_dist, delta_over_beta, processes=1, device='cpu', mult_FWHM=5, mag_fact=1.0, dtype='float')

Linear phase retrieval algorithm for imaging of single-material objects using single propagation distance phase contrast CT that is also popularly known as Paganin’s phase retrieval.

Retrieves the refractive index decrement projections from normalized radiographs acquired at a single propagation distance. These projections are proportional to the negative logarithm of a low-pass filtered version of the input radiographs. As input, it takes normalized radiographs that are acquired at a single object to detector propagation distance and multiple views. It returns the projections of the refractive index decrement using Paganin’s phase retrieval. As additional inputs, this function also requires specification of the pixel width, X-ray energy, propagation distance, and the ratio of delta over beta. The same length units should be used for all input parameters.

Parameters

• rads (numpy.ndarray) – Normalized radiographs.

• pix_wid (float) – Pixel width.

• energy (float) – X-ray energy.

• prop_dist (float) – Propagation distance from object to detector (or scintillator).

• delta_over_beta (float) – Ratio of delta over beta.

• processes (int, default=1) – Number of parallel processes. If running on GPUs, this parameter is the total number of available GPUs.

• device ({'cpu', 'cuda'}, default='cpu') – Target device for running this function. String ‘cuda’ specifies GPU execution.

• mult_FWHM (float, default=5) – Specifies the dominant width of the Fresnel impulse response function in the spatial domain as a multiple of its full width half maximum (FWHM).

• mag_fact (float, default=1) – Magnification factor given by the ratio of the source to detector distance over the source to object distance. This parameter is one for parallel beam X-rays.

• dtype ({'float', 'double'}, default='float') – Specifies the floating point precision during computation. Must be either ‘float’ or ‘double’. ‘float’ uses 32-bit single precision floats and ‘double’ uses 64-bit double precision floats.

Returns

delta_projs – Refractive index decrement projections.

Return type

numpy.ndarray

Notes

The array shape for rads and delta_projs is \((N_{views}, N_y, N_x)\).

Examples

if __name__ == '__main__': #Required for parallel compute
    import numpy as np
    import matplotlib.pyplot as plt
    from phasetorch.sim import simdata
    from phasetorch.pr import paganin

    x, y = np.meshgrid(np.arange(0, 128.0)*0.001, np.arange(0, 128.0)*0.001) #Grid shape 128 x 128. Pixel 0.001mm.
    projs = np.zeros_like(x)[np.newaxis] #Ground-truth projections. Unitary dimension indicates one view.
    mask = (0.032**2 - (x-0.064)**2) > 0 #Mask for non-zero projections of 0.032mm cylinder.
    projs[0][mask] = 2*np.sqrt((0.032**2 - (x-0.064)**2)[mask]) #Generate path lengths.
    delta_projs, beta_projs = projs*1e-7, projs*1e-8 #delta = 1e-7, beta=1e-8

    rads = simdata(delta_projs, beta_projs, 0.001, 20, 100, processes=1, device='cuda') #Radiographs at distances of 100mm and 200mm for 20keV X-rays
    rads = rads.squeeze(1) #Remove the distance dimension

    delta_projs_est = paganin(rads, 0.001, 20, 100, 10, processes=1, device='cuda') #Paganin phase retrieval

    rmse = np.sqrt(np.mean((delta_projs_est - delta_projs)**2)) #Root mean squared error
    print('Percentage RMSE is {}'.format(100*rmse/np.mean(delta_projs)))

    plt.imshow(delta_projs_est[0]) #Estimated projections of delta
    plt.show()

phasetorch.mono.pr.nlpr(rads, pix_wid, energy, prop_dists, delta_projs, beta_projs, weights=None, processes=1, device='cpu', solver_kwargs=None, mult_FWHM=5, dtype='float', ret_pad=False, verbose=False, ret_conp=False)

Non-linear phase retrieval algorithm for imaging of multi-material objects using X-ray phase contrast CT at multiple propagation distances.

Retrieves the projections of the refractive index decrement and absorption index from normalized radiographs acquired at multiple object to detector propagation distances. This function uses an iterative optimization algorithm that inverts a forward measurement model based on the Fresnel transform. As input, this function takes normalized radiographs. Importantly, this algorithm benefits from an intelligent initialization for the projections of the refractive index decrement (delta), \({\delta}\), and/or absorption index \({\beta}\). It returns refined estimates for the projections of the refractive index decrement and absorption index. As additional inputs, it also requires specification of the pixel width, X-ray energy, and propagation distances. Other optional inputs that drives the numerical optimization method used for phase retrieval can also be provided. The same length units should be used for all inputs.

Parameters

• rads (numpy.ndarray) – Normalized radiographs.

• pix_wid (float) – Pixel width.

• energy (float) – X-ray energy.

• prop_dists (list of float) – Propagation distances from object to detector (or scintillator).

• delta_projs (numpy.ndarray) – Initial values for the refractive index decrement projections.

• beta_projs (numpy.ndarray) – Initial values for the absorption index projections.

• weights (numpy.ndarray, optional) – Weights are a relative importance measure or relative information measure for the radiograph pixel values. If not provided, then every pixel in rads is assigned equal importance.

• processes (int, default=1) – Number of parallel processes. If running on GPUs, this parameter is the total number of available GPUs.

• device ({'cpu', 'cuda'}, default='cpu') – Target device for running this function. String ‘cuda’ specifies GPU execution.

• solver_kwargs (dict, optional) – Dictionary of parameters for controlling the LBFGS algorithm used for numerical optimization. Even though optional, it is recommended to intelligently set this parameter. Please see notes for more information.

• mult_FWHM (float, default=5) – Specifies the dominant width of the Fresnel impulse response function in the spatial domain as a multiple of its full width half maximum (FWHM). This is used to determine padding.

• dtype ({'float', 'double'}, default='float') – Specifies the floating point precision during computation. ‘float’ uses 32-bit single precision floats and ‘double’ uses 64-bit double precision floats.

• ret_pad (bool, default=True) – If True, the returned arrays include additional padding.

• verbose (bool, default=False) – If True, print additional messages to aid in debugging.

• ret_conp (bool, default=False) – If True, return parameters related to convergence of optimization algorithm such as final loss function value, etc.

Returns

• delta_projs (numpy.ndarray) – Projections of the refractive index decrement.

• beta_projs (numpy.ndarray) – Projections of the absorption index.

• convg_params (dict (Optional)) – Returns convergence parameters if the parameter ret_conp is True.

Notes

The shape of rads is \((N_{views}, N_{dists}, N_y, N_x)\). If the shape of rads is \((N_{views}, N_y, N_x)\), it is equivalent to \(N_{dists}=1\). If the shape of rads is \((N_y, N_x)\), it is equivalent to \(N_{views}=1\) and \(N_{dists}=1\). The returned arrays delta_projs and beta_projs are of shape \((N_{views}, N_y, N_x)\). It is recommended to set the parameter solver_kwargs intelligently. solver_kwargs is a python dictionary with several key and value pairs. All keys are python strings (str). The value for key ‘solver’ is either ‘native’ or ‘non-native’. ‘native’ refers to the LBFGS solver that is shipped with pytorch. When using this solver, the other solver parameters for LFBGS as defined in torch.optim.LBFGS are specified as key-value pairs. When the value for key ‘solver’ is ‘non-native’, this function uses a different open source implementation for LBFGS with several key benefits over the ‘native’ solver. In this case, the value for key ‘history_size’ is an integer (int) that specifies the history size for LBFGS. Value for key ‘max_iter is an integer (int) that defines the maximum number of LBFGS iterations. Value for key ‘line_search_fn’ must be ‘Wolfe’. Value for key ‘pct_thresh’ is a floating point threshold (float) for percentage change in estimated values. Value for key ‘chk_iters’ is an integer (int) for the number of iterations to check for convergence. LBFGS is stopped when the percentage change in estimates fall below solver_kwargs[‘pct_thresh’] for solver_kwargs[‘chk_iters’] number of consecutive iterations. Value for key ‘cost_convg’ is False.

Examples

if __name__ == '__main__': #Required for parallel compute
    import numpy as np
    import matplotlib.pyplot as plt
    from phasetorch.sim import simdata
    from phasetorch.pr import paganin, nlpr

    #Grid shape 128 x 128. Pixel 0.001mm.
    x, y = np.meshgrid(np.arange(0, 128.0)*0.001,
                       np.arange(0, 128.0)*0.001)
    #Ground-truth projections. Unitary dimension indicates one view.
    projs = np.zeros_like(x)[np.newaxis]
    #Mask for non-zero projections of 0.032mm cylinder.
    mask = (0.032**2 - (x-0.064)**2) > 0
    #Generate path lengths.
    projs[0][mask] = 2*np.sqrt((0.032**2 - (x-0.064)**2)[mask])
    delta_projs, beta_projs = projs*1e-7, projs*1e-8 #delta = 1e-7, beta=1e-8

    #Radiographs at distances of 100mm and 200mm for 20keV X-rays
    rads = simdata(delta_projs, beta_projs, 0.001, 20, [
                   10, 100, 300], processes=1, device='cuda')

    delta_projs_pag = paganin(
        rads[:,1], 0.001, 20, 100, 10, processes=1, device='cuda') #Paganin phase retrieval

    delta_projs_nlpr, beta_projs_nlpr = nlpr(rads, 0.001, 20, [
                                             10, 100, 300], delta_projs_pag, beta_projs=delta_projs_pag/10, processes=1, device='cuda') #Multi distance NLPR

    #Root mean squared error
    rmse = np.sqrt(np.mean((delta_projs_pag - delta_projs)**2))
    print('Percentage RMSE with Paganin is {}'.format(
        100*rmse/np.mean(delta_projs)))
    #Root mean squared error
    rmse = np.sqrt(np.mean((delta_projs_nlpr - delta_projs)**2))
    print('Percentage RMSE for delta with multi-distance NLPR is {}'.format(100*                   rmse/np.mean(delta_projs)))

    #Estimated projections of delta using multi-distance NLPR
    plt.imshow(delta_projs_nlpr[0])
    plt.show()

phasetorch.mono.pr.nlprcon(rads, pix_wid, energy, prop_dists, delta_projs, alpha, gamma, weights=None, processes=1, device='cpu', solver_kwargs=None, mult_FWHM=5, dtype='float', ret_pad=False, verbose=True, ret_conp=False)

Non-linear phase retrieval algorithm for imaging of homogeneous or weakly absorbing objects using X-ray phase contrast CT at a single propagation distance.

Retrieves the projections of the refractive index decrement and absorption index from normalized radiographs acquired at a single object to detector propagation distance. This function uses an iterative optimization algorithm that inverts a forward measurement model based on the Fresnel transform. As input, this function takes normalized radiographs. Importantly, this algorithm benefits from an intelligent initialization for the projections of the refractive index decrement (delta), \({\delta}\). It returns refined estimates for the projections of the refractive index decrement (delta) and absorption index (beta). As additional inputs, it also requires specification of the pixel width, X-ray energy, and propagation distances. Other optional inputs that drives the numerical optimization method used for phase retrieval can also be provided. The same length units should be used for all inputs.

Parameters

• rads (numpy.ndarray) – Normalized radiographs.

• pix_wid (float) – Pixel width.

• energy (float) – X-ray energy.

• prop_dists (list of float) – Propagation distances from object to detector (or scintillator).

• delta_projs (numpy.ndarray) – Initial values for the refractive index decrement projections.

• alpha (float) – Real value of the estimated transmission function’s exponent. Recommended setting is alpha=1 and gamma=delta/beta.

• gamma (float) – Imaginary value of the estimated transmission function’s exponent. Recommended setting is alpha=1 and gamma=delta/beta.

• weights (numpy.ndarray, optional) – Weights are a relative importance measure or relative information measure for the radiograph pixel values. If not provided, then every pixel in rads is assigned equal importance.

• processes (int, default=1) – Number of parallel processes. If running on GPUs, this parameter is the total number of available GPUs.

• device ({'cpu', 'cuda'}, default='cpu') – Target device for running this function. String ‘cuda’ specifies GPU execution.

• solver_kwargs (dict, optional) – Dictionary of parameters for controlling the LBFGS algorithm used for numerical optimization. Even though optional, it is recommended to intelligently set this parameter. Please see notes for more information.

• mult_FWHM (float, default=5) – Specifies the dominant width of the Fresnel impulse response function in the spatial domain as a multiple of its full width half maximum (FWHM). This is used to determine padding.

• dtype ({'float', 'double'}, default='float') – Specifies the floating point precision during computation. ‘float’ uses 32-bit single precision floats and ‘double’ uses 64-bit double precision floats.

• ret_pad (bool, default=True) – If True, the returned arrays include additional padding.

• verbose (bool, default=False) – If True, print additional messages to aid in debugging.

• ret_conp (bool, default=False) – If True, return parameters related to convergence of optimization algorithm such as final loss function value, etc.

Returns

• delta_projs (numpy.ndarray) – Projections of the refractive index decrement.

• beta_projs (numpy.ndarray) – Projections of the absorption index. Redundant since beta_projs can be uniquely determined from delta_projs, alpha, and gamma.

• convg_params (dict (Optional)) – Returns convergence parameters if the parameter ret_conp is True.

Notes

The shape of rads is \((N_{views}, N_{dists}, N_y, N_x)\). If the shape of rads is \((N_{views}, N_y, N_x)\), it is equivalent to \(N_{dists}=1\). If the shape of rads is \((N_y, N_x)\), it is equivalent to \(N_{views}=1\) and \(N_{dists}=1\). The returned arrays delta_projs and beta_projs are of shape \((N_{views}, N_y, N_x)\). It is recommended to set the parameter solver_kwargs intelligently. solver_kwargs is a python dictionary with several key and value pairs. All keys are python strings (str). The value for key ‘solver’ is either ‘native’ or ‘non-native’. ‘native’ refers to the LBFGS solver that is shipped with pytorch. When using this solver, the other solver parameters for LFBGS as defined in torch.optim.LBFGS are specified as key-value pairs. When the value for key ‘solver’ is ‘non-native’, this function uses a different open source implementation for LBFGS with several key benefits over the ‘native’ solver. In this case, the value for key ‘history_size’ is an integer (int) that specifies the history size for LBFGS. Value for key ‘max_iter is an integer (int) that defines the maximum number of LBFGS iterations. Value for key ‘line_search_fn’ must be ‘Wolfe’. Value for key ‘pct_thresh’ is a floating point threshold (float) for percentage change in estimated values. Value for key ‘chk_iters’ is an integer (int) for the number of iterations to check for convergence. LBFGS is stopped when the percentage change in estimates fall below solver_kwargs[‘pct_thresh’] for solver_kwargs[‘chk_iters’] number of consecutive iterations. Value for key ‘cost_convg’ is False.

Examples

if __name__ == '__main__': #Required for parallel compute
    import numpy as np
    import matplotlib.pyplot as plt
    from phasetorch.sim import simdata
    from phasetorch.pr import paganin, nlprcon

    #Grid shape 128 x 128. Pixel 0.001mm.
    x, y = np.meshgrid(np.arange(0, 128.0)*0.001,
                       np.arange(0, 128.0)*0.001)
    #Ground-truth projections. Unitary dimension indicates one view.
    projs = np.zeros_like(x)[np.newaxis]
    #Mask for non-zero projections of 0.032mm cylinder.
    mask = (0.032**2 - (x-0.064)**2) > 0
    #Generate path lengths.
    projs[0][mask] = 2*np.sqrt((0.032**2 - (x-0.064)**2)[mask])
    delta_projs, beta_projs = projs*1e-7, projs*1e-8 #delta = 1e-7, beta=1e-8

    #Radiographs at distances of 100mm and 200mm for 20keV X-rays
    rads = simdata(delta_projs, beta_projs, 0.001, 20,
                   100, processes=1, device='cuda')
    rads = rads.squeeze(1) #Remove unitary dimension

    #Paganin phase retrieval
    delta_projs_pag = paganin(
        rads, 0.001, 20, 100, 10, processes=1, device='cuda')
    delta_projs_nlpr, beta_projs_nlpr = nlprcon(
        rads, 0.001, 20, 100, delta_projs_pag, alpha=1, gamma=10, processes=1, device='cuda') #Single distance NLPR

    #Root mean squared error
    rmse = np.sqrt(np.mean((delta_projs_pag - delta_projs)**2))
    print('Percentage RMSE with Paganin is {}'.format(
        100*rmse/np.mean(delta_projs)))
    #Root mean squared error
    rmse = np.sqrt(np.mean((delta_projs_nlpr - delta_projs)**2))
    print('Percentage RMSE with single-distance NLPR is {}'.format(100*                   rmse/np.mean(delta_projs)))

    #Estimated projections of delta using single-distance NLPR
    plt.imshow(delta_projs_nlpr[0])
    plt.show()

phasetorch.mono.pr.tropt_nlprcon(delta, beta, num_xy, pix_wid, energy, low_trans=0.01)

For constrained NLPR, compute the constraint parameters of alpha and gamma such that the transmission values are approximately scaled between 0 and low_trans.

Parameters

• delta (float) – Refractive index decrement.

• beta (float) – Absorption index.

• num_xy (int) – Maximum number of pixels along a ray path in the xy-plane of object (plane perpendicular to axis of rotation).

• pix_wid (float) – Pixel width.

• low_trans (float) – Lower bound for the real-valued transmission.

Returns

• alpha (float) – Constraint parameter for the real part of the transmission.

• gamma (float) – Constraint parameter for the imaginary part of the transmission.

phasetorch.mono.pr.ctf(rads, pix_wid, energy, prop_dists, regp=None, rel_regp=None, processes=1, device='cpu', mult_FWHM=5, dtype='float')

Phase retrieval algorithm based on the contrast transfer function (CTF) approach.

Parameters

• rads (list of numpy.ndarray) – A list of numpy arrays, each containing the normalized radiograph measured at a given object to detector distance

• pix_wid (float) – Pixel width in mm

• energy (float) – X-ray energy in keV. Assumes monochromatic source.

• prop_dists (list) – A list of all object to detector distances. Must have a one-to-one correspondance with the radiographs passed using above argument.

• regp (float, default=None) – Parameter of Tikhonov regularization.

• rel_regp (float, default=None) – Relative parameter of Tikhonov regularization.

• processes (int, default=1) – Number of parallel processes. If running on GPUs, this parameter is the total number of available GPUs.

• device ({'cpu', 'cuda'}, default='cpu') – Target device for running this function. String ‘cuda’ specifies GPU execution.

• mult_FWHM (float, default=5) – Specifies the dominant width of the Fresnel impulse response function in the spatial domain as a multiple of its full width half maximum (FWHM). Used to determine padding.

• dtype ({'float', 'double'}, default='float') – Specifies the floating point precision during computation. ‘float’ uses 32-bit single precision floats and ‘double’ uses 64-bit double precision floats.

Returns

• delta_projs (numpy.ndarray) – A numpy array of projection of refractive index decrement

• beta_projs (numpy.ndarray) – A numpy array of projection of absorption index

Notes

The array shape for delta_projs and beta_projs is \((N_{views}, N_y, N_x)\), where \(N_{views}\) is the total number of views, \(N_y\) is the number of rows, and \(N_x\) is the number of columns. Optionally, the shape can also be \((N_y, N_x)\), which is equivalent to \(N_{views}=1\). The size of prop_dists is \(N_{dists}\), which is the total number of propagation distances.

phasetorch.mono.pr.tie(rads, pix_wid, energy, prop_dists, regp=None, rel_regp=None, processes=1, device='cuda', mult_FWHM=5, dtype='double')

Phase retrival algorithm based on the Transport of Intensity Equation (TIE).

Parameters

• rads (list of numpy.ndarray) – A list of numpy arrays, each containing the normalized radiograph measured at a given object to detector distance

• pix_wid (float) – Pixel width in mm

• energy (float) – X-ray energy in keV. Assumes monochromatic source.

• prop_dists (list) – A list of two object to detector distances. Must have a one-to-one correspondance with the radiographs passed using above argument.

• regp (float, default=None) – Parameter of Tikhonov regularization.

• rel_regp (float, default=None) – Relative parameter of Tikhonov regularization.

• processes (int, default=1) – Number of parallel processes. If running on GPUs, this parameter is the total number of available GPUs.

• device ({'cpu', 'cuda'}, default='cpu') – Target device for running this function. String ‘cuda’ specifies GPU execution.

• mult_FWHM (float, default=5) – Specifies the dominant width of the Fresnel impulse response function in the spatial domain as a multiple of its full width half maximum (FWHM). Used to determine padding.

• dtype ({'float', 'double'}, default='float') – Specifies the floating point precision during computation. ‘float’ uses 32-bit single precision floats and ‘double’ uses 64-bit double precision floats.

Returns

• delta_projs (numpy.ndarray) – A numpy array of projection of refractive index decrement

• beta_projs (numpy.ndarray) – A numpy array of projection of absorption index

Notes

The array shape for delta_projs and beta_projs is \((N_{views}, N_y, N_x)\), where \(N_{views}\) is the total number of views, \(N_y\) is the number of rows, and \(N_x\) is the number of columns. Optionally, the shape can also be \((N_y, N_x)\), which is equivalent to \(N_{views}=1\). The size of prop_dists is \(N_{dists}\), which is the total number of propagation distances.

phasetorch.mono.pr.mixed(rads, pix_wid, energy, prop_dists, regp=None, rel_regp=None, tol=0, max_iter=3, processes=1, device='cpu', mult_FWHM=5, dtype='double')

Phase retrieval based on the Mixed approach that combines Contrast Transfer Function (CTF) and Transport of Intensity (TIE) phase retrieval methods.

Parameters

• rads (list of numpy.ndarray) – A list of numpy arrays, each containing the normalized radiograph measured at a given object to detector distance

• pix_wid (float) – Pixel width in mm

• energy (float) – X-ray energy in keV. Assumes monochromatic source.

• prop_dists (list) – A list of two object to detector distances. Must have a one-to-one correspondance with the radiographs passed using above argument.

• regp (float, default=None) – Parameter of Tikhonov regularization.

• rel_regp (float, default=None) – Relative parameter of Tikhonov regularization.

• tol (float, default=0) – Stop iterations when percentage change in reconstruction is less than tolerance.

• max_iter (default=3) – Maximum number of iterations. By default, run for 3 iterations irrespective of the tolerance tol.

• processes (int, default=1) – Number of parallel processes. If running on GPUs, this parameter is the total number of available GPUs.

• device ({'cpu', 'cuda'}, default='cpu') – Target device for running this function. String ‘cuda’ specifies GPU execution.

• mult_FWHM (float, default=5) – Specifies the dominant width of the Fresnel impulse response function in the spatial domain as a multiple of its full width half maximum (FWHM). Used to determine padding.

• dtype ({'float', 'double'}, default='float') – Specifies the floating point precision during computation. ‘float’ uses 32-bit single precision floats and ‘double’ uses 64-bit double precision floats.

Returns

• delta_projs (numpy.ndarray) – A numpy array of projection of refractive index decrement

• beta_projs (numpy.ndarray) – A numpy array of projection of absorption index

Notes

The array shape for delta_projs and beta_projs is \((N_{views}, N_y, N_x)\), where \(N_{views}\) is the total number of views, \(N_y\) is the number of rows, and \(N_x\) is the number of columns. Optionally, the shape can also be \((N_y, N_x)\), which is equivalent to \(N_{views}=1\). The size of prop_dists is \(N_{dists}\), which is the total number of propagation distances.

Utility Functions

phasetorch.util.cphase_cproj(arr, energy)

Parameters

• arr (numpy.ndarray) – Complex or real valued phase/absorption images.

• energy (float) – X-ray energy in units of keV.

Returns

proj – Complex or real projection of delta/beta.

Return type

numpy.ndarray

phasetorch.util.cproj_cphase(arr, energy)

Parameters

• arr (numpy.ndarray) – Complex or real valued delta/beta images.

• energy (float) – X-ray energy in units of keV.

Returns

cphase – Complex or real valued phase/absorption images.

Return type

numpy.ndarray

phasetorch.util.get_wavelength(energy)

Compute the wavelength of X-rays from its energy.

Parameters

energy (float) – Energy of X-rays (units of keV).

Returns

wavelength – Wavelength of X-rays in millimeters (mm).

Return type

float

phasetorch.util.path_lengths_cylinder(N_y, N_x, pix_width, center_x, radius)

Return the path length through a cylinder.

Parameters

• N_y (int) – Number of pixels along the vertical direction of a radiograph (y-axis).

• N_x (int) – Number of pixels along the horizontal direction of a radiograph (x-axis).

• pix_width (float) – Pixel width in units of mm.

• center_x (float) – Center of cylinder along the horizontal direction in units of mm.

• radius (float) – Radius of cylinder in units of mm.

Returns

• projs (numpy.ndarray) – Array of X-ray path lengths in units of mm.

• mask (numpy.ndarray) – Mask to indicate the location of non-zero projections.

phasetorch.util.path_lengths_sphere(N_y, N_x, pix_width, center_y, center_x, radius)

Return the path length through a sphere.

Parameters

• N_y (int) – Number of pixels along the vertical direction of a radiograph (y-axis).

• N_x (int) – Number of pixels along the horizontal direction of a radiograph (x-axis).

• pix_width (float) – Pixel width in units of mm.

• center_y (float) – Center of sphere along the vertical direction in units of mm.

• center_x (float) – Center of sphere along the horizontal direction in units of mm.

• radius (float) – Radius of sphere in units of mm.

Returns

• projs (numpy.ndarray) – Array of X-ray path lengths in units of mm.

• mask (numpy.ndarray) – Mask to indicate the location of non-zero projections.