Astrophotography Image Processing Workflows

Hi, my name is Tom Sherlock. I'm a user interface developer at Wolfram Research. Once again, I'm talking about my hobby of astrophotography. I'm not a professional astrophysicist. I just like to go out and take pictures of astronomical objects with my telescope and some equipment. And I thought I would share how I use Mathematica to make this process easier. So let's get started to do astronomy, and in particular, astrophotography involves a fair amount of data processing and you can get commercial packages to do this, but it's always more fun to roll your own and you have more control if you do it yourself. So that's why I use Mathematica. But there's a fair amount of work involved. It's not just a matter of going outside and pointing your camera up at the sky. And I've talked about aspects of this before, so I'm going to continue about this. So this is. I took some pictures here in May of this year of a bright galaxy called M101, also known as the Pinwheel Galaxy. And I also managed to capture a supernova which had just been discovered the previous week. So let's go down. So anyway, it is kind of tricky to do astronomical images of deep sky objects. And the reason why is because all these objects are very dim. If you've ever looked through a telescope with something like the Ring Nebula, unless you've got a very large telescope, what you're looking at is kind of a fuzzy blob. And if it's a little larger telescope, you can start to see the ring structure. But you almost never see any colors unless you've got a really large telescope, like a big Dobsonian scope. But you can capture this stuff with a smaller scope like I've got, and with an astronomical camera. But you have to do some work to get some results. And this isn't just work that amateurs have to do. Everyone has to do it. You know, even the scientists who run the James Webb Space Telescope, they have a fair amount of, of crunching that they have to do in order to get those nice pictures, let alone actual meaningful scientific images. So you always have to. The big battle is maximizing the signal and minimizing the noise. And because the, the signal is so low and the so high, it takes a fair amount of effort. So in order to do this, you have to take multiple images, sometimes through a sequence of filters. If you're trying to get a color picture, you have to take additional calibration images, which are images which capture the amount of noise that's in the image, and you have to subtract that noise out. And then that sequence of images you have to combine in a statistical manner to maximize the signal. And then additionally, you will probably have to do some image processing to sharpen and stretch and remove artifacts. But all that stuff can be done in Mathematica using fairly simple functions. So for this talk, I've tried to recreate as best I can the pipeline from a commercial package like Pix Insight, which is a popular astrophotography system. So let's go on to the next slide. So you can see that the workflow in order to get an image is fairly involved because it needs to be done on a sequence of images. And here I'm talking about monochrome images. If you want to produce a color image, you have to do this same workflow like three or four times through a different filter each time. And that will give you each color channel. And then you combine those at the very end to create a final image. But right now we're just going to keep it simple and talk about the monochrome workflow. And then you can just imagine that multiplied times four. So you start out with a set of what's called FITS data files. FITS is a standard astronomical data format. It's used by all astronomers. Last month I gave a talk about how you can download the FITS images coming from the well produced by the James Webb Space Telescope and create your own color images by combining the scientific data that they encapsulate. So we always start out with the sequence of these, which are a sequence taken of the same object. And each one can be like 15 to 1 minute, 15 seconds to 1 minute long. And the reason why we don't just take one long image is because it's hard to get your telescope to track accurately for a long amount of time, even with automatic guiding. So once you've got your raw images, then you have to calibrate them, which is the process I mentioned before of subtracting out the noise. So you have to take additional images just to capture the noise. And this is actually something that's not too unusual. A lot of cameras that just regular consumer grade cameras, if you take a long exposure, they'll take an additional calibration image and subtract it out kind of behind your back. Once you've got the images calibrated, then you need to align them so they're on top of each other. And that's a little tricky because sometimes they can drift. Telescope alignment isn't perfect. And you have to find the right algorithm to line them up. Sometimes you have to do it by hand, as I talked about last month. Then you normalize them, you stack them, which is the statistical averaging process. You perform additional noise reduction, and then you stretch them out to see the parts of the image that you're really interested in, which is like the nebula or something. And you do final processing, you get the final image. So let's talk a little bit about this calibration process. So calibration means you're capturing additional images which are taken with the shutter closed or a cover over your telescope, and they measure noise characteristics of your camera. And there's basically three types that we'd like to talk about. There's what's called bias frames, and these are just estimate the uncertainty in reading a particular CCD chip cell. And these are typically very short frames, as short as you can get, maybe 5,000ths of a second. And the camera shutters close. So you're just, just estimating what's the uncertainty of reading a particular frame or a particular cell out of the ccd. Then there's dark frames. These estimate the thermal noise from the CCD chip. When you're taking a long exposure, the chip heats up and that produces a glow on the chip. And so you need to, to alleviate this, you typically refrigerate the chip and then you capture this and subtract it out. So the dark frames need to be taken at the same temperature and exposure length as the image frames. So for this reason, it's important to have a regulated cooling system on the CCD camera. And this is my camera, and this whole back part of it is the cooling system. And then you can also take flat frames, which I'm not describing in this particular talk, but this is the reason why you take flat frames. This is A picture of M51, the Whirlpool Galaxy. And you can see these big donuts here. These are tiny bits of dust on the filter that's on the CCD chip. And it produces a shadow in the shape of the telescope aperture. You typically take these after sundown. You get a uniformly illuminated surface, so you can see what the garbage is in your optical system. So the implementation here is all the code I'm going to use. So we're not going to go over that. I'm just going to describe it and then demonstrate it. So to get the bias frames, you have to take a sequence of images, which is what I have already done here. And I'll include that with the notebook. And I Have this code here called Compute Master Bias, which basically just averages out all the bias frames. So I've got 20 bias frames and I exported that into this master bias fit file. Same goes with dark frames. So the dark frames, those are the file names. And here the dark frames also take the name of the master bias file because each dark frame has to have that bias noise subtract from it. So, yeah, so it exported the master bias file and we're going to skip over the flat frame thing. But that would. I'll explain where that comes in the later in the process. So, Okay, so then after you have your noise frames, the bias and dark frames captured, you take what's called the light frames. So for each image subframe, you subtract out the bias frame and you do this using image subtract, and you then subtract out the dark frame and that subtracts out the bias noise and the dark noise. And then the numbers in the fits file are typically the actual counts pulled out of the CCD chips. And we like to work with the range 0 to 00, 0.0 to 1.0, which makes the processing easier. So we have to do image adjust. So if you bring in one of these frames just as it is, it looks all white. It's because it's probably in the range, you know, 100 to 30,000. So we have to do an image adjust and that will automatically normalize the image to this range. Once the frames are calibrated, then they need to be precisely aligned. And this is a pretty tricky process. And I'll describe why. For some images, good images, frankly, there are too many stars. And in order to align the images, you need to pick definite points to compute, alignment, transform. And if you have a bunch of stars which are about the same brightness, it's not going to work out very well because it's hard to find the same set of stars. The best type of algorithm creates basically artificial asterisms of three to five stars. And it computes a hash of for each asterism by computing the relative positions of the stars in each asterism. And it uses that to compare an asterism from one image to the next. And those match up pretty well. But that's a complex algorithm. And last year I demonstrated doing that by calling a Python routine which executes that algorithm. You can, and you can also just manually pick the stars yourself, which works quite well with highly detailed images. Last month I gave a little lecture to about the data science of doing this, and I did it with James Webb Space Telescope images. I produced some nice color images. So this particular notebook does. Since my images are fairly simple and they don't have too many stars, I can do the basically the brightness picking. So the algorithm that I use, which works fairly well for these images, you have to find at least three bright stars in each image. And I do that with max detect. Then I determine the star centroids using component measurements. And then once I have at least three stars in each image, and I do some other checking to make sure that they're pretty much the same stars, I find geometric transform. And that will compute the rigid transform that takes one image into subsequent image, and then I apply that to subsequent images using image perspective transformation. And that brings all the images into alignment with each other. And as I said, it works pretty good as long as there aren't too many stars in the image, and it works approximately as well as the Python routines. You'll see when I demonstrate this that there's some images that I can't find a transform for, but it still works pretty good. Hopefully, if I have enough time, I'll try implementing the full asterism algorithm in Wolfram language in the future. And then once you've got the images all aligned, then you can stack them up. Now, before we can do this, we have to basically adjust the histograms of each image so they match each other. And the reason why is demonstrated by this example here. We're using histogram transform, and PixInsight does the same thing. Of course they call it a little bit different, but we're using histogram transform. You can see these are two images, like an NMR image of someone's head. And if we didn't do, if we try averaging these two images, it would just make the bright parts here darker. And we need to have the images about the same brightness in all parts so that we're actually canceling out the noise and not canceling out the signal. So here this transformed this image so that it's approximately the same brightness values as this one. And then we can take a meaningful statistical average of these two. So for this, we have to apply histogram transforms to all of our aligned images with respect to a reference image, which I take to be the first image. And once that is done, then we can actually average all the images to eliminate the noise. And we use another function in Mathematica called Winsorized mean. And this is a robust averaging process that is used for data that has outlined data values. And this is useful for us because a Lot of times there's hot pixels or dead pixels in an image, and this eliminates them automatically. It will eliminate a certain percentage of these outline values at the tail ends of the of the histogram and replace them with the next closest values. So this is a good way to average out the noise and get the best possible signal. So what we're doing here is I take all the aligned images, extract the image data from it, flatten that out, and then do a winsorized mean on that whole array and that will give us combined data. And then I take that and I reconstruct the average image by partitioning it with the dimensions and constructing an image out of that. So the full pipeline is we do the aligning and stacking. And if we were using flats, we would use brightness equalize to eliminate dust artifacts with our master flat image. And then we'll use, we do a noise reduction pass, and we like to use a total variation filter. And that it just doesn't blur out the image. It still preserves rapid transits in the image, but it doesn't eliminate all of that data because some of it's valid. And then I implemented the IRAF ZSCALE stretch algorithm. And that lets you pick out the interesting parts of the histogram, which is where the, like the dim nebula might be located. And then you can also average out the sky glow by sampling points randomly selected in the image and subtract that out and do a final image adjust. So that's a lot of talking. Let's see if we can make this work here. Here's The Pinwheel Galaxy M101, taken from my backyard in Urbana. This is a performance professional photograph. Now, I didn't take this one, so we can see what it looks like based on our galaxy data. So here I've got my sub file name, the file names of all the light frames. And as you can see, there's about 20 images I took, each one 30 seconds long. So I'm going to run. I got all the stuff I talked about built into this thing called Image Pipeline. So this is doing that whole process I just talked about. So it's aligning the files. Here's our dark frame, our bias frame. This is our target frame. It's writing out these reference images, the reference image, and then the transformed images. You can see in this one, it failed to find a transform, so it's skipping that one. So it's writing out these. These are the aligned images. So here it's bringing back in the aligned files and it's doing the histogram transform pass. And, okay, so it's, it's applied all the histogram transforms. It. It's doing the winsorized mean, it's eliminating the 5% of the outline pixels. We could probably tweak that. And it's constructing the final image. And it's removing the noise using total variation filter and doing the Z scale stretching. So this is the image I got, and it also removed the sky glow. And this, this image picked up the. The supernova that was just discovered the week before, which is, I believe, this star, but I think it may also be this one. It was a bright star that appeared. So that's an example of how you get the Pinwheel Galaxy. And you can do further work on that. Nice, expensive CCD cameras have a debloating feature, but I don't have one of those. So I have a cheap CCD camera. So you can see that the stars in my image are kind of fat. And that is because the longer the exposure you do, the stars will swell up. And it's unaesthetic. So you can fix that, luckily with software. So I just wrote a little routine to debloat the stars. And this uses some additional image processing functionality built into Mathematica. So what it's doing here makes a mask by thresholding the image. So this picks out just the brightest elements of the image, and it inflates the stars by a little bit to make a paint mask. So this is an image which just has like slightly fatter versions of the stars on it. And then it calls our inpaint routine, which, with the image and the paint mask. And this will fill in everywhere in the target image with its surrounding pixels. So basically erases all the stars from the image. Then we can sharpen up the background. And then we take our stars in the. In our mask here, and we shrink them by doing erosion on it with a disk matrix, which does an erosion, a uniform circular erosion, and it shrinks them down and gives them a little bit of a blur so that edge artifacts aren't visible. And it adds that back in. And there's still a few little creaky pieces here, but the stars look sharper. And this is closer to what you'd see in a better picture of the galaxy. So that is about it. There's references here, links to the posting I did a couple months ago about manually aligning James Webb Space Telescope images. And the other posting I did which describes this with an additional example of another galaxy as well. So that's it. And if, if everyone has any questions, this is a good time to ask them. Okay, a couple questions here. I didn't see these. Have you ever tried to do stacking out of a video instead of single photos? Yes, I have. And this works best for planetary photography and I've done some examples of that in previous lectures, but that works better for like planetary photography. You have a video and you take that of like Mars or Jupyter, and by combining thousands of images you can get a fairly sharp image. But you know, it doesn't work miracles. You have to have a fairly steady image to begin with and then, and then it will, it will clean itself up. How long does it take to process one image using the image pipeline? Well, I just processed it, so it took about 30 seconds for 20 subframes. They're not super high resolution subframes, I might add. You can get cameras which are much higher resolution and those would take longer, of course. Did you say what camera you have? Yes, I have the cheapest possible CCD camera. It's an Orion Starshoot G4. It's the cheapest monochrome camera with a filter wheel on it. And it's not a thousands of dollars CCD camera, but it does okay for, for what I like to do. And as, and since the images are not super high resolution, they're pretty easy to, to work with. And it also gives me good opportunities to, to improve things using Mathematica code. And Jeff said, do you have any new targets planned from your darker house? Yes, I'm going to be relocating to a darker sky site with lots more view of the, of the sky instead of that narrow slot that I've got right now in my present place in Urbana. So yes, I'm going to probably take some pictures of the Eagle Nebula and a few other Southern sky objects instead of the few that you keep seeing during these lectures, which are the ones that pass right over my house. So yeah, I've got a lot of plans for that and new projects and all sorts of things that I'm working on in my copious free time. See? Any other questions? Let's see. Yeah, next time Mars comes around, I'll probably be taking some, some of this video of, of Mars. Last time I took a lot of pictures of Mars and I managed to get a kind of a fuzzy picture of it, which is a lot better than the thing that I was looking at through my telescope, which was like something bouncing around like that. But it was. But I think I can probably do better than that. So I'm I've got some plans and maybe a better camera. A better camera that I'm working on. I tend to build my equipment, so I. It's. I'm. I'm super cheap when it comes to this stuff. So that's why I tend to build stuff rather than buy stuff. Okay, I think we are just about done, so I am going to leave here in a minute. If there's any further questions, please ask them now. I will be uploading this notebook and the data so you can access it soon, and then you can play around with this code and hopefully improve on what I have done. But other than that, I will be saying goodbye and have fun with your own astrophotography adventures.