Fix lossless decoding of point transformed and interleaved multi component scans

src/decoder.ts

@@ -228,16 +228,14 @@ export class Decoder {
       this.numComp = this.scan.numComp
       this.selection = this.scan.selection
 
-      if (this.numBytes === 1) {
-        if (this.numComp === 3) {
-          this.getter = this.getValueRGB
-          this.setter = this.setValueRGB
-          this.output = this.outputRGB
-        } else {
-          this.getter = this.getValue8
-          this.setter = this.setValue8
-          this.output = this.outputSingle
-        }
+      // How many components the scan carries decides the layout, not how wide a
+      // sample is: `decodeUnit` dispatches on numComp alone, so a 16 bit colour
+      // scan used to decode three components and write them through the single
+      // component output, which only ever wrote the first
+      if (this.numComp === 3) {
+        this.getter = this.getValueRGB
+        this.setter = this.setValueRGB
+        this.output = this.outputRGB
       } else {
         this.getter = this.getValue8
         this.setter = this.setValue8
@@ -295,7 +293,7 @@ export class Decoder {
         index[0] = 0
 
         for (i = 0; i < 10; i += 1) {
-          pred[i] = 1 << (this.precision - 1)
+          pred[i] = this.initialPrediction()
         }
 
         if (this.restartInterval === 0) {
@@ -340,9 +338,51 @@ export class Decoder {
       }
     } while (current !== 0xffd9 && this.xLoc < this.xDim && this.yLoc < this.yDim && scanNum === 0)
 
+    this.applyPointTransform()
+
     return this.outputData
   }
 
+  /**
+   * The point transform Al (T.81 A.4, the low order bits the encoder discarded
+   * before encoding).
+   */
+  getPointTransform(): number {
+    return this.scan.al
+  }
+
+  /**
+   * Prediction for the first sample of a scan, and for the first sample after a
+   * restart marker: 2^(P - Pt - 1) per T.81 H.1.2.1. Seeding with 2^(P - 1)
+   * instead offsets every sample of the frame, because lossless prediction
+   * chains additively from this value.
+   */
+  initialPrediction(): number {
+    return 1 << (this.frame.precision - this.getPointTransform() - 1)
+  }
+
+  /**
+   * Scales the decoded samples back up by the point transform. The encoder
+   * right shifted them by Al before encoding (T.81 A.4.1), so a frame decoded
+   * without this keeps only 1/2^Al of its dynamic range.
+   *
+   * This runs once the whole frame is decoded rather than in `setValue`,
+   * because prediction reads previously decoded samples back out of
+   * `outputData` and those reads have to stay in the point transformed domain
+   * the arithmetic runs in.
+   */
+  applyPointTransform(): void {
+    const pointTransform = this.getPointTransform()
+
+    if (!pointTransform || !this.outputData) {
+      return
+    }
+
+    for (let i = 0; i < this.outputData.length; i += 1) {
+      this.outputData[i] = this.outputData[i] << pointTransform
+    }
+  }
+
   decodeUnit(prev: number[], temp: number[], index: number[]): number {
     if (this.numComp === 1) {
       return this.decodeSingle(prev, temp, index)
@@ -381,50 +421,48 @@ export class Decoder {
     return (this.getPreviousX(compOffset) + this.getPreviousY(compOffset)) / 2
   }
 
+  /**
+   * Decodes one multi component sample of a lossless scan.
+   *
+   * This is `decodeSingle` applied to each component of an interleaved scan:
+   * every component keeps its own prediction, taken from its own neighbours
+   * through the selector, and contributes one Huffman coded difference per
+   * sample.
+   *
+   * It used to be a copy of the baseline DCT path - it read 63 AC coefficients
+   * per block against the AC Huffman table and multiplied by a quantization
+   * table. A lossless scan carries none of that, so after the first component's
+   * difference the entropy stream was read out of step and the remaining
+   * components never decoded: a colour frame came out with components 2 and 3
+   * stuck at the initial prediction (cornerstone3D#1116).
+   */
   decodeRGB(prev: number[], temp: number[], index: number[]) {
     if (this.selector === null) throw new Error("decode hasn't run yet")
 
-    let actab, dctab, qtab, ctrC, i, k, j
+    let ctrC, i
 
-    prev[0] = this.selector(0)
-    prev[1] = this.selector(1)
-    prev[2] = this.selector(2)
+    for (ctrC = 0; ctrC < this.numComp; ctrC += 1) {
+      prev[ctrC] = this.restarting ? this.initialPrediction() : this.selector(ctrC)
+    }
+
+    this.restarting = false
 
     for (ctrC = 0; ctrC < this.numComp; ctrC += 1) {
-      qtab = this.qTab[ctrC]
-      actab = this.acTab[ctrC]
-      dctab = this.dcTab[ctrC]
       for (i = 0; i < this.nBlock[ctrC]; i += 1) {
-        for (k = 0; k < this.IDCT_Source.length; k += 1) {
-          this.IDCT_Source[k] = 0
-        }
-
-        let value = this.getHuffmanValue(dctab, temp, index)
+        const value = this.getHuffmanValue(this.dcTab[ctrC], temp, index)
 
         if (value >= 0xff00) {
           return value
         }
 
-        prev[ctrC] = this.IDCT_Source[0] = prev[ctrC] + this.getn(index, value, temp, index)
-        this.IDCT_Source[0] *= qtab[0]
+        const n = this.getn(prev, value, temp, index, ctrC)
+        const nRestart = n >> 8
 
-        for (j = 1; j < 64; j += 1) {
-          value = this.getHuffmanValue(actab, temp, index)
-
-          if (value >= 0xff00) {
-            return value
-          }
-
-          j += value >> 4
-
-          if ((value & 0x0f) === 0) {
-            if (value >> 4 === 0) {
-              break
-            }
-          } else {
-            this.IDCT_Source[Decoder.IDCT_P[j]] = this.getn(index, value & 0x0f, temp, index) * qtab[j]
-          }
+        if (nRestart >= Decoder.RESTART_MARKER_BEGIN && nRestart <= Decoder.RESTART_MARKER_END) {
+          return nRestart
         }
+
+        prev[ctrC] += n
       }
     }
 
@@ -438,7 +476,7 @@ export class Decoder {
 
     if (this.restarting) {
       this.restarting = false
-      prev[0] = 1 << (this.frame.precision - 1)
+      prev[0] = this.initialPrediction()
     } else {
       prev[0] = this.selector()
     }
@@ -552,7 +590,7 @@ export class Decoder {
     return code & 0xff
   }
 
-  getn(PRED: number[], n: number, temp: number[], index: number[]) {
+  getn(PRED: number[], n: number, temp: number[], index: number[], compOffset = 0) {
     let result, input
     const one = 1
     const n_one = -1
@@ -565,7 +603,8 @@ export class Decoder {
     }
 
     if (n === 16) {
-      if (PRED[0] >= 0) {
+      // Each component of an interleaved scan predicts from its own value
+      if (PRED[compOffset] >= 0) {
         return -32768
       } else {
         return 32768
@@ -645,7 +684,7 @@ export class Decoder {
     } else if (this.yLoc > 0) {
       return this.getPreviousY(compOffset)
     } else {
-      return 1 << (this.frame.precision - 1)
+      return this.initialPrediction()
     }
   }
 
@@ -694,9 +733,11 @@ export class Decoder {
     const offset = this.yLoc * this.xDim + this.xLoc
 
     if (this.xLoc < this.xDim && this.yLoc < this.yDim) {
-      this.setter(offset, PRED[0], 0)
-      this.setter(offset, PRED[1], 1)
-      this.setter(offset, PRED[2], 2)
+      // Masked like the single component path: lossless prediction is modulo
+      // 2^P, so a sample that wrapped has to be brought back into range
+      this.setter(offset, this.mask & PRED[0], 0)
+      this.setter(offset, this.mask & PRED[1], 1)
+      this.setter(offset, this.mask & PRED[2], 2)
 
       this.xLoc += 1
 

tests/driver-point-transform.test.ts

@@ -0,0 +1,68 @@
+import fs from 'fs'
+import { describe, it, assert } from 'vitest'
+import { Utils, Decoder } from '../src/main.js'
+import { toArrayBuffer } from './utils.js'
+
+// A 16 bit lossless scan with a point transform of 2, from cornerstone3D#757.
+// Every other fixture in this suite has Al = 0, which is why ignoring the point
+// transform went unnoticed: this frame decoded 24576 too high and with a quarter
+// of its dynamic range, so it rendered as a flat white image.
+const jpegDataOffset = 2916
+const jpegDataSize = 7556340
+const buf = fs.readFileSync('./tests/data/jpeg_lossless_sel1-pt2.dcm')
+const data = toArrayBuffer(buf)
+const decoder = new Decoder()
+const output = decoder.decompress(data, jpegDataOffset, jpegDataSize)
+const samples = new Uint16Array(output)
+
+let smallest = Infinity
+let largest = -Infinity
+for (let i = 0; i < samples.length; i += 1) {
+  if (samples[i] < smallest) smallest = samples[i]
+  if (samples[i] > largest) largest = samples[i]
+}
+
+describe('driver-point-transform', function () {
+  it('dimX should equal 3021', function () {
+    assert.equal(3021, decoder.frame.dimX)
+  })
+
+  it('dimY should equal 2467', function () {
+    assert.equal(2467, decoder.frame.dimY)
+  })
+
+  it('number of components should be 1', function () {
+    assert.equal(1, decoder.frame.numComp)
+  })
+
+  it('precision should be 16 and the point transform 2', function () {
+    assert.equal(16, decoder.frame.precision)
+    assert.equal(2, decoder.scan.al)
+  })
+
+  it('decompressed size should be 14905614', function () {
+    assert.equal(14905614, output.byteLength)
+  })
+
+  // The three assertions below are ground truth from outside this library:
+  // the DICOM header of the file, and independent decodes by Weasis 4.7.2 and
+  // by imagecodecs' jpegsof3 (whose raw output matches these values once
+  // shifted left by Al).
+  it('smallest sample should match the file SmallestImagePixelValue of 8192', function () {
+    assert.equal(8192, smallest)
+  })
+
+  it('largest sample should stay within the file LargestImagePixelValue of 24575', function () {
+    assert.equal(24148, largest)
+    assert.isAtMost(largest, 24575)
+  })
+
+  it('sample at (365, 324) should equal 19768, as Weasis reports', function () {
+    assert.equal(19768, samples[324 * 3021 + 365])
+  })
+
+  it('data checksum should equal 2003526094', function () {
+    const checksum = Utils.crc32(output)
+    assert.equal(checksum, 2003526094)
+  })
+})

tests/interleaved-components.test.ts

@@ -0,0 +1,109 @@
+import { describe, it, assert } from 'vitest'
+import { Decoder } from '../src/main.js'
+
+/**
+ * Builds a minimal multi component lossless JPEG, so interleaved scans can be
+ * exercised without shipping a binary fixture.
+ *
+ * The Huffman table holds two DC codes: symbol 0 as '0' and symbol 1 as '10'.
+ * `withAcTable` additionally declares an AC table holding no codes, which is
+ * what the file in cornerstone3D#1116 carries and what used to break decoding.
+ */
+const buildLosslessJPEG = ({
+  precision,
+  width,
+  height,
+  components,
+  entropy,
+  withAcTable = false,
+}: {
+  precision: number
+  width: number
+  height: number
+  components: number
+  entropy: number[]
+  withAcTable?: boolean
+}) => {
+  const bytes: number[] = []
+  const push16 = (v: number) => bytes.push((v >> 8) & 0xff, v & 0xff)
+  const bits = new Array(16).fill(0)
+  bits[0] = 1 // one code of length 1 -> symbol 0
+  bits[1] = 1 // one code of length 2 -> symbol 1
+
+  bytes.push(0xff, 0xd8) // SOI
+
+  bytes.push(0xff, 0xc4) // DHT
+  push16(2 + (1 + 16 + 2) + (withAcTable ? 1 + 16 : 0))
+  bytes.push(0x00, ...bits, 0x00, 0x01) // Tc = 0 (DC), Th = 0, HUFFVAL 0 and 1
+  if (withAcTable) {
+    bytes.push(0x10, ...new Array(16).fill(0)) // Tc = 1 (AC), Th = 0, no codes
+  }
+
+  bytes.push(0xff, 0xc3) // SOF3, lossless huffman
+  push16(8 + 3 * components)
+  bytes.push(precision)
+  push16(height)
+  push16(width)
+  bytes.push(components)
+  for (let c = 0; c < components; c++) {
+    bytes.push(c, 0x11, 0) // Ci, H/V = 1/1, Tq = 0
+  }
+
+  bytes.push(0xff, 0xda) // SOS
+  push16(6 + 2 * components)
+  bytes.push(components)
+  for (let c = 0; c < components; c++) {
+    bytes.push(c, 0x00) // Cs, Td/Ta = 0
+  }
+  bytes.push(1, 0, 0) // Ss = 1 (predictor 1), Se = 0, Ah/Al = 0
+
+  bytes.push(...entropy)
+  bytes.push(0xff, 0xd9) // EOI
+
+  return Uint8Array.from(bytes)
+}
+
+// Two pixels of three components each:
+//   pixel 0: comp 0 SSSS = 0 -> stays on the 128 seed
+//            comp 1 SSSS = 1, magnitude bit 1 -> 129
+//            comp 2 SSSS = 1, magnitude bit 0 -> 127
+//   pixel 1: every component SSSS = 0 -> each keeps its own left neighbour
+// Bits '0' '101' '100' '0' '0' '0', padded with ones.
+const entropy = [0x58, 0x3f]
+
+const decode = (withAcTable: boolean) => {
+  const jpeg = buildLosslessJPEG({
+    precision: 8,
+    width: 2,
+    height: 1,
+    components: 3,
+    entropy,
+    withAcTable,
+  })
+  const decoder = new Decoder()
+  const output = decoder.decompress(jpeg.buffer, 0, jpeg.length)
+
+  return { decoder, samples: Array.from(new Uint8Array(output)) }
+}
+
+describe('interleaved multi component lossless scans', function () {
+  it('decodes every component, each predicting from its own neighbours', function () {
+    const { decoder, samples } = decode(false)
+
+    assert.equal(3, decoder.frame.numComp)
+    assert.equal(3, decoder.scan.numComp)
+    assert.deepEqual(samples, [128, 129, 127, 128, 129, 127])
+  })
+
+  it('decodes the same when the scan declares an empty AC table', function () {
+    // cornerstone3D#1116: decodeRGB used to be a copy of the baseline DCT path,
+    // reading 63 AC coefficients per block against the AC Huffman table. That
+    // stayed harmless only while no AC table existed, so the lookup terminated
+    // at once. A file that declares one - as this ultrasound did - had its
+    // entropy stream read out of step, and components 2 and 3 came back stuck
+    // at the initial prediction of 128.
+    const { samples } = decode(true)
+
+    assert.deepEqual(samples, [128, 129, 127, 128, 129, 127])
+  })
+})