Problem setting up CircularConeFlatGeometry

[sseebbbbee]

I'm trying to use the CircularConeFlatGeometry like below (proj_data is sinogram from the micro-CT), but I'm probably doing something wrong when setting up the geometry as I get an error.

reco_space = FunctionSpace(Cuboid([-20, -20, -20], [20, 20,20]))
discr_reco_space = uniform_discr([-20, -20, -20], [20, 20, 20], [300, 300, 300], dtype='float32')
src_rad=source_to_origin_mm
det_rad=origin_to_detector_mm
angle_intvl = Interval(0, 2*np.pi)
dparams = Rectangle([-50, -50], [50, 50])
agrid = uniform_sampling(angle_intvl, 360, as_midp=False)
dgrid= uniform_sampling(dparams, [558, 558],as_midp=False)
geom=CircularConeFlatGeometry(angle_intvl, dparams, src_rad, det_rad, agrid,dgrid)
proj_data2=np.zeros((360,558,558))
for ii in range(0,360):
    proj_data2[ii,:,:]=proj_data[:,ii,:]
xray_trafo = DiscreteXrayTransform(discr_reco_space, geom,backend='astra')
back_proj = xray_trafo.adjoint(proj_data2)

I get the following error:

    back_proj = xray_trafo.adjoint(proj_data2)
  File "C:\Users\Sebastian\AppData\Roaming\Python\Python27\site-packages\odl\operator\operator.py", line 356, in __call__
    'of {!r}.'.format(x, self.domain, self))
odl.operator.operator.OpDomainError: input array([[[ 0.36181641,  0.36865234,  0.36889648, ...,  0.01251984,
          0.01896667,  0.01416016],
        [ 0.35327148,  0.36499023,  0.37548828, ...,  0.00637054,
          0.00182819,  0.00429535],
        [ 0.36157227,  0.37109375,  0.36987305, ...,  0.01499939,
          0.00553131,  0.00436783],
        ..., 
        [ 0.30029297,  0.3112793 ,  0.31518555, ...,  0.00400162,
          0.00313568,  0.00699997],
        [ 0.31616211,  0.3203125 ,  0.31591797, ...,  0.01052856,
          0.00304031,  0.00230026],
        [ 0.3203125 ,  0.31494141,  0.30273438, ...,  0.00519562,
         -0.00389481,  0.00630951]],

       [[ 0.33837891,  0.33740234,  0.34008789, ..., -0.01937866,
         -0.01716614, -0.01454926],
        [ 0.33642578,  0.33642578,  0.34448242, ..., -0.01482391,
         -0.01667786, -0.01315308],
        [ 0.33740234,  0.33007812,  0.34130859, ..., -0.01901245,
         -0.01864624, -0.01928711],
        ..., 
        [ 0.27612305,  0.29711914,  0.30249023, ..., -0.01986694,
         -0.01424408, -0.0145874 ],
        [ 0.2746582 ,  0.29125977,  0.30200195, ..., -0.02203369,
         -0.01503754, -0.01322174],
        [ 0.2824707 ,  0.29467773,  0.29443359, ..., -0.0171051 ,
         -0.01454163, -0.01670837]],

       [[ 0.32519531,  0.32861328,  0.33911133, ..., -0.01194   ,
         -0.02085876, -0.02671814],
        [ 0.32910156,  0.32666016,  0.33520508, ..., -0.01277161,
         -0.01679993, -0.02377319],
        [ 0.3293457 ,  0.32763672,  0.33520508, ..., -0.0142746 ,
         -0.01361847, -0.01568604],
        ..., 
        [ 0.30004883,  0.29077148,  0.27832031, ..., -0.01185608,
         -0.01499939, -0.01313019],
        [ 0.31323242,  0.29663086,  0.2800293 , ..., -0.01274872,
         -0.01219177, -0.01698303],
        [ 0.29956055,  0.296875  ,  0.28588867, ..., -0.01589966,
         -0.01282501, -0.01919556]],

       ..., 
       [[ 0.28759766,  0.29931641,  0.31616211, ..., -0.03463745,
         -0.03613281, -0.03430176],
        [ 0.296875  ,  0.30908203,  0.31835938, ..., -0.02980042,
         -0.04031372, -0.04525757],
        [ 0.31298828,  0.31591797,  0.31396484, ..., -0.03543091,
         -0.03625488, -0.04318237],
        ..., 
        [ 0.29125977,  0.28540039,  0.28149414, ..., -0.03167725,
         -0.02867126, -0.02870178],
        [ 0.2722168 ,  0.28027344,  0.27587891, ..., -0.02609253,
         -0.0279541 , -0.03292847],
        [ 0.27124023,  0.27954102,  0.28564453, ..., -0.02531433,
         -0.02520752, -0.02482605]],

       [[ 0.32885742,  0.30395508,  0.29223633, ..., -0.02932739,
         -0.03704834, -0.02978516],
        [ 0.32666016,  0.30419922,  0.29125977, ..., -0.03353882,
         -0.03344727, -0.03384399],
        [ 0.32104492,  0.31347656,  0.3034668 , ..., -0.03427124,
         -0.03356934, -0.03326416],
        ..., 
        [ 0.28271484,  0.28540039,  0.28149414, ..., -0.03274536,
         -0.03146362, -0.02841187],
        [ 0.29443359,  0.2902832 ,  0.28613281, ..., -0.0322876 ,
         -0.02912903, -0.034729  ],
        [ 0.28808594,  0.29150391,  0.28686523, ..., -0.03414917,
         -0.03668213, -0.03393555]],

       [[ 0.31811523,  0.30810547,  0.31030273, ..., -0.03118896,
         -0.04290771, -0.04269409],
        [ 0.31811523,  0.30859375,  0.30786133, ..., -0.0302124 ,
         -0.03149414, -0.03735352],
        [ 0.30810547,  0.31030273,  0.30810547, ..., -0.03234863,
         -0.03317261, -0.03659058],
        ..., 
        [ 0.2824707 ,  0.28173828,  0.27807617, ..., -0.03686523,
         -0.03619385, -0.02871704],
        [ 0.27172852,  0.27539062,  0.28173828, ..., -0.0296936 ,
         -0.02912903, -0.03161621],
        [ 0.27172852,  0.26879883,  0.28735352, ..., -0.02696228,
         -0.02687073, -0.0249939 ]]]) not an element of the domain DiscreteLp(
    FunctionSpace(Cuboid([0.0, -50.0, -50.0], [6.2831853071795862, 50.0, 50.0])),
    RegularGrid([0.0, -50.0, -50.0], [6.2831853071795862, 50.0, 50.0], [360, 558, 558]),
    Rn(112091040, 'float32', weight=0.00056412458811)
    ) of Adjoint: DiscreteLp(
    FunctionSpace(Cuboid([0.0, -50.0, -50.0], [6.2831853071795862, 50.0, 50.0])),
    RegularGrid([0.0, -50.0, -50.0], [6.2831853071795862, 50.0, 50.0], [360, 558, 558]),
    Rn(112091040, 'float32', weight=0.00056412458811)
    ) -> uniform_discr([-20.0, -20.0, -20.0], [20.0, 20.0, 20.0], [300, 300, 300]).

EDIT: adler-j formated code.

[adler-j]

The error says that the input data (which is a numpy ndarray), is not an element of the domain of the adjoint. You need to cast the array to a element in the domain by calling

proj_data2 = xray_trafo.range(proj_data2)

This will be improved once we merge issue #100, which will make your code valid.

[adler-j]

Btw some short python hints for your code:

• This is created but never used

  reco_space = FunctionSpace(Cuboid([-20, -20, -20], [20, 20,20]))

• Why as_midp=False here? I think the xray transform expects it to be true

  agrid = uniform_sampling(angle_intvl, 360, as_midp=False)
  dgrid= uniform_sampling(dparams, [558, 558],as_midp=False)

• This

  proj_data2=np.zeros((360,558,558))
  for ii in range(0,360):
    proj_data2[ii,:,:]=proj_data[:,ii,:]

  Can be written more pythonically

  proj_data2 = proj_data.swapaxes(0, 1).copy() #Need to copy since swapaxes returns a view

[sseebbbbee]

I guessed that it could be that but didn't know how to solve it. Now when writing

proj_data2 =xray_trafo.range(proj_data2)

I get:

TypeError: 'DiscreteLp' object is not callable

About as_midp=False, I just copied from https://github.com/odlgroup/odl/blob/issue-87__conebeam/examples/astra_parallel_2d_fp.py where you have False for angles, but yes if set True the angle apcing will be correct.

Thanks for the third tip, didn't know about that function.

[sseebbbbee]

I might also ask here, for CircularConeFlatGeometry you write:

The source moves on a circle with radius r, and the detector reference point is opposite to the source on a circle with radius R and aligned tangential to the circle.

If I want the detector to not be exactly opposite to the source and maybe translated 50units in x-direction. How do I set that up?

[kohr-h]

You can make the detector itself not centered around zero. Apart from units, you can do

det_params = odl.Rectangle([-width / 2, -height / 2] - (shift_x, 0), [width / 2, height / 2]  - (shift_x, 0))

and then use that guy to initialize the geometry. Note that the interface is still a bit unwieldy since this functionality was added just lately (on the branch). We're going to write convenience functions for the creation of geometries.

Edit: syntax corrected

[kohr-h]

proj_data2 =xray_trafo.range(proj_data2)

That one was a typo, it should be xray_trafo.range.element(proj_data2)

About as_midp=False, I just copied from https://github.com/odlgroup/odl/blob/issue-87__conebeam/examples/astra_parallel_2d_fp.py where you have False for angles, but yes if set True the angle apcing will be correct.

Here's what happens:

For as_midp=True, the sampling points (=angles in this case) are interpreted as midpoints of small subintervals, so when sampling [-1, 1] with 2 points, you will get -0.5 and 0.5, which is probably not what you want (but this makes most sense for the detector). If you set as_midp=False you'll get the points -1 and 1 instead, which is correct for the angles, but wrong for the detector.

The problem is that we're not properly subdividing and allowing different as_midp in different axes. There is an issue about that, it will be solved, but for now you have to live with the current situation, sorry.

[sseebbbbee]

Great, then I assume i can do the same with the reconstruction space? I mean, that it also don't have to be centered around zero?

I guess you are starting to hate me, but now I get the following error:

Configuration Error in ReconstructionGeometry2D: No GridColCount tag specified.Traceback (most recent call last):
  File "C:\Users\Sebastian\Dropbox\Masterarbete\kod\optimerad kod\ODL_real_data.py", line 86, in <module>

    back_proj = xray_trafo.adjoint(proj_data2)
  File "C:\Users\Sebastian\AppData\Roaming\Python\Python27\site-packages\odl\operator\operator.py", line 372, in __call__
    result = self._call(x, *args, **kwargs)
  File "C:\Users\Sebastian\AppData\Roaming\Python\Python27\site-packages\odl\tomo\operators\xray_trafo.py", line 181, in _call
    return self.forward._call_adjoint(inp)
  File "C:\Users\Sebastian\AppData\Roaming\Python\Python27\site-packages\odl\tomo\operators\xray_trafo.py", line 163, in _call_adjoint
    self.domain)
  File "C:\Users\Sebastian\AppData\Roaming\Python\Python27\site-packages\odl\tomo\backends\astra_cpu.py", line 185, in astra_cpu_backward_projector_call
    ndim=reco_space.grid.ndim)
  File "C:\Users\Sebastian\AppData\Roaming\Python\Python27\site-packages\odl\tomo\backends\astra_setup.py", line 433, in astra_data
    return create(astra_dtype_str, astra_geom)
  File "C:\Users\Sebastian\Anaconda\lib\site-packages\astra\data3d.py", line 41, in create
    return d.create(datatype,geometry,data)
  File "astra\data3d_c.pyx", line 81, in astra.data3d_c.create (astra\data3d_c.cpp:2256)
Exception: Geometry class not initialized.

Edit @kohr-h: added fences to console output

[kohr-h]

Great, then I assume i can do the same with the reconstruction space? I mean, that it also don't have to be centered around zero?

I realized that I wrote it in the wrong order. The Rectangle takes minimum point as first and maximum point as second argument. So the correct syntax is

det_params = odl.Rectangle([-width / 2, -height / 2] - (shift_x, 0), [width / 2, height / 2]  - (shift_x, 0))

And yes, for the volume you can do just the same, and it should work. Let's see :-)

[kohr-h]

Regarding your error: @moosmann can you check that one? Maybe it's already fixed upstream.

[kohr-h]

@sseebbbbee another GitHub hint: if you paste console output use "fences" so that the text will be displayed verbatim, i.e.

My fenced text It can stretch over multiple lines.

[moosmann]

I'm going to check the error.

[adler-j]

For as_midp=True, the sampling points (=angles in this case) are interpreted as midpoints of small subintervals, so when sampling [-1, 1] with 2 points, you will get -0.5 and 0.5, which is probably not what you want (but this makes most sense for the detector). If you set as_midp=False you'll get the points -1 and 1 instead, which is correct for the angles, but wrong for the detector.

I disagree. He currently gives angles as:

angle_intvl = Interval(0, 2*np.pi)
agrid = uniform_sampling(angle_intvl, n, as_midp=False)

with n=360. If we set n=2 we get:

>>> angle_intvl = Interval(0, 2*np.pi)
>>> agrid = uniform_sampling(angle_intvl, 2, as_midp=False)
>>> agrid.points()
array([[ 0.        ],
       [ 6.28318531]])

This obviously makes no sense (both angles are the same modulo 2pi). Hence, unless the user gives explicit angles, as_midp=True would be correct.

[wjp]

I'm jumping into the discussion a bit late, but I think here the problem with the last error is that it's still trying to use a CPU projector for a 3D geometry. (Which Astra doesn't support.)

[kohr-h]

I'm jumping into the discussion a bit late, but I think here the problem with the last error is that it's still trying to use a CPU projector for a 3D geometry. (Which Astra doesn't support.)

Right @wjp. We should raise in that case.

This obviously makes no sense (both angles are the same modulo 2pi). Hence, unless the user gives explicit angles, as_midp=True would be correct.

No @adler-j. That one of them is false doesn't mean that the other one is correct. Basically both are wrong. The False case produces 0, 2 * pi which is nonsense as you write, but the True case gives you the angles pi / 4, 3 * pi /4, which makes no real sense either. The correct way to do this would be to recognize that we have a periodic structure here and default to 0, pi when as_midp=False.

In conclusion, we need something which tells that a discretization is periodic.

[adler-j]

The correct way to do this would be to recognize that we have a periodic structure here and default to 0, pi when as_midp=False.

This still doesn't make much sense, you would get basically the same projection from both sides, at least in parallel beam that really makes no sense. My opinion is that we use as_midp=True as the default since this causes the "sensible" answer of two projections "optimally spaced", aka at pi/2 difference.

Anyway, going from 0 to pi gives you a sub-sampled problem in fan/cone beam CT. You would need to go from 0 to pi + fan beam angle. It is probably simpler to default to full angle (2 pi) scans in this case. It is also the default in many applications such as in linear accelerators and in most (all?) CT machines.

In conclusion, we need something which tells that a discretization is periodic.

Agree, but it is only the case in parallel beam geometries that it is pi periodic, else it is 2pi periodic.

[sseebbbbee]

I'm jumping into the discussion a bit late, but I think here the problem with the last error is that it's still >trying to use a CPU projector for a 3D geometry. (Which Astra doesn't support

Is this something I can change?

[kohr-h]

Is this something I can change?

If you find the issue and send a patch, sure :-). But I think we can find the source of the error quite quickly.

[sseebbbbee]

I think I'll let you guys do it :)

[moosmann]

I'll soon make a push where the problem is fixed including your code in the examples folder.

[moosmann]

Just pushed. It should work now. The example is in odl/examples/tomo/xray_trafo.py

[kohr-h]

This still doesn't make much sense, you would get basically the same projection from both sides, at least in parallel beam that really makes no sense. My opinion is that we use as_midp=True as the default since this causes the "sensible" answer of two projections "optimally spaced", aka at pi/2 difference.

Anyway, going from 0 to pi gives you a sub-sampled problem in fan/cone beam CT. You would need to go from 0 to pi + fan beam angle. It is probably simpler to default to full angle (2 pi) scans in this case. It is also the default in many applications such as in linear accelerators and in most (all?) CT machines.

This is correct, but you indicate yourself that "correct" vs. "incorrect" behavior is very problem-dependent. I think we need a broader scope than just "makes sense for CT" or not. I'll put down some thoughts in the coming time, but my current take on the whole thing is the following:

• Grids are just grids, i.e. a sparse representation of points in space
• Tessellations come on top and interpret the points in terms of a partition of a spatial domain (as a generalization of the as_midp flag).
• Those are (for now) very simplistic, rectangular, but have two axes of parametrization: midpoint yes/no and periodic yes/no.
• Regarding where to start and where to end, I think that the as_midp rule for periodic=False should be handled as it is now, and just the periodic case for as_midp=False should be adapted that you get the grid x0 + j * dx, j=0,...,n-1, i.e. everything except the last point. Selection of sensible default start end end angles should then be up to the geometry.

I'll make an issue on the periodic tessellations.

[adler-j]

Sounds good on all points.

[adler-j]

Is this issue done? @sseebbbbee.

[sseebbbbee]

Sorry for the late reply, I had no opportunity to test until now. I still have problems. Now when using https://github.com/odlgroup/odl/blob/issue-87__conebeam/examples/tomo/xray_trafo.py (and alsowith my code) then I get the following error:

    discr_proj_data = xray_trafo(discr_vol_data)
  File "c:\users\sebastian\odl\odl\operator\operator.py", line 372, in __call__
    result = self._call(x, *args, **kwargs)
  File "c:\users\sebastian\odl\odl\tomo\operators\xray_trafo.py", line 165, in _call
    self.range)
  File "c:\users\sebastian\odl\odl\tomo\backends\astra_cuda.py", line 92, in astra_gpu_forward_projector_call
    vol_id = astra_data(vol_geom, datatype='volume', data=vol_data)
  File "c:\users\sebastian\odl\odl\tomo\backends\astra_setup.py", line 427, in astra_data
    return link(astra_dtype_str, astra_geom, data.asarray())
  File "C:\Users\Sebastian\Anaconda\lib\site-packages\astra\data3d.py", line 61, in link
    return d.create(datatype,geometry,data,True)
  File "astra\data3d_c.pyx", line 72, in astra.data3d_c.create (astra\data3d_c.cpp:2142)
  File "C:\Users\Sebastian\Anaconda\lib\site-packages\astra\pythonutils.py", line 49, in geom_size
    elif geom['type'] == 'parallel' or geom['type'] == 'fanflat':
KeyError: 'type'

[moosmann]

I'm not sure about this, but maybe you have to use the latest ASTRA version from github.

[sseebbbbee]

I thought I could update to that version just by copying this folder https://github.com/astra-toolbox/astra-toolbox/tree/master/python/astra into the site-packages directory but apperently not. I get the error:

    import astra
  File "C:\Users\Sebastian\Anaconda\lib\site-packages\astra\__init__.py", line 26, in <module>
    from . import matlab as m
  File "C:\Users\Sebastian\Anaconda\lib\site-packages\astra\matlab.py", line 38, in <module>
    from . import astra_c
ImportError: cannot import name astra_c

There is no astra_c.pyd file only astra_c.pyx in the folder, maybe that's the reason but I don't know.

I think i might give up on using ODL for now, but thanks for the help, hopefully you found some error that is helpful for you.

[adler-j]

Astra is a compiled package, so you need to install it properly, the instructions are here:

https://github.com/astra-toolbox/astra-toolbox/

[sseebbbbee]

But the instructions there is just for creating .mex files for matlab? I did like above for version 1.6

[adler-j]

You are actually correct, it doesn't mention installing python on windows. Perhaps @wjp has a word on this?

[wjp]

Building python extensions on Windows is sadly a rather painful process. But you're right, we should really write it up.

[adler-j]

Is there pre-built 1.7 for him to use?

[wjp]

Sure: http://sourceforge.net/projects/astra-toolbox/files/astra-1.7beta/

[adler-j]

Can you try that one @sseebbbbee?

[sseebbbbee]

Thanks! Now it's working!

[adler-j]

Great! Please ask if there is anything else.

[sseebbbbee]

I will, don't worry ;) Should I ask here or make new issues? Astra has an already implemented feldkamp algorithm, is there anyway to use that one together with ODL?

[adler-j]

Best make new issues so we dont get enormous discussions that are hard to follow.

FDK i dont know, but i think not. Perhaps @moosmann knows? We could need to implement it as an inverse.

[kohr-h]

Oh god, no, not as an inverse. Far from it. Convenience has its limits.

[adler-j]

Well make it a free-standing class, but I don't see why not make it an inverse property. Most standard libraries (matlab and skimage.transform at least) have the radon <-> iradon "pair" and in many applications, having a FBP initial guess is cool and simply calling op.inverse would do well.

If we dont do inverse in this case I really don't see any cases where we will.

[adler-j]

I'll make a new issue on FBP, closing this one.