On Feb 26, 3:31 pm, "Michael48" <m...@mirametrics.com> wrote:
> Yes Steve, I did read the manual and I don't "guess". I have a PhD from a
> top university and 20 years in the business. The point is that I made a
> hueristic inference and it turned out wrong.

When you write "heuristic inference," I read "guess".  The question
is, why did you make a "heuristic inference" from the prototype,
rather than following the explicit documentation?

You read the manual and you missed something that is clearly spelled
out several times.  This happens to the best of us.  But I'm not sure
the blame lies entirely with the manual.

> And speaking of hueristics,
> why even use the "x" in "nx", or similarly for the other axes?

Because 2d and 3d FFTs most commonly correspond to spatial dimensions,
and that is how I'm used to thinking of them.

> Are you not
> thinking of Cartesian x, y, and z axes? If so, then which direction is
> conventional for the x axis?

In row-major order, the first (discontiguous) dimension is
conventional for the x axis, since the dimensions are usually listed
from left to right alphabetically.   How the data is stored in memory
has nothing to do with how the data is displayed on the screen.

> And then which direction would "nx" reference
> the count of? I, for one, don't think it is unambiguous.

I never claimed it was unambiguous.   I don't think that there *is* an
unambiguous naming convention.  I think you have no choice but to read
the manual and see what it defines the data layout to be.

> You describe the
> array as "n1 x n2 x n3... in row major order". But the prototypes use
> nx,ny,nz, and not n1,n2,n3. Where's the association of nx to n1, or n3 to
> nx, or whatever? Is it right alongside?

Yes, it's right alongside.  If you look in the reference manual, it
states that the array is "nx x ny x nz" with "row-major order".  In
the tutorial section, it is also right alongside, stating, "All of
these transforms operate on contiguous arrays in the C-standard row-
major order, so that the last dimension has the fastest-varying index
in the array."

> Here's my suggestion: How about scrapping "nx,ny,nz" altogether for the
> function prototypes and using instead "n1,n2,n3" ? Then state simply that
> n1 refers to the most slowly varying subscript.

That's a reasonable, constructive suggestion.  Thank you.  I don't
think there's anything wrong with calling them "nx, ny, nz", but "n0,
n1, n2" has the advantage of using a nomenclature consistent with the
arbitrary-dimensional case.

Regards,
Steven G. Johnson

Yes Steve, I did read the manual and I don't "guess". I have a PhD from a
top university and 20 years in the business. The point is that I made a
hueristic inference and it turned out wrong. And speaking of hueristics,
why even use the "x" in "nx", or similarly for the other axes? Are you not
thinking of Cartesian x, y, and z axes? If so, then which direction is
conventional for the x axis? And then which direction would "nx" reference
the count of? I, for one, don't think it is unambiguous. You describe the
array as "n1 x n2 x n3... in row major order". But the prototypes use
nx,ny,nz, and not n1,n2,n3. Where's the association of nx to n1, or n3 to
nx, or whatever? Is it right alongside?

Here's my suggestion: How about scrapping "nx,ny,nz" altogether for the
function prototypes and using instead "n1,n2,n3" ? Then state simply that
n1 refers to the most slowly varying subscript.

Michael

>On Feb 26, 12:37 pm, "Michael48" <m...@mirametrics.com> wrote:
>> I really do appreciate your work of FFTW, and I understand how array
>> ordering works in Fortran vs C. But that wasn't my point. My point was
>> that your doc explicltly states the parameter passing order as "foo(
nx,
>> ny, nz, ..."  So my mistake was to associate (nx) with the number of
>> columns and (ny) with number of rows. If you are looking at a 2-D image
on
>> a screen, then don't you think that is somewhat hueristic? If I have a
C
>> program, with its multi-D arrays in a native row-major order, and my
code
>> calls your C-library and the docs show the parameter order as nx, ny,
nz,
>> then can you see that there's room for confusion at the implementation
>> end?
>
>There is room for confusion whether we list them as (nx,ny,nz) or
>(nz,ny,nx): that's why you have to read the manual instead of
>guessing.
>
>In English, it is conventional to name the dimensions "x", "y", "z"
>listed from left to right.  It is also conventional in discussing row-
>major order to read the dimensions from left to right, with the right-
>most being the contiguous dimension.
>
>Conversely, it's not reasonable to expect FFTW's data formats to be
>determined by how you happen to display your 2d data on a screen.
>
>The fact is, when you use any library that processes multidimensional
>data, you have to find out how the data must be stored---it is foolish
>to guess.   And we say *over and over* in the manual that the data is
>in row-major order.  In the tutorial.  In the reference section.  In
>the section on multi-dimensional array formats.  As far as I can tell,
>in every place the dimensions are discussed.
>
>I'm not sure what you expect us to do.  If we label the dimensions
>(nz,ny,nx), then at *least* as many people will be confused in the
>opposite direction.  And for people who just read the function
>prototype and don't read the accompanying explanations in the manual,
>nothing we write will make any difference...
>
>Regards,
>Steven G. Johnson
>
>

On Feb 26, 12:37 pm, "Michael48" <m...@mirametrics.com> wrote:
> I really do appreciate your work of FFTW, and I understand how array
> ordering works in Fortran vs C. But that wasn't my point. My point was
> that your doc explicltly states the parameter passing order as "foo( nx,
> ny, nz, ..."  So my mistake was to associate (nx) with the number of
> columns and (ny) with number of rows. If you are looking at a 2-D image on
> a screen, then don't you think that is somewhat hueristic? If I have a C
> program, with its multi-D arrays in a native row-major order, and my code
> calls your C-library and the docs show the parameter order as nx, ny, nz,
> then can you see that there's room for confusion at the implementation
> end?

There is room for confusion whether we list them as (nx,ny,nz) or
(nz,ny,nx): that's why you have to read the manual instead of
guessing.

In English, it is conventional to name the dimensions "x", "y", "z"
listed from left to right.  It is also conventional in discussing row-
major order to read the dimensions from left to right, with the right-
most being the contiguous dimension.

Conversely, it's not reasonable to expect FFTW's data formats to be
determined by how you happen to display your 2d data on a screen.

The fact is, when you use any library that processes multidimensional
data, you have to find out how the data must be stored---it is foolish
to guess.   And we say *over and over* in the manual that the data is
in row-major order.  In the tutorial.  In the reference section.  In
the section on multi-dimensional array formats.  As far as I can tell,
in every place the dimensions are discussed.

I'm not sure what you expect us to do.  If we label the dimensions
(nz,ny,nx), then at *least* as many people will be confused in the
opposite direction.  And for people who just read the function
prototype and don't read the accompanying explanations in the manual,
nothing we write will make any difference...

Regards,
Steven G. Johnson

Steve,

I really do appreciate your work of FFTW, and I understand how array
ordering works in Fortran vs C. But that wasn't my point. My point was
that your doc explicltly states the parameter passing order as "foo( nx,
ny, nz, ..."  So my mistake was to associate (nx) with the number of
columns and (ny) with number of rows. If you are looking at a 2-D image on
a screen, then don't you think that is somewhat hueristic? If I have a C
program, with its multi-D arrays in a native row-major order, and my code
calls your C-library and the docs show the parameter order as nx, ny, nz,
then can you see that there's room for confusion at the implementation
end? 

Michael

>On Feb 26, 10:10 am, "Michael48" <m...@mirametrics.com> wrote:
>> I do appreciate the difference between the order of parameter passing
and
>> whether or not to transpose twice. I know you can't speak for the FFTW
>> authors, but why on earth would the documentation specify "nx,ny, nz"
>> *everywhere*, when it really should always be "nz, ny, nx" ?? That is
just
>> so arcane! Which is why it cost me most of 2 days of time and
frustration
>> before I realized there was something I was simply not getting about
using
>> this library. I am grateful to you for pointing me the right way!
>
>I am one of the FFTW authors, and I can point out that FFTW's
>dimension ordering conforms to the standard for C programs and for
>most other modern programming languages other than Fortran.
>
>In a C program, if you have an array[nx][ny], it is stored with the ny
>dimension contiguous, and this is specified as an "nx, ny" array in
>FFTW.  This is called "row-major order" for an nx by ny array.
>
>The other convention is "column-major" order, in which an nx by ny
>array is stored with the nx dimension contiguous.  This is used in
>Fortran for an array(nx,ny), and in various other programs derived
>from Fortran (such as Matlab).  A column-major array can be specified
>in FFTW merely by passing the dimensions in reverse order.
>
>We can't help it that there are two common conventions for laying out
>multidimensional arrays in memory; we have to pick one ordering, and
>it is a C library so we pick the C convention.  I will note, however,
>that there is an *entire section* of the FFTW manual discussing
>multidimensional array formats in great detail.  I'm not sure what
>more you can expect.
>
>Regards,
>Steven G. Johnson
>
>

On Feb 26, 10:10 am, "Michael48" <m...@mirametrics.com> wrote:
> I do appreciate the difference between the order of parameter passing and
> whether or not to transpose twice. I know you can't speak for the FFTW
> authors, but why on earth would the documentation specify "nx,ny, nz"
> *everywhere*, when it really should always be "nz, ny, nx" ?? That is just
> so arcane! Which is why it cost me most of 2 days of time and frustration
> before I realized there was something I was simply not getting about using
> this library. I am grateful to you for pointing me the right way!

I am one of the FFTW authors, and I can point out that FFTW's
dimension ordering conforms to the standard for C programs and for
most other modern programming languages other than Fortran.

In a C program, if you have an array[nx][ny], it is stored with the ny
dimension contiguous, and this is specified as an "nx, ny" array in
FFTW.  This is called "row-major order" for an nx by ny array.

The other convention is "column-major" order, in which an nx by ny
array is stored with the nx dimension contiguous.  This is used in
Fortran for an array(nx,ny), and in various other programs derived
from Fortran (such as Matlab).  A column-major array can be specified
in FFTW merely by passing the dimensions in reverse order.

We can't help it that there are two common conventions for laying out
multidimensional arrays in memory; we have to pick one ordering, and
it is a C library so we pick the C convention.  I will note, however,
that there is an *entire section* of the FFTW manual discussing
multidimensional array formats in great detail.  I'm not sure what
more you can expect.

Regards,
Steven G. Johnson

Hendrik,

Thanks for your help. Now I get the forward/inverse pair for the 1042x1032
array in about 0.6 seconds and the result is clean to about 1 part in 10^14
or better. It is a significant improvement over my old FFT code.

I do appreciate the difference between the order of parameter passing and
whether or not to transpose twice. I know you can't speak for the FFTW
authors, but why on earth would the documentation specify "nx,ny, nz"
*everywhere*, when it really should always be "nz, ny, nx" ?? That is just
so arcane! Which is why it cost me most of 2 days of time and frustration
before I realized there was something I was simply not getting about using
this library. I am grateful to you for pointing me the right way!

Michael

>Michael48 wrote:
> > There's clearly omething I am missing. The docs say, for every
example,
> > that the arguments are passed as nx, ny, nz. So, I guess that nx does
not
> > mean "number of columns" and ny does not mean "number of rows" ? I am
> > confused: Should pass the args in the order nx, ny, nz but transpose
the
> > data buffer before r2c and after c2r?  Or should I reverse the order
of
> > parameter passing and leave the data array intact?
>
>Transposing a 1MP image takes some time.
>Passing the args in a different order doesn't cost anything.
>
> > For the in-place case, I get great precision:
>
> > Here's the remaining wrinkle: If I use an auxiliary buffer rather
than
> > in-place, then there's problem in how the data are packed and
unpacked,
> > from the point of view of the column count.
>
>If you do a R2C FFT out-of-place, the memory strides (# of bytes
>between two line starts) of input and output image are different,
>see
>
>http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html
>
>
>Hendrik vdH
>

Michael48 wrote:
 > There's clearly omething I am missing. The docs say, for every example,
 > that the arguments are passed as nx, ny, nz. So, I guess that nx does not
 > mean "number of columns" and ny does not mean "number of rows" ? I am
 > confused: Should pass the args in the order nx, ny, nz but transpose the
 > data buffer before r2c and after c2r?  Or should I reverse the order of
 > parameter passing and leave the data array intact?

Transposing a 1MP image takes some time.
Passing the args in a different order doesn't cost anything.

 > For the in-place case, I get great precision:

 > Here's the remaining wrinkle: If I use an auxiliary buffer rather than
 > in-place, then there's problem in how the data are packed and unpacked,
 > from the point of view of the column count.

If you do a R2C FFT out-of-place, the memory strides (# of bytes
between two line starts) of input and output image are different,
see

http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html


Hendrik vdH

Michael48 wrote:
> There's clearly omething I am missing. The docs say, for every example,
> that the arguments are passed as nx, ny, nz. So, I guess that nx does not
> mean "number of columns" and ny does not mean "number of rows" ? I am
> confused: Should pass the args in the order nx, ny, nz but transpose the
> data buffer before r2c and after c2r?  Or should I reverse the order of
> parameter passing and leave the data array intact?

Transposing a 1MP image takes some time.
Passing the args in a different order doesn't cost anything.

> For the in-place case, I get great precision:

> Here's the remaining wrinkle: If I use an auxiliary buffer rather than
> in-place, then there's problem in how the data are packed and unpacked,
> from the point of view of the column count.

If you do a R2C FFT out-of-place, the memory strides (# of bytes
between two line starts) of input and output image are different,
see

http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html

Another thing is that in out-of-place transforms, the input array
data will generally not be preserved unless you don't specify the
FFTW_PRESERVE_INPUT planner flag.
Maybe you just compared to you overwritten input array.


Hendrik vdH

That fixed the major problem but there is one remaining wrinkle to work
out. Thank you!!!

There's clearly omething I am missing. The docs say, for every example,
that the arguments are passed as nx, ny, nz. So, I guess that nx does not
mean "number of columns" and ny does not mean "number of rows" ? I am
confused: Should pass the args in the order nx, ny, nz but transpose the
data buffer before r2c and after c2r? Or should I reverse the order of
parameter passing and leave the data array intact? The 2D image is stored
as a single buffer with the column number (nx) varying most rapidly. And I
am calling the C libraries from C code. Something is wrong here!!

For the in-place case, I get great precision:

   (forward >> inverse DFT) - orignal < ~10^-18.

Here's the remaining wrinkle: If I use an auxiliary buffer rather than
in-place, then there's problem in how the data are packed and unpacked,
from the point of view of the column count. Maybe this is a transpose/arg
ordering issue like I mentioned above. I need to investigate this some
more.

Thank you again, for getting me over this first hurdle.


Michael

>Michael48 wrote:
>> I have been struggling to figure out some strange results from FFTW. I
get
>> a angled banding in the DFT and the inverse gives an image that looks
>> generally OK except that the signal varies periodically, consistent
with
>> those bands, by up to 1% above and below the original image signal.
>> 
>> I have put some picture and description here: 
>> http://www.mirametrics.com/pub/fftw_strange_results.htm
>
>Looks as if you got your memory management/addressing wrong.
>My first idea is that you gave FFTW the dimensions in the wrong
>order. FFTW wants the dimension with consecutive memory as
>the last dimension. So if your image is 1042x1032 and the
>1042 pixels in a line are in consecutive memory locations,
>you've got to plan a 1032x1042 FFT.
>
>Maybe you could show some source?
>
>
>Hendrik vdH
>
>

"Hendrik van der Heijden" <hvdh@gmx.de> wrote in message 
news:45da252d$0$6392$9b4e6d93@newsspool2.arcor-online.net...
> Michael48 wrote:
>> I have been struggling to figure out some strange results from FFTW. I 
>> get
>> a angled banding in the DFT and the inverse gives an image that looks
>> generally OK except that the signal varies periodically, consistent with
>> those bands, by up to 1% above and below the original image signal.
>>
>> I have put some picture and description here: 
>> http://www.mirametrics.com/pub/fftw_strange_results.htm
>
> Looks as if you got your memory management/addressing wrong.
> My first idea is that you gave FFTW the dimensions in the wrong
> order. FFTW wants the dimension with consecutive memory as
> the last dimension. So if your image is 1042x1032 and the
> 1042 pixels in a line are in consecutive memory locations,
> you've got to plan a 1032x1042 FFT.
>
> Maybe you could show some source?
>
>
> Hendrik vdH
>

It might also help to see the original image.  Starting from the DFTs is 
always questionable.

Fred