Using ImageJ / Fiji for astro image processing

[vlaiv]

I planned to this thread a long time ago, but somehow never got to it.

Envisioned like sort of tutorial (or collection there of), I'll post here some techniques on how to process astro images, or better call it data reduction stage (prior to actual processing step and up to stacking phase). I invite all interested members to give their suggestions on methods they use with this software package. Any questions or need for clarification are welcome of course.

I'll also include some additional plugins that I wrote for this software package for everyone to use.

So what is ImageJ / Fiji / AstroImageJ? It is open source / free software package for scientific image manipulation. It's written in Java (so you can run it on OS of choice), and it features nice plugin repository and scripting capabilities (you can write your own scripts or plugins for it). It has a lot of features out of the box. ImageJ is basic / "naked" distribution of software. Fiji is just ImageJ (as acronym says) but loaded with bunch of most used and very useful plugins - this version I recommend for installation.

AstroImageJ is separate project, I believe - it is also ImageJ but comes loaded with astronomy useful plugins and macros and its modified for astro use. For anyone looking to do a bit of photometry / astrometry and plate solving / measuring exo planets and such - this is package that I would recommend - it should support all features present in Fiji as well - but I have not tried it (probably requires manual install). So for my initial "tutorial-lets", I will focus on Fiji. You can install them all side by side, as installation is merely unzipping to a user chosen directory / folder (it comes with JRE so I don't think that you need to install Java separately - but not sure about that one, I had Java Runtime already on my computer).

Here is website that you can download software from:

https://imagej.net/Fiji/Downloads

(browse website for pure ImageJ download, or search google for AstroImageJ)

And here is (as promised in another thread), first of tutlets (can we name short tutorials like that - I kind of like the name ...)

Basic image calibration with Fiji

(note: tutorial, or tutlet :D, expects that you use fits files, have set point cooling and such, for best results. You can use method outlined with different acquisition methods / data, if you convert to fits format, but but be aware that that will have its own peculiarities and drawbacks)

1. Creating master dark

This is basic version of master dark file with no fancy algorithms involved - it is ok for basic usage, but for any serious work I would recommend sigma rejection on darks (due to cosmic rays / radiation hits - I get those regularly in some subs, particularly with longer exposure times).

- Start up Fiji
- Use File / Import / Image Sequence

- Select folder with your dark subs / any sub from the list
- Set import options

Here you can choose to load certain number of files, starting from image number, with increments of and file name containing (so you can basically load mix / subset of files if you choose to do so). It will tell you total size of loaded data, and here is one
important thing that you need to decide.

Using virtual stacks - if you check this option, files will be opened but not really loaded in memory until you perform operation on them - this is useful if total data set size is larger than your computer RAM and you want to avoid swapping to disk. You don't have to worry about this option if your data set is smaller and you have enough computer RAM for all operations.

For time being, let's leave that check box unmarked.

- Convert to 32 bit format via Image / Type / 32 bit

(note: this operation is going to double memory consumption of your stack, so you can run out of RAM and start swapping at this point if your data set is large. Using virtual stack will not help here as there is manipulation of data and whole virtual stack will still be loaded in memory when converting to 32 bit, if you are bound by ram and have a large set - there are ways to go about it when creating simple average stack - will explain in later tutorials).

- Finally produce master dark as a stack of slices by - Image / Stacks / Z project

And select Average method

This will create new image that is your master dark, use File / Save as / Fits to save it to disk and close all opened images.

2. Creating master flat

- repeat process for flat darks - and save that as well (close all)

- repeat process for master flat up to point where you should save master flat - don't save it yet

- close stack of flats that you used to create master flat

- open saved flat dark via File / Open ...

- Now select following action: Process / Image calculator

And set following options - Image1 should be your stacked Flat (one that you did not save), operation is subtraction, Image2 should be your loaded dark_flat and two check boxes ticked

After you do this - new image will be created - save that image as your master flat, and close all ...

3. Calibrating your subs

- Open stack of your subs
- Convert stack to 32 bit precision
- Load master dark
- Load master flat
- Subtract master dark via Process / Image calculator, select subs stack as image1, select operation subtract, select image2 to be master dark, check 32 bit precision, leave create new image unchecked (we will do calculations "in place" for this). When prompted perform operation on whole stack
- Divide with master flat via Process / Image calculator, select subs stack as image1, select operation divide, select image2 to be master flat, again check 32 bit precision, and uncheck new image. Perform on whole stack.

- Finally perform File / Save as / Image sequence on stack of subs, enter name that you want and numbering scheme (or use slice names), and select output directory to save your calibrated subs

done.

You can open these calibrated subs in software of choice and align and stack them - just make sure that you don't repeat any calibration steps as these are calibrated frames.

For quick software bin if your camera is oversampling, and to gain some SNR and reduce file size, use:

Image / Transform / Bin as last step prior to saving result subs

Select X and Y shrink values to be same (either 2 for x2 bin, or 3 for x3 bin) and average as your method. This will do the binning on your stack and just save it after.

[vlaiv]

Let's continue with these tutlets, in last tutorial I mentioned Sigma Reject stacking as a method, and want to expand on that.

Sigma Reject Z Project

When taking darks - because they are usually long exposure (to match lights) unlike flats and flat darks - it can happen that stray particles can hit sensor and cause anomalous values in one or few adjacent pixels. This also happens when we take lights - and we often use sigma reject method of stacking to deal with such cases.

ImageJ does not come with suitable plugin for this stacking method, so I wrote one.

You can download it from here: http://serve.trimacka.net/astro/imagej_plugins/Sigma_Reject_Z_Project.class

Installation is simple - open folder where you installed ImageJ/Fiji and look for plugins sub folder. In that folder create another folder with meaningful name where you will be placing plugins related to astronomy - I call my Astro. Just copy above file in that folder and restart ImageJ/Fiji. That's it - it should show up under Plugins menu - Sub menu item will be called by the same name as folder you created, like this:

How does it work?

Review calibration process, and creation of master dark from previous tutlet. Instead of using Image / Stacks / Z project command - use Plugins / Astro / Sigma Reject Z Project command.

(note: make sure you call command on stack of images, and that stack is 32bit float format - otherwise I'm not sure it will work properly)

It will bring up following dialog:

What do all these options do?

Sigma represent clip value of standard deviation. Each pixel in stack of images is averaged in regular stacking (all pixels in exact same position across subs), and if you look at them the will form "population", or sequence of pixel values. Such sequence of values forms certain distribution, and it should be close to Gaussian distribution. When we find mean value of Gaussian distribution and standard deviation, we can assert certain things. One of those assertions: what is the probability that certain value falls in range Mean +/- Standard deviation of given distribution. This graph sums it nicely:

So any particular value that belongs to this distribution has 99.7% chance that it will be larger than mean - 3*stddev and smaller than mean+3*stddev.

This is something that we can exploit to remove anomalous values from distribution - outliers. Above sigma reject works like this: calculates mean and stddev of pixel from all subs. Then it goes thru each pixel and looks if particular value is higher or lower than mean value + or - stddev that you specified as Sigma parameter (choosing correct sigma depends on number of subs, but Sigma 3 is pretty good value, as most of the time we will not stack 1000 subs - where on average we can expect 3 pixels to be anomalous). After first round of rejection - we calculate stddev and mean of remaining pixels (not excluded) and we repeat process again - comparing all pixels with these new values and doing further rejections.

How many rounds are we performing is specified by Iterations parameter.

Last two options are helpful for further analysis of data, but you don't need them for calibration.

Produce rejection map - creates additional image - one which pixels hold value of how many values are rejected for that particular pixel - ideally you want dark image full of 0 values with sporadic valued of 1 (often) and 2 (very rarely).

Produce stddev map - creates additional image - this one contains standard deviation of remaining (not rejected pixel values) - handy if you want to estimate noise of your distribution, calculate SNR per pixel for stack of lights, etc ...

[vlaiv]

Tips and tricks volume 1

Here is small list of useful tricks related to above section. All trick are related to the fact that ImageJ needs to load whole stack in memory to do some operations, and if you are using CMOS sensor and short (ish) exposures, you can end up with heaps of data to process. I run my instance of Windows with minimal swap (almost disabled) as I like responsive system, and don't like idea of windows swapping out parts of RAM to disk at is leisure (although I have snappy SSD for quite some time now, it's been habit of mine since HDD days). Disabling swap is not something that I would recommend unless you know what you are doing and what are implications and limitations of that.

Be as it may, I often find myself dealing with hundreds of subs from ASI1600. Single 16bit sub has about 32MB in size. Convert it to 32bit (as you really should for any sort of processing) and it grows double in size to ~64MB. If you need to deal with ~240 subs (and that would be 4h worth of 1 minute exposures) - we are looking at 15 and a half GB of data. My astro laptop is packing 4GB RAM and desktop 8GB, so both are lacking memory resources for such large volume of data. This brings us to list of tips - all are dealing with issues of memory usage and the way to get above accomplished if you lack RAM and don't want your computer to swap to infinity and back.

1. Creating master flat when you have large number of subs

This is straight forward average method of creating flats (no sigma reject - that will be covered in subsequent tips).

- Open your image sequence by File / Import / Image sequence - be sure to check Use Virtual Stack at this point

- Don't convert to 32bit at this point, but rather do Z Project by Image / Stack / Z Project, and select Addition as method. This is important, because if you select average - you will get 16bit result and we don't want that. If you select add, ImageJ will for some reason produce 32bit result (it might produce 16 bit result if number of subs is small and sum of every pixel fits in 16bit range - this is just guess since I have not yet had that case, anyway in that case, convert result to 32bit before proceeding).

- Take your 32bit Z project result, and perform Process / Math / Divide, and there enter number of subs that you just added together (number of flats that you stacked). Leave this image open, and close stack of flats.

- Open your flat dark subs the same way as flats and perform same steps to get stacked image (be sure to divide with number of flat darks rather than flats like in first example)

- In the end perform Process / Image calculator, select Flat as image1, and Flat dark as image2, operation is subtract, uncheck create new image, and check 32bit result. After this save flat image as your master_flat ...

I regularly use this method to produce master flat directly on laptop (has only 4GB RAM) as I usually shoot 256 flats / flat darks per filter - saves quite a bit of time copying files over to desktop for processing (64MB copied instead of ~30GB)

2. Creating master dark out of large number of subs with sigma reject method.

This one is a bit more involved, and you need to judge how large "chunk" of data that you work on will be. Depending on size of "chunk" you will need to repeat this number of times (total size / chunk size).

- Open your dark stack via File / Import / Image sequence - be sure to use virtual stack

- Now perform Image / New ...

Use 32bit, Fill with black, enter width and height of your dark subs (basically number of pixels your camera produces), and set Slices to 1

This will create new "empty" image that we will use to "put our results in"

- Select back stack and select Rectangle selection tool

- Draw rectangular region of full width starting at top of your image and of certain height - this will be your first "chunk", so height should be approximately height / number of chunks. Size of chunk should be such that it comfortably "sits" in RAM when converted to 32 bit. In above example of 256 dark subs, I would use at least 4 chunks (~16GB / 4 = ~4GB, suitable for 8GB machine).

- Press ctrl+T - this will bring up ROI manager window and add your selection to the list (we need this to accurately replicate selection later)

- Now use right button click on selection and choose Duplicate

When prompted, choose duplicate stack option

This will make a copy of selection as new stack of images - this time in memory.

- Perform Image / Type / 32bit on this new stack

- Perform Plugins / Astro (or however you named your astro plugins folder) / Sigma reject z project on this

- With resulting image perform Edit / Selection / Select All

- Then copy it to clipboard via Edit / Copy

- Select empty image from the start (one that will hold our result), and then in ROI Manager window click on selection that we added - it should mark exactly the same zone as was used in stack to copy section that we have just been working on

- After that do Edit / Paste - and result of Sigma reject stacking should be copied in our result image in correct place

(Note: you might be surprised to see that pasted area is white - something that you did not expect, but this is because display range of new image is 0-1, and your dark values are in ADU units - larger than 1. If you want to adjust display range use - Image / Adjust / Brightness Contrast option

Here you can set minimum and maximum display values. You can adjust here to your liking, but don't use buttons on the bottom and that way you will just change display of image and not actual pixel values)

- Close all images except original stack and result image

- Now you need to select another part / chunk of the original image, and you can do that by slightly overlapping with previous selection and going "down" in horizontal direction, or there is simpler method. "Grab" top central anchor of selection (small rectangle) with mouse+left click (hold) - and move it almost all the way to bottom edge of selection. Then grab bottom center anchor and move it further down until selection is high enough to represent single chunk.

- Again press ctrl + T to add new selection region to ROI manager, and repeat above steps until you have covered whole image. Note, new selection will be added to ROI manager list, so don't confuse them when pasting into result image (select proper one).

Calibrating to proper electron count values

Sometimes we want to measure electron values in our subs, and we want our subs to be calibrated so they represent true electron counts rather than arbitrary ADU values, or values in range 0-1.

To achieve this we need to include couple of steps.

1. Removing bit offset with cameras with lower bit ADC.

Some sensors have 12bit or 14bit ADC but produce images in 16bit mode. In this case actual pixel values in image are multiplied with 16 for 12bit ADC, or with 4 for 14bit ADC (there is general formula - multiplier = 2^(16bit - Xbit) where X is number of bits in ADC - it even works for 16 bit ADC: 2^(16-16) = 2^0 = 1 - hence multiplier is 1 )

This means that we need to divide subs with multiplier at some points in workflow. We do this by Image / Math / Divide - and enter above multiplier.

- Flats don't need to be handled like this, we will handle master flat in different manner.

- Darks need to be handled like this, so divide "whole stack" after you convert it to 32 bit (you can even do it while still 16 bit, since all pixel values will be divisible with multiplier without remainder - you can check this by examining pixel values in image - just move your mouse cursor over image and value under it will be displayed in ImageJ main window).

- Lights need to be handled the same way as darks - so divide with multiplier after you convert them to 32bit.

2. Normalizing master flat

Since master flat represents percentage of light reaching sensor in particular place, values for master flat need to be in 0-1 range. Here are couple methods how to scale master flat to that range.

2.1 Using max value

- Open your master flat and run Analyze / Measure

Note Max value in results. In case you don't see Max being listed in results, you can edit what sort of measurements appear when you run above command via Analyze / Set measurements

- Divide your master flat with Max value via Process / Math / Divide, and save your master flat

2.2 Using median / mean of selection

Above example will give you slightly higher values than real because there is some noise in master flat - and Max value is not true peak intensity, but true peak intensity + noise. If we want more precision we should use following approach:

- Open master flat, and use selection tool to make selection in central / brightest part of our master. You can use circular or rectangle selection, but try to include in selection only bright part of flat - avoid any dust circles or artifacts.

- Again use Analyze / Measure and note either mean or median value (use median if you can't completely avoid dust circles)

- Divide your master flat (but before you divide, deselect previous selection as divide command will work only on selected region - not what you want) with mean or median value (which ever you chose).

2.3 Using histogram peak

I've discovered that there is checkerboard pattern in my flats with ASI1600, so above methods are not quite precise. To deal with this particular problem, best way to determine scaling factor for your master flat would be:

- Open master flat and same as in previous example - make selection in central / brightest part of image.

- Run Analyze / Histogram command

As you see in the image, there are two peaks in histogram - signature of checkerboard pattern present. We want to select "middle / highest point" of right peak - one representing higher values in the image.

- Place mouse at center/top of rightmost peak and read value displayed in histogram window (bottom right - value=1670.728 in image example)

- Use this value to divide your master flat with.

3. Converting to electron count.

After you finish calibration of your lights, only thing left to do is convert from ADU values to electron values. For this, you simply need to multiply with e/ADU value for your camera (or gain setting that you've used) - Process / Math / Multiply.

You can look up this value in specs for CCD camera, or gain graph of your CMOS camera. Some software manages to read this value and puts it in fits header when capturing. I use SGP, and value is written in FITS header.

To access this fits header information - select Image / Show Info

You can read it from EGAIN value in result list.

 

 

[Gasman]

Hi Vlaiv

I was trying out Astroimagej last night and came across the Spectroscopy plugin but can’t seem to find how to use it?. The help file says Its easy to plot a spectrum from a fits file once the file is calibrated  using the Astro tools but doesn’t say how to do this would you have any idea please?.

thanks

Steve

[Merlin66]

Steve,

I appears that there is a spectroscopic plug in for ImageJ, haven't tried it but may be interesting....

http://www.astro.physik.uni-goettingen.de/~hessman/ImageJ/Book/Image-based Spectroscopy/install.html

[Gasman]

Thanks Ken I managed to install the Spectro plugin ok but a bit confusing how to proceed with it. Imagej simply says once the fits is calibrated then it’s easy to process the spectrum but no info on how to calibrate the fits to start with. Must be some help files somewhere in the mist I’ll delve further ?.

Steve

[vlaiv]

There are couple of types of calibration that you can perform on image, that I'm aware of.

1. Regular calibration that we use in image processing as a first step - this is needed to get only light signal out of the image without impact of other factors (bias / dark current, any shadows / vignetting, pixel sensitivity etc).

2. Spatial calibration of the image - this converts pixel "distances/widths" to some other units - this is useful for astrometric measurements, or in microscopy to determine size of objects - lengths.

You can do that by Analyze / Set scale menu item.

See following section for further details: https://imagej.net/Spatial_Calibration

AstroImageJ probably has plate solving and will set spatial scale when you plate solve your image, so any distance or size measurement that you do will be in appropriate units (like arc seconds).

3. Calibrate pixel values - this is what is described above when converting ADU to electron count (linear function), but you can use different function to perform value calibration - this is useful if you want your image for example to contain pixel values not in number of electrons but magnitude per arc second squared, or similar.

Again details are found here:

https://imagej.nih.gov/ij/docs/examples/calibration/

Now from reading

http://www.astro.physik.uni-goettingen.de/~hessman/ImageJ/Book/Image-based Spectroscopy/index.html

Section on calibration, only normal calibration is mentioned as a step (above point 1., or see the first post in this thread for calibration).

Spatial calibration in imageJ is simple one - using regular distance - this is not suitable to convert to wavelength for spectral curve because most spectrographs have some function that defines relationship between distance from 0 order image (or other known marker) to particular wavelength - that is not trivial (not constant like above mentioned spatial calibration).

As far as I can tell - Spectroscopy plugin is capable of doing this measurements, as text on given link states:

"Once a spectral image has been calibrated, the spectra can be displayed as traditional graphs. ImageJ normally can only produce profile averaged plots of image intensity versus pixel distance along the profile or some other distance quantity linearly related to pixels, so a special spectrogram plotter was created in the Spectroscopy package."

So with normal ImageJ, you can do 1. - calibrate image for light response, 2 - do spatial calibration, 3 - plot profile (just select distance measure tool and draw a line over spectrum and select Analyze/Plot profile), but such graph will be "distorted" - X axis will provide distance in imaging plane to zero order image but not actual wavelengths. So Spectroscopy plugin has a proper function to extract such graph and make X axis be actual wavelengths.

Have a read on above link for Spectroscopy plugin, and pay attention to this section:

http://www.astro.physik.uni-goettingen.de/~hessman/ImageJ/Book/Image-based Spectroscopy/index.html#imagej

It describes process quite thoroughly - and in first sentence there is clue of what sort of calibration they expect:

" After the raw images have been calibrated for the normal image effects (bias, dark-current, flatfield), " - so we are talking about point 1. calibration - regular light response calibration of the image.