Sun Photos

Testing my new AVX mount. With a level bubble, compass and GPS coordinates, I was able to get it aligned fairly accurately, enough to capture the sun and have it remain in the field of view on the order of minutes.

The setup. A 127mm Maksutov-Cass on a Celestron AVX mount, with a Baden solar filter. The camera used is a DSLR Canon Rebel T3.

When I took image, I was a little disappointed, there were some defects (which at first I thought were sunspots :-( ). This is due to dust motes.

Sun. You can see defects in image

Unfortuntaley, I didn't take a flat field. However, I did have a few images of the sun that were shifted, since my tracking wasn't perfect. So I masked out the dust motes:

Masking defects

And chose two images decently shifted apart, and re-shifted them. There are sophisticated methods out there for this but for this application, I decided just to shift manually and look at the absolute difference of the image and the second image reshifted. When they align, you should see a uniform noisy sun. (Else, you'll see some bright regions which represent the overlap)

Taking two images, top left and top right are two images where the sun moved in field of view. If I shift them back by the correct amount and take difference, we see a nice noisy image as lower left figure.

After creating a mask and knowing the shift, you can then write some quick code to average the two together, as below:

def combinephotos(imgs, shifts, mask):
    imgresults = np.zeros_like(imgs[0],dtype=float)
    imgcts = np.zeros_like(imgs[0],dtype=float)

    for i in range(len(imgs)):
        imgresults += np.roll(np.roll(imgs[i]*mask,shifts[i][0], axis=0), shifts[i][1],
                axis=1)
        imgcts += np.roll(np.roll(mask,shifts[i][0], axis=0), shifts[i][1],
                axis=1)

    imgresults /= imgcts

    # parts that average more than one image, add poisson noise
    #w = np.where(imgcts > 1)
    #imgresults[w] = np.random.poisson(imgresults[w])
    return imgresults, imgcts

shifts_array = [
        [0,0],
        [-28, -51],
        ]

img_final, imgcts = combinephotos([Grey1, Grey2], shifts_array, mask)

There are much more efficient methods of course, but this allows more control over your data, and better understand exactly what it is you want. For example, I don't like inpainting. This method does not require any inpainting or blurring whatsoever. This is the raw data. The only potential issue is the loss of pixel resolution from the error in the estimation of your shifts. However, there are ways to actually get subpixel resolution, something I'll ignore here because this image is already quite large!

The resultant number of images averaged perpixel. White represents two while black represents one. There are no zero regions.

The final image after combining.

Image saved to a JPEG.

Getting Hands Dirty in Python

I learned python about a year and a half ago and it's been a huge time saver. I've moved from a language where you had to write a library for everything to a language where a library for everything already exists. Here's something I wrote quickly this morning. Very easy, and you can do it too. If you're interested in better understanding your astronomy data for educational purposes, this post will help get you started.
Before we begin, you'll need to have installed these libraries:
Python
pylab
rawkit
Pillow
LibRaw (Without this library, we couldn't do this thanks to developers for all this hard work!)
I'm happy to provide more instructions if need be.
I finally bought a mount (Celestron AVX) and I'm beginning to take my own images.



Here is some quick code to do this. Note, some steps are not very 'pythonic'. I'm a strong believer in breaking conventions to better understand data, and worry about conforming to them when packaging code.
This allows for much quicker analysis and gaining a deeper understanding of the data, which is *much* more important than the code!!! (coding is easy, making sense of data is the challenge...)

In this code, you'll learn to:
1. Extract the raw data from your image
2. Extract each color from your image (if needed), and make a grey scale image
3. Rescale your image in different ways to try to better emphasize features without brining out the noise

This code should be fairly straightforward and is meant just as a motivational exercise. There are better ways to do these things, and nice statistics you can use to better understand your data (like histograms). I hope this motivates you in the same way it excites me to have so much control over my data!

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99

''' Astronomy grey image analysis ''' # This imports functions into global namespace # useful for quick hacking, NOT recommended for writing programs from pylab import * # makes plots appear if in interactive mode (comment out if not) ion() # used for saving images from PIL import Image # used for reading Raw images from rawkit.raw import Raw # choose your filename here import os.path filename = os.path.expanduser("~/pictures/17-03-18-custer-whirlpool/IMG_3787.CR2") # Grab the data, really easy im = Raw(filename) data = np.array(im.raw_image()) # The data is in a Bayer Color Matrix here for Canon, so let's # just quickly extract each color by choosing the right pixels # quick way to get from Bayer matrix. Note there are way more # sophisticated algorithms out there than this method... But this # yields the rawest untouched data possible # # On Rebel T3, pixels are: # [ R G ] # [ G B ] # so index every other row/column from the right start (0 or 1) # start:end:skip means start at start, end to end and skip every other row # if "end" not specified, then just go to end of array # so 0::2 means start at 0, go to the end, and skip every other row/column # data is indexed as data[row, colum]a (row goes up/down in image, column goes # left/right) R = data[0::2,0::2] # for green, you average the two together G = (data[0::2,1::2]+data[1::2,0::2])*.5 B = data[1::2,1::2] # sum colors to get greyscale (note this isn't 100% correct, you need # to weight colors appropriately, but this is just messing around first) Grey = R + G + B # take logarithm, sometimes this could look better LGrey = np.log10(Grey) # try other combinations Grey_squared = Grey**2 # so now, let's plot the image, plot it first without "vmin" and "vmax" # and mouse over the plot to see what values you see # when happy, choose a good "vmin" and "vmax", this worked for this image set_cmap("Greys_r") figure(4);clf(); imshow(LGrey[200:900,600:1300],vmin=3.8,vmax=3.9) figure(5);clf(); imshow(Grey[200:900,600:1300],vmin=10**3.8,vmax=10**3.9) figure(6);clf(); imshow(Grey_squared[200:900,600:1300], vmin=4e7, vmax=6e7) # now, if we want to save image, it needs to have its pixel values go from # 0 to 2**8-1 (255). You don't want a larger range since your eye can't really tell # the difference beyond that, roughly # normalizing function: def normalize(img, mn, mx, logimg=True): ''' normalize to a uint8 This is also known as byte scaling. ''' dynamic_range = 2**8-1 img = img-mn img /= (mx-mn) img *= dynamic_range img = np.minimum(np.maximum(img, 0), dynamic_range) img = img.astype(np.uint8) return img img_grey = normalize(Grey[200:900, 600:1300], 10**3.8, 10**3.9) limg_grey = normalize(LGrey[200:900, 600:1300], 3.8, 3.9) img_grey_squared = normalize(LGrey[200:900, 600:1300], 4e7, 6e7) # Image comes from the PIL library (imported at top) img = Image.fromarray(img_grey) limg = Image.fromarray(limg_grey) img_grey_squared = Image.fromarray(img_grey_squared) # saving img.save("../storage/01-whirlpool-oneimage.png") limg.save("../storage/01-whirlpool-log-oneimage.png") limg.save("../storage/01-whirlpool-squared-oneimage.png")

Finally, you can run the code by running it in ipython (interactive mode):

$ ipython
[1] : %run mycode.py

or just python:

$ python run mycode.py

The whirpool galaxy. 5 min exposure taken with a Canon Rebel T3 on an Orion Maksutov-Cass 127mm on a Celestron AVX mount, autoguided.

The whirpool galaxy, logarithm of image. 5 min exposure taken with a Canon Rebel T3 on an Orion Maksutov-Cass 127mm on a Celestron AVX mount, autoguided.

The whirpool galaxy, square of image. 5 min exposure taken with a Canon Rebel T3 on an Orion Maksutov-Cass 127mm on a Celestron AVX mount, autoguided.

Right now, all images look the same, but there will be cases where different transformations help (this is similar to playing with the linearization curve). We'll mess around with colors later.