Use sparse matrices for dose actor?

[djboersma]

I think that a large part of the memory foot print in radiotherapy simulations comes from the dose actor. Often we want the dose distribution with a high resolution but only a small number of voxels actually has a nonzero dose. But we also want to be able run multiple instances of Gate on a multicore computer, preferably on Gate instance per core, without having to require enormous amount of RAM. On a typical computer cluster a single core job is supposed to consume at most a few GiB. But with a big CT image and maximum dose actor information the Gate process easily blows up to much more than 5 GiB. There are other RAM hogs in Gate (GateSourceTPSPencilBeam source is one, but yesterday I submitted a pull request that hopefully solves that), but DoseActor is one of the biggest I think.

Has the Gate community every considered using "sparse matrices" instead of a data vector with entries for all voxels in a dose distribution? Am I too optimistic in thinking that we can save lots of RAM with sparse matrices?

(A sparse matrix would not work for the patient CT image, obviously...)

[dsarrut]

Hi, Sure that sparse matrices could be useful, but it is already possible to reduce the size of the dose actor by simple options. It is not necessary that the nb of dosels (scoring pixel for dose) is the same than the nb of pixels in the CT image. Use the setResolution macro to reduce it to the minimum required size and setPosition to put it in the right position inside the CT image. As example, you may have a 512x512x300 CT images, with a small 128x128x128 dose matrix included in the CT and centered at the right position. The two matrices (CT and dose) are independent. When using mhd file format (recommended) the origin of the dose matrix will be adapted from the CT origin and the setPosition such that VV (for example) displays everything correctly.

I think spare matrix could be a large work not only in Gate, but also for the analysis of the output (how to display ?). It could be useful and I can help if someone decide to go to it, but IMHO it is not a priority.

[djboersma]

So far I have used dose matrices with the same voxelization as the patient images because I want to do dose tracking, i.e. deform the dose distributions with deformation vector fields that link a reference image to the simulated "moving" image. This dose tracking may in fact be possible with a smaller dose matrix as well, I don't know. Already with the same size dose matrix some of the deformation vectors point outside of the image volume, so in principle the dose tracking code should already deal with that.

The purpose of using sparse matrices in Gate would be to reduce the memory usage while running the simulation. At the end of the simulation we could convert them to a normal format, so that there are no extra complications for analysis and visualization.

Developing sparse matrices from scratch would be time consuming. There are libraries/toolkits with mature and efficient implementations already available (e.g. SuiteSparse: http://faculty.cse.tamu.edu/davis/suitesparse.html). We could (c)make the use of this optional, so that normal users are not forced to install those extra dependencies.

[djboersma]

Hm. SuiteSparse actually only works with 2D matrices. Of course we can "flatten" a 3D matrix into a 2D matrix. A bit of a hack... I'll look a bit further.

[djboersma]

Of course ITK has a "sparse image" class, but then this would be available only to ITK-enabled builds. And moreover it would be nicer to have the "sparse image" class be more similar to the existing GateImage* classes. I am now experimenting with a very simple idea: a GateSparseImageT class that is very similar to GateImageT, except that the std::vector<PixelType> data member is replaced with a std::map<size_t,PixelType> and an PixelType mInsideValue. The latter contains the value of "boring" pixels (usually zero), the former contains pairs of (index,value) pairs for the "interesting" pixels.

[BishopWolf]

I agree Dose Image object shall be sparse matrix on memory, we already bundled a MetaIO and a dicom images from ITK, so having also sparse class shall not be difficult.

We shall start thinking in having ITK as a requirement, I dont think it is too hard to do that, we can avoid then to bundle its code, just put a link to the already tested version and thats it, teach users how to compile it just like geant4.

Sparse matrixes there, and matrixes in general are thread safe so we can rid out some classes we have that are not thread safe. This can be the way to follow for multithreading gate.

If you plan to have our own GateSparseImageT, then please make the reading and writing thread safe with flags already on them.

[BishopWolf]

You shall also have a private boolean minmaxset variable, a private PixelValue min and max and a private getminmax function. inside the public getmin and getmax, watch the state of minmaxset if true just return min or max, if false call getminmax function which assigns, this shall be the getminmax function

template<class PixelType>
void
GateSparseImageT<PixelType>::GetMinMax() {
    max = (data.empty() ? mInsideValue : std::max(mInsideValue,std::max_element(begin(), end(), second_less<PixelType>)->second));
    min =(data.empty() ? mInsideValue : std::min(mInsideValue,std::min_element(begin(), end(), second_less<PixelType>)->second));
    minmaxset = true;
}

template<class PixelType>
PixelType
GateSparseImageT<PixelType>::GetMaxValue() {
  if !(minmaxset) GetMinMax();
  return max;
}

template<class PixelType>
PixelType
GateSparseImageT<PixelType>::GetMinValue() {
  if !(minmaxset) GetMinMax();
  return min;
}

Also in all setValues you must put minmaxset to false.

[brenthuisman]

Just to pile on: ROOT has a THnSparse class, but I haven't had the time (or expertise) to work on that. I work with 4D images, so memory size is of interest to me.

[djboersma]

It is nice to have some options, thanks for all these suggestions. I think I can alter the dose actor such that it can be used with any of the sparse image classes (the sparse image from ITK, THnSparse from ROOT or my own hack with std::map), depending on availability of the needed libraries. I don't have time do this this week, maybe next week.