Allow Cimg::get_warp output CImg to be provided

[rtobar]

While doing some profiling on our R imager package that uses CImg underneath, I found that the fundamental CImg::get_warp operation allocates the output CImg object internally, then returns it. This prevents us from saving an allocation when returning the data back to R, which needs to be done on a different buffer.

Locally I have the following patch that allows me to avoid this, and I was wondering whether this would be something that would be considered as a desired feature upstream. In other words: would such a PR be welcome?

diff --git a/inst/include/CImg.h b/inst/include/CImg.h
index 9bb8997..59162fa 100644
--- a/inst/include/CImg.h
+++ b/inst/include/CImg.h
@@ -38609,8 +38609,26 @@ namespace cimg_library {
                                     "have different XYZ dimensions.",
                                     cimg_instance,
                                     p_warp._width,p_warp._height,p_warp._depth,p_warp._spectrum,p_warp._data);
 
       CImg<T> res(p_warp._width,p_warp._height,p_warp._depth,_spectrum);
+      do_warp(p_warp, res, mode, interpolation, boundary_conditions);
+      return res;
+    }
+
+    //! Warp image content by a warping field, with the user providing the output image \newinstance
+    template<typename t>
+    void do_warp(const CImg<t>& p_warp, CImg<t> &res, const unsigned int mode=0,
+                     const unsigned int interpolation=1, const unsigned int boundary_conditions=0) const {
+      if (is_empty() || !p_warp) return;
+      if (mode && !is_sameXYZ(p_warp))
+        throw CImgArgumentException(_cimg_instance
+                                    "warp(): Instance and specified relative warping field (%u,%u,%u,%u,%p) "
+                                    "have different XYZ dimensions.",
+                                    cimg_instance,
+                                    p_warp._width,p_warp._height,p_warp._depth,p_warp._spectrum,p_warp._data);
+      if (!res.is_sameXYZ(p_warp) || res._spectrum != _spectrum)
+        throw CImgArgumentException("warp(): output buffer has XYZ dimensions different warping field (%u,%u,%u) "
+                                    "or different spectrum than source image (%u)",
+                                    p_warp._width, p_warp._height, p_warp._depth, _spectrum);
 
       if (p_warp._spectrum==1) { // 1D warping
         if (mode>=3) { // Forward-relative warp
@@ -39422,7 +39440,6 @@ namespace cimg_library {
             }
         }
       }
-      return res;
     }
 
     //! Generate a 2D representation of a 3D image, with XY,XZ and YZ views.

The way I'm then using this is (very roughly) like this:

// inputs
unsigned int mode = ...;
unsigned int interpolation = ...;
unsigned int boundary_conditions = ...;
CImg<double> img = ...;
CImg<double> wrp = ...;
double *my_r_buffer = ...;

// Create a shared CImg object with the R buffer, and the `warp` result there
CImg<double> output {my_r_buffer, wrp.width(), wrp.height(), wrp.depth(), img.spectrum(), true};
img.do_warp(wrp, output, mode, interpolation, boundary_conditions);
return my_r_buffer;

[dtschump]

Well, if I understand well, you will get this kind of issues with all the CImg<T>::get_*() method, or is it specific to warp() ?

[rtobar]

I haven't really dug that deep to see if this is a general issue with all get_ methods, so I'd assume you're correct. I focused on warp because that's the one in particular we were more worried about in terms of performance.

The question still stands: would such a patch be welcomed? I wouldn't be able to provide a patch for all get_ methods, in case you worry about maintaining full consistency, but maybe ths could be a start. Maybe I could have also named do_warp as get_warp and make it an overload instead of bringing a new name into the mix.

[dtschump]

Considering that it is not really a issue with the get_warp() method, I doubt such a patch is the solution.
Maybe we can think about a more general (and elegant) solution that could apply to all get_*() methods instead.

[rtobar]

Thanks @dtschump for that feedback. I agree that this is a wider problem and not specific to the warp function.

OTOH I don't think there's a solution that will allow you to offer this functionality in a generic way for all functions, given how their outputs' sizes depend so heavily on the input values of each operation. If anything, I think the "generic" solution would be something along the lines of implementing three variants for each method, just like there are two at the moment:

CImg<T> &op(a, b, c); // inline
CImg<T> get_op(a, c, b); // new instance
void get_op(CImg<T> &out, a, c, b); // raw operation

This is obviously a big effort that I wouldn't be able to undertake, and probably neither would you without some clear incentive. If that's the situation, and you'd want to have an all-or-nothing solution, then I'll stick to my patch on our local copy; otherwise if you would consider adding this functionality on a function-by-function basis then I'd be happy to upstream this patch, and others that might come.

[dtschump]

I do actually think there is a generic solution.
Maybe this can be a good start. The example below apply some CImg methods on buffers manually allocated.
All temporary allocations that may happen will be freed by CImg destructors.

#include "CImg.h"
using namespace cimg_library;

int main(int argc, char **argv) {

  // Allocate memory for input image.
  const int width = 512, height = 512, depth = 1, spectrum = 3;
  const float *p_in = new float[width*height*depth*spectrum];

  // Allocate memory for output image.
  const float *p_out = new float[width*height*depth*spectrum];

  // Call a CImg method to fill values in 'p_in'.
  CImg<float>(p_in,width,height,depth,spectrum,true).rand(0,255);

  // Call a method so that p_out contains the resulting image.
  // All temporary allocations will be freed by CImg.
  CImg<float>(p_out,width,height,depth,spectrum,true) =
    CImg<float>(p_in,width,height,depth,spectrum,true).get_blur(10);

  // Display result.
  CImg<float>(p_out,width,height,depth,spectrum,true).display();

  // Free memory.
  delete[] p_in;
  delete[] p_out;

  return 0;
}

Wouldn't it be something like this that you are looking for?