🚀 go-pugleaf

In article <9209011949.AA19843@hp-and.an.hp.com>, susanm@hp-and.an.hp.com (Susan Margulies) writes:
> In recent posts I have seen some discussion about the phase
> information which results from a two-dimensional discrete Fourier
> Transform.  I am new to 2D Signal and Image Processing, and as such
> have my own question to throw into the discussion.
> 
> Basically, I am trying to come up with a way of visualizing the
> phase information as a function of position.  I am already
> displaying the magnitude information by mapping transform
> coefficients into grey scale information.  But since it is the phase
> information which is necessary (and sufficient) for image
> reconstruction, I would prefer to be studying the phase information.

One possibility would be to map phase information to hue in the
HSV color model.  The magnitude could be mapped to value in the HSV
model as well.  The result could be displayed as an image.

I haven't tried this myself, but I think the result could be interesting.

Pat

> From: susanm@hp-and.an.hp.com (Susan Margulies)

>Basically, I am trying to come up with a way of visualizing the
>phase information as a function of position. ...

Hello!

    I was interested in the same problem, and have heard other people
mension that displaying the phase as colour was effective.

    I have tried this, and it can produce some very pretty pictures.
A good colour mapping seems to be

  0 degrees -> red
 90 degrees -> blue
180 degrees -> green
270 degrees -> yellow

    As well as producing pretty pictures, this mapping has some nice
perceptual properties: colours that are 180 degrees apart are also
complementary as far as the human visual system is concerned, because
colour information gets into the brain as a luminance, blue/yellow
balance and green/red balance (I think!).

    I've put a pretty picture generated by this method in the
anonymous ftp area on our machine here; the directory is
/anonymous@csis.dit.csiro.au:pub/fft , and the two files are fft.ras.Z
(an eight bit rasterfile with phase/amplitude colurmap) and
input.ras.Z, which was the input file used to generate the pretty
picture.

    Note that to actually produce the fft image here, it is actually
necessary to translate the input image so that its centre is on the
origin, not in the middle of the space; if you perform an fft on the
untranslated image, you get a phase modulation of 180 degrees per
pixel, which smears out the colours and makes them impossible to see.

    Unfortunately, the phase examined in this way in any real (as in
real-life) image looks like gobbledygook.

-peter
--
Peter Fletcher,            Internet : peterf@csis.dit.csiro.au
PhD Student                Physical : CSIRO Division of Information Technology,
Phone    : +61-6-2750914              ANU, Acton, Canberra ACT AUSTRALIA
-------------------------------------------------------------------------------

Susan Margulies writes
> Basically, I am trying to come up with a way of visualizing the
> phase information as a function of position.  I am already
> displaying the magnitude information by mapping transform
> coefficients into grey scale information.

Map phase into color, or specifically hue. That is, map magnitude into L,
map phase into H, set S to 1.0 (or something close) and then convert
from HSL into RGB.

We did this and it sometimes worked. And sometimes didn't. Strongly
depends upon the sort of structures you are looking for.

Another to try is map magnitude into red and map phase into green. Some
shape information might come out better. Please note that the phase
map is circular to linear --- provide a rotation knob to alter the map
interactively to shift around the discontinuous point.

Good luck.

Mark R. Thomsen

In <9209011949.AA19843@hp-and.an.hp.com> susanm@hp-and.an.hp.com (Susan Margulies) writes:

>Basically, I am trying to come up with a way of visualizing the
>phase information as a function of position.  I am already
>displaying the magnitude information by mapping transform
>coefficients into grey scale information.  But since it is the phase
>information which is necessary (and sufficient) for image
>reconstruction, I would prefer to be studying the phase information.
>(I am doing graduate thesis research on image compression.)

One thing that you might want to take into account is that the Fourier
transform of an ordinary image is hermitic, ie if you do a transform
AxB real image -> AxB complex Fourier transform, you got quite some
redundancies.
Doing instead AxB image -> -(A/2+1)xB phase and (A/2+1)xB modulus
will give you a nicely viewable (and displayable) Fourier transform.
Also pseudocolor is no mean idea to display intensities ;-)

regards,
	spz
-- 
        spz is S.Petra Zeidler  Copyright 1992 S.Petra Zeidler
spz@specklec.mpifr-bonn.mpg.de  
 MPI fuer Radioastronomie Bonn  
    Auf dem Huegel 69, Germany  reading prescription: cum kilo salis

In article <9209011949.AA19843@hp-and.an.hp.com>, susanm@hp-and.an.hp.com (Susan Margulies) writes:
...
> Basically, I am trying to come up with a way of visualizing the
> phase information as a function of position.  I am already
> displaying the magnitude information by mapping transform
> coefficients into grey scale information.  But since it is the phase
> --
> Susan E. Margulies       
> HP Imaging Systems       
> 508-681-2153 (work)      
> susanm@hp-and.an.hp.com  


A couple of folks have responded with suggestions regarding using hue and
luminance (or value if you talk HSV) to show the two parameters.  We have
been working on multiparameter visualizations as applied to medical image
display (for instance showing metabolic and anatomic information in a single
image as opposed to two separate images).  The primary method I have
investigated for doing this is similar to what the other folks have
suggested, that being, using hue to encode one dimension and luminance the
other dimension.  This can be fairly easily done on standard unix
workstations (with 8 bit color, or better yet deeper framebuffers).  The
following is a summary of advice based on our current (limited) experience.

1) Do the lower frequencey signal with hue and the higher frequency signal with
luminance (humans don't perceive high frequency hue changes very
quantatively).  

2) Don't use low level luminance values (because you can't distinguish hues at
near black)

3) Provide interactive control.  For instance we allow the users to threshold
from 0-100% the hue mapping and 0-100% of the luminance mapping.  This
allows them to "fade" either dimension in and out.

4) A hue range from Red to Green works fairly well for a lot of visual
system reasons.  Avoid crossing over several different major colors (ie. no
rainbow hue attempts).  We are fairly picky on our choices (I use a
perceptually linear greyscale base on the luminance side, and at each a
luminance scale step a peceptually linear hue scale between two hues on the
hue side.  Note that you want to be careful in your hue choice offerings and
to offer some altervatives (as we were reminded by a couple of our RG
"colorblind" observers).

-Brad Hemminger
bmh@rad.unc.edu
Dept of Radiology
Medical Image Display Research Lab
University of North Carolina at Chapel Hill

In the referenced article spz@specklec.mpifr-bonn.mpg.de (S.P.Zeidler) writes:
>In <9209011949.AA19843@hp-and.an.hp.com> susanm@hp-and.an.hp.com (Susan Margulies) writes:
>
>One thing that you might want to take into account is that the Fourier
>transform of an ordinary image is hermitic, ie if you do a transform
>AxB real image -> AxB complex Fourier transform, you got quite some
>redundancies.

Serious question: so why not use the Hartley transform which does not
introduce these redundancies in the first place?

>Doing instead AxB image -> -(A/2+1)xB phase and (A/2+1)xB modulus
>will give you a nicely viewable (and displayable) Fourier transform.

... or use your favourite Fourier package and display the Hartley transform?

Regards, Peter

PS The Hartley transform can be regarded as a `tidier' Fourier transform.

   a) Uses cos(x)+sin(x) rather than cos(x)+i*sin(x).
   b) The inverse is the same operation as the transform.
   c) Takes real => real.
   d) Does not introduce redundancies into the representation.
   e) Even clearer relationship to the sine and cosine transforms.
   d) Easy to convert between Hartley and Fourier transforms.

In article <1992Sep7.110102.25536@bnr.co.uk> pgh@bnr.co.uk (Peter Hamer) writes:

   In the referenced article spz@specklec.mpifr-bonn.mpg.de (S.P.Zeidler) writes:
   >In <9209011949.AA19843@hp-and.an.hp.com> susanm@hp-and.an.hp.com (Susan Margulies) writes:
   >One thing that you might want to take into account is that the Fourier
   >transform of an ordinary image is hermitic, ie if you do a transform
   >AxB real image -> AxB complex Fourier transform, you got quite some
   >redundancies.

To be precise, half of your transform is redundant. If you pack things
carefully, you can even transform an N x M real array into an N/2 x M
complex array. If you are less aggressive in packing and waste a
little space, you can trivially get an (N+2)/2 x (M+2) complex output
array for your N x M real input.

   Serious question: so why not use the Hartley transform which does not
   introduce these redundancies in the first place?

I have tried both. It's relatively easy not to generate the half of
the Fourier transform that you don't need, and the result is the
familar and easy to use complex representation. It's messy and error
prone trying to rewrite the simple formulas (e.g., for filtering) in
the Fourier domain for the Hartley transform, and, as far as I can
tell, the result isn't any faster (if anything, it's likely to be
slower because you need more array indexing).

   >Doing instead AxB image -> -(A/2+1)xB phase and (A/2+1)xB modulus
   >will give you a nicely viewable (and displayable) Fourier transform.

   ... or use your favourite Fourier package and display the Hartley transform?

If you look at the Hartley transform of an image, you are looking at
both phase and magnitude information, interwoven in some visually
complicated way. The nice thing about the Fourier transform is the
separation into phase and magnitude.

   Regards, Peter

   PS The Hartley transform can be regarded as a `tidier' Fourier transform.

      a) Uses cos(x)+sin(x) rather than cos(x)+i*sin(x).
      b) The inverse is the same operation as the transform.
      c) Takes real => real.
      d) Does not introduce redundancies into the representation.
      e) Even clearer relationship to the sine and cosine transforms.

It looks neat on paper, but I don't think it's tidier in practice.

      d) Easy to convert between Hartley and Fourier transforms.

If you are going to compute the FHT and then reconstitute Fourier
coefficients (or phase/magnitude information) from the result, you are
probably better off using the Fourier transform in the first place.

					Thomas.

In article <1992Sep7.110102.25536@bnr.co.uk>
  pgh@bnr.co.uk (Peter Hamer) writes:
>Serious question: so why not use the Hartley transform which does not
>introduce these redundancies in the first place?

Well, in at least some fields, (mine, for instance), the fundamental observable
actually starts in the Fourier plane.  We measure Fourier coefficients
directly, and reconstruct the image from them.  And, as someone else already
mentioned, it's really not all that painful to go from a real full plane to
a complex half plane.  There is a lot of bookkeeping that must go into our
programs, no matter which transform we use.  Keeping track of transform
direction just isn't that much added complexity.  Finally, is or isn't the
FHT kosher to use these days?  When first it came out, Stanford was trying
to milk it for all it was worth.  I actually sent for one of their information
packets.  They wanted something like 5% of gross sales per copy or some hefty
minimum, for any program which used the transform.  That all by itself was
a good reason not to use the transform at the time, when it doesn't really
provide much more in the way of functionality than the FFT.

-- 
| Daniel Briggs  (dbriggs@nrao.edu)                      |        USPA B-14993
| New Mexico Tech / National Radio Astronomy Observatory |            DoD #387
| P.O. Box O / Socorro, NM 87801   (505) 835-7391        |
Support the League for Programming Freedom  (info from league@prep.ai.mit.edu)

In article <1992Sep8.133449.26552@zia.aoc.nrao.edu> dbriggs@zia.aoc.nrao.edu (Daniel Briggs) writes:

   direction just isn't that much added complexity.  Finally, is or isn't the
   FHT kosher to use these days?  When first it came out, Stanford was trying
   to milk it for all it was worth.  I actually sent for one of their information
   packets.  They wanted something like 5% of gross sales per copy or some hefty
   minimum, for any program which used the transform.  That all by itself was
   a good reason not to use the transform at the time, when it doesn't really
   provide much more in the way of functionality than the FFT.

GOOD QUESTION! I've been trying to find this out for a while, too. I'm not
sure if Stanford can lay claim to much other than the original algorithm
implementations that Bracewell published, and there have been other
equivalent implementations published (IEEE ASSP had one a few years back, I
can get the specifics if need be). BTW, there is a serious bug in the
BASIC FHT code Bracewell published in his book, anyway, so I always wondered
how that would affect any claims they tried to make on copyrights.

Has anybody got the lowdown on this stuff? It seems a shame for everybody
doing audio and other real-valued stuff to have to deal with ffts when
there are some advantages to fhts. I've got a lot of stuff I could
publish or even release pd if I knew for sure that Stanford wasn't going
to come after my butt later on...
--
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
* Eric Jacobsen        Honeywell Inc.                           *
|     1-602-436-3941   Commercial Flight Systems, Air Transport |
* HVN       436-3941   US Mail: P.O. Box 21111, MS 2Q34C1       *
| FAX 1-602-436-1310   Phoenix AZ 85036-1111                    |
* INet: ericj@hwcae.honeywell.com                               *
| UUCP: uunet!hwcae.honeywell.com!ericj@uunet.uu.net            |
* The opinions expressed are my own and therefore better        *
| thought out than my employer's.                               |
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

In article <TMB.92Sep7153447@arolla.idiap.ch> tmb@arolla.idiap.ch (Thomas M. Breuel) writes:

      Serious question: so why not use the Hartley transform which does not
      introduce these redundancies in the first place?

   I have tried both. It's relatively easy not to generate the half of
   the Fourier transform that you don't need, and the result is the
   familar and easy to use complex representation. It's messy and error
   prone trying to rewrite the simple formulas (e.g., for filtering) in
   the Fourier domain for the Hartley transform, and, as far as I can
   tell, the result isn't any faster (if anything, it's likely to be
   slower because you need more array indexing).

If done right you can get speedups with fhts, but they are not dramatic.
If there are symmetries in the inputs (not likely with most images), you
can get even more benefit with the fht.

      >Doing instead AxB image -> -(A/2+1)xB phase and (A/2+1)xB modulus
      >will give you a nicely viewable (and displayable) Fourier transform.

      ... or use your favourite Fourier package and display the Hartley transform?

   If you look at the Hartley transform of an image, you are looking at
   both phase and magnitude information, interwoven in some visually
   complicated way. The nice thing about the Fourier transform is the
   separation into phase and magnitude.

Ick. If you look at the raw output of an fht it is indeed gross. The
even and odd components of the output are interwoven and I've never
seen anybody do anything practical with it in this form. It is easy,
though to extract just the magnitude or just the phase from this
"Hartley vector", with no more computation than that required from
an fft output.

   If you are going to compute the FHT and then reconstitute Fourier
   coefficients (or phase/magnitude information) from the result, you are
   probably better off using the Fourier transform in the first place.

It may depend on what you're doing or your specfic implementation, but
I've had a lot of good experiences with the fht. Generally there is
a noticeable reduction in computational requirements with the fht.

BTW, I don't want to start any "transform wars", but the fft is
definitely not as universal a tool as it used to be.
--
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
* Eric Jacobsen        Honeywell Inc.                           *
|     1-602-436-3941   Commercial Flight Systems, Air Transport |
* HVN       436-3941   US Mail: P.O. Box 21111, MS 2Q34C1       *
| FAX 1-602-436-1310   Phoenix AZ 85036-1111                    |
* INet: ericj@hwcae.honeywell.com                               *
| UUCP: uunet!hwcae.honeywell.com!ericj@uunet.uu.net            |
* The opinions expressed are my own and therefore better        *
| thought out than my employer's.                               |
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*

In <1992Sep7.110102.25536@bnr.co.uk> pgh@bnr.co.uk (Peter Hamer) writes:

>Serious question: so why not use the Hartley transform which does not
>introduce these redundancies in the first place?
>... or use your favourite Fourier package and display the Hartley transform?

My favourite Fourier transform package is a Radix-8 FFT that does
512x512 float -> 256x512 complex in 4.3 sec on a DECStation 5000/200,
and was written by myself (that's why it is my favourite ;-), and
doesn't contain a Hartley transform.

Seriously, though: I need phase and modulus seperately, and it wouldn't
make sense to first compute 10,000 Hartley transforms followed by the
same number of phase/modulus separations, because that would take even
longer than shoving the stuff thru the optimized FFT.

regards,
	spz
-- 
        spz is S.Petra Zeidler  Copyright 1992 S.Petra Zeidler
spz@specklec.mpifr-bonn.mpg.de  
 MPI fuer Radioastronomie Bonn  
    Auf dem Huegel 69, Germany  reading prescription: cum kilo salis

> 
> My favourite Fourier transform package is a Radix-8 FFT that does
> 512x512 float -> 256x512 complex in 4.3 sec on a DECStation 5000/200,
> and was written by myself (that's why it is my favourite ;-), and
> doesn't contain a Hartley transform.

Now that we are on to 'favourite FFTs', mine is Lynn ten Eykes multi-radix FFT
which does FFTs of any sized multi-dimensioned array with dimensions that have 
prime factors less than 49 (X must be even, however) the output of a 200x200 
array of reals is 101x200 complex.  It is very fast, due to special code to 
handle factors of 2,3,4,5 & 8, in addition to the primes to 49.  We have 
ported it to IBM, VMS, Ultix and UNICOS (Cray).  Its good... (less than 3 
secs for 512x512 reals on a DECstation 5000/240)

+----------------------------------------------------------------------------+
|Ross Smith,  Research Computing Resource, Department of Cell Biology, NYU-MC|
|550 First Ave., NYC, 10016-6402.  Phone: (212) 263-5356: FAX: (212) 263-8139|
|E-Mail:  SMITH@NYUMED.BITNET (BITNET),   SMITH@MCCLB0.MED.NYU.EDU (Internet)|
+----------------------------------------------------------------------------+

In <1992Sep9.083621.7140@mcclb0.med.nyu.edu> smith@mcclb0.med.nyu.edu writes:
>> 512x512 float -> 256x512 complex in 4.3 sec on a DECStation 5000/200,
>Its good... (less than 3 secs for 512x512 reals on a DECstation 5000/240)

Seems that both routines keep the FP pipeline busy and use about the same
number of operations (the 5000/240 runs at 40MHZ, the 5000/200 runs at 25MHZ).

Regards,
-- 
Michael van Elst
UUCP:     universe!local-cluster!milky-way!sol!earth!uunet!unido!mpirbn!p554mve
Internet: p554mve@mpirbn.mpifr-bonn.mpg.de
                                "A potential Snark may lurk in every tree."

smith@mcclb0.med.nyu.edu writes:

>> 
>> My favourite Fourier transform package is a Radix-8 FFT that does
>> 512x512 float -> 256x512 complex in 4.3 sec on a DECStation 5000/200,
>> and was written by myself (that's why it is my favourite ;-), and
>> doesn't contain a Hartley transform.

>Now that we are on to 'favourite FFTs', mine is Lynn ten Eykes multi-radix FFT
>which does FFTs of any sized multi-dimensioned array with dimensions that have 
>prime factors less than 49 (X must be even, however) the output of a 200x200 
>array of reals is 101x200 complex.  It is very fast, due to special code to 
>handle factors of 2,3,4,5 & 8, in addition to the primes to 49.  We have 
>ported it to IBM, VMS, Ultix and UNICOS (Cray).  Its good... (less than 3 
>secs for 512x512 reals on a DECstation 5000/240)

Well, while we're at it.
My favorite FFT is a number theoretic transform which works modulo some
big prime number p at the end of your number range, where p-1 is dividable
by N. Uses modulo arithmetic instead of complex arithmetic, and makes
use of the fact that in these number systems Nth primitive roots of 1 exist
which are simple integers. (Base 2 FFTs work quite ok).

Digital filters implemented with these FFTs have no calculation rounding
errors whatsoever, and 16 bit are sufficient for almost all conceivable
picture filters. A pity signal processors are that bad at modulo arithmetic.

I have recieved several requests for the FFT code.  I have made a 
concatinated file called 'ten_eyke.txt' and put it in out anonymous area for 
people to FTP if they want.  For VMS people, take 'ten_eyke_pak.com' which is 
a VMS_SHARE file.  The concatinated file should be just fine, however, for 
people cursed with UNI*X: the code compiles without problem under Ultrix, and 
also on all the other platforms we have tried.

This routine, and a _WHOLE_LOT_MORE_ is a part of our image processing 
package called MDPP (Micrograph Data Processing Program) which runs on DEC 
workstations (either VAX/VMS or Ultrix).  It is a general purpose image 
processing package, originally targetted at problems in electron microscopy, 
but now very much more general in scope.  It requires X-windows for image
display, but there are versions for VWS and a small number of other displays, 
like the Lexidata and the Gould IP8500 and IP9000.  If you are interested 
in this you can order MDPP from the DECUS software Library which charges a 
nominal fee for handling and media.  From memory, the part number is V00535.

+----------------------------------------------------------------------------+
|Ross Smith,  Research Computing Resource, Department of Cell Biology, NYU-MC|
|E-Mail:  SMITH@NYUMED.BITNET (BITNET),   SMITH@MCCLB0.MED.NYU.EDU (Internet)|
+----------------------------------------------------------------------------+