SyneRBI / SIRF

Main repository for the CCP SynerBI software
http://www.ccpsynerbi.ac.uk
Other
58 stars 29 forks source link

Sptr functinality for sirf::AcquisitionsVector #243

Closed johannesmayer closed 5 years ago

johannesmayer commented 5 years ago

Currently the sirf::AcquisitionsVector carries a data member of type std::vector<gadgetron::shared_ptr<ISMRMRD::Acquisition> > acqs_;

However, there is no function directly appending a gadgetron::shared_ptr<ISMRMRD::Acquisition> with call-by-value argument to have a memory-efficient adding of shared pointers. Currently adding an ISMRMRD::Acquisition involves gadgetron::make_shared which generates a duplication of memory (c.f. #204 ) .

Functions get_sptr_acquisition(gadgetron::shared_ptr<ISMRMRD::Acquisition>) and append_sptr_acquisition(gadgetron::shared_ptr<ISMRMRD::Acquisition>) are needed to allow constructions such as the following without memory duplication:

sirf::AcquisitionsVector all_acquis = some filled AcquisitionsVetor;

size_t num_acq = all_acquis.items();
size_t num_iterations = 1000;

std::vector< sirf::AcquisitionsVector > acq_vec_vector;
for(size_t i=0; i<num_iterations; i++)
    {
    std::cout << " NUM IT "<<i<<std::endl;

        sirf::AcquisitionsVector temp_vec;
    temp_vec.copy_acquisitions_info( all_acquis ) ;

    for(size_t i_acq=0; i_acq<num_acq; i_acq++)
    {   
        auto sptr_acq = all_acquis.get_sptr_acquisition(i_acq);
        temp_vec.append_sptr_acquisition(sptr_acq);
    }
    acq_vec_vector.push_back( temp_vec );
}

Also the index function MRAcquisitionData::get_acquisition() calling this->index() is not save as

int MRAcquisitionData::index(int i)
        {
            if (index_ && i >= 0 && i < (int)number())
                return index_[i];
            else
                return i;
        }

does not protect from i being out of range.

KrisThielemans commented 5 years ago

I don't quite understand this example I'm afraid. As far as I can see, for all practical purposes, temp_vec will now be identical to all_acquis (ok, it's a different vector, but it points to all the same objects). So, that seems to mean that you'd now have 1000 copies in acq_vec_vector.

I'm afraid I fail to see why this would be useful. As an example, suppose I now modify one of them via acq_vec_vector[5].get_sptr_acquisition(3), as far as I can see actually all 1000 of them will be modified (as they all point to the same object). That is probably not what you want, and incredibly dangerous.

Can you maybe elaborate a bit more what you want? Either you need 1000 different objects, in which case cloning seems unavoidable, or you need only 1, or I suppose you need to pretend you have 1000 objects, but it's fine that they're all the same (which is what you seem to have here).

johannesmayer commented 5 years ago

I want to pretend to have a 1000 objects. In my development branch I have a class sirf::MRDynamic which carries two things. 1.) all readout lines of the simulation process in the form of a sirf::AcquisitionsVector 2.) a signal (such as an ECG or a breathing navigator signal) Now we can do this:

sirf::AcquisitionsVector all_acquisitions;
all_acquisitions.read( SOME_ISMRMRD_FILENAME.h5 );

MRDynamic mr_dyn;
mr_dyn.set_signal( some_signal_in_correct_format );
mr_dyn.set_acquisitions( sirf::AcquisitionsVector& all_acquisitions );

mr_dyn.bin_acquisitions();

The last call to bin_acquisitions() actually then bins the data based on the supplied signal, i.e. fills an std::vector < sirf::AcquisitionsVector > inside MRDynamic with the readout lines sorted such that for 10 bins, one has 10 disjoint (!) sets of ISMRMRD::Acquisitions each strored in it's own AcquisitionsVector.

If you have multiple dynamics, e.g. cardiac and respiratory motion, and each is simulated in 10 different states, then one needs to simulate the phantom in 10x10 = 100 different (cardiac, respiratory) dynamic states, and only acquire the readout lines for which the patient would have been in the respective ( cardiac, respiratory ) state by intersecting the binned datasets dynamic with the ones of every other dynamic.

But if passing the data, and binning the data, intersecting the data, and then handling the data inside the simulation class, and passing it to the MRAcquisitionModel, and then storing it afterwards always duplicates the data then with a 5GB dataset this becomes very memory consuming. So I decided it's more efficient to let everything point to the one place in memory it actually exists despite the danger that I can now manipulate it from many points.

KrisThielemans commented 5 years ago

but if you just need 1000 identical copies, wouldn't it be better to make that explicit in the code? i.e. have only 1 object instead of 1000. If that cannot be done, why not have

std::vector< std::shared_ptr<sirf::AcquisitionsVector > > acq_vec_sptr_vector(1000, all_acquis_sptr);

But I'm afraid I still fail to see how this is supposed to work. The intersecting will mean that each of the 1000 acquisitions could have different read-out lines, and therefore need to be different. ok, I guess it's possible that some of them would be the same read-out lines, but then if there's contrast-kinetics, the data would still be different.

sorry, I'm slightly confused.

johannesmayer commented 5 years ago

The terminology is such that, one ISMRMRD::Acquisition is one single readout line.

Currently the only way to access a readout line in SIRF is to create a temporary object of type ISMRMRD::Acquisition and call the accessor of the AcquisitionsVector

sirf::AcquisitionsVector acquis_vec = all_ro_lines;

sirf::AcquisitionsVector other_acquis_vec;
other_acquis_vec.copy_acquisitions_info( acquis_vec );

for( size_t i_acq=0; i_acq<all_ro_lines.number(); i_acq++)
{
    ISMRMRD::Acquisition temp_acq;
    acquis_vec.get_acquisition( num_acquisition, temp_acq);

    // do whatever you like with the acquisition

    other_acquis_vec.append_acquisition( temp_acq );
}

But the in gadgetron_data_containers.h the appending function calls make shared on the temporary object:

void sirf::AcquisitionsVector::append_acquisition(ISMRMRD::Acquisition& acq)
{
    acqs_.push_back(gadgetron::shared_ptr<ISMRMRD::Acquisition>
        (new ISMRMRD::Acquisition(acq)));
}

Since the shared pointer is created on the temporary acquisition, its memory is duplicated, and the object is not deleted as a shared pointer is pointing to it.

That is I cannot access do something with them and append them again somewhere else without a duplication in memory. And this does not change if I access them via a shared pointer of the AcquisitionsVector since at some point I need to access the readout lines.

johannesmayer commented 5 years ago

The ISMRMRD::Acquisition::resize() method allows to shrink the data portion of an Acquisition object down to 1 element. [0 elements is not possible since the realloc() behaviour is not well-defined and hence the ISMRMRD team excluded this possibility.]

Downsizing however, will reduce memory use by a factor of (#readout points) * (#coil channels) which is usually O(1000). Duplication of an empty dataset is therefore not a problem memory-wise. I have to confirm this behaviour though, when reverting my code to use the AcquisitionsVector::get_acquisition() method and let the testjobs run.