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

SIRFReg #211

Merged
merged 391 commits into from Jan 30, 2019

Conversation

Projects
None yet
8 participants
@rijobro
Copy link
Contributor

commented Sep 27, 2018

This wraps NiftyReg's registration and resampling functionality.

@casperdcl

This comment has been minimized.

Copy link
Member

commented Sep 27, 2018

surely would be nicer to

@KrisThielemans

This comment has been minimized.

Copy link
Member

commented Sep 27, 2018

the dependency is via the SuperBuild.

There's another at https://github.com/KCL-BMEIS/niftyreg. We were thinking to point to that one. Seems that the sourceforge version is currently a clone of that. Maybe you can ask Marc which one to use...

@KrisThielemans

This comment has been minimized.

Copy link
Member

commented Sep 27, 2018

@rijobro , can you have a look at Travis? All failing.

@rijobro

This comment has been minimized.

Copy link
Contributor Author

commented Sep 27, 2018

Travis fails because I use a different data commit (it's a PR in SIRF data, here. It's only 1.1MB.

Got the images from Slicer example data (BRAINS module), which as far as I can see is under the Apache licence (link). Is that ok?

@KrisThielemans

This comment has been minimized.

Copy link
Member

commented Sep 28, 2018

ok. I missed that one. sorry.

once merged, you'll have to update the data submodule in this PR which is not as trivial as it should be.

@casperdcl

This comment has been minimized.

@rijobro

This comment has been minimized.

Copy link
Contributor Author

commented Sep 28, 2018

Comments from @DANAJK :

I spotted:
src/Registration/mReg/libload_sirfreg.m
error message mentions mutil… but it seems to be miutil…

src/Registration/mReg/+mSIRFReg/setParameter.m
has mutilities, mUtilities and miutilities – in either code or comments – are these all correct?

src/Registration/mReg/tests/test_mSIRFReg.m
You should use fullfile to compose filenames with paths – this provides robustness against the / and \ file separators on different systems and I think fullfile acts sensibly when there is one missing, or one too many file separators.

Variable output_path might be confusingly named as it appears to be both a path and part of the filename.

In src/Registration/cReg/SIRFRegNiftyAladinSym.cpp
You have the American spelling of ‘neighbor’, but elsewhere you use the English ‘neighbour’

parser.add_key ( "SetInterpolationToNearestNeighbor", reg_aladin<T>::SetInterpolationToNearestNeighbor

src/Registration/cReg/SIRFReg.h
comment about need to give path on line 41

@rijobro

This comment has been minimized.

Copy link
Contributor Author

commented Sep 28, 2018

Comments from @casperdcl :

That's why I'd assumed it was a copy-paste job rather than a proper wrapping. 14k indeed. If there's some pattern to the wrapping maybe we could make use of cmake's configure_file?

@rijobro

This comment has been minimized.

Copy link
Contributor Author

commented Sep 28, 2018

@DANAJK ,

  • Good spot for the spelling in the first point
  • The get/set functions were written by Evgueni and they exist in mSTIR and mGadgetron, so I copied them over to mSIRFReg, too. I don’t know what the difference is between mutilities, miutilities, etc., but everything seems to work.
  • OK for fullfile in Matlab
  • OK for output_path -> output_prefix
  • The American spelling is only used when wrapping NiftyReg. In the case of SetInterpolationToNearestNeighbor, it’s a callback function, so I have to match their spelling (it would only be used in a parameter file). When I’m adding my own functionality, however, I spell it the correct way ;)
  • Ok for the last point on missing documentation.

@casperdcl ,
I’m not sure which part you’re referring to as copy-paste. The wrapping into Matlab and Python? Or the wrapping of NiftyReg itself at the C++ level?

  • If it was the first, it’s not copy-paste. To get the same functionality, you have to write the C-interface and then the relevant python and Matlab classes (as discussed here). So obviously the goal is to get duplicated functionality in the different languages, but it’s certainly (and unfortunately) not a copy-paste job!
  • If it was the second, I haven’t copied any source code from NiftyReg. These are wrappers around the NiftyReg classes (with a fair bit of new functionality), so no copy-paste there either.

@rijobro rijobro force-pushed the SIRFReg branch 2 times, most recently from b4b65d1 to 0572d9e Oct 3, 2018

@rijobro

This comment has been minimized.

Copy link
Contributor Author

commented Oct 3, 2018

Ahead of the code review in two days, I thought it would be useful to give a summary of the class hierarchy of this PR. I'll also circulate an email with the doxygen.

SIRFReg

There is an abstract base class, SIRFReg, from which SIRFRegNiftyAladinSym and SIRFRegNiftyF3dSym inherit. Setup methods in SIRFReg include set_floating_image, set_reference_image and set_parameter_file. There is a virtual void method, update. Once the registration is done, return the warped image with get_output and get the forwards and backwards deformation and displacements with get_deformation_field_fwrd etc.

With aladin, you can also get the forward and backwards transformation matrices. With f3d, you can set the initial transformation, as well as floating and reference time points.

Nifti images

Nifti images are tricky because their data is stored as void*. The datatype is stored, but this means that when you want to cast the data back to perform any operation, you need a bunch of if statements. To simplify that, I only deal internally with float (change the datatype if necessary).

There are 5 classes for images: NiftiImage, NiftiImage3D, NiftiImage3DTensor, NiftiImage3DDisplacement and NiftiImage3DDeformation. NiftiImage is the base class, from which NiftiImage3D and NiftiImage3DTensor inherit. NiftiImage3DDisplacement and NiftiImage3DDeformation then inherit from NiftiImage3DTensor.

NiftiImage has methods such as IO, ==, !=, +, - (add and subtract either a single value or an image), *, get_min, get_max, get_sum, fill, get_norm, crop and deep_copy.

NiftiImage3D can also copies data to and from SIRF's PETImageData, and will do the same for MRImageData at some point in the future.

NiftiImage3DTensor is the general case for displacement and deformations. It has methods create_from_3D_image, save_to_file_split_xyz_components and flip_component (which flips the x-, y- or z-component).

NiftiImage3DDisplacement has a method create_from_def, and NiftiImage3DDeformation has a method create_from_disp. They also both use multiple inheritance; they inherit from SIRFRegTransformation as well as NiftiImage3DTensor.

NiftiImage3DDeformation also has methods create_from_cpp and compose_single_deformation (in which multiple deformations are joined together into one).

SIRFRegTransformation

This is an abstract base class. It has the virtual void method get_as_deformation_field. As long as we can get any transformation as a deformation field, we can compose multiple transformations together. The derived classes are SIRFRegMat44, NiftiImage3DDisplacement and NiftiImage3DDeformation.

SIRFRegMat44 stores a 4x4 matrix for affine transformations. Its methods include reading, writing ==, !=, *, get_as_deformation_field (for which a reference image has to be supplied to get image size, voxel size etc.), deep_copy, fill, get_determinant and print.

SIRFRegImageWeightedMean

To do the weighted mean of a series of images, set the input images with add_image (giving the image and the weighting), perform the mean with update and get the result with get_output.

SIRFRegNiftyResample

Set the input with set_reference_image and set_floating_image. Add transformations with add_transformation_affine, add_transformation_displacement and add_transformation_deformation. Other methods are set_interpolation_type, update and get_output.

Parser

There is a simple parser for supplying the parameter files for the NiftyReg registrations. Classes include SIRFRegParser and SIRFRegParserKey. Set the callback object, set the filename, add in any keys and then parse.

Misc

Namespace SIRFRegMisc that contains all the long functions that I didn't want to put in the classes. Functions include open_nifti_image, save_nifti_image, copy_nifti_image, do_nifti_image_metadata_match, do_nifti_image_metadata_elements_match, dump_headers, dump_header_element and change_datatype.

Executables

Executables include sirfreg_affine_to_disp, sirfreg_dump_nifti_info, sirfreg_aladin, sirfreg_tensor_split_join, sirfreg_change_datatype and sirfreg_crop_image.

C-wrapper, Matlab, Python

I won't get into that here, but so far all methods that exist in C++ are also in Python and Matlab.

@ashgillman

This comment has been minimized.

Copy link
Member

commented Oct 4, 2018

This is a rather large PR, 14k lines. Lots of boilerplate. I wonder about the sustainability of the current architecture if this is what's required to add new functionality. Unrelated, but how much work is in the changes to SWIG needed for MATLAB?

@KrisThielemans

This comment has been minimized.

Copy link
Member

commented Oct 4, 2018

Let's concentrate on the C++ code tomorrow.

The wrapping is a problem. we can briefly discuss tomorrow as well.

@evgueni-ovtchinnikov

This comment has been minimized.

Copy link
Contributor

commented Oct 4, 2018

@ashgillman, SWIG: have alook at STIR/src/swig/stir.i - 16k lines of interface (in Richard's code hardly a 1k), and Matlab is still not officially supported

@rijobro

This comment has been minimized.

Copy link
Contributor Author

commented Oct 4, 2018

As @KrisThielemans says, best not to get bogged down in wrapping. Just for interest, here's the output from cloc run in the SIRFReg folder.

Language files blank comment code
C++ 24 710 927 3146
MATLAB 16 163 379 1186
Python 2 279 248 1181
C/C++ Header 19 473 766 891
CMake 9 55 198 171
C 1 3 17 157
SWIG 1 1 0 6
--- --- --- --- ---
SUM: 72 1684 2535 6738

The first row encompasses both the C++ code and the C-interface to Python and Matlab, which when separated out look like this:

Language files blank comment code
C++ 22 670 822 2385
C-interface 4 40 105 761
@ashgillman
Copy link
Member

left a comment

A number of smallish comments.

I think the main things are just:

  • .par files in SIRF?
  • lots of std::cerr <<, perhaps we should use some logging module so these can be disabled, say while using Python/MATLAB interactively.

# NiftiImage
def try_niftiimage():
time.sleep(0.5)

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

Why all the time.sleep()s? If they aren't necessary they will just slow down the testing process.

This comment has been minimized.

Copy link
@rijobro

rijobro Oct 5, 2018

Author Contributor

I put in the sleep because C++ output isn't necessarily in sync with python. So I waited at the start and end of each test. It's used 48 times, so 24 seconds in total. Was useful during debugging, could potentially remove it now (although output risks not making much sense anymore).

time.sleep(0.5)
sys.stderr.write('\n# --------------------------------------------------------------------------------- #\n')
sys.stderr.write('# Starting NiftiImage test...\n')
sys.stderr.write('# --------------------------------------------------------------------------------- #\n')

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

I'd suggest print('message', file=sys.stderr)

flo_f3d = pSIRFReg.NiftiImage3D(flo_f3d_filename)

# NiftiImage
def try_niftiimage():

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

if you name this test_niftiimage() nose should automatically run the tests.

Then you don't need the test() function.

Actually, this is not just easier but also gives more fine-grained output from nose.

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

Actually, for these tests, a more pythonic layout would be

from unittest import TestCase

class TestNIfTIImage(TestCase):
    def setUp(self):
        self.a = pSIRFReg.NiftiImage()
        self.b = pSIRFReg.NiftiImage(ref_aladin_filename)

    def test_default_constructor(self):
        if self.a.handle is None:
            raise AssertionError()

    def test_save_to_file(self):
        self.b.save_to_file(save_nifti_image)

    def test_addition(self):
        results = self.b - self.b
        if abs(results.get_max()) > 0.0001:
            raise AssertionError('NiftiImage __sub__ failed.')

etc.

Should give better output from Nose and test might even run in parallel, not sure.

This comment has been minimized.

Copy link
@rijobro

rijobro Oct 5, 2018

Author Contributor

I intentionally called them something other than test_ so that I could choose which tests I wanted to run during the debugging. You're right, I should do this before merging.

epsilon *= required_accuracy_compared_to_max;

if (norm > epsilon) {
cout << "\nImages are not equal (norm > epsilon).\n";

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

Probably best not to print any output here.

return norm <= epsilon

This comment has been minimized.

Copy link
@rijobro

rijobro Oct 5, 2018

Author Contributor

I'll put in #ifndef NDEBUG

namespace sirf {
/// SIRFReg parser
template<class Z>
class SIRFRegParser

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

Is this for .par files? Do we want these in SIRF?

{
// Check the nifti exists
if (!this->is_initialised()) {
std::cout << "\nWarning: Nifti image not initialised, can't fill image.\n";

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

Does SIRF have a warning() function like STIR?

This comment has been minimized.

Copy link
@evgueni-ovtchinnikov

evgueni-ovtchinnikov Oct 8, 2018

Contributor

no, it does not, it uses STIR's TextWriter and TextWriterHandle to handle output to stderr, stdout and files

if (!image.is_initialised())
throw runtime_error("Cannot save image to file.");

cout << "\nSaving image to file (" << filename << ")..." << flush;

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

Too much output?

const shared_ptr<nifti_image> im1_sptr = im1.get_raw_nifti_sptr();
const shared_ptr<nifti_image> im2_sptr = im2.get_raw_nifti_sptr();

bool images_match = true;

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

Probably not important, but using

images_match = do_nifti_image_metadata_elements_match("analyze75_orient",im1_sptr->analyze75_orient,im2_sptr->analyze75_orient))
            && do_nifti_image_metadata_elements_match("byteorder",       im1_sptr->byteorder,       im2_sptr->byteorder       )
...

allows short circuiting

}

/// Dump info of multiple nifti images
void dump_headers(const vector<NiftiImage> &ims)

This comment has been minimized.

Copy link
@ashgillman

ashgillman Oct 5, 2018

Member

"print" terminology used elsewhere, rather than dump

free(im->data);

// Set the datatype
if (typeid(newType) == typeid(bool)) im->datatype = DT_BINARY;

This comment has been minimized.

Copy link
@ashgillman
@ashgillman

This comment has been minimized.

Copy link
Member

commented Oct 5, 2018

Also, perhaps the PLS stuff can be cherry picked into a separate PR?

@KrisThielemans
Copy link
Member

left a comment

some relatively small comments. A problem is the close tie-up with Nifti.

void set_floating_image(const NiftiImage3D &floating_image) { _floating_image = floating_image; }

/// Update
virtual void update() = 0;

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

standard SIRF terminology is process()

virtual ~SIRFReg() {}

/// Set parameter file
void set_parameter_file(const std::string parameter_filename) { _parameter_filename = parameter_filename; }

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

prefer const std::string&

/// Get forward deformation field image
NiftiImage3DDeformation get_deformation_field_fwrd() const { return _def_image_fwrd; }
/// Get backward deformation field image
NiftiImage3DDeformation get_deformation_field_back() const { return _def_image_back; }

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

Potential scope for confusion between AcquisitionModel::forward and backward, where it is the transpose. I think the deformation field here is supposed to be the inverse.

Default implementation should throw error

virtual void parse_parameter_file() = 0;

/// Check parameters
virtual void check_parameters();

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

probably should be const.
Document what it does if there's something wrong.

NiftiImage3D _warped_image;

/// Forward displacement field image
NiftiImage3DDisplacement _disp_image_fwrd;

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

storing both displacements and deformations?

void SIRFRegMat44::save_to_file(const string &filename) const
{
// Check that input isn't blank
if (filename.size() == 0)

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

recommended isempty()

/// Save to file
void save_to_file(const std::string &filename) const;

/// Fill

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

can't see why we'd ever want to fill with anything like this. only thing that would make sense is to set it to idenity matrix.

virtual SIRFRegMat44 deep_copy() const;

/// Save to file
void save_to_file(const std::string &filename) const;

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

make virtual such that it can be overriden

/// Multiplication operator
SIRFRegMat44 operator* (const SIRFRegMat44 &other) const;

/// Overload [] operator (const)

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

not so sure that we should expose this. in any case, should return a float[4].

float *operator [](int i) { return _tm.m[i]; }

/// Get raw mat44
const mat44 &get_raw_mat44() const { return _tm; }

This comment has been minimized.

Copy link
@KrisThielemans

KrisThielemans Oct 5, 2018

Member

what's mat44? NiftyReg type?

#include <_reg_tools.h>
#include "stir_data_containers.h"

using namespace std;

This comment has been minimized.

Copy link
@bathomas

bathomas Oct 5, 2018

Contributor

Do you need to pull in all of std just for cout and flush?

@rijobro

This comment has been minimized.

Copy link
Contributor Author

commented Oct 5, 2018

Following the code review today, the following changes have been suggested:

  • Superclasses should exist for resampling and registration which no nothing about NiftyReg/nifti. This will allow for inclusion of different registration packages in the future.
  • Registration classes should return SIRFRegTransformation. You should be able to get this as an affine/deformation/displacement transformation. If a non-rigid registration algorithm was used, an error should be thrown upon trying to return an affine transformation matrix.
  • Interpolation types in the registration class should use enum
  • Namespaces should be used in python/Matlab. Should be pSIRF.Reg.Nifty.RigidRegistration instead of pSIRFRegNiftyAladin. This will allow for generic changing of registration engines.
  • Instead of relying solely on parameter files, why not have a set_parameter(parameter_name, val) method (this saves exposing all methods of all registration packages, which is not sustainable).
  • Return `shared_ptr to avoid copying where possible.
  • Check SimpleITK python registration to look at similarities/differences, and decide if anything needs to change to make them any closer (might make it easier for incorporating other registration packages in the future).
  • A new hierarchy for images has been proposed. More details can be found in issue #212.

@rijobro rijobro referenced this pull request Oct 5, 2018

Closed

SIRFReg code review #213

38 of 41 tasks complete

@rijobro rijobro force-pushed the SIRFReg branch from 25b3f71 to cab6fbf Oct 5, 2018

end
function set_only_2D(self, name)
%***SIRF*** Sets only 2D.
mSTIR.setParameter(self.handle_, 'PLSPrior', 'only_2D', name, 'f')

This comment has been minimized.

Copy link
@evgueni-ovtchinnikov

evgueni-ovtchinnikov Oct 9, 2018

Contributor

the last argument must be 'i'

This comment has been minimized.

Copy link
@rijobro

rijobro Oct 9, 2018

Author Contributor

There are a few files that I didn't mean to make changes to, including this one. Almost all the changes should be in src/Registration, so I'll just reset this file.

This comment has been minimized.

Copy link
@rijobro

rijobro Oct 9, 2018

Author Contributor

@evgueni-ovtchinnikov I was just reseting this, but it seems it doesn't exist on the master. Did you only wrap PLSPrior to Python?

rijobro added some commits Jan 28, 2019

Construct tensor image from three SIRFImageData
Use SIRFImageData instead of NiftiImageData3D

@rijobro rijobro merged commit 51d131e into master Jan 30, 2019

0 of 3 checks passed

Codacy/PR Quality Review Hang in there, Codacy is reviewing your Pull request.
Details
continuous-integration/travis-ci/pr The Travis CI build is in progress
Details
continuous-integration/travis-ci/push The Travis CI build is in progress
Details

@rijobro rijobro deleted the SIRFReg branch Jan 30, 2019

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.