Simulation von Beugungsbildern mit ImageJ

Hallo zusammen,

nachdem nun mein anderer Thread (Simulation Einfluss Fangspiegel/Spinne/etc auf Beugungsbild) langsam heiss läuft und das Interesse an verschiedenen Obstruktionen und Geometrien gross und grösser wird ist es Zeit, das "wie" aus diesem Thread auszulagern.

Angefangen hatte es ja mit meiner Frage, welchen Einfluss diverse Geometrien auf das Beugungsbild haben und wie man das simulieren und visualisieren kann. Erste Versuche waren kläglich gescheitert, ein Ausflug zu Maskulator (http://www.njnoordhoek.com/?p=376) schien erfolgversprechend. Nach ausgiebiger Lektüre des Threads https://www.cloudynights.com/t…what-to-do-what-to-avoid/ in Cloudynights und einigen Papers habe ich weiter mit Fiji (ImageJ) gespielt. Maskulator hat ein paar Probleme, die die Genauigkeit der Resultate begrenzt. So ist die FFT auf 2kx2k begrenzt, die Dateien haben grundsätzlich nur 8-bit Tiefe. Vergleiche mit verschiedenen Geometrien sind so nur bedingt möglich. Leider ist Niels, der Autor von Maskulator vor ein paar Jahren verstorben, so dass die Software meines Wissens nicht mehr weiterentwickelt wird.

Als dann der Astrotreff kurzzeitig gegroundet war habe ich die Zeit genutzt, meine Versuche in Makros zu giessen, so dass die Auswertung automatisch abgearbeitet werden kann. Die Makros sind weder professionell noch optimiert, tun aber, was ich von ihnen erwarte.

So wie es aussieht, scheint das Thema auch andere zu interessieren, deshalb möchte ich euch das ganze hier so vorstellen, dass es (hoffentlich) auch ohne meine Hilfe möglich ist, selber solche Auswertungen zu machen.

Wie gesagt arbeite ich (nicht nur für diese Aufgabe) mit Fiji (ImageJ) (https://imagej.net/software/fiji/). Fiji steht für Fiji ist just ImageJ und ist eine sofort lauffähige Version von ImageJ. Dieses wiederum ist ein in Java geschriebenes Framework zur Bildverarbeitung und -analyse im wissenschaftlichen Bereich. Wie in diesem Bereich üblich wird davon ausgegangen, dass der Benutzer mit einer mehr oder weniger ungeordneten Werkzeugkiste umgehen kann. Wer magische Knöpfe wie in Photoshop erwartet, hat eine steinige Lernkurve vor sich, die sich aber lohnt.

Ich spare mir eine Einführung in ImageJ, erst mal einfach etwas spielen und das System dahinter verstehen. Für uns wichtig ist hier: die Software kann vollständig skriptbasiert verwendet werden. Theoretisch können auch sehr viele Dateien als Stack verarbeitet werden, was aber hier nicht benötigt wird.

Also: ImageJ (Fiji) installieren, dann seit ihr (fast) bereit. Einzig der Plugin Polar Transformer (https://imagej.nih.gov/ij/plugins/polar-transformer.html) wird noch benötigt und muss installiert werden.

Das Vorgehen ist folgendermassen:

1. Als Ausgangspunkt brauchen wir eine Datei mit der gewünschten Maske. Um Artefakte bei der Fouriertransformation zu minimieren, muss die eigentliche Maske auf idealerweise <10% der Bildgrösse begrenzt werden. Das heisst eine Maske von 1500x1500 px (Pixeln) braucht einen schwarzen Trauerrand, der das Bild auf >15kx15k aufbläst. Ich habe mich für 2^14x2^14 (16384x16384)px entschieden. Da es nur transparent/opak gibt, ist die Dateigrösse bei PNG Format trotzdem nur ca. 45kB gross.

Die Masken habe ich in Freecad erstellt, als SVG exportiert, mit Inkscape geladen und als 1bit PNG mit 16kx16k Pixeln exportiert. Vorteil: in Freecad kann die gewünschte Geometrie exakt konstruiert werden, Inkscape arbeitet nativ mit SVG und kann Bitmaps exportieren.

Leider können solche Dateien nicht hier hochgeladen werden. Ich werde sie auf meinen SFTP Server laden und dann den Zugang hier posten.

2. ImageJ starten:

3. Maske(n) (PNG Datei(en)) damit öffnen. Das kann entweder über den Reiter "File" geschehen oder indem man die Datei(en) auf das oben gezeigte Fenster zieht.

4. Makro "calc_mask.ijm" und "calc_diff.ijm" öffnen.

5. Maskendatei selektieren

6. "calc_mask.ijm" selektieren und mit RUN starten:

Dieser Skript extrahiert den zentralen Teil der Maske für Dokumentationszwecke.

7. "calc_diff.ijm" starten. Nun rattert es eine Weile.

- Der Skript macht erst mal eine FFT (Fouriertransformation) der Maske, was auf meinem Rechner ca, 25s dauert. Je nach Rechnerleistung kann das auch länger gehen.

- Der zentrale Ausschnitt der so gerechneten PSF (Point Spread Function) wird als 32bit TIFF abgespeichert.

- Anschliessend wird logarithmiert und auf -mag (Minus Magnitude) skaliert. Damit das ganze visuell besser erkennbar ist, wird eine farbige LUT (Look Up Table) verwendet.

- dieses Bild wird als 8bit PNG abgespeichert.

8. Jetzt gehts rund. Aus diesem Bild wird ein 3D Plot erstellt und ebenfalls als PNG abgespeichert. Im Prinzip kann man in diesem Plot live rumfahren, wer das will, führt den Skript nur bis dahin aus (Im Skript auswählen, was man ausführen will und mit (CTRL-Shift-R) ausführen.

9. Anschliessend wird eine Linie vom Zentrum zum rechten Rand markiert. Entlang dieser wird ein Plot erzeugt, der die Intensität vom Zentrum her abbildet.

Dieser Plot wird optisch etwas aufgemöbelt und als PNG abgespeichert. Die hinterlegten Daten werden zusätzlich als CSV Datei abgeschpeichert.

Wer genug hat, kann hier aufhören. Aber: ausser bei rotationssymmetrischen Beugungsbildern schwankt die Intensität über den Winkel. Und deshalb machen wir weiter.

10. Nun wird das Beugungsbild abgewickelt. Mit einer Polartransformation werden vom Zentrum her Kreise über 360° als Geraden dargestellt. Das gibt eine Art Landkarte.

Da dieses Plugin nur mit 8-bit oder 16-bit Daten arbeitet, müssen einige Informationen zwischengespeichert und die Bilder entsprechend skaliert werden.

11. Nun kommt eines der mächigsten Tools von ImageJ zur Anwendung der ROI-Manager. Dieser erlaubt es, entweder mit geeigneten Operation automatisch oder, wie in diesem Fall, manuell ROI (Regions of Interest) zu definieren und anschliessend darüber weitere Operationen auszuführen. Wir sind interessiert, den polaren Verlauf der Beugungsringe darzustellen. Der Skript erstellt dazu Linien im Abstand von 17px, was dem Abstand der Maxima entspricht.

12. Nun rattert der ROI Manager und erstellt 11 Plots der ersten 11 Beugungsringe. Der Plot wird ebenfalls als PNG Datei, die Daten dahinter als CSV abgespeichert

13. Aufräumen. Alle Dateien werden geschlossen.

Das ganze ist ein wildes Gezuckel und am Ende liegen ein paar neue Dateien auf der Disk, hier für eine obstruktionsfreie Apertur:

<Maskendateiname>_diff.tif (PSF, linear)

Zu gross zum Hochladen (ca. 16MB)...

<Maskendateiname>_diff.png (PSF, logarithmisch in -mag)

<Maskendateiname>_sp.png (PSF Plot in 3D)

<Maskendateiname>_plot.tif (PSF Plot entlang einer Linie vom Zentrum nach aussen)

<Maskendateiname>_values.csv (Plotdaten dazu)

<Maskendateiname>_pol.png (PSF als Polardarstellung)

<Maskendateiname>_angular_psf.png (Plot dazu)

<Maskendateiname>_apsf.csv (Plotdaten dazu)

Wie schon geschrieben, Skripte können auch partiell ausgeführt werden, allerdings muss bedacht werden, dass Variablen nur während der Laufzeit bekannt sind, Schrittweises Ausführen (z.B. zum Debuggen) kann da etwas mühsam sein.

Die Forensoftware teilt mir gerade mit, dass ich zu geschwätzig bin und die Nachrichtenlänge auf 10k Zeilen begrenzt ist. Die Skipte folgen deshalb als eigener Beitrag.

Wenn Fragen sind (was ich annehme, das wird nicht auf Anhieb funktionieren), einfach fragen. Ich werde den Beitrag entsprechend anpassen und ergänzen.

Have fun

Herzliche Grüsse Robert

Hier die Skripte.

So, hier noch die Skripte in voller Hässlichkeit:

// extract mask

// get relevant data for output file name 
title = getTitle();
path = getInfo("image.directory");
file_name_no_extension = File.getNameWithoutExtension(title);

// calculate coordinates for central selection rectangle
image_width = getWidth();
image_center = image_width / 2;
mask_width = image_width / 8;

// select central part of image and copy to clipboard
//setTool("rectangle");
makeRectangle(image_center - mask_width / 2, image_center - mask_width / 2, mask_width, mask_width);
setBackgroundColor(255, 255, 255);
run("Copy");

// paste from Internal Clipboard
run("Internal Clipboard");

// save file
saveAs("PNG", path + file_name_no_extension+"_mask.png");

// close both files
close();

// calculate diffraction pattern of file
// output is a series of files

// get relevant data for output file name 
title = getTitle();
path = getInfo("image.directory");
file_name_no_extension = File.getNameWithoutExtension(title);

// different maximal intensity can be selected. Scale is in -mag (astronomical magnitude)
max_intensity = 0;
//max_intensity = 2;
//max_intensity = 5;

// create output directory from intensity
output_path = path + max_intensity + File.separator;
File.makeDirectory(output_path);

// calculate coordinates for central selection rectangle
image_width = getWidth();
image_center = image_width / 2;
PSF_width = image_width / 8;

// calculate PSF from input file.
run("FFT Options...", "raw do");

selectWindow("PS of "+title);

// select central part of image and copy to clipboard
//setTool("rectangle");
makeRectangle(image_center - PSF_width / 2, image_center - PSF_width / 2, PSF_width, PSF_width);
setBackgroundColor(255, 255, 255);
run("Copy");

// paste from Internal Clipboard
run("Internal Clipboard");
// get minimum and maximum values
resetMinAndMax();
getMinAndMax(min, max);
// scale image so max value is 1.0
run("Divide...", "value="+max);
// save image as 32bit TIFF
saveAs("TIFF", output_path+file_name_no_extension+"_diff.tif");

// calculate log (log_base_e) of image
run("Log");
// scale for -mag (astronomical magnitude)
scale_factor = 5 / log(100);
run("Multiply...", "value="+scale_factor);
// again scale intensity so max value is 0mag
resetMinAndMax();
getMinAndMax(min, max);
run("Subtract...", "value="+max);
// only show intensity values from 15mag to 0mag
setMinAndMax(-15 - max_intensity, 0 - max_intensity);
// use "Thermal" LUT for better visualisation
//run("Apply LUT");
run("Thermal");
// save as 8 bit PNG file
saveAs("PNG", output_path+file_name_no_extension+"_diff.png");

// create 3D Surface Plot from file and save as PNG file
run("3D Surface Plot", "plotType=3 smooth=0 colorType=0 scale=0.9 snapshot=1 rotationZ=45");
saveAs("PNG", output_path+file_name_no_extension+"_sp.png");

// select line from center to right edge
selectWindow(file_name_no_extension+"_diff.png");
makeLine(PSF_width / 2, PSF_width / 2, PSF_width - 1, PSF_width / 2);
// create intensity plot across line
run("Plot Profile");
Plot.setLogScaleY(false);
// only show values from 15mag to 0mag and only 11th maximum (ring n to ring n+1 = 17px)
Plot.setLimits(0.0,191.0,-15.0,0.0);
Plot.setFontSize(12.0);
Plot.setStyle(0, "blue,none,2.0,Line");
Plot.setXYLabels("Distance (pixels)", "Relative mag");
Plot.setFormatFlags("11001100111111");
// save plot as PNG
saveAs("PNG", output_path+file_name_no_extension+"_plot.png");

// export plot data as CSV
run("Clear Results");
Plot.getValues(xpoints, ypoints); 
for (i=0; i<xpoints.length; i++)
  setResult("Value", i, ypoints[i]);
updateResults();
saveAs("Results", output_path+file_name_no_extension+"_values.csv");
run("Close");
selectWindow(file_name_no_extension+"_plot.png"); 
close();

// now create polar visualisation (needs Polar Transformer Plugin)
// select central part of diffraction pattern
// calculate coordinates for central selection rectangle
selectWindow(file_name_no_extension+"_diff.png");
image_width = getWidth();
image_center = image_width / 2;
PSF_width = image_width / 4;

// select central part of image and copy to clipboard
//setTool("rectangle");
makeRectangle(image_center - PSF_width / 2, image_center - PSF_width / 2, PSF_width, PSF_width);
setBackgroundColor(255, 255, 255);
run("Copy");

run("Internal Clipboard");
run("Select All");

// calculate min value for later rescaling intensity
getMinAndMax(min, max);
minmin=-255/min;
setOption("ScaleConversions", true);
// polar transformation only works in 8-bit or 16-bit mode
run("8-bit");
run("Polar Transformer", "method=Polar degrees=360 default_center for_polar_transforms,");
run("Thermal");
run("Rotate 90 Degrees Right");
// change back to 32-bit mode and rescale intensity to -mag
run("32-bit");
run("Subtract...", "value="+255);
run("Divide...", "value="+minmin);
resetMinAndMax();

// only use data within circle of 180px
makeRectangle(0, 0, 359, 179);
setBackgroundColor(255, 255, 255);
run("Copy");
run("Internal Clipboard");
// save as PNG
saveAs("PNG", output_path+file_name_no_extension+"_pol.png");

// create ROI (regions of interest) table. We are interested in the first 11 diffraction rings
run("ROI Manager...");
for(i=1; i < 12; i++) {
    makeLine(0, 10*i + 7, 359, 10*i + 7);
    roiManager("Add"); }
// selectWindow("Profiles");
roiManager("Multi Plot");
Plot.setStyle(0, "gray,none,2.0,Line");
Plot.setStyle(1, "red,none,2.0,Line");
Plot.setStyle(2, "green,none,2.0,Line");
Plot.setStyle(3, "blue,none,2.0,Line");
Plot.setStyle(4, "magenta,none,2.0,Line");
Plot.setStyle(5, "cyan,none,2.0,Line");
Plot.setStyle(6, "yellow,none,2.0,Line");
Plot.setStyle(7, "orange,none,2.0,Line");
Plot.setStyle(8, "Red,none,2.0,Line");
Plot.setStyle(9, "Green,none,2.0,Line");
Plot.setStyle(10, "Blue,none,2.0,Line");
Plot.setLogScaleY(false);
Plot.setLimits(0.0,360.0,-15.0,0.0);
Plot.setFontSize(12.0);
Plot.setXYLabels("Distance (pixels)", "Value");
Plot.setFormatFlags("11001100111111");

// save polar visualisation as PNG
saveAs("PNG", output_path+file_name_no_extension+"_angular_psf.png");

// create plot from ROI data and save as CSV
Plot.showValues();
saveAs("Results", output_path+file_name_no_extension+"_apsf.csv");

// clean up and close as open files
run("Close");
close("ROI Manager");
//close("*");
close();
close();
close();
close();
close();
close();
close();
close();