> I'm currently working on using numpy arrays with routines in hsi. I
> have managed to write code which succeeds in doing so, but it was
> quite involved insofar as I modified libpano, a bit of hugin code
> (PanoToolsInterface.c and -.h) and the interface definition file
> hsi.i, and I had to incorporate a bit of obscure (but free and
> documented) helper code. Interested? read on!

> Why numpy? If any of you have been using python for numeric computing,
> you are bound to have stumbled upon numpy:

> http://www.numpy.org/

> Numpy defines a data type called ndarray (short for n-dimensional
> array) which provides infrastructure for homogenous arrays of
> arbitrary dimensionality. To process data in such arrays, the most
> efficient method is to use vectorized C routines - the C code does all
> the looping, offsetting and adressing fast and efficiently, and the
> python code merely passes in the ndarray objects. Numpy merely
> provides the ndarray code, but it's the foundation for other modules
> which require this data type to operate.

> The module I'm intersted in using is scipy, short for scientific
> python. Scipy provides a collection of routines working on ndimages
> (you guessed it, it's n-dimensional images, so it's good for higher-D
> images as well, like videos or tomographs, but it will process plain
> old 2D images just as happily). This package has interesting stuff
> like B-spline interpolation, filtering, etc:

> http://docs.scipy.org/doc/scipy/reference/ndimage.html

> For an initial trial I contained data in numpy arrays and used python
> adressing and looping to feed the data to hsi (specifically to
> HuginBase::PTools::Transform::transform). Performance wasn't very
> good: on a test set of 9000000 coordinates, the procedure took some
> two minutes. So I tried to speed things up. First I created a
> vectorized version of execute_stack in libpano, like with this
> prototype:

> int v_execute_stack_new ( double* x_dest, // x coordinates destination
>                           double* y_dest, // y coordinates destination
>                           double* x_src,  // x coordinates source
>                           double* y_src,  // y coordinates source
>                           int count,      // number of coordinates
>                           void* params ); // pointer to vector of
> fDesc

> then I wrote a method for hugin's Transform object using this
> vectorized execute_stack routine, so I had an overloaded
> Transform::transform

>         bool transform( double* x_dest,
>                         double* y_dest,
>                         double* x_src,
>                         double* y_src,
>                         int count ) const;

> But how to get the numpy arrays in there? I found some magic code
> described in

> http://docs.scipy.org/doc/numpy/reference/swig.interface-file.html

> and the actual code at

> https://github.com/numpy/numpy/blob/master/doc/swig/numpy.i

> With these building blocks, I managed to SWIG-%extend hsi's interface
> to class Transform so I could pass in numpy arrays and the glue code
> would unpack the raw data from them and pass them on to the vectorized
> transform routine. Now when I ran the equivalent of my previous test,
> the code performed in under 18 seconds, so the speedup was
> considerable (like, sevenfold).

> I basically followed the method outlined in above-mentioned interface-
> file.htm in the section

> Beyond the Provided Typemaps

> and the solution it recommends for situations like mine when there's
> several ndarrays to pass in, ending up with this swig code in hsi.i:

> // we use the readymade numpy.i typemap for double arrays
> // on the signature of the %extended transform routine below:

> %apply (double* IN_ARRAY1, int DIM1) {(double* xd, int xd_size),
>                                       (double* yd, int yd_size),
>                                       (double* xs, int xs_size),
>                                       (double* ys, int ys_size)} ;

> %include <panotools/PanoToolsInterface.h>

> // KFJ 2011-11-26 extension to Tranform::transform to handle
> // data in numpy arrays with the vectorized version of execute_stack

> %extend HuginBase::PTools::Transform
> {
>     // vectorized version of transform, to transform
>     // many coordinates at once. note this signature only
>     // serves to allow the typemaps from numpy.i to take effect.

>     bool transform(double* xd, int xd_size,
>                    double* yd, int yd_size,
>                    double* xs, int xs_size,
>                    double* ys, int ys_size) const
>         {
>             if (    xd_size != ys_size
>                  || yd_size != ys_size
>                  || xs_size != ys_size )
>                {
>                    PyErr_Format(PyExc_ValueError,
>                                 "transform: Arrays of different
> lengths given");
>                    return 0;
>                }
>             // all arrays same size. proceed.
>             return $self->transform ( xd, yd, xs, ys, ys_size ) ;
>         } ;

> }

> I'm happy I managed to get this coded (it took me days of reading and
> experimenting - once I had the solution it was straightforward, but
> the SWIG programming is still kind of alien to me), so now I'm
> wondering how to go on about my discoveries.

> As I've mentioned, I modified libpano. These modifications were done
> because I found the vectorization was best done at the lowest level
> possible. Using an unmodified libpano and doing the vectorization in
> hugin would cost performance, but not too much - I'd estimate 50%, but
> I haven't timed it precisely.

> The most significant gain is by using vectorized C/C++ code at all, as
> compared to looping in pthon. Whichever way, the vectorized code hurts
> noone since it merely adds a (small) amount of C/C++ code on top of
> what is there already, so I reckon neither libpano nor hugin would
> object (it's even feasible that either project might find it useful).
> Pulling numpy.i into hsi produces a bit more bloat, but once it's
> there, interfacing to numpy arrays becomes quite straightforward,
> thereby providing an interface between hsi and numpy. This could be
> made optional somehow, but I'd even propose to consider making it
> integral and distributing numpy.i with hugin source.

> The code I wrote is currently only dealing most efficiently with
> contiguous data; modifying it to cater for strided data (which would
> be appropriate in handling numpy ndarrays) would be feasible, but
> require a bit more involved programming.

> Any comments/questions?

> Kay

> Given that Hugin now does Python, it is natural to extend some of the
> core functions of libpano to support numpy, and I'm glad to hear you
> are doing that.  This should open the Hugin platform to lots of
> interesting new developments.

> > Given that Hugin now does Python, it is natural to extend some of the
> > core functions of libpano to support numpy, and I'm glad to hear you
> > are doing that.  This should open the Hugin platform to lots of
> > interesting new developments.

> I've at least made a start. I suppose there might be a good deal of
> code in libpano which would lend itself to vectorization. I know very
> little of libpano yet; the place where I did my bit of vectorization
> resulted from my current need of mass-transforming coordinates, and
> when I dug into hugin's wrapping of the relevant panotools functions I
> found that execute_stack was my ideal target routine. The
> vectorization could have been done in the hugin code, but with a
> slight performance penalty. I've posted to the panotools-dev list to
> ask what they think of the idea.

> The vectorization of the stack execution and the interface to numpy
> are done at different levels, though. My modifications of libpano and
> hugin's panotools interface code merely introduce vectorized versions
> of execute_stack() and Transform::transform(), respectively, and are
> totally unaware of and independent from numpy. The numpy compatibility
> is introduced in the swig interface and is therefore only relevant to
> and usable from python code which imports hsi.

> > Since libpano13 is supposed to be an independent library, such
> > modifications should not depend on any supporting code in Hugin.  Of
> > course there would be no harm in making parts of Hugin depend on the
> > new libpano code.

> The libpano code would be extended, but no existing code in libpano
> would be modified. On my (Linux) system, after a recompile of libpano,
> the new routines were instantly available to hugin, which would now
> depend on the new routines in libpano to do it's vectorized
> Transform::transform, but, again, no existing hugin code would be
> modified, only some new code added. The changes to libpano and hugin
> would not introduce dependencies to new libraries or such.

> The changes to hsi.i (the swig interface file defineing the python
> interface) introduce the ability to deal with numpy arrays. There is
> one new dependency here, which is from numpy.i, which is a bit of
> source code, so can be distributed with the source. I think the
> resulting python module should still be independent of numpy - and if
> it is passed numpy arrays as arguments, numpy would have to have been
> imported in the first place.

> So far, my approach is quite ad-hoc; it's really only for the task at
> hand, but I feel that having found a way here is relevant and should
> make further investigation easier. I also hope this will open up new
> avenues. Scipy provides a lot of very advanced code, and it even has a
> fair bit of overlap with hugin/libpano - apart from the scipy.ndimage
> module there are other bits which might be exploitable in panorama
> programming. And a lot of other interesting python modules also use
> numpy arrays to handle data.

> On the other hand I found that the panotools code might be useful for
> other python code - it does some things really well and may even do
> some things no other tool can do so easily. If this takes off, it
> should be a win-win situation. I sometimes toyed with the idea to just
> wrap libpano with swig, but I reckon that huign has added a great lot
> of useful code on top of libpano, and therefore I stuck with the 'wrap
> libhuginbase' approach.

> Kay