Startrails automatisch stacken mit ImageMagick oder Python

inklusive dark, bias and flat frames

Im Sommerurlaub habe ich zum ersten Mal versucht eine Startrail-Aufnahme zu machen, d.h. nicht eine, sondern 26 á 60s mit einem 7.5mm Weitwinkel-Objektiv, um diese anschließen zu Sternenspuren übereinander zu legen (stacken). Da meine Kamera einen Fleck auf dem Sensor, ein hohes Rauschen und das Objektiv eine ausgeprägte Vignette hat, habe ich noch ein paar dark, bias und flat Aufnahmen analog zum Stacken von Astrofotos mit Siril gemacht. Diese habe ich dann in einem ersten Versuch mühsamer Handarbeit in Gimp miteinader kombiniert und mit dem Modus “Aufhellen” übereinander gelegt. Da dachte ich mir, dass muss aber einfacher gehen, wenn ich das nochmal machen will!

1. Versuch mit ImageMagick

Also habe ich ein Skript geschrieben, dass die RAW-Aufnahmen erstmal über das Kommandozeilentool von darktable in 16-bit tif umwandelt und anschließend mit ImageMagick die gemittelten Dark, Bias und Flat-Frames berechnet und diese dann über die Formel (light - dark) / (flat - bias) mit den einzelnen Startrail-Aufnahmen (light) verechnet. Anschließen werden die so korrigierten Bilder durch “Aufhellen” übereinander gelegt.

#!/bin/bash

set -euo pipefail

## function to flatten bias darks and flats
function flatten_bdf_raws() {
  dir=$1

  cd $dir
  ls *.RW2 | xargs -n 1 sh -c 'darktable-cli $0 ../general.xmp ${0%.*}.tif --bpp 16'
  
  imgs=( $(ls *.tif) )
  transpc=$(echo 1/${#imgs[*]} | bc -l)

  first_img=${imgs[0]}
  unset imgs[0]

  rotation=$(exiftool -Rotation $first_img)

  if [[ "$rotation" =~ "90 CW" ]]; then
    convert -quiet $first_img -rotate 270 -alpha set -background none -channel A -evaluate multiply 1 +channel ${dir}.tif
  elif [[ "$rotation" =~ "270 CW" ]]; then
    convert -quiet $first_img -rotate 90 -alpha set -background none -channel A -evaluate multiply 1 +channel ${dir}.tif
  else
    convert -quiet $first_img -alpha set -background none -channel A -evaluate multiply 1 +channel ${dir}.tif
  fi

  for img in ${imgs[*]}; do
    fbase=$(basename $img .tif)
    rotation=$(exiftool -Rotation $img)
    
    if [[ "$rotation" =~ "90 CW" ]]; then
      convert -quiet $img -rotate 270 -alpha set -background none -channel A -evaluate multiply $transpc +channel ${fbase}_alpha.tif
    elif [[ "$rotation" =~ "270 CW" ]]; then
      convert -quiet $img -rotate 90 -alpha set -background none -channel A -evaluate multiply $transpc +channel ${fbase}_alpha.tif
    else
      convert -quiet $img -alpha set -background none -channel A -evaluate multiply $transpc +channel ${fbase}_alpha.tif
    fi
    
    composite -depth 16 -quiet -gravity center ${fbase}_alpha.tif ${dir}.tif ${dir}_tmp.tif
    mv ${dir}_tmp.tif ${dir}.tif
  done
}

base=$(pwd)

with_darks=0
if [ -d darks ]; then
  flatten_bdf_raws darks
  with_darks=1
  cd $base
fi

with_biases=0
if [ -d biases ]; then
  flatten_bdf_raws biases
  with_biases=1
  cd $base
fi

with_flats=0
if [ -d flats ]; then
  flatten_bdf_raws flats
  with_flats=1  
  cd $base
fi

## flat - bias
if [[ $with_biases -eq 1 && $with_flats -eq 1 ]];then
  convert -quiet biases/biases.tif flats/flats.tif -compose minus -composite divisor_tmp.tif
  convert divisor_tmp.tif -alpha off divisor.tif
  rm divisor_tmp.tif
## flats only
elif [[ $with_biases -eq 0 && $with_flats -eq 1 ]];then
  convert flats/flats.tif -alpha off divisor.tif
fi

cd lights
ls *.RW2 | xargs -n 1 sh -c 'darktable-cli $0 ../general.xmp ${0%.*}.tif --bpp 16'

imgs=( $(ls *.tif) )

for img in ${imgs[*]}; do
  fbase=$(basename $img .tif)
  rotation=$(exiftool -Rotation $img)
  
  if [[ "$rotation" =~ "90 CW" ]]; then
    convert -quiet $img -rotate 270 ${fbase}_land.tif
  elif [[ "$rotation" =~ "270 CW" ]]; then
    convert -quiet $img -rotate 90 ${fbase}_land.tif
  else
    cp $img ${fbase}_land.tif
  fi


  ## (light - dark) / (flat - bias) bzw.  (light - dark) / flat
  if [[ $with_darks -eq 1 && $with_flats -eq 1 ]];then
    convert -quiet ../darks/darks.tif ${fbase}_land.tif -compose minus -composite ${fbase}_land_tmp.tif
    convert -quiet ../divisor.tif ${fbase}_land_tmp.tif -compose divide -composite ${fbase}_corr.tif
    rm ${fbase}_land_tmp.tif
  ## light - dark
  elif [[ $with_darks -eq 1 ]];then
    convert -quiet ../darks/darks.tif ${fbase}_land.tif -compose minus -composite ${fbase}_corr.tif
  ## light / (flat - bias) bzw. light / flat
  elif [[ $with_darks -eq 0 && $with_flats -eq 1 ]];then
    convert -quiet ../divisor.tif ${fbase}_land.tif -compose divide -composite ${fbase}_corr.tif
  ## No correction at all
  else
    mv ${fbase}_land.tif ${fbase}_corr.tif
  fi
  
  ## combine the lightest pixels from each picture
  if [ $img == ${imgs[0]} ]; then
    mv ${fbase}_corr.tif lights.tif
  else 
    convert ${fbase}_corr.tif lights.tif -compose lighten -composite lights_tmp.tif
    mv lights_tmp.tif  lights.tif
  fi
done

cd $base

Mit dem Ergebnis war ich aber überhaupt nicht zufrieden:

  • die Leuchtspur des Schiffes unten links ist völlig ausgebrannt
  • rotes Rauschen, besonders in den Ecken
  • der Hintergrund ist viel zu hell

Also habe ich verschiedene Kombinationen durch weglassen, der dark, bias und/oder flat Frames versucht (deshalb auch die vielen ‘if’-Verzweigungen). Das beste Ergenis lieferte nur das dark von den lights abgezogen und dann aufhellen:

Das sieht schon deutlich besser aus. So wirklich zufrieden war ich aber immer noch nicht. Vor allem durch das ständige rausschreiben der Zwischenschritte und wieder einlesen, ist das Skript nicht besondes schnell.

2. Versuch in Python

Zum Beschleunigen muss also alles im Arbeitsspeicher verechnet werden, ohne ständigen Input/Output von der Festplatte. Meine Wahl fiel auf Python mit rawpy und numpy. Das hat den weiteren Vorteil, dass die interne Verrechnung in Fließkommazahlen erfolgt anstatt mit Ganzzahlen wie oben.

#!/usr/bin/python3

import os
import argparse
import rawpy as rp
import numpy as np
import tifffile
from glob import glob

## for random file name generation
import random
import string

def getImageData(path, extension='', verbose=0):
  files = sorted(glob("%s/*%s" %(path, extension))) # key=os.path.getmtime
  for f in files:
    with rp.imread(f) as raw:
      if verbose > 0:
        print("Read '%s'" %f)
      ## "no_auto_scale=True" macht es zu dunkel
      default_kwargs = dict(output_bps=16,
                            use_camera_wb=False,
                            use_auto_wb=False,
                            no_auto_bright=True,
                            user_flip=0)
      img = raw.postprocess(**default_kwargs)
      if 'data' in locals():
        data = np.append(data, np.expand_dims(img, axis=3), axis=3)
      else:
        data = np.expand_dims(img, axis=3)
  return(data)

## https://e2eml.school/convert_rgb_to_grayscale.html
def luminance(img, depth=16, perceptual=True, gamma=False):
  r = img[:,:,0].copy()
  g = img[:,:,1].copy()
  b = img[:,:,2].copy()

  if gamma == True:
    r = r / (2 ** depth - 1)
    g = g / (2 ** depth - 1)
    b = b / (2 ** depth - 1)
    r[r <= 0.04045] = r[r <= 0.04045] / 12.92
    r[r > 0.04045]  = ((r[r > 0.04045] + 0.055) / 1.055) ** 2.4
    g[g <= 0.04045] = g[g <= 0.04045] / 12.92
    g[g > 0.04045]  = ((g[g > 0.04045] + 0.055) / 1.055) ** 2.4
    b[b <= 0.04045] = b[b <= 0.04045] / 12.92
    b[b > 0.04045]  = ((b[b > 0.04045] + 0.055) / 1.055) ** 2.4
    r = r * (2 ** depth - 1)
    g = g * (2 ** depth - 1)
    b = b * (2 ** depth - 1)
    
  if perceptual == True:
    l = 0.2126 * r +  0.7152 * g +  0.0722 * b
  else:
    l = (r + g + b) / 3
  return(l)

def main():
  parser = argparse.ArgumentParser(description='Stack a bunch of startrail images')
  parser.add_argument('-v', dest='verbose', action='count', default=0, help='verbosity level')
  parser.add_argument('-g', dest='gamma', action='store_true', default=False, help='use gamma correction to find brighter pixels (default=False)')
  parser.add_argument('-l', dest='dir_lights', metavar='DIR', default='lights', help='Folder containing lights (default="lights")')
  parser.add_argument('-d', dest='dir_darks', metavar='DIR', default='darks', help='Folder containing darks (default="darks")')
  parser.add_argument('-f', dest='dir_flats', metavar='DIR', default='flats', help='Folder containing flats (default="flats")')
  parser.add_argument('-b', dest='dir_biases', metavar='DIR', default='biases', help='Folder containing biases (default="biases")')
  parser.add_argument('-e', dest='extension', metavar='EXTENSION', default='.RW2', help='Input filename extension (default=".RW2")')
  parser.add_argument('-o', dest='outfile', metavar='FILE', default='result.tif', help='Output filename (default="result.tif")')

  args = parser.parse_args()

  if args.verbose > 1:
    prefix_out=''.join(random.choices(string.ascii_uppercase + string.digits, k=8))

  images = getImageData(args.dir_biases, args.extension, args.verbose)
  bias = np.average(images, 3)
  images = getImageData(args.dir_flats, args.extension, args.verbose)
  flat = np.average(images, 3)
  images = getImageData(args.dir_darks, args.extension, args.verbose)
  dark = np.average(images, 3)

  ## Calculate the divisor for flat/bias correction
  ## was ist richtig?
  ## https://forum.astronomie.de/threads/bias-darks-und-flats-ein-ueberblick.305382/post-1570605
  ## (light - dark) / (flat - bias)

  divisor = flat - bias
  ## set pixels with higher bias values to "black"
  divisor[divisor < 0] = 65535
  divisor = divisor + 1

  lights = getImageData(args.dir_lights, args.extension, args.verbose)

  for i in range(lights.shape[3]):
    if args.verbose > 0:
      print("Process %i/%i frames" % (i+1, lights.shape[3]))
    corr = (lights[:, :, :, i] - dark ) / divisor
    corr[corr < 0] = 0
    corr_lum = np.repeat(luminance(corr, gamma=args.gamma)[:, :, np.newaxis], 3, axis=2)
    if 'final' in locals():
      final[corr_lum > final_lum] = corr[corr_lum > final_lum]
      final_lum = np.repeat(luminance(final, gamma=args.gamma)[:, :, np.newaxis], 3, axis=2)
    else:
      final = corr.copy()
      final_lum = corr_lum

    if args.verbose > 1:
      fout = "%s_%04i.tif" % (prefix_out, i)
      corr = (corr - np.min(corr)) / (np.max(corr) - np.min(corr)) * 65535
      tifffile.imsave(fout, np.array(corr, dtype="uint16"))
      print("Saved corrected image as '%s'" % fout)
      
  final = (final - np.min(final)) / (np.max(final) - np.min(final)) * 65535

  tifffile.imsave(args.outfile, np.array(final, dtype="uint16"))
  if args.verbose > 0:
    print("Saved final image as '%s'" % args.outfile)
  
if __name__ == '__main__':
  main()

Das Skript ist viel schneller und das Ergebnis ist meiner Meinung nach deutlich besser, um es in Gimp weiter zu bearbeiten.

Weiterverarbeitung mit Gimp

Mein Ablauf danach in Gimp sah dann aus wie folgt:

  1. Despakle (Flecken entfernen), um Hot/dead Pixel zu entfernen.
  2. Layer 4x kopieren (den Vordergrund ausschneiden), in jede Richtung um 2-5 Pixel verschieben und mit “Lighten” (Nur Aufhellen) wieder zusammenfügen. Das vergrößert die Sternenspuren und entfernt die Lücken zwischen den einzelnen Aufnahmen:
  1. Lichter und Schatten anpassen.
  2. Unschärfe maskieren
  3. Selektiver Weichzeichner
  4. Kurven anpassen

Ich denke das Ergebnis kann sich trotz Dämmerung und aufgehendem Vollmond sehen lassen. Wenn jemand eine bessere Methode mit Opensource-Software unter Linux hat, freue ich mich über eine Nachricht.

Ergänzung 5.11.

Im Forum von astronomie.de hab ich den Tip bekommen es mal mit Siril zu versuchen. Wenn man Bilder ohne Sternenregistierung und der Methode “Pixel Maximum Stacking” stacked, bekommt man das beste Ergebnis. Also alle obige Mühe war nur zum Testen und Verständnis ;-)

weitere Quellen

Anregungen zur Astrofoto-Bearbeitung mit Gimp hatte ich von focustoinfinity.de


See also