Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

How to treat mixed uniform and non-uniform discretization in tomography? #286

Open
kohr-h opened this issue Feb 21, 2016 · 1 comment
Open

Comments

@kohr-h
Copy link
Member

kohr-h commented Feb 21, 2016

Problem

In tomography, one often has the case that the angles are not uniformly sampled, while the detector has a uniform pixelization. Currently these two partitions are thrown together to one big partition, which results in a TensorGrid based partition, i.e. for later code, the compound partition is irregular as a whole.

This is a bit unfortunate since we don't have an integrator for general tensor grid based function spaces. An approximation to it could be implemented with 1d multiplications, but that's yet to be done. (For CPU-based spaces this wouldn't be hard, but CUDA is nowhere near.) So instead, we pretend the space is uniformly discretized and take the average angle spacing and the detector pixel size as basis for the weighting, which is okay for now, I guess.

Possible solutions

The question is how we want to treat this in the long run. We could

  1. make a product space with one component per angle or
  2. make DiscreteLp (and RectPartition) more intelligent in the handling of regularity per axis.

Pros and cons

The product space approach would probably be easier to implement, especially regarding CUDA. However, for a large number of angles (> 1000), this can become slow. Perhaps.

The second approach with a more intelligent way of handling discretizations would be the cleaner solution, and also nice to have generally. But as mentioned above, the implementation for CUDA is AFAIK all but trivial (you need to multiply along the first axis with the cell_sizes array), and piping this computation through Numpy altogether will be even slower than the product space version, probably.

Related issues: #156, #201

@kohr-h kohr-h changed the title How to treat mixed uniform and non-uniform discretization? How to treat mixed uniform and non-uniform discretization in tomography? Feb 21, 2016
@kohr-h kohr-h added this to the 0.3 milestone Mar 11, 2016
@adler-j adler-j removed the question label May 30, 2016
@kohr-h kohr-h removed this from the 0.4 milestone Aug 26, 2016
@kohr-h
Copy link
Member Author

kohr-h commented Jan 23, 2017

Infrastructure-wise this is solved by #841, only the operators need to understand grids that are uniform only in certain axes.

kohr-h pushed a commit to kohr-h/odl that referenced this issue Dec 11, 2017
Changes in detail:
- Add dtype with shape to DiscreteLp (mostly __repr__, factory
  functions and some downstream methods). As a consequence,
  `shape_[in,out]` and `ndim_[in,out]` are added for the
  different types of axes, as well as `scalar_dtype`.
- Add `PerAxisWeighting` and make it the default for
  `DiscreteLp` type spaces. Reason: this way the `tspace`
  knows how to deal with removed axes etc. This is important
  for a smooth experience with indexing and reductions over
  axes.
  Helpers for slicing weightings help structure this task.
- Implement `__getitem__` for `TensorSpace` and `DiscreteLp`,
  including (hopefully) reasonable propagation of weights.
  The new `simulate_slicing` utility function simplifies
  this task.
- Allow indexing with ODL tensors of boolean or integer dtype.
- Implement correct weighting for backprojections with
  non-uniform angles, using per-axis weighting and a new
  helper `adjoint_weightings` to apply the weightings in an
  efficient way.
  The correct weighting from the range of `RayTransform`
  is determined by the new `proj_space_weighting` helper.
- Change the space `_*_impl` methods to always expect and
  return Numpy arrays, and adapt the calling code.
- Change behavior of `norm` and `dist` to ignoring weights
  for `exponent=inf`, in accordance with math.
- Improve speed of `all_equal` for comparison of arrays.
- Account for `None` entries in indices in the
  `normalized_index_expression` helper, thus allowing
  creation of new axes.
- Remove `dicsr_sequence_space`, it was largely unused and
  just a maintenance burden. Use a regular `uniform-discr`
  from zero to `shape` instead.
- Remove `Weighting.equiv()` mehtods, never used and hard
  to maintain (n^2 possibilities).
- Remove the (largely useless) `_weighting` helper to create
  weighting instances since it would have been ambiguous
  with sequences of scalars (array or per axis?). Also remove
  the `npy_weighted_*` functions, they were useless, too.
- Remove some dead code from tomo/util.
- A bunch of minor fixes, as usual.

Closes: odlgroup#908, odlgroup#907, odlgroup#1113, odlgroup#965, odlgroup#286, odlgroup#267, odlgroup#1001
kohr-h pushed a commit to kohr-h/odl that referenced this issue Dec 11, 2017
Changes in detail:
- Add dtype with shape to DiscreteLp (mostly __repr__, factory
  functions and some downstream methods). As a consequence,
  `shape_[in,out]` and `ndim_[in,out]` are added for the
  different types of axes, as well as `scalar_dtype`.
- Add `PerAxisWeighting` and make it the default for
  `DiscreteLp` type spaces. Reason: this way the `tspace`
  knows how to deal with removed axes etc. This is important
  for a smooth experience with indexing and reductions over
  axes.
  Helpers for slicing weightings help structure this task.
- Implement `__getitem__` for `TensorSpace` and `DiscreteLp`,
  including (hopefully) reasonable propagation of weights.
  The new `simulate_slicing` utility function simplifies
  this task.
- Allow indexing with ODL tensors of boolean or integer dtype.
- Implement correct weighting for backprojections with
  non-uniform angles, using per-axis weighting and a new
  helper `adjoint_weightings` to apply the weightings in an
  efficient way.
  The correct weighting from the range of `RayTransform`
  is determined by the new `proj_space_weighting` helper.
- Change the space `_*_impl` methods to always expect and
  return Numpy arrays, and adapt the calling code.
- Change behavior of `norm` and `dist` to ignoring weights
  for `exponent=inf`, in accordance with math.
- Improve speed of `all_equal` for comparison of arrays.
- Account for `None` entries in indices in the
  `normalized_index_expression` helper, thus allowing
  creation of new axes.
- Remove `dicsr_sequence_space`, it was largely unused and
  just a maintenance burden. Use a regular `uniform-discr`
  from zero to `shape` instead.
- Remove `Weighting.equiv()` mehtods, never used and hard
  to maintain (n^2 possibilities).
- Remove the (largely useless) `_weighting` helper to create
  weighting instances since it would have been ambiguous
  with sequences of scalars (array or per axis?). Also remove
  the `npy_weighted_*` functions, they were useless, too.
- Remove some dead code from tomo/util.
- A bunch of minor fixes, as usual.

Closes: odlgroup#908
Closes: odlgroup#907
Closes: odlgroup#1113
Closes: odlgroup#965
Closes: odlgroup#286
Closes: odlgroup#267
Closes: odlgroup#1001
kohr-h pushed a commit to kohr-h/odl that referenced this issue Jun 30, 2018
Changes in detail:
- Add dtype with shape to DiscreteLp (mostly __repr__, factory
  functions and some downstream methods). As a consequence,
  `shape_[in,out]` and `ndim_[in,out]` are added for the
  different types of axes, as well as `scalar_dtype`.
- Add `PerAxisWeighting` and make it the default for
  `DiscreteLp` type spaces. Reason: this way the `tspace`
  knows how to deal with removed axes etc. This is important
  for a smooth experience with indexing and reductions over
  axes.
  Helpers for slicing weightings help structure this task.
- Implement `__getitem__` for `TensorSpace` and `DiscreteLp`,
  including (hopefully) reasonable propagation of weights.
  The new `simulate_slicing` utility function simplifies
  this task.
- Allow indexing with ODL tensors of boolean or integer dtype.
- Implement correct weighting for backprojections with
  non-uniform angles, using per-axis weighting and a new
  helper `adjoint_weightings` to apply the weightings in an
  efficient way.
  The correct weighting from the range of `RayTransform`
  is determined by the new `proj_space_weighting` helper.
- Change the space `_*_impl` methods to always expect and
  return Numpy arrays, and adapt the calling code.
- Change behavior of `norm` and `dist` to ignoring weights
  for `exponent=inf`, in accordance with math.
- Improve speed of `all_equal` for comparison of arrays.
- Account for `None` entries in indices in the
  `normalized_index_expression` helper, thus allowing
  creation of new axes.
- Remove `dicsr_sequence_space`, it was largely unused and
  just a maintenance burden. Use a regular `uniform-discr`
  from zero to `shape` instead.
- Remove `Weighting.equiv()` mehtods, never used and hard
  to maintain (n^2 possibilities).
- Remove the (largely useless) `_weighting` helper to create
  weighting instances since it would have been ambiguous
  with sequences of scalars (array or per axis?). Also remove
  the `npy_weighted_*` functions, they were useless, too.
- Remove some dead code from tomo/util.
- A bunch of minor fixes, as usual.

Closes: odlgroup#908
Closes: odlgroup#907
Closes: odlgroup#1113
Closes: odlgroup#965
Closes: odlgroup#286
Closes: odlgroup#267
Closes: odlgroup#1001
kohr-h pushed a commit to kohr-h/odl that referenced this issue Sep 12, 2018
Changes in detail:
- Add dtype with shape to DiscreteLp (mostly __repr__, factory
  functions and some downstream methods). As a consequence,
  `shape_[in,out]` and `ndim_[in,out]` are added for the
  different types of axes, as well as `scalar_dtype`.
- Add `PerAxisWeighting` and make it the default for
  `DiscreteLp` type spaces. Reason: this way the `tspace`
  knows how to deal with removed axes etc. This is important
  for a smooth experience with indexing and reductions over
  axes.
  Helpers for slicing weightings help structure this task.
- Implement `__getitem__` for `TensorSpace` and `DiscreteLp`,
  including (hopefully) reasonable propagation of weights.
  The new `simulate_slicing` utility function simplifies
  this task.
- Allow indexing with ODL tensors of boolean or integer dtype.
- Implement correct weighting for backprojections with
  non-uniform angles, using per-axis weighting and a new
  helper `adjoint_weightings` to apply the weightings in an
  efficient way.
  The correct weighting from the range of `RayTransform`
  is determined by the new `proj_space_weighting` helper.
- Change the space `_*_impl` methods to always expect and
  return Numpy arrays, and adapt the calling code.
- Change behavior of `norm` and `dist` to ignoring weights
  for `exponent=inf`, in accordance with math.
- Improve speed of `all_equal` for comparison of arrays.
- Account for `None` entries in indices in the
  `normalized_index_expression` helper, thus allowing
  creation of new axes.
- Remove `dicsr_sequence_space`, it was largely unused and
  just a maintenance burden. Use a regular `uniform-discr`
  from zero to `shape` instead.
- Remove `Weighting.equiv()` mehtods, never used and hard
  to maintain (n^2 possibilities).
- Remove the (largely useless) `_weighting` helper to create
  weighting instances since it would have been ambiguous
  with sequences of scalars (array or per axis?). Also remove
  the `npy_weighted_*` functions, they were useless, too.
- Remove some dead code from tomo/util.
- A bunch of minor fixes, as usual.

Closes: odlgroup#908
Closes: odlgroup#907
Closes: odlgroup#1113
Closes: odlgroup#965
Closes: odlgroup#286
Closes: odlgroup#267
Closes: odlgroup#1001
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants