Hi,
this comment might be counterproductive to this discussion:
However.
I'm also using numpy as a basis for putting together an image analysis
"platform"/library.
This includes both 2D and higher dimensional images.
Since all my code, if not n Python, is written in C or C++ and not
Fortran, I decided early on that I had to get used to "invese
indexing", as in
image[y,x] or image[z,y,x]
"Inverse" only refers to the common (alphabetecally oriented)
intuition to say "at point x,y,z" rather than "z,y,x".
I also argueed that mathematical matrices (mostly / often) use an
index order as "row , column" which essentially is "y, x" (if the
matrix is seen as an image)
In any case, getting used to a "z,y,x" order saved me lots of
technical troubles, and after all, most code is or will be C/C++ and
not Fortran.
Hope I was not entirey misguided,
Sebastian Haase
On Nov 9, 2007 12:37 PM, Hans Meine <[EMAIL PROTECTED]> wrote:
> Am Freitag, 09. November 2007 00:16:12 schrieb Travis E. Oliphant:
> > C-order is "special" in NumPy due to the history. I agree that it
> > doesn't need to be and we have taken significant steps in that
> > direction.
>
> Thanks for this clarifying statement.
>
> > Right now, the fundamental problem is probably due to the
> > fact that the output arrays are being created as C-order arrays when the
> > input is a Fortran order array. Once that is changed then we can cause
> > Fortran-order inputs to use the simplest path through the code.
>
> That looks like a plan. I have read about __array_wrap__ in your book, and it
> sounds exactly like what we need for our imaging library; this way we can
> ensure that the result of operations on images is again an image. Here, it
> is crucial that the resulting image has Fortran order, too (for C/C++
> algorithms to agree with the Python idea of the images' orientation).
>
> > Fortran order arrays can be preserved but it takes a little extra work
> > because backward compatible expectations had to be met. See for example
> > the order argument to the copy method of arrays.
>
> What do you mean exactly (if you have something specific in mind at all)?
> Should the order be configurable for all operations, and default to "C" for
> backwards compatibility?
> (BTW: copy() has no arguments yet in my numpybook, page 57.)
>
> At last, I had a look at the relevant code for adding arrays; this is what I
> found: In core/src/umathmodule.c.src, the code templates for the relevant
> ufuncs (e.g. add) are defined, which will get expanded to e.g. DOUBLE_add.
> This will get wrapped into an ufunc in 'InitOperators' using
>
> f = PyUFunc_FromFuncAndData(add_functions, add_data, add_signatures, 18,
> 2, 1, PyUFunc_Zero, "add",
> "adds the arguments elementwise.", 0);
>
> and inserted as "add" into a dictionary which is later passed to
> PyArray_SetNumericOps. This will install the ufunc in the static
> NumericOps struct 'n_ops', which is used in the 'array_add' function
> installed in the PyTypeObject.
>
> Thus, it is the ufunc code (not suprisingly), which seems to call
> DOUBLE_add for the inner loop:
>
> | static void
> | DOUBLE_add(char **args, intp *dimensions, intp *steps, void *func)
> | {
> | register intp i;
> | intp is1=steps[0],is2=steps[1],os=steps[2], n=dimensions[0];
> | char *i1=args[0], *i2=args[1], *op=args[2];
> | for(i=0; i<n; i++, i1+=is1, i2+=is2, op+=os) {
> | *((double *)op)=*((double *)i1) + *((double *)i2);
> | }
> | }
>
> If I understood David correctly, he suggested an unstrided variant of
> this inner loop, which simply uses pointer increments. While this is
> a good idea (also probably quite some work), the real thing bugging me is
> that the above DOUBLE_add could (and should!) be called by the ufunc
> framework in such a way that it is equally efficient for C and Fortran
> arrays. (I.e. with steps set to sizeof(double), without prior
> copying/conversion.) It could even be used with negative step sizes when the
> input and output arrays have different orders. Disclaimer: So far, I did not
> look at the ufunc code yet, so I do not know which code paths exist.
>
> BTW: What are the "void *func" in core/src/umathmodule.c.src arguments
> for? Shouldn't that be "void * /*func*/", since the parameters are
> never used?
>
> --
> Ciao, / /
> /--/
> / / ANS
>
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>
_______________________________________________
Numpy-discussion mailing list
Numpy-discussion@scipy.org
http://projects.scipy.org/mailman/listinfo/numpy-discussion
