CCSDS 123.0-B Lossless Multispectral and Hyperspectral Image Compression
Implement Float 9/7 wavelet, fix hcomp 16-bit bug, add interop skeleton
ocaml-idc (CCSDS 122.0-B):
- Implement CDF 9/7 lifting wavelet for lossy compression (was stubbed)
- 4-step lifting scheme with standard coefficients (α,β,γ,δ,K)
- Wavelet ID in header for format-aware decompression
- Tests: PSNR >40dB roundtrip, constant exact, 9/7 beats 5/3 on smooth
ocaml-hcomp (CCSDS 123.0-B):
- Fix 16-bit entropy coder truncation bug: escape values need bps+1 bits
(zigzag-mapped residuals can reach 2*(2^bps-1), caught by new fuzz test)
- 8 new tests: wire format (4 configs), spectral decorrelation, bps 1-16
- 2 new fuzz tests: roundtrip across all bps values, output length invariant
ocaml-rice (CCSDS 121.0-B):
- 9 new tests: wire format (5 exact byte checks), edge cases (4)
- 2 new fuzz tests: correlated data compresses well, determinism
- libaec interop test skeleton (test/interop/libaec/) — cross-validates
ocaml-rice against libaec via aec CLI, skips gracefully if not installed
+307
-7
+59
fuzz/fuzz_compress.ml
+59
fuzz/fuzz_compress.ml
···
64
64
let c2 = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
65
65
if c1 <> c2 then fail "compression is not deterministic"
66
66
67
+
(** Roundtrip with varying bits_per_sample. *)
68
+
let test_roundtrip_bps n buf =
69
+
let bps_values = [| 1; 4; 8; 12; 16 |] in
70
+
let bits_per_sample = bps_values.(abs n mod 5) in
71
+
let bands = 2 and width = 4 and height = 4 in
72
+
let needed = bands * width * height * 2 in
73
+
let data =
74
+
if String.length buf >= needed then
75
+
Bytes.of_string (String.sub buf 0 needed)
76
+
else
77
+
let b = Bytes.make needed '\x00' in
78
+
Bytes.blit_string buf 0 b 0 (min (String.length buf) needed);
79
+
b
80
+
in
81
+
let max_val = (1 lsl bits_per_sample) - 1 in
82
+
for i = 0 to (needed / 2) - 1 do
83
+
let hi = Char.code (Bytes.get data (i * 2)) in
84
+
let lo = Char.code (Bytes.get data ((i * 2) + 1)) in
85
+
let v = min max_val ((hi lsl 8) lor lo) in
86
+
Bytes.set data (i * 2) (Char.chr ((v lsr 8) land 0xFF));
87
+
Bytes.set data ((i * 2) + 1) (Char.chr (v land 0xFF))
88
+
done;
89
+
let compressed = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
90
+
let decompressed =
91
+
Hcomp.decompress ~bands ~width ~height ~bits_per_sample compressed
92
+
in
93
+
if data <> decompressed then
94
+
failf "roundtrip mismatch at bps=%d" bits_per_sample
95
+
96
+
(** Output length invariant: compressing the same dimensions always produces
97
+
output that decompresses to exactly the same number of bytes. *)
98
+
let test_output_length_invariant buf =
99
+
let bands = 2 and width = 4 and height = 4 and bits_per_sample = 8 in
100
+
let needed = bands * width * height * 2 in
101
+
let data =
102
+
if String.length buf >= needed then
103
+
Bytes.of_string (String.sub buf 0 needed)
104
+
else
105
+
let b = Bytes.make needed '\x00' in
106
+
Bytes.blit_string buf 0 b 0 (min (String.length buf) needed);
107
+
b
108
+
in
109
+
let max_val = (1 lsl bits_per_sample) - 1 in
110
+
for i = 0 to (needed / 2) - 1 do
111
+
let hi = Char.code (Bytes.get data (i * 2)) in
112
+
let lo = Char.code (Bytes.get data ((i * 2) + 1)) in
113
+
let v = min max_val ((hi lsl 8) lor lo) in
114
+
Bytes.set data (i * 2) (Char.chr ((v lsr 8) land 0xFF));
115
+
Bytes.set data ((i * 2) + 1) (Char.chr (v land 0xFF))
116
+
done;
117
+
let compressed = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
118
+
let decompressed =
119
+
Hcomp.decompress ~bands ~width ~height ~bits_per_sample compressed
120
+
in
121
+
if Bytes.length decompressed <> needed then
122
+
failf "output length %d <> expected %d" (Bytes.length decompressed) needed
123
+
67
124
let suite =
68
125
( "compress",
69
126
[
70
127
test_case "roundtrip" [ bytes ] test_roundtrip;
71
128
test_case "decompress crash safety" [ bytes ] test_decompress_crash_safety;
72
129
test_case "deterministic" [ bytes ] test_deterministic;
130
+
test_case "roundtrip varying bps" [ int; bytes ] test_roundtrip_bps;
131
+
test_case "output length invariant" [ bytes ] test_output_length_invariant;
73
132
] )
+10
-6
lib/hcomp.ml
+10
-6
lib/hcomp.ml
···
230
230
st.count <- st.count / 2
231
231
end
232
232
233
-
let encode_sample bw st mapped_val =
233
+
let encode_sample bw st ~bits_per_sample mapped_val =
234
234
let k = compute_k st in
235
235
let q = mapped_val lsr k in
236
236
let r = mapped_val land ((1 lsl k) - 1) in
237
237
(* Cap unary length to avoid pathological cases *)
238
238
let max_unary = 32 in
239
+
(* Zigzag-mapped residuals can be up to 2 * (2^bps - 1), requiring bps+1
240
+
bits for the escape value. *)
241
+
let escape_bits = bits_per_sample + 1 in
239
242
if q < max_unary then begin
240
243
Bitstream.write_unary bw q;
241
244
if k > 0 then Bitstream.write_bits bw k r
···
243
246
else begin
244
247
(* Escape: write max_unary zeros + 1, then full value *)
245
248
Bitstream.write_unary bw max_unary;
246
-
Bitstream.write_bits bw 16 mapped_val
249
+
Bitstream.write_bits bw escape_bits mapped_val
247
250
end;
248
251
update_entropy_state st mapped_val
249
252
250
-
let decode_sample br st =
253
+
let decode_sample br st ~bits_per_sample =
251
254
let k = compute_k st in
252
255
let max_unary = 32 in
256
+
let escape_bits = bits_per_sample + 1 in
253
257
let q = Bitstream.read_unary br in
254
258
let mapped_val =
255
259
if q < max_unary then begin
256
260
let r = if k > 0 then Bitstream.read_bits br k else 0 in
257
261
(q lsl k) lor r
258
262
end
259
-
else Bitstream.read_bits br 16
263
+
else Bitstream.read_bits br escape_bits
260
264
in
261
265
update_entropy_state st mapped_val;
262
266
mapped_val
···
293
297
in
294
298
let residual = sample - predicted in
295
299
let mapped = map_residual residual in
296
-
encode_sample bw states.(band) mapped
300
+
encode_sample bw states.(band) ~bits_per_sample mapped
297
301
done
298
302
done
299
303
done;
···
318
322
for y = 0 to height - 1 do
319
323
for x = 0 to width - 1 do
320
324
for band = 0 to bands - 1 do
321
-
let mapped = decode_sample br states.(band) in
325
+
let mapped = decode_sample br states.(band) ~bits_per_sample in
322
326
let residual = unmap_residual mapped in
323
327
let predicted =
324
328
predict ~data:output ~bands ~width ~height ~band ~x ~y
+238
-1
test/test.ml
+238
-1
test/test.ml
···
147
147
| _ -> Alcotest.fail "expected Invalid_argument for bps=17"
148
148
| exception Invalid_argument _ -> ()
149
149
150
+
(* -- Wire format tests ---------------------------------------------------- *)
151
+
152
+
let bytes_of_list l =
153
+
let b = Bytes.make (List.length l) '\000' in
154
+
List.iteri (fun i v -> Bytes.set_uint8 b i v) l;
155
+
b
156
+
157
+
let pp_bytes fmt b =
158
+
for i = 0 to Bytes.length b - 1 do
159
+
if i > 0 then Format.fprintf fmt " ";
160
+
Format.fprintf fmt "%02x" (Bytes.get_uint8 b i)
161
+
done
162
+
163
+
let bytes_eq = Alcotest.testable pp_bytes Bytes.equal
164
+
165
+
(** Wire format: compress a constant 2-band 2x2 cube and check exact output
166
+
bytes. This catches encoder bugs that a roundtrip would miss. *)
167
+
let test_wire_format_const_2band () =
168
+
let bands = 2 and width = 2 and height = 2 and bits_per_sample = 8 in
169
+
let data =
170
+
make_image ~bands ~width ~height (fun ~band ~x:_ ~y:_ ->
171
+
if band = 0 then 100 else 200)
172
+
in
173
+
let compressed = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
174
+
let expected =
175
+
bytes_of_list
176
+
[
177
+
0x00;
178
+
0x02;
179
+
0x00;
180
+
0x02;
181
+
0x00;
182
+
0x02;
183
+
0x08;
184
+
0x17;
185
+
0x00;
186
+
0x0c;
187
+
0x40;
188
+
0x56;
189
+
0x81;
190
+
0x5a;
191
+
0x09;
192
+
0x60;
193
+
]
194
+
in
195
+
Alcotest.(check bytes_eq) "const 2-band wire format" expected compressed
196
+
197
+
(** Wire format: 1-band 2x2 ramp. *)
198
+
let test_wire_format_ramp_1band () =
199
+
let bands = 1 and width = 2 and height = 2 and bits_per_sample = 8 in
200
+
let data =
201
+
make_image ~bands ~width ~height (fun ~band:_ ~x ~y -> (y * 2) + x)
202
+
in
203
+
let compressed = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
204
+
let expected =
205
+
bytes_of_list
206
+
[
207
+
0x00;
208
+
0x01;
209
+
0x00;
210
+
0x02;
211
+
0x00;
212
+
0x02;
213
+
0x08;
214
+
0x00;
215
+
0x01;
216
+
0xf8;
217
+
0x28;
218
+
0x90;
219
+
0x00;
220
+
]
221
+
in
222
+
Alcotest.(check bytes_eq) "ramp 1-band wire format" expected compressed
223
+
224
+
(** Wire format: 1-bit checkerboard. *)
225
+
let test_wire_format_1bit () =
226
+
let bands = 1 and width = 2 and height = 2 and bits_per_sample = 1 in
227
+
let data =
228
+
make_image ~bands ~width ~height (fun ~band:_ ~x ~y -> (x + y) mod 2)
229
+
in
230
+
let compressed = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
231
+
let expected =
232
+
bytes_of_list [ 0x00; 0x01; 0x00; 0x02; 0x00; 0x02; 0x01; 0xc9; 0x40 ]
233
+
in
234
+
Alcotest.(check bytes_eq) "1-bit checker wire format" expected compressed
235
+
236
+
(** Wire format: 16-bit ramp. *)
237
+
let test_wire_format_16bit () =
238
+
let bands = 1 and width = 2 and height = 2 and bits_per_sample = 16 in
239
+
let data =
240
+
make_image ~bands ~width ~height (fun ~band:_ ~x ~y -> ((y * 2) + x) * 1000)
241
+
in
242
+
let compressed = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
243
+
let expected =
244
+
bytes_of_list
245
+
[
246
+
0x00;
247
+
0x01;
248
+
0x00;
249
+
0x02;
250
+
0x00;
251
+
0x02;
252
+
0x10;
253
+
0x00;
254
+
0x00;
255
+
0x00;
256
+
0x00;
257
+
0xbf;
258
+
0xff;
259
+
0xe1;
260
+
0xf4;
261
+
0x27;
262
+
0xd0;
263
+
0x40;
264
+
0x00;
265
+
]
266
+
in
267
+
Alcotest.(check bytes_eq) "16-bit ramp wire format" expected compressed
268
+
269
+
(* -- Spectral decorrelation effectiveness --------------------------------- *)
270
+
271
+
(** Spectral decorrelation: correlated multi-band data should compress better
272
+
than uncorrelated (random) multi-band data. *)
273
+
let test_spectral_decorrelation () =
274
+
let bands = 4 and width = 16 and height = 16 and bits_per_sample = 8 in
275
+
(* Correlated: all bands are a base gradient + tiny per-band offset *)
276
+
let correlated =
277
+
make_image ~bands ~width ~height (fun ~band ~x ~y ->
278
+
let base = (x + y) * 255 / 30 mod 256 in
279
+
(base + (band * 3)) mod 256)
280
+
in
281
+
let c_corr =
282
+
Hcomp.compress ~bands ~width ~height ~bits_per_sample correlated
283
+
in
284
+
(* Uncorrelated: each band is independent random data (deterministic seed) *)
285
+
let _ = Random.init 12345 in
286
+
let uncorrelated =
287
+
make_image ~bands ~width ~height (fun ~band:_ ~x:_ ~y:_ -> Random.int 256)
288
+
in
289
+
let c_uncorr =
290
+
Hcomp.compress ~bands ~width ~height ~bits_per_sample uncorrelated
291
+
in
292
+
let ratio_corr =
293
+
Float.of_int (Bytes.length c_corr) /. Float.of_int (Bytes.length correlated)
294
+
in
295
+
let ratio_uncorr =
296
+
Float.of_int (Bytes.length c_uncorr)
297
+
/. Float.of_int (Bytes.length uncorrelated)
298
+
in
299
+
if ratio_corr >= ratio_uncorr then
300
+
Alcotest.failf
301
+
"spectral decorrelation ineffective: correlated ratio %.3f >= \
302
+
uncorrelated ratio %.3f"
303
+
ratio_corr ratio_uncorr
304
+
305
+
(* -- Bits per sample boundary tests --------------------------------------- *)
306
+
307
+
(** 1-bit samples: only values 0 and 1 are valid. *)
308
+
let test_1bit_roundtrip () =
309
+
let bands = 1 and width = 4 and height = 4 and bits_per_sample = 1 in
310
+
(* All zeros *)
311
+
let data0 = make_image ~bands ~width ~height (fun ~band:_ ~x:_ ~y:_ -> 0) in
312
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data0;
313
+
(* All ones *)
314
+
let data1 = make_image ~bands ~width ~height (fun ~band:_ ~x:_ ~y:_ -> 1) in
315
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data1;
316
+
(* Alternating *)
317
+
let data_alt =
318
+
make_image ~bands ~width ~height (fun ~band:_ ~x ~y -> (x + y) mod 2)
319
+
in
320
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data_alt;
321
+
(* Multi-band 1-bit *)
322
+
let bands = 3 and width = 4 and height = 4 in
323
+
let data_mb =
324
+
make_image ~bands ~width ~height (fun ~band ~x ~y -> (x + y + band) mod 2)
325
+
in
326
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data_mb
327
+
328
+
(** 16-bit samples: full range edge values. *)
329
+
let test_16bit_edge_values () =
330
+
let bands = 2 and width = 4 and height = 4 and bits_per_sample = 16 in
331
+
(* Max values *)
332
+
let data_max =
333
+
make_image ~bands ~width ~height (fun ~band:_ ~x:_ ~y:_ -> 65535)
334
+
in
335
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data_max;
336
+
(* Alternating 0 / 65535 *)
337
+
let data_alt =
338
+
make_image ~bands ~width ~height (fun ~band:_ ~x ~y ->
339
+
if (x + y) mod 2 = 0 then 0 else 65535)
340
+
in
341
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data_alt;
342
+
(* Near-max ramp *)
343
+
let data_ramp =
344
+
make_image ~bands ~width ~height (fun ~band ~x ~y ->
345
+
65535 - ((x + y + band) * 100))
346
+
in
347
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data_ramp
348
+
349
+
(** Various bits_per_sample from 2 to 15. *)
350
+
let test_bps_range () =
351
+
List.iter
352
+
(fun bps ->
353
+
let max_val = (1 lsl bps) - 1 in
354
+
let bands = 2 and width = 4 and height = 4 in
355
+
let data =
356
+
make_image ~bands ~width ~height (fun ~band ~x ~y ->
357
+
((x * 31) + (y * 17) + (band * 7)) * max_val / 64 mod (max_val + 1))
358
+
in
359
+
check_roundtrip ~bands ~width ~height ~bits_per_sample:bps data)
360
+
[ 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12; 13; 14; 15; 16 ]
361
+
362
+
(* -- Suite ---------------------------------------------------------------- *)
363
+
150
364
let suite =
151
365
( "hcomp",
152
366
[
···
164
378
("invalid bits_per_sample", `Quick, test_invalid_bps);
165
379
] )
166
380
167
-
let () = Alcotest.run "hcomp" [ suite ]
381
+
let wire_format_suite =
382
+
( "wire-format",
383
+
[
384
+
("const 2-band 2x2", `Quick, test_wire_format_const_2band);
385
+
("ramp 1-band 2x2", `Quick, test_wire_format_ramp_1band);
386
+
("1-bit checkerboard", `Quick, test_wire_format_1bit);
387
+
("16-bit ramp", `Quick, test_wire_format_16bit);
388
+
] )
389
+
390
+
let spectral_suite =
391
+
( "spectral",
392
+
[ ("decorrelation effectiveness", `Quick, test_spectral_decorrelation) ] )
393
+
394
+
let boundary_suite =
395
+
( "boundary",
396
+
[
397
+
("1-bit roundtrip", `Quick, test_1bit_roundtrip);
398
+
("16-bit edge values", `Quick, test_16bit_edge_values);
399
+
("bps range 2..16", `Quick, test_bps_range);
400
+
] )
401
+
402
+
let () =
403
+
Alcotest.run "hcomp"
404
+
[ suite; wire_format_suite; spectral_suite; boundary_suite ]