Read FITS 16-bit integers as signed and expose BZERO/BSCALE in im.info

Tests/test_file_fits.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import struct
 from io import BytesIO
 
 import pytest
@@ -29,6 +30,20 @@ def test_gzip1() -> None:
         assert_image_equal_tofile(im, "Tests/images/m13.fits")
 
 
+def test_gzip_16bit_consecutive_samples() -> None:
+    # The decompressed payload of a GZIP_1 image may hold consecutive 2-byte
+    # samples rather than the widened 4-byte groups of Tests/images/m13_gzip.fits
+    with Image.open("Tests/images/gzip_1_16bit.fits") as im:
+        assert im.mode == "I"
+        assert [im.getpixel((x, 0)) for x in range(im.width)] == [
+            -32768,
+            -1,
+            0,
+            1,
+            32767,
+        ]
+
+
 def test_gzip_decode() -> None:
     with open("Tests/images/m13_gzip.fits", "rb") as fp:
         data = bytearray(fp.read())
@@ -93,3 +108,68 @@ def test_comment() -> None:
     image_data = b"SIMPLE  =                    T / comment string"
     with pytest.raises(OSError):
         FitsImagePlugin.FitsImageFile(BytesIO(image_data))
+
+
+def fits_image(bitpix: int, values: list[int], extra: list[str] = []) -> bytes:
+    cards = [
+        "SIMPLE  =                    T",
+        f"BITPIX  = {bitpix:20d}",
+        "NAXIS   =                    2",
+        f"NAXIS1  = {len(values):20d}",
+        "NAXIS2  =                    1",
+        *extra,
+        "END",
+    ]
+    header = b"".join(c.ljust(80).encode() for c in cards)
+    header += b" " * (-len(header) % 2880)
+    data_fmt = {8: "B", 16: "h", 32: "i"}[bitpix]
+    data = struct.pack(f">{len(values)}{data_fmt}", *values)
+    return header + data + b"\0" * (-len(data) % 2880)
+
+
+def test_16bit_is_signed() -> None:
+    # FITS 16-bit integers are signed
+    with Image.open(BytesIO(fits_image(16, [-32768, -1, 0, 1, 32767]))) as im:
+        assert im.mode == "I"
+        assert [im.getpixel((x, 0)) for x in range(im.width)] == [
+            -32768,
+            -1,
+            0,
+            1,
+            32767,
+        ]
+
+
+def test_16bit_scaling_keywords_in_info() -> None:
+    # Unsigned 16-bit data is stored as signed with BZERO=32768;
+    # the scaling keywords are exposed without being applied
+    with Image.open(
+        BytesIO(
+            fits_image(
+                16, [-32768, -32767, -1, 0, 32767], ["BZERO   =                32768"]
+            )
+        )
+    ) as im:
+        assert im.mode == "I"
+        assert im.info["BZERO"] == 32768
+        assert [im.getpixel((x, 0)) for x in range(im.width)] == [
+            -32768,
+            -32767,
+            -1,
+            0,
+            32767,
+        ]
+
+
+def test_bscale_keyword_in_info() -> None:
+    with Image.open(
+        BytesIO(
+            fits_image(
+                16,
+                [-32768, -1, 0, 1, 32767],
+                ["BSCALE  = 0.5", "BZERO   =                100"],
+            )
+        )
+    ) as im:
+        assert im.info["BZERO"] == 100
+        assert im.info["BSCALE"] == 0.5

src/PIL/FitsImagePlugin.py

@@ -112,8 +112,10 @@ def _parse_headers(
         if number_of_bits == 8:
             self._mode = rawmode = "L"
         elif number_of_bits == 16:
-            self._mode = "I;16"
-            rawmode = "I;16B"
+            # FITS 16-bit integers are signed (FITS Standard 4.0, Section 5.2);
+            # unsigned data is stored as signed with BZERO=32768
+            self._mode = "I"
+            rawmode = "I;16BS"
         elif number_of_bits == 32:
             self._mode = "I"
             rawmode = "I;32B"
@@ -124,6 +126,21 @@ def _parse_headers(
             msg = "Unsupported number of bits"
             raise OSError(msg)
 
+        # The physical value is BZERO + BSCALE * array value (FITS Standard 4.0,
+        # Section 4.4.2.5). Scaling is left to the caller; the keywords are
+        # exposed in im.info.
+        for keyword in (b"BZERO", b"BSCALE"):
+            value = headers.get(keyword)
+            if value is None:
+                continue
+            try:
+                number = float(value.replace(b"D", b"E"))
+            except ValueError:
+                continue
+            self.info[keyword.decode("ascii")] = (
+                int(number) if number.is_integer() else number
+            )
+
         args: tuple[str | int, ...]
         if decoder_name == "raw":
             args = (rawmode, 0, -1)
@@ -150,15 +167,37 @@ def decode(self, buffer: Image.DecoderInput) -> tuple[int, int]:
             rows = []
             offset = 0
             number_of_bytes = abs(self.args[0]) // 8
+            # The decompressed payload is either widened to 4-byte groups, as in
+            # existing files, or holds consecutive samples (see #10309)
+            widened = len(value) == self.state.xsize * self.state.ysize * 4
             for y in range(self.state.ysize):
                 row = bytearray()
                 for x in range(self.state.xsize):
-                    row += value[
-                        offset + 4 - 1 : offset + (4 - number_of_bytes) - 1 : -1
-                    ]
-                    offset += 4
+                    if number_of_bytes == 2:
+                        # 16-bit integers are signed; expand to the 4-byte "I" rawmode
+                        sample = (
+                            value[offset + 2 : offset + 4]
+                            if widened
+                            else value[offset : offset + 2]
+                        )
+                        pixel = int.from_bytes(sample, "big", signed=True)
+                        row += pixel.to_bytes(4, "little", signed=True)
+                    elif widened:
+                        row += value[offset + 4 - number_of_bytes : offset + 4][::-1]
+                    else:
+                        row += value[offset : offset + number_of_bytes][::-1]
+                    offset += 4 if widened else number_of_bytes
                 rows.append(row)
-        self.set_as_raw(bytes([pixel for row in rows[::-1] for pixel in row]))
+        data = bytes([pixel for row in rows[::-1] for pixel in row])
+        # The 4-byte rows are little-endian; the default rawmodes for "I" and "F"
+        # are native-endian, so pass explicit little-endian ones for big-endian hosts
+        if number_of_bytes == 2:
+            # 16-bit samples are expanded to signed 32-bit
+            self.set_as_raw(data, "I;32S")
+        elif self.mode == "F":
+            self.set_as_raw(data, "F;32F")
+        else:
+            self.set_as_raw(data)
         return -1, 0