Aligning sun images using Python

Fix it in post. It’s a common phrase in photography and one which is widely recognised as landing you in hot water if you’re not careful. Fix it in post-production, fix it after the fact. “Hm, I’ll just fix that in post.”

Well when I spent a day taking a few hundred photos of Venus crossing the sun last year with the intention of creating a time-lapse but without a tracking system, “fix it post” was my mantra by necessity. For 7 hours I followed the sun across the sky just by nudging my camera mount this way or that. As a result the sun was in a different place in every photo. The fact I wasn’t standing on the equator at the time meant it rotated too.

How to align all the photos? I puzzled over this for quite a while. My first tactic was to seek out existing software to do the job: RegiStax and DeepSkyStacker both offer alignment algorithms for astro images and both are freely available. But DeepSkyStacker is geared towards star fields, and after a couple of hours fruitlessly faffing about with RegiStax (I don’t think it likes big movements between frames), I decided to write my own bit of code to align all the sun images.

Some things to consider then:

• The sun is in a different place in every image;
• it’s also arbitrarily rotated throughout all the images; and
• over 7 hours, the sun itself will have apparently spun as we orbit it.

Sunspots are great reference features for an alignment task like this, but the third factor above needs to be remembered.  We can’t simply pile up the sunspots on top of each other because they move slightly over time.  Luckily there’s a way around this though.

A PYTHON SCRIPT FOR ALIGNING SUN IMAGES

Python is an excellent language for stuff like this. I’d actually never used it prior to this occasion, but I think its testament to the sheer usability of Python that I managed to come up with the following scripts in a day, and I’m hardly an accomplished programmer. Another of Python’s great assets is the fantastic Python Imaging Library (PIL), which makes opening and manipulating images a breeze.

THE TECHNIQUE

The general procedure is broadly divided into three stages:

1. Pre-processing
   • Rotate by 90 degrees if needed
   • Set white balance, exposure, contrast etc.
   • Export the files with a numerical index
2. A Python script to crop the sun out of the original picture
   • Open the image
   • Find the sun and crop around it
   • Save the image
3. A Python script to rotate the sun
   • Open the image
   • Find the centre coordinates of two sunspots
   • Find the angle of a line between the two sunspots
   • Rotate the sun to the desired orientation
   • Save the image
4. A Python script to scale the sun and create the final output
   • Scale the sun to the desired final size
   • Expand the black space if necessary
   • Save the final image

Cropping tightly around the sun in stage 2 means that the sun is automatically centred. Best of all, the gradual rotation of the sun is accounted for in stage 3 by the fact we are setting the *relative* position of the sunspots to each other rather than their absolute positions in the image frame.

Below is a more detailed description of each step. You’ll find the Python code at the bottom of the page.

Stage 1.  Pre-processing

First we need to rotate all images to the nearest 90 degree orientation to each other. This is because the code is only going to locate sunspots that are within a certain neighbourhood; for example if the camera was changed from a horizontal to a vertical orientation while shooting, then some of the images should be rotated by 90 degrees so as to approximately match the orientation of the sun in other frames. The actual location of the sun within the frame doesn’t matter.

Next we set all the usual image processing parameters such as white balance, exposure, contrast etc. so that the image looks as good as it can. Of course it’s important that all the images share the same settings for consistency.

All the photos then need to be exported to a directory somewhere, with a numerical index appended to the file name.

Stage 2.  A Python script to crop the sun out of the original picture (code below – “suncrop.py”)

Opening the image is quite straightforward. The PIL module can open loads of file types (mine were 8-bit TIFFs).

The sun image is zero-thresholded so that the surrounding sky is made all black, then another nifty PIL function returns the coordinates of the smallest bounding box which contains non-zero pixels. In other words it draws a box around the sun and cuts it out of the original image.

The image is saved and the process repeats for the remaining images in the sequence.

Stage 3.  A Python script to rotate the sun (code below – “sunrot.py”)

Take a look at the picture below. It’s a stack-up of all the images produced by my first script, and it shows how the sunspots all share a common neighbourhood. The size and position of this neighbourhood is what we give to this next script.

Composite ‘difference image’ of all the time-lapse frames. The rotation of the two main sunspots can clearly be seen. Each red box marks a neighbourhood in which we know the sunspot will definitely be.

Each image is opened as before. An inverted version is also generated, so that sunspots appear as white on black rather than the other way round.

Now for the main trick: finding the centre coordinates of two sunspots. The program has already been told where look for each spot in terms of a central pixel and a surrounding box of a given size (350×350 in my case). The inverted sun image is cropped to this small area and the result is zero-thresholded so that all pixels below a given value (~60) are set to zero. We now have a black square with a white spot somewhere inside it.

Next is to calculate the centre of mass of this small section of the image. Finding the centre of mass of the cropped pixels is a simple matrix operation (found in the SciPy Python module) which will return the pixel coordinates of where the image is ‘balanced’. Brighter pixels have more ‘mass’, so the centre of mass naturally moves towards regions with a higher overall brightness. Remember that the sunspot appears as a ‘heavy’ white spot on a background which is almost totally black. Assuming there aren’t two major spots on opposite sides of the cropped image, the centre of mass will be somewhere quite near the sunspot we’re trying to locate.

The image is re-cropped around the previously calculated centre of mass and then the centre of mass is calculated again; this process iterates until the centre of mass exactly matches the centre pixel of the cropped image. Finally the cropped image is shrunk down to a more sensible size (e.g. 64×64 pixels) which is likely to contain only the sunspot we’re interested in, and the centre of mass is calculated one last time. This ensures than any other spots, not dominant but influential, are cropped away and no longer shift the ‘balance’ of the pixels. Eventually we end up with exact coordinates for the centre of mass of our two chosen sunspots!

Given the location of the spots, some trigonometry (inverse tangent) can be used to find the angle of an imaginary line drawn between them.

The ‘sunspot angle’ is used to rotate the sun images to the same angular position

All that remains is to rotate the sun to the desired orientation. PIL has a nice built-in function to rotate an image by a specified number of degrees, which in this case is the difference between the sunspot angle calculated above and the angle we want the spots at. I shot my Transit of Venus time-lapse in Australia, so I chose to represent the sun in “southern hemisphere mode” (i.e. as we’d see it standing at the South Pole). I worked out that this meant a sunspot angle of about 50 degrees for the spots present on that day.

The rotated image is saved and the above is repeated for the remaining frames.

Stage 4.  A Python script to scale the sun and create the final output (code below – “sunscale.py”)

There are a couple more steps if you want to end up with frames of a certain size, such as full HD (1920×1080 pixels) for example. This involves scaling the rotated sun images and placing them on a black rectangle of the desired size. Again, PIL makes this nice and easy.

The final output image is saved and the above is repeated for the remaining frames.

An that’s it!  And here’s the proof that it works!

THE CODE

Here’s the code for the scripts described above. The first is a simple module used by the other files. They’re unlikely to work for a completely different set of images on first try, but with some tinkering it should. Feel free to copy and modify this for your own use :)

sun.py

import os, sys
from PIL import Image, ImageOps

def thresh1(x): #zero threshold for finding sunspots
 if x < 60:
 x = 0
 return x

def thresh2(x): #zero threshold for cropping around sun
 if x < 128:
 x = 0
 return x

def bbcrop(image, border):
 box1 = image.point(thresh2).getbbox()
 box2 = (box1[0]-border, box1[1]-border, box1[2]+border, box1[3]+border)
 return image.crop(box2)

suncrop.py

import os, sys, sun
from PIL import Image

for i in range (1,400):		#image index range

	infile = "img/in/sun-%04d.tif" % i
	outfile = "img/crop/sun-%04d.tif" % i

	im = Image.open(infile)
	print "\nFile %s (%s, %s, %d x %d)" % (infile, im.format, im.mode, im.size[0], im.size[1])

	print "Finding sun..."
	imtrim = sun.bbcrop(im, 200)

	print "Saving to %s" % outfile
	imtrim.save(outfile)

print "Done!\n"

sunrot.py

import os, sys, math, sun
from PIL import Image, ImageOps
from scipy import ndimage
import numpy as np

def findspot(cboxctr, cboxsize, shrinkboxsize, maxiter):

	cboxtl = ((cboxctr[0]-int(cboxsize[0]/2)) , (cboxctr[1]-int(cboxsize[1]/2)))
	cboxbr = ((cboxctr[0]+int(cboxsize[0]/2)) , (cboxctr[1]+int(cboxsize[1]/2)))

	cropbox = (cboxtl + cboxbr)
	ccrop = imthresh.crop(cropbox)

	#ccrop.show()

	done = 0
	shrunk = 0
	iteration = 0
	prescom = (-1,-1)
	while done == 0:

		carr = np.array(ccrop)
		com = ndimage.measurements.center_of_mass(carr)
		abscom = ( int(cboxtl[0]+com[1]) , int(cboxtl[1]+com[0]) )

		if abscom != cboxctr:

			cboxctr = (abscom[0], abscom[1])
			cboxtl = ((cboxctr[0]-int(cboxsize[0]/2)) , (cboxctr[1]-int(cboxsize[1]/2)))
			cboxbr = ((cboxctr[0]+int(cboxsize[0]/2)) , (cboxctr[1]+int(cboxsize[1]/2)))

			cropbox = (cboxtl + cboxbr)
			ccrop = imthresh.crop(cropbox)

			#ccrop.show()	#for debugging

			iteration += 1
			if iteration >= maxiter:
				print "%s: Max iterations reached!" % infile
				sys.exit()
				done = 1

		elif shrunk == 0:

			cboxsize = shrinkboxsize
			cboxctr = (abscom[0], abscom[1])
			cboxtl = ((cboxctr[0]-int(cboxsize[0]/2)) , (cboxctr[1]-int(cboxsize[1]/2)))
			cboxbr = ((cboxctr[0]+int(cboxsize[0]/2)) , (cboxctr[1]+int(cboxsize[1]/2)))

			cropbox = (cboxtl + cboxbr)
			ccrop = imthresh.crop(cropbox)

			shrunk = 1

		else:
			prescom = ( cboxtl[0]+com[1] , cboxtl[1]+com[0] )	#precise CoM (absolute)
			done = 1

	ccrop = imtrim.crop(cropbox)
	#ccrop.show()
	return prescom

def findangle(spot1, spot2):
	angle = math.degrees(math.atan((abs(spot1[1]-spot2[1]))/(abs(spot1[0]-spot2[0]))))
	return angle

######

desiredspotangle = 49+180	#desired sunspot angle in degrees

for i in range (1,400):		#image index range

	infile = "img/crop/sun-%04d.tif" % i
	outfile = "img/rot/sun-%04d.tif" % i

	im = Image.open(infile)
	print "\nFile %s (%s, %s, %d x %d)" % (infile, im.format, im.mode, im.size[0], im.size[1])

	imtrim = im
	immono = ImageOps.equalize(ImageOps.grayscale(imtrim))
	imneg = ImageOps.invert(imtrim)
	imthresh = imneg.point(sun.thresh1)

	spot1 = findspot( (1075,845), (350,350), (64,64), 10 )
	print "Sunspot 1 is at (%f, %f)" % (spot1[0], spot1[1])
	spot2 = findspot( (845,1105), (350,350), (64,64), 10 )
	print "Sunspot 2 is at (%f, %f)" % (spot2[0], spot2[1])

	spotangle = findangle(spot1,spot2)
	print "Angle between sunspots is %f degs" % spotangle

	print "Rotating..."
	im = imtrim.rotate((desiredspotangle-spotangle), Image.BICUBIC, 0)

	print "Trimming..."
	im = sun.bbcrop(im, 160)

	print "Saving to %s" % outfile
	im.save(outfile)

print "Done!\n"

sunscale.py

import os, sys, sun
from PIL import Image, ImageOps

for i in range (1,400):		#image index range

	infile = "img/rot/sun-%04d.tif" % i
	outfile = "img/scale/sun-%04d.tif" % i

	outdims = (1920,1080)
	imout = Image.new("RGB",outdims,0)

	im = Image.open(infile)
	print "\nFile %s (%s, %s, %d x %d)" % (infile, im.format, im.mode, im.size[0], im.size[1])

	print "Scaling..."
	mindim = min(outdims)
	im = sun.bbcrop(im, 100)
	im = ImageOps.fit(im, (mindim,mindim), Image.ANTIALIAS, 0, (0.5,0.5))

	print "Framing..."
	imout.paste(im, ( ((imout.size[0]-im.size[0])/2) ,0))

	print "Saving to %s" % outfile
	imout.save(outfile)

print "Done!\n"