Added optional argument out in all gradient-computing methods

[evgueni-ovtchinnikov]

Changes in this pull request

Added optional argument out in all gradient-computing methods for compatibility with CIL.

Testing performed

Related issues

Fixes #1242.

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
• I have implemented unit tests that cover any new or modified functionality
• 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:

• 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.

[KrisThielemans]

This would work, and we can merge as-is.

However, it doesn't really provide an optimisation for the user. For instance, checking

SIRF/src/xSTIR/cSTIR/cstir.cpp

Lines 1246 to 1252 in 2c66faf

	Prior3DF& prior = objectFromHandle<Prior3DF>(ptr_p);
	STIRImageData& id = objectFromHandle<STIRImageData>(ptr_i);
	Image3DF& image = id.data();
	shared_ptr<STIRImageData> sptr(new STIRImageData(image));
	Image3DF& grad = sptr->data();
	prior.compute_gradient(grad, image);
	return newObjectHandle(sptr);

we will be creating a new image, then filling out with its content. Better would be to directly call STIR's compute_gradient with the out image. This could possibly be done like this

        if out is None:
            out = ImageData()
        grad.handle = pystir.cSTIR_priorGradient(self.handle, image.handle, out.handle)
        check_status(grad.handle)
        return out

accompanied by a check in the wrapper if the handle points to an "empty" image and if so, create a new one. Easy?

By the way, shared_ptr<STIRImageData> sptr(new STIRImageData(image)); creates a copy, while actually we only need a zero image, which would be faster to construct. Do we have get_uniform_copy(0) or similar on the C++ side?

[evgueni-ovtchinnikov]

I fixed it this way:

        if out is None:
            grad = ImageData() # does next to nothing
        else:
            grad = out
        grad.handle = pystir.cSTIR_priorGradient(self.handle, image.handle)
        check_status(grad.handle)
        return grad

[KrisThielemans]

How does this remove the copy? The C function will still do the fill into a temporary object (as it doesn't get the grad handle)?

[evgueni-ovtchinnikov]

I got rid of out.fill(grad). Replacing the image copy constructor with something like get_uniform_copy(0) (which we do not have in SIRF at C++ level - does STIR have it?), is a separate issue.

[KrisThielemans]

I might misunderstand the handle details, but as far as I can see

        grad.handle = pystir.cSTIR_priorGradient(self.handle, image.handle)

gets essentially rid of any allocated memory that out was pointing to (I guess the handle destructor will be called). In any case,

SIRF/src/xSTIR/cSTIR/cstir.cpp

Line 1249 in 2c66faf

shared_ptr<STIRImageData> sptr(new STIRImageData(image));

still has allocation of a new object, as opposed to reusing whatever was in out.

So, while this doesn't have the fill anymore, it still doesn't reuse memory out. I think that can only be done by passing the out.handle to the C function.

[evgueni-ovtchinnikov]

ok, see your point now, thanks!

trying this now:

        if out is None:
            out = ImageData()
        if out.handle is None:
            out.handle = pystir.cSTIR_priorGradient(self.handle, image.handle)
        else:
            assert_validities(image, out)
            pystir.cSTIR_computePriorGradient(self.handle, image.handle, out.handle)
        check_status(out.handle)
        return out

not sure I could manage to avoid copies with PLSPriorGradient though

[KrisThielemans]

can't check right now, but unfortunate to have the duplication. Cannot pass 0 or an empty handle (or default constructed one)?

we do need a test though.

[KrisThielemans]

Something like

void*
cSTIR_computeObjectiveFunctionGradient(void* ptr_f, void* ptr_i, int subset, void* ptr_g)
{
	try {
		ObjectiveFunction3DF& fun = objectFromHandle< ObjectiveFunction3DF>(ptr_f);
		STIRImageData& id = objectFromHandle<STIRImageData>(ptr_i);
               shared_ptr<STIRImageData> sptr;
		if (ptr_g==0)
		{
			Image3DF& image = id.data();
			STIRImageData* ptr_id = new STIRImageData(image);
			sptr.reset(ptr_id);
			ptr_g = newObjectHandle(sptr);
		}
		else
		{
			sptr = getObjectSptrFromHandle(ptr_g);
		}
		Image3DF& gd = sptr->data();
                ...
		return ptr_g;
}

Of course, there's large scope for using auto here.

[KrisThielemans]

not sure I could manage to avoid copies with PLSPriorGradient though

It's a bug Edit: it's a bad name. See #1248

[KrisThielemans]

Thanks! Please create an issue such that we don't forget to add tests, as well as reduce code duplication.