PIL does not visualise FITS files

Question

I will start this post by saying that although I can code a little, I'm not super good at it yet so apologies if the answer is really obvious and I do not see it.

I have a series of .fits files that I need to visualise through PIL. I wanted to do it in an easy manner by just using the fact that PIL has support for FITS.

The images used are from the WISE database for the M1 nebulae.

When I tried to run the following code, all I get is a grainy image:

from PIL import Image

fits_file = Image.open('image.fits')
fits_file.show()

If I try to read it through the Astropy fits module instead, it sometimes returns the right image and sometimes returns a white image (which I believe is due to the min/max being over 255).

from astropy.io import fits
from PIL import Image

fits_file = fits.open('image.fits')
fits_img = fits_file[0].data
pil_img = Image.fromarray(fits_img)
pil_img.show()

Am I missing something as to why I cannot get all of the images to visualise correctly?

EDIT: I have decided to provide the images I get with my codes.

The first block of code shows the following image:

The second block of code shows the following image:

The original image as visualised by astropy + matplotlib:

EDIT 2: Following a link provided in the comments I changed the second code to:

from astropy.io import fits
from PIL import Image

fits_file = fits.open('image.fits')
fits_img = fits_file[0].data
pil_img = Image.fromarray(fits_img, mode="F")
pil_img.show()

This returns the same image as the first code.

Answer 1

I think PIL has a bug and reads the FITS file with the wrong byte order. If I open the file, convert to Numpy array, swap the byte order to accommodate the bug and then scale and save as PNG it works fine:

#!/usr/bin/env python3

import numpy as np
from PIL import Image

# Oopen sample FITS image and convert to Numpy array
im = Image.open('sample.fits')
na = np.array(im)

# BYTESWAP FOR SOME UNKNOWN REASON
nb = na.byteswap()

# Scale data to range 0..65535 so we can see it
nc = (nb - nb.min()) * 65535/(nb.max() - nb.min())

# Save as 16-bit PNG
Image.fromarray(nc.astype(np.uint16)).save('result.png')