diff --git a/src/decoder.ts b/src/decoder.ts index 798e8fb..69d1f2e 100644 --- a/src/decoder.ts +++ b/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 diff --git a/tests/data/jpeg_lossless_sel1-pt2.dcm b/tests/data/jpeg_lossless_sel1-pt2.dcm new file mode 100644 index 0000000..300c648 Binary files /dev/null and b/tests/data/jpeg_lossless_sel1-pt2.dcm differ diff --git a/tests/driver-point-transform.test.ts b/tests/driver-point-transform.test.ts new file mode 100644 index 0000000..e60ee06 --- /dev/null +++ b/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) + }) +}) diff --git a/tests/interleaved-components.test.ts b/tests/interleaved-components.test.ts new file mode 100644 index 0000000..a0b37ff --- /dev/null +++ b/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]) + }) +}) diff --git a/tests/point-transform.test.ts b/tests/point-transform.test.ts new file mode 100644 index 0000000..d919fc5 --- /dev/null +++ b/tests/point-transform.test.ts @@ -0,0 +1,115 @@ +import { describe, it, assert } from 'vitest' +import { Decoder } from '../src/main.js' + +/** + * Builds a minimal single component lossless JPEG, so a point transform can be + * exercised without shipping a binary fixture. + * + * The Huffman table holds two codes: symbol 0 as '0' and symbol 1 as '10', + * which is enough to encode the SSSS categories these cases need. + */ +const buildLosslessJPEG = ({ + precision, + width, + height, + pointTransform, + entropy, +}: { + precision: number + width: number + height: number + pointTransform: number + entropy: number[] +}) => { + const bytes: number[] = [] + const push16 = (v: number) => bytes.push((v >> 8) & 0xff, v & 0xff) + + bytes.push(0xff, 0xd8) // SOI + + bytes.push(0xff, 0xc4) // DHT + push16(2 + 1 + 16 + 2) + bytes.push(0x00) // Tc = 0 (DC), Th = 0 + 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(...bits) + bytes.push(0x00, 0x01) // HUFFVAL + + bytes.push(0xff, 0xc3) // SOF3, lossless huffman + push16(8 + 3) + bytes.push(precision) + push16(height) + push16(width) + bytes.push(1) // Nf + bytes.push(1, 0x11, 0) // C = 1, H/V = 1/1, Tq = 0 + + bytes.push(0xff, 0xda) // SOS + push16(6 + 2) + bytes.push(1) // Ns + bytes.push(1, 0x00) // Cs = 1, Td/Ta = 0 + bytes.push(1, 0, pointTransform & 0x0f) // Ss = 1 (predictor 1), Se = 0, Ah/Al + + bytes.push(...entropy) + bytes.push(0xff, 0xd9) // EOI + + return Uint8Array.from(bytes) +} + +// Two samples: the first codes SSSS = 0 ('0', difference 0) so it lands on the +// initial prediction, the second codes SSSS = 1 ('10' plus magnitude bit '1', +// difference +1). Bits '0' '10' '1' padded with ones is 0x5f. +const entropy = [0x5f] + +const decode = (precision: number, pointTransform: number) => { + const jpeg = buildLosslessJPEG({ + precision, + width: 2, + height: 1, + pointTransform, + entropy, + }) + const decoder = new Decoder() + const output = decoder.decompress(jpeg.buffer, 0, jpeg.length) + + return { + decoder, + samples: Array.from(precision > 8 ? new Uint16Array(output) : new Uint8Array(output)), + } +} + +describe('point transform (T.81 A.4.1)', function () { + it('parses Al out of the scan header', function () { + const { decoder } = decode(8, 1) + assert.equal(decoder.scan.al, 1) + assert.equal(decoder.getPointTransform(), 1) + }) + + it('seeds the prediction with 2^(P - Pt - 1) and scales the output back up', function () { + // Seed is 2^(8 - 1 - 1) = 64, second sample 65, both shifted left by 1. + // Ignoring Al gave [128, 129]: right for a flat frame, but the differences + // keep only half their amplitude, so the frame loses dynamic range. + const { samples } = decode(8, 1) + assert.deepEqual(samples, [128, 130]) + }) + + it('handles a larger point transform', function () { + // Seed 2^(8 - 2 - 1) = 32, second sample 33, both shifted left by 2 + const { samples } = decode(8, 2) + assert.deepEqual(samples, [128, 132]) + }) + + it('is a no-op when there is no point transform', function () { + // Seed 2^(8 - 1) = 128, no shift + const { samples } = decode(8, 0) + assert.deepEqual(samples, [128, 129]) + }) + + it('applies to 16 bit frames', function () { + // Seed 2^(16 - 2 - 1) = 8192, second sample 8193, both shifted left by 2. + // This is the shape of the DICOM frame in cornerstone3D#757: a 16 bit + // lossless scan with Al = 2, which decoded 24576 high with a quarter of its + // dynamic range. + const { samples } = decode(16, 2) + assert.deepEqual(samples, [32768, 32772]) + }) +})