Skip to content

Commit

Permalink
BUG: Fix problems with ordered-subsets not passing CI tests (#557)
Browse files Browse the repository at this point in the history 
Removes the tiny constants added to division as part of #545 because they do not seem to actually prevent the instability in these expectation maximization algorithms. Instead, if instability is observed, the number of update blocks should be set to 1.

This PR also fixes a bug where ospml_quad was receiving an integer pointer cast from a float pointer which caused the blocks to be indexed incorrectly.

This PR also adds more complicated testing to the ordered-subset algorithms so that they are actually tested with more than 1 ordered subset. This seems to be stable for the test dataset used in the CI.
  • Loading branch information
@carterbox
carterbox committed Sep 18, 2021
1 parent 3bc0e9f commit bf34307
Show file tree
Hide file tree
Showing 9 changed files with 90 additions and 59 deletions.
3 changes: 1 addition & 2 deletions source/libtomo/recon/osem.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,6 @@
// ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.

#include <float.h>
#include <string.h>

#include "recon.h"
Expand Down Expand Up @@ -163,7 +162,7 @@ osem(const float* data, int dy, int dt, int dx, const float* center, const float
if(sum_dist2 != 0.0f)
{
ind_data = d + p * dx + s * dt * dx;
upd = data[ind_data] / (simdata[ind_data] + FLT_MIN);
upd = data[ind_data] / simdata[ind_data];
for(n = 0; n < csize - 1; n++)
{
update[indi[n]] += upd * dist[n];
Expand Down
29 changes: 15 additions & 14 deletions source/tomopy/util/dtype.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,6 @@
# ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE #
# POSSIBILITY OF SUCH DAMAGE. #
# #########################################################################

"""
Module for internal utility functions.
"""
Expand All @@ -61,22 +60,23 @@

logger = logging.getLogger(__name__)


__author__ = "Doga Gursoy"
__copyright__ = "Copyright (c) 2015, UChicago Argonne, LLC."
__docformat__ = 'restructuredtext en'
__all__ = ['as_ndarray',
'as_dtype',
'as_float32',
'as_int32',
'as_uint8',
'as_uint16',
'as_c_float_p',
'as_c_int',
'as_c_int_p',
'as_c_float',
'as_c_char_p',
'as_c_void_p']
__all__ = [
'as_ndarray',
'as_dtype',
'as_float32',
'as_int32',
'as_uint8',
'as_uint16',
'as_c_float_p',
'as_c_int',
'as_c_int_p',
'as_c_float',
'as_c_char_p',
'as_c_void_p',
]


def as_ndarray(arr, dtype=None, copy=False):
Expand Down Expand Up @@ -121,6 +121,7 @@ def as_c_int(arr):


def as_c_int_p(arr):
arr = arr.astype(np.intc, copy=False)
c_int_p = ctypes.POINTER(ctypes.c_int)
return arr.ctypes.data_as(c_int_p)

Expand Down
2 changes: 1 addition & 1 deletion source/tomopy/util/extern/recon.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -363,7 +363,7 @@ def c_ospml_quad(tomo, center, recon, theta, **kwargs):
dtype.as_c_int(kwargs['num_iter']),
dtype.as_c_float_p(kwargs['reg_par']),
dtype.as_c_int(kwargs['num_block']),
dtype.as_c_float_p(kwargs['ind_block'])) # TODO: should be int?
dtype.as_c_int_p(kwargs['ind_block']))


def c_pml_hybrid(tomo, center, recon, theta, **kwargs):
Expand Down
Binary file modified test/test_tomopy/test_data/bart.npy
View file
Binary file not shown.
Binary file modified test/test_tomopy/test_data/osem.npy
View file
Binary file not shown.
Binary file modified test/test_tomopy/test_data/ospml_hybrid.npy
View file
Binary file not shown.
Binary file modified test/test_tomopy/test_data/ospml_quad.npy
View file
Binary file not shown.
111 changes: 71 additions & 40 deletions test/test_tomopy/test_recon/test_algorithm.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@
import unittest
from ..util import read_file
from tomopy.recon.algorithm import recon
from numpy.random import default_rng
from numpy.testing import assert_allclose
import numpy as np

Expand Down Expand Up @@ -85,13 +86,21 @@ def test_art(self):
rtol=1e-2)

def test_bart(self):
assert_allclose(recon(self.prj,
self.ang,
algorithm='bart',
num_iter=4,
num_block=3),
read_file('bart.npy'),
rtol=1e-2)
rng = default_rng(0)
ind_block = np.arange(len(self.ang))
rng.shuffle(ind_block)
assert_allclose(
recon(
self.prj,
self.ang,
algorithm='bart',
num_iter=4,
num_block=3,
ind_block=ind_block,
),
read_file('bart.npy'),
rtol=1e-2,
)

def test_fbp(self):
assert_allclose(recon(self.prj, self.ang, algorithm='fbp'),
Expand Down Expand Up @@ -188,47 +197,69 @@ def test_mlem_gpu(self):
assert_allclose(result, read_file('mlem_accel_gpu.npy'), rtol=1e-2)

def test_osem(self):
assert_allclose(recon(self.prj,
self.ang,
algorithm='osem',
num_iter=4,
num_block=3),
read_file('osem.npy'),
rtol=1e-2)
rng = default_rng(0)
ind_block = np.arange(len(self.ang))
rng.shuffle(ind_block)
assert_allclose(
recon(
self.prj,
self.ang,
algorithm='osem',
num_iter=4,
num_block=3,
ind_block=ind_block,
),
read_file('osem.npy'),
rtol=1e-2,
)

def test_ospml_hybrid(self):
assert_allclose(recon(self.prj,
self.ang,
algorithm='ospml_hybrid',
num_iter=4,
num_block=3),
read_file('ospml_hybrid.npy'),
atol=1e-6)
rng = default_rng(0)
ind_block = np.arange(len(self.ang))
rng.shuffle(ind_block)
assert_allclose(
recon(
self.prj,
self.ang,
algorithm='ospml_hybrid',
num_iter=4,
num_block=3,
ind_block=ind_block,
),
read_file('ospml_hybrid.npy'),
rtol=1e-2,
)

def test_ospml_quad(self):
assert_allclose(recon(self.prj,
self.ang,
algorithm='ospml_quad',
num_iter=4,
num_block=3),
read_file('ospml_quad.npy'),
atol=1e-6)
rng = default_rng(0)
ind_block = np.arange(len(self.ang))
rng.shuffle(ind_block)
assert_allclose(
recon(
self.prj,
self.ang,
algorithm='ospml_quad',
num_iter=4,
num_block=3,
ind_block=ind_block,
),
read_file('ospml_quad.npy'),
rtol=1e-2,
)

def test_pml_hybrid(self):
assert_allclose(recon(self.prj,
self.ang,
algorithm='pml_hybrid',
num_iter=4),
read_file('pml_hybrid.npy'),
rtol=1e-2)
assert_allclose(
recon(self.prj, self.ang, algorithm='pml_hybrid', num_iter=4),
read_file('pml_hybrid.npy'),
rtol=1e-2,
)

def test_pml_quad(self):
assert_allclose(recon(self.prj,
self.ang,
algorithm='pml_quad',
num_iter=4),
read_file('pml_quad.npy'),
rtol=1e-2)
assert_allclose(
recon(self.prj, self.ang, algorithm='pml_quad', num_iter=4),
read_file('pml_quad.npy'),
rtol=1e-2,
)

def test_sirt(self):
result = recon(self.prj, self.ang, algorithm='sirt', num_iter=4)
Expand Down
4 changes: 2 additions & 2 deletions test/test_tomopy/test_recon/test_rotation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -104,14 +104,14 @@ def test_find_center(self):

def test_find_center_vo(self):
sim = read_file('sinogram.npy')
cen = find_center_vo(sim, smin=-10, smax=10)
cen = find_center_vo(sim, smin=-10, smax=10, ncore=2)
assert_allclose(cen, 44.75, rtol=0.25)

def test_find_center_vo_with_downsampling(self):
sim = read_file('sinogram.npy')
sim = zoom(sim[:, 0, :], (45, 22), order=3, mode='reflect')
sim = np.expand_dims(sim, 1)
cen = find_center_vo(sim, smin=-10, smax=10)
cen = find_center_vo(sim, smin=-10, smax=10, ncore=2)
assert_allclose(cen, 1002.0, rtol=0.25)

def test_find_center_pc(self):
Expand Down

0 comments on commit bf34307

Please sign in to comment.