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

Support channeled dicom slices #16

Merged
merged 6 commits into from
May 14, 2020
Merged

Conversation

jond01
Copy link
Contributor

@jond01 jond01 commented May 11, 2020

Refactored _merge_slice_pixel_arrays to support channeled slices, e.g. RGB.
Discussion in PR #14.
Currently, there is no special test for RGB slices case, I will probably add it a bit later.

if rescale:
voxels = np.empty((num_columns, num_rows, num_slices), dtype=np.float32, order='F')
for k, dataset in enumerate(sorted_slice_datasets):
first_dataset = sorted_slice_datasets[0]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I won't have a chance to read through this in detail until later, but two comments:

  1. We should add a test
  2. What happens if some slices are RGB and others aren't? Would the error message make the problem clear to a user?

Copy link
Contributor

@johndgiese johndgiese May 11, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

BTW thanks for making a PR for this---it will definitely make dicom-numpy more useful and generic!

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

My pleasure :)

  1. Yes, will add a bit later
  2. I do not think it is a legitimate situation (?). My guess is that along with Rows, Columns, BitsStored - the attribute SamplesPerPixel is necessary for calculating pixel_array (and its shape) from PixelData. If there are some slices with different shape, they cannot be combined.
    The current error will be something related to numpy's incompatible shapes.
    In order to avoid it, maybe we should add SamplesPerPixel to the invariant_properties list? (in _validate_slices_form_uniform_grid)

@jond01
Copy link
Contributor Author

jond01 commented May 12, 2020

After inspecting in detail the combine_slices code, I found some optimization possibilities:

  • _sort_by_slice_position is called twice:
    First in the voxels (_merge_slice_pixel_arrays),
    Second in the affine (_ijk_to_patient_xyz_transform_matrix).
    Obviously, it can be done once in the outer combine_slices function.
  • _slice_positions is called 4 times:
    _validate_slices_form_uniform_grid - 1 time
    _slice_spacing - 1 time
    _sort_by_slice_position - 1 time, but _sort_by_slice_position is called twice --> 2 times
  • _extract_cosines is called 6 times:
    _validate_image_orientation - 1 time
    _ijk_to_patient_xyz_transform_matrix - 1 time
    _slice_positions - 1 time, but _slice_positions is called 4 times --> 4 times

All these calls better be done once (and can be).
However,

  1. this may incorporate substantial redesign of the code,
  2. it will change the functionality of _merge_slice_pixel_arrays and _ijk_to_patient_xyz_transform_matrix, which should not be used alone probably (?), but I saw some special tests for _merge_slice_pixel_arrays.

@johndgiese - what do you think?

@johndgiese
Copy link
Contributor

@jond01 it would be great if we could remove some or all of these duplicate calls. When this library was written we were optimizing for clarity of the code over performance. I think performance matters, so removing redundant calls would be welcome. The only situation I can see where we would not want to do this is if the refactor made the code significantly more complicated.

Regarding:

it will change the functionality of _merge_slice_pixel_arrays and _ijk_to_patient_xyz_transform_matrix, which should not be used alone probably (?), but I saw some special tests for _merge_slice_pixel_arrays.

We probably have some tests for private functions as a convenience, but if the python functions are private they are not meant to be used directly.

@jond01
Copy link
Contributor Author

jond01 commented May 13, 2020

Ok,

  • RGB slices test added
  • SamplesPerPixel added to invariant_properties (and MockSlice for tests), in order to avoid different values between slices

Regarding the optimizations - I do not know when I'll have time for that, so let's keep it as an optional improvement issue.

The shape of pixel_array is:
pixel_array.T.shape == (rows, cols)
Or (RGB):
pixel_array.T.shape == (channels, rows, cols)
Copy link
Contributor

@johndgiese johndgiese left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for adding the test---this looks great! I left a few suggestions for naming improvements, including for the existing names we had used (na and nb, which weren't very helpful).

shape = pixel_array.shape
if len(shape) == 2:
na, nb = shape
SamplesPerPixel = 1
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
SamplesPerPixel = 1
samples_per_pixel = 1

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍

na, nb = shape
SamplesPerPixel = 1
else:
na, nb, SamplesPerPixel = shape
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
na, nb, SamplesPerPixel = shape
num_rows, num_columns, samples_per_pixel = shape

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍 except of the order:
num_columns, num_rows, samples_per_pixel = shape

self.Columns = nb
self.Columns = na
self.Rows = nb
self.SamplesPerPixel = SamplesPerPixel
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
self.SamplesPerPixel = SamplesPerPixel
self.SamplesPerPixel = samples_per_pixel

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍

na, nb = pixel_array.shape
shape = pixel_array.shape
if len(shape) == 2:
na, nb = shape
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
na, nb = shape
num_rows, num_columns = shape

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

num_columns, num_rows = shape
(because the pixel_array is not transposed:
pixel_array[j, i] == pixel_array.T[i, j])


self.pixel_array = pixel_array

self.SeriesInstanceUID = 'arbitrary uid'
self.SOPClassUID = 'arbitrary sopclass uid'
self.PixelSpacing = [1.0, 1.0]
self.Rows = na
self.Columns = nb
self.Columns = na
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that your code appears to swap na and nb. My suggestions here swap them back---I am assuming your swap was unintentional.

Suggested change
self.Columns = na
self.Columns = num_columns

Copy link
Contributor Author

@jond01 jond01 May 14, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍 for naming, but the swap is intentional:
#16 (comment)

self.Rows = na
self.Columns = nb
self.Columns = na
self.Rows = nb
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
self.Rows = nb
self.Rows = num_rows

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍

@@ -63,5 +70,16 @@ def axial_slices():
]


@pytest.fixture
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍

@@ -63,5 +70,16 @@ def axial_slices():
]


@pytest.fixture
def axial_rgb_slices():
# SamplesPerPixel = 3
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The purpose and meaning of this comment wasn't entirely clear to me

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yep, unnecessary - removed

def test_rgb_axial_set(self, axial_rgb_slices):
combined, _ = combine_slices(axial_rgb_slices)

manually_combined = np.stack([ds.pixel_array for ds in axial_rgb_slices], axis=0).T
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the transpose expected? Perhaps it was necessary because of the na and nb swap?

Copy link
Contributor Author

@jond01 jond01 May 14, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is expected because we return the transposed pixel array:

voxels[:, :, k] = dataset.pixel_array.T

With more convenient shape:
voxels.shape == (num_rows, num_columns, num_slices)
Or:
voxels.shape == (num_channels, num_rows, num_columns, num_slices)

@johndgiese johndgiese merged commit ee570d7 into innolitics:master May 14, 2020
@johndgiese
Copy link
Contributor

Nice work @jond01 ! Thank you :)

@jond01
Copy link
Contributor Author

jond01 commented May 14, 2020

You're welcome, but it seems like the last commit is missing - can you reopen this PR for a moment?

@johndgiese
Copy link
Contributor

@jond01 sorry about that; we have not had too many outside contributors to our open source projects, so I am still getting used to the typical process.

It doesn't appear that I can re-open the PR. Would you mind making a new one with just your commit? Sorry again about that. I should have waited until you said you were ready to merge.

@jond01
Copy link
Contributor Author

jond01 commented May 14, 2020

Sure, no problem at all. I'm also getting used to it

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants