getting hsi to work with numpy
kfj
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:
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
Tom Sharpless
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.
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.
Cheers, Tom
kfj
>
> 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
JohnPW
Am I correct in thinking that HSI is not functional on the Mac
platform?
--John
kfj
> This sounds interesting.
> Am I correct in thinking that HSI is not functional on the Mac
> platform?
I suppose this is still so. Noone has got it to run on MacOS (on
neither of the various versions), let alone built whatever Mac users
need to install it :(
Kay