PRO IMAGE_POWER, fileName, spectrum, job, $
                             XSIZE=xSize,     YSIZE=ySize, $
                             PS=PS_Print,     HOLD=hold,        EPS=eps, $ 
                             NOREBIN=norebin, PRIM=Print_Image, AXIS=axis, $
                             NEGATIV=negativ, NONOISE=noiseLevel, $
                             INFO=Get_Info,   BATCH=batch,      NOTRIM=notrim, $
                             PHASE=phaseSp,   ADD=addphase,     CLEAN=cleanSp, $
                             SHOW3D=show3D,   SAVE=save,        WINDOW=filter
;===============================================================================
; (c) 1990 - khd c/o FU BERLIN.      A l l    r i g h t s    r e s e r v e d. ||
;      --    No part of this software package may be reproduced, transmitted, ||
;     1991   transcribed,  stored in a retrieval system,   or translated into ||
;            any form by any means without the   w r i t t e n  permission of ||
;            Karl-Heinz Dittberner-FU BERLIN, Arnimallee 22, D-1000 Berlin 33 ||
;===============================================================================
; PRO:IMAGE_POWER.PRO                             K.-H. Dittberner - 29.JUL.1990
; V 01.26a                                        P 476/16         - 24.JUL.1992
;
; IDL(V2.0)-Routine:   Logarithmisches Powerspektrum und Phasenspektrum eines 
;                      (Pixel-)Bildes, das als Bilddaten-File im TIFF-Format 
;                      (Industrienorm) vorliegen muss.
;
;          HINWEISE:   Um eine weitere Bearbeitung des Spektrums (z.B. Filtern)
;                      zu erleichtern, wird die VAR spectrum (Feld mit komplexen
;                      Zahlen) sowie ueber die VAR job der Name des lfd. Jobs
;                      an das rufende Programm (evtl. auch eine interaktive
;                      Sitzung) uebergeben. Durch Angabe des Keywords /SAVE kann
;                      die VAR spectrum zusaetzlich in einer File gesichert wer-
;                      den. Ausserdem werden die VAR power (Pixel-Bild des 
;                      Powerspektrums) und die VAR image (Ausschnitt des Bildes)
;                      in neuen TIFF-Files gespeichert.
;
;          KEYWORDS:   
;           /ADD     = Kann zusaetzlich zum Keyword /PHASE angegeben werden;
;                      dann wird das Phasenspektrum auf einem 2. Blatt ausgege-
;                      ben; enthaelt implizit /PHASE.
;           /AXIS    = Im Spektrum wird zusaetzlich ein Achsenkreuz durch den 
;                      Nullpunkt ausgegeben.
;           /BATCH   = Muss angegeben werden, wenn dieses Programm im IDL Batch-
;                      Mode laufen soll.
;           /CLEAN   = Die VAR spectrum wird gemaess dem gewaehlten NONOISE-
;                      Pegel bereinigt, d.h. bei den Frequenzen fuer die 
;                      power < NONOISE ist, werden die komplexen Werte => 0
;                      gesetzt. Das Phasenspektrum wird ggfs. entsprechend
;                      bereinigt.
;           /EPS     = Die PostScript-Ausgabe erfolgt als Encapsulated Post-
;                      Script-File (EPS-File) unter dem Namen "job.EPS"; keine
;                      Ausgabe auf dem Postscript-Drucker. EPS-Files werden 
;                      ausschliesslich fuer den Import in andere Applikationen 
;                      verwendet; z.B. fuer den PageMaker auf einem MAC.
;           /HOLD    = Die PostScript-File wird zunaechst nicht gedruckt.
;                      Sie wird unter dem Namen "job.PS" abgelegt.
;           /INFO    = Es werden Informationen ueber die gelesene Bilddaten-File
;                      ausgegeben (zur Fehlersuche).
;           /NEGATIV = Die Ausgabe des Powerspektrums erfolgt in Negativ-Dar-
;                      stellung; enthaelt implizit NONOISE.
;            NONOISE = Das Powerspektrum soll fuer die Darstellung aufbereitet
;                      werden. Implementiert ist z.Zt. (siehe Pkt.6) ein ein-
;                      faches Abschneiden des Rauschhintergrunds. Angegeben 
;                      wird der Rauschpegel (noiseLevel) in dB (Dezibel);
;                      z.B. NONOISE=-40. Defaultwert ist -60 dB.
;           /NOREBIN = Im Pkt.3 wird die Bildgroesse nicht halbiert.
;           /NOTRIM  = Muss/sollte angegeben werden wenn das Bild bereits mit
;                      einem anderen Programm beschnitten worden ist; z.B. mit
;                      TRIM_IMAGE oder CONVERT_TO_TIFF. Das Bild wird dann in
;                      diesem Programm nicht mehr beschnitten.
;           /PHASE   = Statt des Powerspektrums wird das Phasenspektrum des 
;                      Bildes ausgegeben mit phi=-Pi => SW und phi=+Pi => WE.
;           /PRIM    = PostScript-Ausgabe des Spektrums  m i t  dem Bild/Foto-
;                      (ausschnitt).
;           /PS      = Nur PostScript-Ausgabe des Spektrums; keine Ausgabe des 
;                      Spektrums auf dem Grafik-Terminal.
;           /SAVE    = Das unbearbeitete komplexe Spektrum (VAR spectrum) wird
;                      in einer File "job.SPECTRUM" im IDL Save-Format fuer
;                      eine spaetere IDL Restore-Operation gesichert.
;           /SHOW3D  = Das Powerspektrum wird als 3D-Oberflaeche ausgegeben. 
;                      (Z.Zt. noch nicht implementiert).
;            WINDOW  = Hiermit kann die Art der Vorfilterung des Bildes unter
;                      Pkt.4 beeinflusst werden. Zugelassene Window-Filter:
;                      o  WINDOW = "BH" :  BLACKMANN-HARRIS Window mit -74 dB.
;                      o  WINDOW = "HAN":  HANNING Window mit -32 dB.
;                      o  WINDOW = "NO" :  Es erfolgt keine Vorfilterung.
;                      Defaultwert ist "BH".
;            XSIZE   = Festlegen der Groesse des Bildausschnitts in
;            YSIZE     pixel. Default: 256 pixel x 256 pixel. 
;===============================================================================
; Liste der verwendeten Programme der Abt. Wiss. Datenverarbeitung am IfP-FUB:
; ----------------------------------------------------------------------------
; o  P476/2   PRINT_PLOT      PostScript-Ausgabe der aktuellen Grafik-File.
; o  P476/9   BH_WINDOW       BLACKMANN-HARRIS Windows fuer 1D und 2D-Anwendung.
; o  P476/13  TIFF_READ       Lesen von Bilddaten im TIFF-Format.
; o  P476/14  TIFF_WRITE      Schreiben von Bilddaten im TIFF-Format.
; o  P476/15  TRIM_IMAGE      Interaktives Zuschneiden eines Pixel-Bildes.
; o  P476/21  HAN_WINDOW      HANN(ing) Windows fuer 1D und 2D-Anwendung.
;===============================================================================
; Quelle:  Do-it-yourself!      Merkposten: o  Keyword /SHOW3D behandeln ?
;          C90 234                          o  Menu-Technik ?
;                                           o  2. Y-Achse mit i=0..ny/2 ?
; Aufruf:  z.B. IMAGE_POWER, "FOTO.TIF",    o  
;               spektrum, job [, ...]       o  Bei /AXIS Zahlen weg; wie ?
;                                           o  Plot_Shift, Plot_VSize allgem. ?
;                                           o  
;                                           o  Keyword /PCFILE einfuehren ?
;===============================================================================
; Aenderungen: [10/8/90-khd] => V 01.3:
;                 - Ueberarbeitet und Kommentierung zugefuegt.
;              [13/8/90-khd] => V 01.4:
;                 - Achsenbeschriftung fuer power in Pkt.8 zugefuegt.
;                 - REBIN des Bildes in Pkt.3 eingefuehrt. Diese Option kann
;                   durch das Keyword /NOREBIN abgeschaltet werden.
;              [15/8/90-khd] => V 01.5:
;                 - Bei der Ausgabe des Powerspektrums mit Achsen (Pkt.7+8) ist
;                   zur Erzielung der Device-Unabhaengigkeit eine besondere
;                   Technik des Ueberlagerns des Bildes mit der Ausgabe der
;                   Achsen (inkl. Beschriftung) erforderlich!
;              [16/8/90-khd] => V 01.6:
;                 - Keyword /PRIM zum optionalen Ausdruck des Bildausschnittes
;                   unter dem Powerspektrum eingefuehrt, noch nicht implemen- 
;                   tiert; das funktioniert nur bei einer PostScript-Ausgabe.
;              [17/8/90-khd] => V 01.7:
;                 - Frage nach PS-Druck im Pkt.10 zugefuegt.
;                 - Ausgabe von Titel + SubTitel mit XYOUTS auf DATA-Koordinaten
;                   umgestellt.
;                 - Bei PS-Ausgabe in Pkt.8 mit DEVICE auch XSIZE + YSIZE ge- 
;                   setzt, sodass das origin. Seitenverhaeltnis erhalten bleibt.
;              [20/8/90-khd] => V 01.8:
;                 - In Pkt.8 Berechnung von XSIZE + YSIZE geaendert wg. Bild-
;                   groesse durch Normierung auf y1/Plot_VSize; bisher auf
;                   (y1-y0)/Plot_VSize normiert.
;                 - CONST rand=100 => 50 reduziert, was zur Beschriftung reicht.
;              [21/8/90-khd] => V 01.9:                           (Copy => PRZY)
;                 - Ausdruck des Bildausschnitts in Pkt.9 voll implementiert.
;                 - Bildausschnitt wird als neue TIFF-File ausgegeben.
;                 - Aufbereitung des Powerspektrums (Pkt.6) durch Einfuehren des
;                   Keywords /NONOISE vorbereitet.
;                 - ON_ERROR,2 aktiviert und Kommentierung ergaenzt.
;              [22/8/90-khd] => V 01.10:                 (Copy => SEIDLER, PRZY)
;                 - Pkt.8+9 zur Erzielung einer noch besseren Device-Unabhaen-
;                   gigkeit (Normierung, Seitenverhaeltnisse, Sizing, etc.) gem.
;                   IMAGE_CONT.PRO aus User-Library ueberarbeitet.
;                 - Keyword /BATCH ausgetestet.
;              [23/8/90-khd] => V 01.11:              (Copy => AG GRUE + AG KOE)
;                 - In Pkt.6 ein einfaches Abschneiden (optional) des Rausch-
;                   hintergrunds implementiert.
;                 - Ausgabe des Creation-Dates unterm Titel zugefuegt (Pkt.8).
;                 - Keyword /NEGATIV eingefuehrt.
;              [28/8/90-khd] => V 01.12:                  (Copy => SEIDLER + DK)
;                 - Modifikation des Abschneidens des Rauschhintergrunds; dazu:
;                 - In Pkt.5 Umstellung der Berechnung des Powerspektrums auf
;                   Dezibel (dB); Bezugswert ist Wert bei f=0 cycles/pixel und
;                 - das Keyword NONOISE mit VAR noiseLevel definiert.
;                 - Keyword /HOLD eingefuehrt (Pkt.10).
;                 - Vorbereitung der Implementation der Ausgabe des Powerspek-
;                   trums als EPS-File.
;              [28/8/90-khd] => V 01.13:
;                 - Default-Wert von NONOISE=-80 dB => -60 dB geaendert.
;              [11/9/90-khd] => V 01.14:                       (Copy => SEIDLER)
;                 - Notiz zum C-Right zugefuegt.
;                 - Keyword /SHOW3D eingefuehrt, aber noch nichts programmiert.
;                 - Dazu neuen Pkt.9 "Ausgabe des Powerspektrums als 3D-Ober-
;                   flaeche" vorbereitend eingefuehrt.
;                 - In Pkt.1 Fehlermoeglichkeit durch /NONOISE (=+1dB!) abge-
;                   fangen; wird nunmehr als  noiseLevel = -60 dB gewertet.
;              [30/10/90-khd] => V 01.15:
;                 - Fuer Aufruf des UP TRIM_IMAGE (das geht erst ab V01.7) Key-
;                   words XSIZE und YSIZE zur Voreinstellung der Groesse des 
;                   Bildausschnitts eingefuehrt (Default: 256 pixel x 256 pixel).
;                 - Bei Ausgabe des Powerspektrums x_Titel + y_Titel eingefuehrt.
;                 - Beschriftung der Achsen auf internat. Norm => f/... (Groes- 
;                   sengleichung!) umgestellt.
;                 - In Pkt.8 Normierung fuer XSIZE bzw. YSIZE geaendert; dazu
;                   CONST Plot_HSize neu eingefuehrt.
;                 - In Pkt.8 beim 2. XYOUTS -0.67 => -0.73 geaendert.
;              [1/11/90-khd] => V 01.16:
;                 - In Pkt.3 Aufruf von TRIM_IMAGE erfolgt nunmehr mit XSIZE=...
;                   und YSIZE=...; /POTENZ entfaellt damit.
;                 - In Pkt.8 Normierung fuer XSIZE bzw. YSIZE nochmals verbes-
;                   sert wg. besserer Platzausnutzung auf dem Papier..
;                 - Keyword /AXIS zur Ausgabe von 2 Center-Lines eingefuehrt;
;                   unter Pkt.8 implementiert.
;              [6/11/90-khd] => V 01.17:                  (Copy => SEIDLER + DK)
;                 - Im Pkt.8 bei AXIS ... CHARSIZE=0. zugefuegt wg. kl. Zahlen.
;                 - In Pkt.5 wird nunmehr zusaetzlich das Phasenspektrum 
;                   berechnet.
;                 - Keyword /PHASESP eingefuehrt sowie in Pkt.7+8 dazu Ausgabe-
;                   moeglichkeit des Phasenspektrums zugefuegt.
;              [12/11/90-khd] => V 01.18:        (Copy => AG GRUESSER + SEIDLER)
;                 - Keyword /PHASESP => /PHASE umbenannt.
;                 - Keyword /CLEAN eingefuehrt; entsprechende Ergaenzungen im
;                   Pkt.1 + Pkt.6 vorgenommen.
;                 - Vorne Liste der verwendeten WDV-Programme zugefuegt.
;                 - Keyword /ADD eingefuehrt; noch nicht implementiert.
;              [19/11/90-khd] => V 01.19:                         (Copy => PRZY)
;                 - Aufruf geaendert durch Zufuegen der neuen Ausgangs-VAR 
;                   job; daher im ges. Programm VAR file => job umgestellt.
;                 - In Pkt.1 NOT KEYWORD_SET(addphase)... getilgt und ...
;                 - unter Pkt.11 /ADD implementiert.
;              [20/11/90-khd] => V 01.20:
;                 - Im Pkt.5 + Pkt.11 Boolean-VAR PhaseSp_Printed eingefuehrt
;                   wg. Abbruch der Schleife.
;              [27/11/90-khd] => V 01.21:             (Copy => AG GRUE + AG KOE)
;                 - Keyword /NOTRIM eingefuehrt zum Ueberspringen 
;                   der Ausfuehrung von Pkt.3.
;              [17/12/90-khd] => V 01.22:
;                 - Bedeutung von /HOLD geaendert; dazu unter Pkt.11 bei
;                   PRINT_PLOT,...  /NOPRINT zugefuegt.
;              [17/1/91-khd] => V 01.23:
;                 - Bedeutung von /EPS geaendert, so dass mit zusaetzlichem
;                   /PRIM auch das Bild/der Bildausschnitt in der EPS-File 
;                   ausgegeben werden kann; dazu im Pkt.10 in der Keyword-
;                   Abfrage  AND NOT KEYWORD_SET(eps)  getilgt.
;              [29/1/91-khd] => V 01.24:                       (Copy => SEIDLER)
;                 - Dazu (bei /EPS) im Pkt.8 die BoundingBox bzgl. YSIZE ver-
;                   groessert sowie ein Shift auf den oberen Plot eingefuehrt.
;              [25/2/91-khd] => V 01.25:                (Copy => SEIDLER + PRZY)
;                 - Neues Keyword /SAVE eingefuehrt zur Sicherung
;                   der VAR spectrum.
;                 - Dieses in Pkt.5 realisiert.
;              [11/3/91-khd] => V 01.26:                       (Copy => SEIDLER)
;                 - Neues Keyword WINDOW eingefuehrt.
;                 - Dieses im Pkt.1 + Pkt.4 implementiert.
;                 - Im Pkt.11 Abfrage der Ausgabeart (EPS oder PS) ergaenzt.
;                 - Im Pkt.13 Info ueber Spektrum-File zugefuegt.
;              [*/*/92-khd] => V 01.27:
;                 - Ausgabe eines Titels mit Versionsbezeichnung zugefuegt.
;                 - Check der Existenz der Input-File zugefuegt.
;===============================================================================
; Dieses Programm berechnet das logarithmische Powerspektrum und das Phasen-
; spektrum eines Bildes/Fotos in der ueblichen (normierten) Form, d.h. die 
; Frequenz f = 0 cycles/pixel liegt in der Mitte und an den Raendern gilt
; |f|=0,5 cycles/pixel (Nyquist-Frequenz).
;
; Die Spektralwerte bei negativen Frequenzen (diese resultieren aus den komple-
; xen Drehzeigern der FFT, die sich im mathematisch negativen Sinn drehen) sind
; zwar redundant, muessen aber fuer eine sich anschliessende Verarbeitung (z.B.
; eine Filterung im Frequenzbereich mit anschliessender FFT-Ruecktransformation)
; erhalten bleiben. Wird dieses nicht beachtet, was in der Literatur durchaus 
; zu bestaunen ist, werden total falsche Ergebnisse erzielt!
;
; Die Groessen an den Achsen des (log.) Powerspektrums sind:
;
;     o  X-Achse:  Frequenz   fx = (i/nx) * fabt_x    mit   i = 0..|nx/2|
;     o  Y-Achse:  Frequenz   fy = (j/ny) * fabt_y    mit   j = 0..|ny/2|
;     o  Z-Achse:  Amplitude  p  = 10*lg(|FFT(image)|^2)    (siehe Pkt.5)
;
; Die Groessen an den Achsen des Phasenspektrums sind:
;
;     o  X-Achse:  Frequenz   fx = (i/nx) * fabt_x    mit   i = 0..|nx/2|
;     o  Y-Achse:  Frequenz   fy = (j/ny) * fabt_y    mit   j = 0..|ny/2|
;     o  Z-Achse:  Phase     phi = arctan(Im{FFT(image)}/ Re{FFT(image)})
;
; Ueblicherweise wird die Abtastfrequenz fabt_x = fabt_y = fabt = 1 cycle/pixel
; sein, da das Bild mit Scannern in beiden Richtungen mit der gleichen Pixel-
; Aufloesung d (pixel/mm) abgetastet wird. Auch bei Bildern gilt fuer die 
; Nyquist-Frequenz:  fnyq = fabt/2,  woraus sich fuer  i = nx/2  bzw.  j = ny/2
; ergibt. Die Nyquist-Frequenz ist also die maximal im Powerspektrum bzw. im
; Phasenspektrum dargestellte 
; Frequenz (am Bildrand):
;                                   fnyq = (nx/2*nx) * fabt
;                                   fnyq = 1 cycle/2 pixel = 0,5 cycles/pixel.
;
; Wenn die Pixel-Anzahl in X- und Y-Richtung gleich ist, also  nx = ny = n ist,
; dann laesst sich die Frequenz auch in cycles/image angeben.
; Beispiel:   n = 512,                      f = (i/n) * fabt
;             d.h.                          f = (i/n) * 2 * fnyq
;             512 Pixel = 1 image           f = (i/n) * 2 * 0,5 cycles/pixel
;                                           f =  i * cycles/(n * pixel)
;                                           f =  i * cycles/image
;                                           f/(cycles/image) = i   mit  i=0..n/2
;             Bei i=n/2=256:                fnyq = 256 cycles/image
;===============================================================================
; 0. Prolog, Variablen und Konstanten:
;===============================================================================
;
  On_Error, 2                               ; Return => Ruf. Programm.

  Print, ""
  Print, "**************  IMAGE_POWER: Spektren eines Bildes  V 1.0.26a  **************"
  Print, ""

  rand  = 50                                ; Rand fuer (DEC-)Window in Pixel.
    w0  = 10                                ; Fuer Bezugskoo. des (DEC-)Windows.
  Mask  = "177776                           ; Maske fuer gerade Zahl.
  Norm  = 256                               ; Vereinbarte normale Bildgroesse
  ok    = ""                                ; in Pixel.
 answer = ""                                ; Typ-Deklaration.
  PSres = 2540.                             ; point/inch = ipt/inch werden von
                                            ; IDL in PostScript-Ausgaben verwen-
                                            ; det; PS-Standard = 72 point/inch !

  answer   = FindFile(fileName)             ; Ist die Datei vorhanden ?
  If (answer(0) eq "") $
     then  Message, "Die Datei existiert nicht."
  fileName = STRUPCASE(fileName)            ; Ermittlung des File-Namens (Job)
  cutLen   = STRPOS(fileName, ".TIF")       ; ohne Suffix; z.B. wird aus
  cutLen   = LONG(cutLen(0))                ; fileName = "MAC:FOTO.TIF"
  job      = STRMID(fileName,0,cutLen)      ; =>         "MAC:FOTO"
;;PRINT, "o  File:  ", job 
;
;===============================================================================
; 1. Behandlung der Keywords (Default-Werte):
;===============================================================================
;
  IF NOT KEYWORD_SET(filter)  $
     THEN   filter = "BH"     $
     ELSE   filter = STRUPCASE(filter)

  IF NOT KEYWORD_SET(xSize)    AND  NOT KEYWORD_SET(ySize) $
     THEN    BEGIN
             xSize = Norm
             ySize = Norm
  ENDIF

  IF KEYWORD_SET(Print_Image)  AND  NOT KEYWORD_SET(PS_Print) $
     THEN    PS_Print = 1                   ; Print_Image  n u r  bei PostScript.

  IF KEYWORD_SET(PS_Print) $
     THEN    PS_Asked = 1  $
     ELSE    PS_Asked = 0

  IF KEYWORD_SET(cleanSp)      AND  NOT KEYWORD_SET(noiseLevel) $
     THEN    noiseLevel = 1

  IF KEYWORD_SET(negativ)      AND  NOT KEYWORD_SET(noiseLevel) $
     THEN    noiseLevel = -60               ; /NEGATIV nur mit NONOISE sinnvoll.

  IF KEYWORD_SET(noiseLevel) $
     THEN    BEGIN
             IF (noiseLevel EQ 1) $
                 THEN  noiseLevel = -60     ; Sonst wird +1 dB genommen!
  ENDIF

  IF KEYWORD_SET(eps)          AND  NOT KEYWORD_SET(PS_Print) $
     THEN    PS_Print = 1

  IF KEYWORD_SET(hold)         AND  NOT KEYWORD_SET(PS_Print) $
     THEN    PS_Print = 1

  IF KEYWORD_SET(show3D) $
     THEN    negativ = 0                    ; Sonst wird's ein Tal!
;
;===============================================================================
; 2. Lesen der Bilddaten aus einer File im TIFF-Format:
;===============================================================================
;
  IF KEYWORD_SET(Get_Info) $
     THEN  image = TIFF_READ(fileName, $    ; => Programm P 476/13 der Abt.
                            /INFO,/BLACK) $ ;    Wiss.Datenverarbtg. des IfP.
     ELSE  image = TIFF_READ(fileName)     
;
;===============================================================================
; 3. Bildausschnitt interaktiv festlegen:
;===============================================================================
;                           
  IF KEYWORD_SET(notrim)  $
     THEN   GOTO, Pre_Filter                ; Bild ist bereits beschnitten.

  IF NOT KEYWORD_SET(batch) $
     THEN BEGIN
;;        TRIM_IMAGE, image, $              ; => Programm P 476/15 der Abt.
;;                    /POTENZ               ;    Wiss.Datenverarbtg. des IfP.
          TRIM_IMAGE, image, $              ; => Programm P 476/15 der Abt.
                      XSIZE=xSize, $        ;    Wiss.Datenverarbtg. des IfP.
                      YSIZE=ySize
  ENDIF

    dim = SIZE(image)
     nx = dim(1)                            ; Pixel-Anzahl in X-Richtung.
     ny = dim(2)                            ; Pixel-Anzahl in Y-Richtung.
     nx = FIX(nx AND Mask)                  ; Gradzahligmachen der
     ny = FIX(ny AND Mask)                  ; Spalten und Zeilen
  image = image(0:nx-1, 0:ny-1)             ; des Bildes.

  IF KEYWORD_SET(norebin) OR (nx LE Norm) OR (ny LE Norm) THEN  GOTO, Pre_Filter

  image = REBIN(image, nx/2, ny/2)          ; Bildgroesse halbieren.
    dim = SIZE(image)
     nx = dim(1)                            ; Pixel-Anzahl in X-Richtung.
     ny = dim(2)                            ; Pixel-Anzahl in Y-Richtung.
     nx = FIX(nx AND Mask)                  ; Gradzahligmachen der
     ny = FIX(ny AND Mask)                  ; Spalten und Zeilen
  image = image(0:nx-1, 0:ny-1)             ; des Bildes.
;
;===============================================================================
; 4. Vorfilterung des Bildes mit einem Tiefpass-Window: 
;===============================================================================
;
;    HINWEIS:  Um im Spektrum von vornherein Artefakte durch ein "Wrap-around" 
;              hoher Frequenzanteile zu vermeiden, muss das Bild/Foto  v o r 
;              der Berechnung des Powerspektrums geeignet tiefpassbegrenzt 
;              werden. Falls mit dem Keyword WINDOW="..." nichts anderes ver-
;              langt ist, wird hierzu das 4-gliedr. BLACKMAN-HARRIS Window 
;              - ein WDV-Programm - verwendet, dessen (stoerende) Seitenbaender
;              einen Abstand von -74 dB (!) haben.

Pre_Filter:
  CASE filter OF
   "BH" : lim_image = BH_WINDOW(image,LEVEL=74)  ; => Programm P 476/9 der Abt.
                                                 ; Wiss.Datenverarbtg. des IfP.
   "HAN": lim_image = HAN_WINDOW(image,LEVEL=32) ; => Programm P 476/21 der Abt.
                                                 ; Wiss.Datenverarbtg. des IfP.
   "NO" : lim_image = image
    ELSE: MESSAGE, "Window " + filter + " ist nicht definiert."
  ENDCASE
;
;===============================================================================
; 5. Berechnung des Powerspektrums und des Phasenspektrums:
;===============================================================================
;
  dim = SIZE(lim_image)
   nx = dim(1)                              ; Pixel-Anzahl in X-Richtung.
   ny = dim(2)                              ; Pixel-Anzahl in Y-Richtung.

  spectrum = FFT(lim_image, -1)             ; Diskrete Fourier-Transformation
                                            ; => f-Bereich.
  IF KEYWORD_SET(save) $
     THEN  BEGIN
           spectr_File = job + ".SPECTRUM"
           SAVE, spectrum, FILENAME=spectr_File   ; Spektrum archivieren.
  ENDIF

   Null_dB = ABS(spectrum(0,0))

     power = SHIFT(20.*ALOG10(ABS(spectrum)/Null_dB), nx/2, ny/2)  ; In Dezibel.

     phase = SHIFT(ATAN(IMAGINARY(spectrum),FLOAT(spectrum)), nx/2, ny/2)
                                            ; In Radian: [-Pi .. +Pi].

      time = SYSTIME(0)                     ; Zeitpunkt der Berechnung 
                                            ; festhalten (Creation-Date).
  PhaseSp_Printed = 0
;
;===============================================================================
; 6. Aufbereitung des Powerspektrums und des Phasenspektrums (optional!):
;===============================================================================
;
;    HINWEIS:   An dieser Stelle besteht die Moeglichkeit, das Programm 
;               um eine Bearbeitung bzw. Aufbereitung des Powerspektrums
;               nach eigenen Wuenschen zu ergaenzen; z.B. um eine Behand-
;               lung des Rauschhintergrunds, die z.Zt. im folgenden im-
;               plementiert ist.
;
  IF KEYWORD_SET(noiseLevel) $
     THEN BEGIN
          maxPower = MAX(power, MIN=minPower)                 ; Alle Spektral-
          power = BYTSCL(power, MIN=maxPower + noiseLevel, $  ; Amplituden, die
                                MAX=maxPower, TOP=255B)       ; -noiseLevel dB
                                            ; unter der max. Amplitude liegen,
                                            ; werden => 0 => Schwarz gesetzt.

          IF KEYWORD_SET(cleanSp) $         ; Spektren von obsoleten Spektral-
             THEN BEGIN                     ; Anteilen befreien:
                  is_obsolete = WHERE( SHIFT(power, -nx/2, -ny/2) EQ 0B)
                  spectrum(is_obsolete) = COMPLEX(0.,0.)
                  is_obsolete = WHERE( power EQ 0B)
                  phase(is_obsolete)    = !PI                 ; d.h. => WE.
          ENDIF         

          IF KEYWORD_SET(negativ) $         
             THEN  power = 255B - power     ; => Negativ-Darstellung (0 => Weiss).
  ENDIF
;
;===============================================================================
; 7. Setzen von allgemeinen Parametern fuer die folgenden Ausgaben:
;===============================================================================
;
Power_Output:
  wx = nx + 2 * rand                        ; Groesse des (DEC-)Windows
  wy = ny + 2 * rand                        ; = Fenster auf dem Schirm.

  x0 = rand   &   x1 = x0 + nx - 1          ; In diesem Fenster: Eck-Koo.
  y0 = rand   &   y1 = y0 + ny - 1          ; des Plot-Windows in Pixel.

  !X.CHARSIZE =  0.9     &     !Y.CHARSIZE =  0.9
  !P.TICKLEN  = -0.03   
  !X.TICKS    =  2       &     !Y.TICKS    =  2
  !X.MINOR    =  5       &     !Y.MINOR    =  5
  Plot_Oben   =  6.0     &     Plot_Unten  =  0.5        ; Nur fuer PostScript;
  Plot_Rand   =  1.5                                     ; alle Zahlenangaben
  Plot_VSize  =  5.0     &     Plot_HSize  =  6.5        ; in Inches.
  Plot_Shift  =  Plot_Oben - Plot_Unten
  Shift       =  FIX(Plot_Shift * PSres)                 ; Fuer direkten PS-Code
  Shift       =  STRCOMPRESS(Shift, /REMOVE)             ; zum Shiften des Plot-
                                                         ; Windows.
;_______________________________________________________________________________

  IF NOT KEYWORD_SET(phaseSp) $
     THEN  BEGIN
           Titel        = "Log. Powerspektrum von " + fileName 
           Window_Titel = "Log. Powerspektrum"
           ENDIF $
     ELSE  BEGIN
           Titel        = "Phasenspektrum von " + fileName
           Window_Titel = "Phasenspektrum"
           ENDELSE
;_______________________________________________________________________________

  IF NOT KEYWORD_SET(noiseLevel)  $
     THEN  noiseLevel = -120                ; Nur fuer die Ausgabe!
;_______________________________________________________________________________

  IF KEYWORD_SET(PS_Print) $                ; Je nach Ausgabe-Art unterschied- 
     THEN  BEGIN                            ; lich zu setzende Parameter:
           x_Titel  = "f!Ix!N  / [cycles/pixel]"
           y_Titel  = "f!Iy!N  / [cycles/pixel]"
           IF NOT KEYWORD_SET(phaseSp) $
              THEN  SubTitel = "n!Ix!N = "   + STRCOMPRESS(nx,/REMOVE) + " pixel;" + $
                               "  n!Iy!N = " + STRCOMPRESS(ny,/REMOVE) + " pixel;" + $
                               "  noise < "  + STRCOMPRESS(noiseLevel,/REMOVE) + " dB" $
              ELSE  SubTitel = "n!Ix!N = "   + STRCOMPRESS(nx,/REMOVE) + " pixel;" + $
                               "  n!Iy!N = " + STRCOMPRESS(ny,/REMOVE) + " pixel;" + $
                               "  black: phi = - Pi;  white: phi = + Pi"
           ENDIF $
     ELSE  BEGIN                            ; Keine PostScript-Ausgabe:
           x_Titel  = "fx  / [cycles/pixel]"
           y_Titel  = "fy  / [cycles/pixel]"
           IF NOT KEYWORD_SET(phaseSp) $
              THEN  SubTitel = "nx = "       + STRCOMPRESS(nx,/REMOVE) + " pixel;" + $
                               "  ny = "     + STRCOMPRESS(ny,/REMOVE) + " pixel;" + $
                               "  noise < "  + STRCOMPRESS(noiseLevel,/REMOVE) + " dB" $
              ELSE  SubTitel = "nx = "       + STRCOMPRESS(nx,/REMOVE) + " pixel;" + $
                               "  ny = "     + STRCOMPRESS(ny,/REMOVE) + " pixel;" + $
                               "  black: phi = - Pi;  white: phi = + Pi"
           ENDELSE
;_______________________________________________________________________________

  simx    = FLOAT(nx)                       ; Groesse des Bildes in
  simy    = FLOAT(ny)                       ; DEVICE-Koordinaten.
  x_Achse = FLOAT([0, nx]-nx/2)/simx
  y_Achse = FLOAT([0, ny]-ny/2)/simy
;
;===============================================================================
; 8. Ausgabe des Power- oder Phasenspektrums als Pixel-Bild mit Achsen:
;===============================================================================
;
;    HINWEIS:  Die Ausgabe erfolgt in sogenannter Overlay-Technik!
;
  IF KEYWORD_SET(batch)  OR  (PS_Asked EQ 1)  THEN   GOTO, Get_PlotWindow

  WINDOW, XPOS=w0,  YPOS=w0, $              ; Neues Fenster erzeugen.
          XSIZE=wx, YSIZE=wy, $
          RETAIN=2, /FREE, $
          TITLE=Window_Titel
  LOADCT, 0, /SILENT                        ; Laden der Color-Table: Lin. Grau.
;_______________________________________________________________________________

Get_PlotWindow:
  IF KEYWORD_SET(PS_Print) $                ; Bei PostScript-Ausgabe:
     THEN BEGIN
          OldDevice = !D.NAME
          SET_PLOT, "PS", /COPY
          IF (y1 GE x1) $
             THEN  d = y1/Plot_VSize $      ; Bezugswert: Aufloesung in 
             ELSE  BEGIN                    ; X/Y-Richtung in pixel/inch.
                   d = x1/Plot_HSize        
                   IF (y1/d GT Plot_VSize) $
                      THEN  d = y1/Plot_VSize
                   ENDELSE
          IF KEYWORD_SET(eps) $
             THEN  DEVICE, /ENCAPSULATED, FILENAME=job + ".EPS", $
                            XSIZE  =x1/d, YSIZE  = Plot_Shift + y1/d, /INCHES, $
                           /TIMES, /ITALIC, FONT_SIZE=12 $
             ELSE  DEVICE,  XOFFSET=Plot_Rand, YOFFSET=Plot_Oben, $
                            XSIZE  =x1/d,      YSIZE  =y1/d, /INCHES, $
                           /TIMES, /ITALIC, FONT_SIZE=12
  ENDIF   ;KW_S(PS_Print)
;_______________________________________________________________________________

  PLOT, [x0,x1], [y0,y1], /NODATA, $        ; Hiermit nur Definition der Ein-
         XSTYLE=4, YSTYLE=4                 ; stellungen des Plot-Windows;
                                            ; wg. der Device-Unabhaengigkeit.

  pwx   = !X.WINDOW                         ; Jeweils 2-elem. Vektor
  pwy   = !Y.WINDOW                         ; in NORM-Koordinaten.
  px    = pwx * !D.X_VSIZE                  ; => DEVICE-Koordinaten.
  py    = pwy * !D.Y_VSIZE                  ; => DEVICE-Koordinaten.
  spwx  = px(1) - px(0)                     ; Groesse des Plot-Windows
  spwy  = py(1) - py(0)                     ; in DEVICE-Koordinaten.
  aspim = simx/simy                         ; Seitenverhaeltnis des Bildes.
  asppw = spwx/spwy                         ; Seitenverhaeltnis des Plot-Windows.
  a     = aspim/asppw                       ; Verhaeltnis der Seitenverhaeltnisse.
;_______________________________________________________________________________

  IF KEYWORD_SET(PS_Print) $                ; Ausgabe des Pixel-Bildes
     THEN BEGIN                             ; von power (Powerspektrum) bzw.
          IF (a GE 1.0) $                   ; von phase (Phasenspektrum).
             THEN  spwy = spwy/a $          ; Groesse des Plot-Windows an
             ELSE  spwx = spwx*a            ; das Pixel-Bild anpassen.
          IF KEYWORD_SET(eps) $                          ; Nur bei EPS-File:
             THEN  DEVICE, OUTPUT='0 ' + Shift + $       ; Plot-Window nach
                                  ' translate'           ; oben; keine "!!
          IF NOT KEYWORD_SET(phaseSp) $
             THEN  TVSCL, power, px(0), py(0), /DEVICE, $
                                 XSIZE= spwx, $
                                 YSIZE= spwy  $
             ELSE  TVSCL, phase, px(0), py(0), /DEVICE, $
                                 XSIZE= spwx, $
                                 YSIZE= spwy 
          ENDIF $
     ELSE BEGIN                             ; Keine PostScript-Ausgabe.
          IF NOT KEYWORD_SET(phaseSp) $
             THEN  TVSCL, power, px(0), py(0), /DEVICE $
             ELSE  TVSCL, phase, px(0), py(0), /DEVICE
          spwx = simx                       ; Groesse des Plot-Windows auf die
          spwy = simy                       ; Groesse des Pixel-Bildes setzen.
          ENDELSE
;_______________________________________________________________________________

  x0 = px(0)     &     x1 = px(0) + spwx    ; Eckpunkte des aktuellen Plot-
  y0 = py(0)     &     y1 = py(0) + spwy    ; Windows in DEVICE-Koordinaten.

  Plot_Window = [x0, y0, x1, y1]            ; Definition des aktuellen Plot-
                                            ; Windows in DEVICE-Koordinaten.

  PLOT, [x0,x1], [y0,y1], /NODATA, /DEVICE,$; Ausgabe der Achsen-Box in DEVICE-
         POSITION=Plot_Window, FONT=0, $    ; Koordinaten mit Beschriftung.
         XRANGE=x_Achse, YRANGE=y_Achse, $  
         XSTYLE=1, YSTYLE=1, /NOERASE, $
         XTITLE= x_Titel, $
         YTITLE= y_Titel 
;_______________________________________________________________________________

  IF KEYWORD_SET(PS_Print) $                ; Titel ueber Plot-Mitte zentriert:
     THEN   DEVICE, /TIMES, /BOLD, /ITALIC, FONT_SIZE=12
  XYOUTS, 0.0, 0.62,  /DATA, Titel, ALIGNMENT=0.5, FONT=0

  IF KEYWORD_SET(PS_Print) $                ; Untertitel und Creation-Date 
     THEN   BEGIN                           ; (nur bei PostScript-Ausgabe):
            DEVICE, /TIMES, FONT_SIZE=9
            XYOUTS, 0.0,  0.56, /DATA, "Creation-Date:  " + time, $
                                                 ALIGNMENT=0.5, FONT=0
            XYOUTS, 0.0, -0.73, /DATA, SubTitel, ALIGNMENT=0.5, FONT=0
  ENDIF
;_______________________________________________________________________________

  IF KEYWORD_SET(axis) $                    ; Achsenkreuz durch den Nullpunkt:
     THEN  BEGIN
           IF KEYWORD_SET(negativ) $
              THEN  farbe = 0      $                     ; Linie => SW.
              ELSE  farbe = 255B                         ; Linie => WE.
           AXIS, 0, 0, XAXIS=0,  /DATA,  COLOR=farbe, $  ; x-Linie d. Nullpkt.
                       THICK=1.0, TICKLEN=0., CHARSIZE=0.
           AXIS, 0, 0, YAXIS=0,  /DATA,  COLOR=farbe, $  ; y-Linie d. Nullpkt.
                       THICK=1.0, TICKLEN=0., CHARSIZE=0.
  ENDIF
;_______________________________________________________________________________
  GOTO, Image_Output
;
;===============================================================================
; 9. Ausgabe des Powerspektrums als 3D-Oberflaeche (optional!):
;===============================================================================
;
Power_3D:
  MESSAGE, "Dreidimensionale Ausgabe ist noch nicht implementiert."
;
;===============================================================================
; 10. Ausgabe des Bildausschnitts (optional!):
;===============================================================================
;
Image_Output:
  IF KEYWORD_SET(Print_Image) $
     THEN BEGIN
          DEVICE, OUTPUT='0 -' + Shift + ' translate', $ ; Plot-Window nach
                  /TIMES, /ITALIC, FONT_SIZE=12          ; unten; keine "!!

          TVSCL, image, px(0), py(0), /DEVICE, $         ; Ausgabe des Pixel-
                 XSIZE= spwx,  YSIZE= spwy               ; Bildausschnitts.

          PLOT, [x0,x1], [y0,y1], /NODATA, /DEVICE, $    ; Ausgabe eines 
                 POSITION=Plot_Window, FONT=0, $         ; Bildrandes mit
                 XRANGE=[0,nx], YRANGE=[0,ny], $         ; Beschriftung 
                 XSTYLE=1, YSTYLE=1, /NOERASE, $         ; durch den
                 XTITLE="i!Ix!N  / pixel", $             ; Pixel-Index.
                 YTITLE="i!Iy!N  / pixel"
  ENDIF
;_______________________________________________________________________________

  fileName_I = job + "_PART.TIF"            ; Ausgabe des Bildausschnitts als
                                            ; Bilddaten-File im TIFF-Format:

  TIFF_WRITE, image, fileName_I, $          ; => Programm P 476/14 der Abt.
                     ORIENTATION=1          ;    Wiss.Datenverarbtg. des IfP.
;
;===============================================================================
; 11. Ausgabe der PostScript-File des Powerspektrums bzw. des Phasenspektrums:
;===============================================================================
;
Do_Output:
  IF NOT PS_Asked $
     THEN BEGIN
          PS_Asked = 1
          PRINT, FORMAT='($, "o  Soll das Spektrum ausgegeben werden ? [N]")'
          READ, ok    &    ok = STRUPCASE(ok)
          IF (ok EQ "Y") OR (ok EQ "YES") OR (ok EQ "J") OR (ok EQ "JA") $
              THEN BEGIN
                   PS_Print = 1
                   PRINT, FORMAT='($,"o  In welcher Art (EPS oder PS) ? [EPS]")'
                   READ, answer    &       answer = STRUPCASE(answer)
                   IF (answer EQ "") THEN  answer = "EPS"
                   CASE answer OF
                     "EPS": BEGIN
                            eps = 1
                            GOTO, Power_Output       ; Am Anfang von Pkt.7.
                            END
                      "PS": GOTO, Power_Output       ; Am Anfang von Pkt.7.
                      ELSE: BEGIN
                            MESSAGE, "Falsche Ausgabe-Art:  " + $
                                     answer, /CONTINUE
                            PS_Asked = 0
                            GOTO, Do_Output
                            END
                   ENDCASE
          ENDIF
  ENDIF

  IF NOT PhaseSp_Printed AND $
     KEYWORD_SET(addphase) $                ; => Auf 2. Blatt 
     THEN BEGIN                             ;    das Phasenspektrum ausgeben.
          PhaseSp_Printed = 1     
          phaseSp = 1
          GOTO, Power_Output                ; Am Anfang von Pkt.7.
     ENDIF

  IF KEYWORD_SET(PS_Print) $
     THEN BEGIN
          IF KEYWORD_SET (eps) $
             THEN  DEVICE, /CLOSE_FILE $    ; EPS-File schliessen.
             ELSE  IF KEYWORD_SET(hold) $
                      THEN  PRINT_PLOT, job + ".PS", /NOPRINT  $
                      ELSE  PRINT_PLOT      ; PS-Ausgabe des Spektrums
          SET_PLOT, OldDevice               ; und ggfs. des Bildausschnittes.
  ENDIF   ;KW_S(PS_Print)
;
;===============================================================================
; 12. Ausgabe des Powerspektrums als Bilddaten-File im TIFF-Format(ohne Achsen):
;===============================================================================
;
  fileName_P = job + "_POWER.TIF"           ; Bildung des File-Namens fuer das
                                            ; Powerspektrum: Der File-Name muss
                                            ; den Suffix ".TIF" erhalten, da
                                            ; manche (MAC-)Applikationsprogramme
;;PRINT, fileName_P                         ; dieses erwarten.
  
  TIFF_WRITE, power, fileName_P, $          ; => Programm P 476/14 der Abt.
                     ORIENTATION=1          ;    Wiss.Datenverarbtg. des IfP.
;
;===============================================================================
; 13. Schluss:
;===============================================================================
;
Ende:
  IF NOT KEYWORD_SET(batch) $
     THEN BEGIN
          PRINT, ""
          PRINT, "o  Die Daten des komplexen 2D Spektrums stehen fuer"
          PRINT, "   eine (interaktive) Weiterverarbeitung bereit."
          IF KEYWORD_SET(save) $
             THEN  PRINT, "   Ausserdem in File:  " + spectr_File
  ENDIF
  RETURN
  END
;===============================================================================