Re-implementing Python DataContainer methods in C++

evgueni-ovtchinnikov commented 1 year ago

Changes in this pull request

A large number of DataContainer methods are currently implemented in Python with numerous duplications of data by as_array/fill copies. This PR aims at reducing unnecessary copying by moving as much computation as possible to C++.

Testing performed

Related issues

Closes #283

Checklist before requesting a review

[ ] I have performed a self-review of my code
[ ] I have added docstrings/doxygen in line with the guidance in the developer guide
[x] I have implemented unit tests that cover any new or modified functionality
[x] The code builds and runs on my machine
[ ] CHANGES.md has been updated with any functionality change

Contribution Notes

Please read and adhere to the contribution guidelines.

Please tick the following:

[x] The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in SIRF (the Work) under the terms and conditions of the Apache-2.0 License.

paskino commented 1 year ago

@evgueni-ovtchinnikov I merged #1180 removing the Ubuntu 20.04 build. If you merge master into here you should get a green tick from the other builds on the GHA

KrisThielemans commented 1 year ago

getting close?

I would suggest renaming test5.py and tests_six.py to a name that says what they do. Much more useful when they fail 😀

evgueni-ovtchinnikov commented 1 year ago

@KrisThielemans One of the builds keeps failing with

/usr/bin/ld: /home/runner/install/lib/libIO.a(GEHDF5ListmodeInputFileFormat.cxx.o): in function stir::GE::RDF_HDF5::GEHDF5ListmodeInputFileFormat::read_from_file(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&) const: GEHDF5ListmodeInputFileFormat.cxx:(.text+0x83): undefined reference to stir::GE::RDF_HDF5::CListModeDataGEHDF5::CListModeDataGEHDF5(std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> > const&)

Looks like STIR I/O issue, does it not?

KrisThielemans commented 1 year ago

@KrisThielemans One of the builds keeps failing with
Looks like STIR I/O issue, does it not?

yes, unfortunately. This seems to have been there quite a while. It's the DEVEL_BUILD=on version, so using STIR master. Not sure why it's there, especially as STIR itself is building fine. I cannot fix this now, so I suppose we have to ignore this here for the moment. Apologies.

evgueni-ovtchinnikov commented 1 year ago

I have re-designed DataContainer methods dot, sum and max. Here is how it looks for sum now:

SIRF.py:

    def sum(self):
        z = numpy.zeros((2,), dtype=numpy.float32)
        try_calling(pysirf.cSIRF_compute_sum(self.handle, z.ctypes.data))
        if z[1] == 0:
            return z[0]
        return z[0] + 1j*z[1]

csirf.cpp:

void*
cSIRF_compute_sum(const void* ptr_x, void* ptr_z)
{
    try {
        auto const& x = objectFromHandle<DataContainer>(ptr_x);
        x.sum(ptr_z);
        return new DataHandle;
    }
    CATCH;
}

(so, no complex variables anywhere in this interface any more)

gadgetron_data_containers.cpp:

void
MRAcquisitionData::sum(void* ptr) const
{
    int n = number();
    complex_float_t z = 0;
    ISMRMRD::Acquisition a;
    for (int i = 0; i < n;) {
        if (!get_acquisition(i, a)) {
            i++;
            continue;
        }
        z += MRAcquisitionData::sum(a);
        i++;
    }
    complex_float_t* ptr_z = static_cast<complex_float_t*>(ptr);
    *ptr_z = z;
}

(with the same MRAcquisitionData::sum(const ISMRMRD::Acquisition& acq_a) as before)

stir_data_containers.cpp:

void
STIRAcquisitionData::sum(void* ptr) const
{
    int n = get_max_segment_num();
    double t = 0;
    for (int s = 0; s <= n; ++s)
    {
        SegmentBySinogram<float> seg = get_segment_by_sinogram(s);
        SegmentBySinogram<float>::full_iterator seg_iter;
        for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();)
            t += *seg_iter++;
        if (s != 0) {
            seg = get_segment_by_sinogram(-s);
            for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();)
                t += *seg_iter++;
        }
    }
    float* ptr_t = static_cast<float*>(ptr);
    *ptr_t = (float)t;
}

If this is ok, then I believe I am done with this PR.

SyneRBI / SIRF