Masking a FITS image

Hi,

sharing this because (1) I’ve found a lot of open questions on Google, (2) I want to hear feedback about how to do this better :slight_smile:

If I have a lot of FITS files, and I need to mask something out (e.g. a foreground whose stacking would cause annoyance), how do I do it? Well… by the power of Python!

Step #1: Create your masking file. Make sure the mask is the same size as the FITS file[1]. White is in, black is out, grey is partway. The file must be a grayscale PNG, i.e. single color channel.

Step #2: Run the following code:

import time
import numpy as np
from astropy.io import fits
from PIL import Image

start = time.perf_counter()

# Light
image_data = fits.getdata('/tmp/input.fit')

# Mask
#
# Assumes single grayscale channel. An RGB PNG will raise
# an error later, because the numpy array would have one more axis.
mask = np.array(Image.open('/tmp/star_mask.png'))
print('Data loaded after ', time.perf_counter() - start)

# Normalize it to [0..1], to permit useful multiplication.
mask = ((mask - np.amin(mask)) / np.amax(mask))

# Also reverse it, because FITS files and PNG files have different
# vertical order. Whoopsie.
mask = np.flip(mask, axis=0)

print('Mask normalized and reversed after ', time.perf_counter() - start)

# Mask the image.
masked_image = np.array(
  (np.multiply(image_data[0], mask),
   np.multiply(image_data[1], mask),
   np.multiply(image_data[2], mask)
  ),  
  dtype=np.uint16
)
print('Mask applied after ', time.perf_counter() - start)

fits.CompImageHDU(masked_image).writeto('/tmp/output.fit', overwrite=True)
print('File written after ', time.perf_counter() - start)

As a result, you’ll have a uint16 compressed fits file which will only include the parts you masked in.

Having this capability built into Siril directly would be wonderful, and presumably a lot more performant!

Hope this helps, and looking forward to any feedback.


[1] For some reason, the FITS file generated by Siril was 8 pixels wider and taller than my original raws. I’d love an explanation!:slight_smile: But you can quickly get the canonical number:

>>> from astropy.io import fits
>>> fits.open('your.fit').info()

Filename: /tmp/output.fit
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       4   ()      
  1  COMPRESSED_IMAGE    1 CompImageHDU     10   (6032, 4032, 3)   int16  

Hi,
Masking an image is to multiply an original image and a mask defined in the range [0,1], giving a new image (name it R). Then you transforms the R affecting to areas of R where R is not 0. Finally you must delete the influence of mask. And this is not a direct task: if mask was >0 you didn’t need to do anything, but in those pixels with 0 in mask you have to copy de original pixels of the original image. You can do this adding to R the product or the original image and the inverted mask.

I commented this just yesterday with an algorithm to do it by hand. Of course, compute this algorithm would be better :slight_smile: