Python Wang weights don't work appropriately.

[rupikaraj]

Expected Behavior

I have been trying to reconstruct CBCT 670 projections, with a detector shift of 3 pixels. The expected output is given below.
https://drive.google.com/file/d/1VfDFkuLNrmzMu2Yst8B4_Q09Yymoc_FI/view?usp=drive_link

Actual Behavior

By testing various values, I found that offDetector = [-3, -3] gives the best output.

The output I obtained is in the link below.
https://drive.google.com/file/d/1xf-21ZAjK-m9QC6srX1m6VRd7t3ss2S4/view?usp=drive_link

The geometry parameters are available in this link:
https://drive.google.com/file/d/1OA8CvD5uBea98snv2_GKwvEg0q2l40FP/view?usp=drive_link

Code to reproduce the problem (If applicable)

% I am using FDK for reconstruction. 
recon= FDK(projections, geo, angles);

Can anyone help me improve this reconstruction output? Are there any other factors which I should consider?

Specifications

• MATLAB version: 2021b
• OS: windows 11
• CUDA version: 11.1

[AnderBiguri]

hi @rupikaraj
I wonder if your offDetector is correct, as its in TIGRE its in mm, not pixels. I can see your detector pixels are 0.139, so if you know its 3 pixels of offset, it should be 0.139*3, not 3 that you have offset. Also, is your detector offset in both axis, or just one? You are setting it in both axis, but its often only offseted in horizonal axis in reality.

[Mint41]

hi @rupikaraj I wonder if your offDetector is correct, as its in TIGRE its in mm, not pixels. I can see your detector pixels are 0.139, so if you know its 3 pixels of offset, it should be 0.139*3, not 3 that you have offset. Also, is your detector offset in both axis, or just one? You are setting it in both axis, but its often only offseted in horizonal axis in reality.

I have encountered a similar issue in the reconstruction of CBCT。 I possess a set of multi-shadow data, which I intend to reconstruct using the FDK algorithm. The parameters for this reconstruction are as follows:

geo.DSD = 1500
geo.nDetector = np.array([768, 1024])
geo.sDetector = np.array([768 * 0.388, 1024 * 0.388])
geo.dDetector = geo.sDetector / geo.nDetector
geo.offOrigin = np.array([0, 0, 0])
geo.offDetector = np.array([0, -160])
geo.nVoxel = np.array([256, 512, 512])
geo.sVoxel = np.array([256, 512, 512])
geo.dVoxel = geo.sVoxel / geo.nVoxel
geo.rotDetector = np.array([0, 0, 0])
geo.mode = "cone"

The corresponding scanning parameters are as follows:

<ImagerLat Name="Imager Lateral  [IEC Continuous, mm]" Description="Lateral displacement of the imager in half fan mode">-160</ImagerLat>
      <ImagerLng Name="Imager Longitudinal  [IEC Continuous, mm]" Description="Longitudinal displacement of the imager">0</ImagerLng>
      <ImagerSizeX Name="Imager Size X [pixels]" Description="Size of the image detector in X direction">1024</ImagerSizeX>
      <ImagerSizeY Name="Imager Size Y [pixels]" Description="Size of the image detector in Y direction">768</ImagerSizeY>
      <ImagerResX Name="Imager Resolution X [mm/pixel]" Description="Resolution of the image detector in X direction ">0.388</ImagerResX>
      <ImagerResY Name="Imager Resolution Y [mm/pixel]" Description="Resolution of the image detector in Y direction">0.388</ImagerResY>

The result is：

I am certain that there are no issues with the angle and the Source-to-Detector Distance (DSD). I suspect that the problem lies within the offset. How might I verify the value of this offset?

[AnderBiguri]

@Mint41 as I said to the other user, almost surely your lateral displacement is in pixels. However Tigre takes mm, so you need to convert.

[Mint41]

@Mint41 as I said to the other user, almost surely your lateral displacement is in pixels. However Tigre takes mm, so you need to convert.
-160
I'm sorry, but within this parameters, it appears that the unit of measurement is millimeters, not pixels.

[Mint41]

<ImagerLat Name="Imager Lateral [IEC Continuous, mm]" Description="Lateral displacement of the imager in half fan mode">-160</ImagerLat>

[AnderBiguri]

@Mint41 You are correct! I think you error is no coming from the gemetry, but from the FDK algorithm.
If you use e.g. os_sart, do you get the same rings?

[Mint41]

Indeed, after employing the OS-SART, I obtained more satisfactory reconstruction results. Therefore, perhaps the reconstruction parameters are not at fault. Could you please enlighten me as to what aspect within the FDK algorithm might have led to these outcomes?

[rupikaraj]

Hi @AnderBiguri , I tried changing the offdetector values as you suggested and the output is still not good. Just like @Mint41 , I am getting better outputs in OSSART. Could you suggest how to improve the quality of output using FDK ?

[AnderBiguri]

@Mint41 @rupikaraj there is likely some issue in the Wang weigths of the filtering of FDK. Its not a user-parameter, but a mistake on my side, I think.

There is a function that uses "Wang weights" to apply a weight to correct for detector offset, the relevant code is this:

TIGRE/Python/tigre/algorithms/single_pass_algorithms.py

Lines 103 to 144 in 7279e01

	def preweighting2(proj, geo, theta):
	"""
	Preweighting using Wang function
	:param proj: np.array(dtype=float32),
	Data input in the form of 3d
	:param geo: tigre.utilities.geometry.Geometry
	Geometry of detector and image (see examples/Demo code)
	:param theta: np.array(dtype=float32)
	Angles of projection, shape = (nangles,3) or (nangles,)
	:return: np.array(dtype=float32), np.array(dtype=float32)
	"""
	offset = geo.offDetector[1]
	offset = offset + (geo.DSD / geo.DSO) #* geo.COR
	us = np.arange(-geo.nDetector[1]/2 + 0.5, geo.nDetector[1] /
	2 - 0.5 + 1) * geo.dDetector[1] + abs(offset)
	us = us * geo.DSO / geo.DSD
	abstheta = np.abs(theta * geo.DSO / geo.DSD)
	w = np.ones(proj[0, :, :].shape)
	for ii in range(geo.nDetector[1]):
	t = us[ii]
	if np.abs(t) <= abstheta:
	w[:,ii] = 0.5 * (np.sin((np.pi / 2) * np.arctan(t /
	geo.DSO) / (np.arctan(abstheta / geo.DSO))) + 1)
	if t < -abstheta:
	w[:,ii] = 0
	if theta < 0:
	w = np.fliplr(w)
	proj_w = proj.copy() # preallocation
	for ii in range(proj.shape[0]):
	proj_w[ii, :, :,] = proj[ii, :, :] * w * 2
	return proj_w, w

I believe that the MATLAB side does not have this issue. The relevant function is this: https://github.com/CERN/TIGRE/blob/master/MATLAB/Utilities/redundancy_weighting.m

If you can find differences, I suggest you try to fix them yourself, and test if that fixes the issue. I will try to investigate myself as soon as I have some time, but it won't be very soon, apologies for that, as I am very busy. If you find the issue, let me know and I can update TIGRE to get it fixed.

[Mint41]

thanks for your patient elucidation

[Mint41]

When I endeavoured to proceed without the utilization of Wang weighting
by change the code as
# dowang = kwargs["dowang"] if "dowang" in kwargs else False
dowang = False
the reconstructed images still manifested the same rings. Hence, I think that the issue does not reside within preweighting. Could you suggest where else the problem might potentially originate?

[AnderBiguri]

You may have misunderstood. Want weights are required to get a good image with detector offset. I think that my version of the Wang weigths may have a bug, some small mistake that makes them not perfect, and therefore you get the errors. If you remove Wang weigths then you are guaranteed to get a bad image.

[Medphys-KM]

In MATLAB, I always meet similar problems. Recently, I find if the 'parker' is set to be 'flase', the artifact can be reduced. This is confusing. It seems that other user meet the problem before in issues 382.
#382

[AnderBiguri]

@Medphys-KM I believe the the MATLAB version is now correct, but the python version indeed seems to have a bug