Default center of rotation might be 0.5 pixels off

[dmpelt]

Working on some code for generating mathematical phantoms, I noticed that I was getting center of rotation artifacts when reconstructing the phantoms with TomoPy. I first thought that I was wrongly generating projections in my code, but now I think it might be due to the default center of rotation in TomoPy. To test this, I created the following minimal working example:

import numpy as np

# Create single projection of a disk
x = np.linspace(-1,1,1024)
y = np.zeros(1024)
y[np.abs(x)<=0.1] = 2*np.sqrt(0.1*0.1-x[np.abs(x)<=0.1]**2)

# Put into array with 1024 projections
projs = np.zeros((1024,1,1024), dtype=np.float32)
projs[:] = y
angles = np.linspace(0, np.pi, 1024, False)

# Reconstruct using TomoPy
import tomopy
rec_wrong_cor = tomopy.recon(projs, angles, algorithm='gridrec')

rec_correct_cor = tomopy.recon(projs, angles, algorithm='gridrec', center=(projs.shape[2]-1)/2)

The code generates the projection of a disk, centered on the detector, and uses the same projection for all angles (i.e. the COR is in the middle of the disk). When reconstructing with the default COR that tomopy uses, I get COR artifacts. The last line results in a reconstruction without any COR artifacts. Reconstructing with astra also results in an image that is similar to rec_correct_cor. Images of both reconstructions (gray scale chosen to show the problem better) are shown below. Left is rec_wrong_cor, right is rec_correct_cor.

In practice, it probably doesn't really matter, since for experimental data you typically have to pick a COR that is different from the default in any case.

[dgursoy]

The default rotation axis is assumed to be at the center of the detector's field-of-view (i.e. projs.shape[2]/2), which I believe makes sense. But you should create your projections likewise centered with a False flag like this:

x = np.linspace(-1, 1, 1024, False)

[dmpelt]

I agree that the center of the detector's field-of-view is the best default rotation axis. However, isn't the position of this center in pixel coordinates equal to (projs.shape[2]-1)/2? For example, if you have a detector with two elements, the pixel coordinates that represent the centers of these pixels are located at 0 and 1, respectively. The center of the detector itself is then positioned at location 0.5, which is precisely (projs.shape[2]-1)/2. This is the same for larger detectors.

For this example, I have a disk located at the origin, and a detector that is precisely centered at the origin as well. In this case, I would expect the projection of the disk also be centered on the detector, which would not be the case when using x = np.linspace(-1, 1, 1024, False) instead of x = np.linspace(-1, 1, 1024).

[altaskier]

But for Fourier transforms we go from -N/2 to 0 to +N/2-1 on arrays, so that the center for 256 is at 128 not 127.5.

…

On May 24, 2019, at 09:39, Daniel M. Pelt ***@***.***> wrote: I agree that the center of the detector's field-of-view is the best default rotation axis. However, isn't the position of this center in pixel coordinates equal to (projs.shape[2]-1)/2? For example, if you have a detector with two elements, the pixel coordinates that represent the centers of these pixels are located at 0 and 1, respectively. The center of the detector itself is then positioned at location 0.5, which is precisely (projs.shape[2]-1)/2. This is the same for larger detectors. For this example, I have a disk located at the origin, and a detector that is precisely centered at the origin as well. In this case, I would expect the projection of the disk also be centered on the detector, which would not be the case when using x = np.linspace(-1, 1, 1024, False) instead of x = np.linspace(-1, 1, 1024). — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or mute the thread.

[dgursoy]

@dmpelt Now I see your point but the center range is from the left edge of your detector's FOV to the right edge of your detector's FOV. Assume we have 2 pixels. Then the range for the center is from 0 to 2 (in the FOV), where 0 means the axis is at the left edge of your detector and 2 means it is the right edge of the detector. So, 1 is at the center.

[carterbox]

I believe this is resolved, so I am closing it.