CCSDS 123.0-B Lossless Multispectral and Hyperspectral Image Compression
Rename ccsds-122 → idc, ccsds-123 → hcomp
- ocaml-idc: Image Data Compression (CCSDS 122.0-B)
- ocaml-hcomp: Hyperspectral Compression (CCSDS 123.0-B)
These are the names used by the space community. Updated index.mld.
+733
+42
README.md
+42
README.md
···
1
+
# hcomp
2
+
3
+
CCSDS 123.0-B Lossless Multispectral and Hyperspectral Image Compression in
4
+
pure OCaml.
5
+
6
+
## Overview
7
+
8
+
Predictor-based lossless compression for multi-band (hyperspectral) images
9
+
following the CCSDS 123.0-B standard, used for space science data systems.
10
+
The algorithm uses a local difference predictor with spatial and spectral
11
+
neighbors followed by an adaptive entropy coder.
12
+
13
+
## Features
14
+
15
+
- Local difference predictor using spatial and spectral neighbors
16
+
- Adaptive entropy coder (Rice-like variable-length coding)
17
+
- Band-interleaved-by-pixel (BIP) sample ordering
18
+
- Support for 1--16 bits per sample
19
+
20
+
## Installation
21
+
22
+
```
23
+
opam install hcomp
24
+
```
25
+
26
+
## Usage
27
+
28
+
```ocaml
29
+
(* Compress multi-band image data *)
30
+
let compressed =
31
+
Hcomp.compress ~bands:4 ~width:256 ~height:256
32
+
~bits_per_sample:12 image_data
33
+
34
+
(* Decompress *)
35
+
let restored =
36
+
Hcomp.decompress ~bands:4 ~width:256 ~height:256
37
+
~bits_per_sample:12 compressed
38
+
```
39
+
40
+
## Licence
41
+
42
+
ISC License. See [LICENSE.md](LICENSE.md) for details.
+20
dune-project
+20
dune-project
···
1
+
(lang dune 3.21)
2
+
(name hcomp)
3
+
4
+
(generate_opam_files true)
5
+
(implicit_transitive_deps false)
6
+
7
+
(source (tangled gazagnaire.org/ocaml-hcomp))
8
+
9
+
(maintainers "Thomas Gazagnaire")
10
+
(authors "Thomas Gazagnaire")
11
+
12
+
(package
13
+
(name hcomp)
14
+
(synopsis "CCSDS 123.0-B Lossless Multispectral and Hyperspectral Image Compression")
15
+
(description
16
+
"Predictor-based lossless compression for multi-band (hyperspectral) images following the CCSDS 123.0-B standard. Uses a local difference predictor with spatial and spectral neighbors and an adaptive entropy coder.")
17
+
(depends
18
+
(ocaml (>= 5.1))
19
+
(alcotest (and (>= 1.0) :with-test))
20
+
(crowbar (and (>= 0.2) :with-test))))
+22
fuzz/dune
+22
fuzz/dune
···
1
+
(executable
2
+
(name fuzz)
3
+
(modules fuzz fuzz_compress)
4
+
(libraries hcomp alcobar))
5
+
6
+
(rule
7
+
(alias runtest)
8
+
(enabled_if
9
+
(<> %{profile} afl))
10
+
(deps fuzz.exe)
11
+
(action
12
+
(run %{exe:fuzz.exe})))
13
+
14
+
(rule
15
+
(alias fuzz)
16
+
(enabled_if
17
+
(= %{profile} afl))
18
+
(deps fuzz.exe)
19
+
(action
20
+
(progn
21
+
(run %{exe:fuzz.exe} --gen-corpus corpus)
22
+
(run afl-fuzz -V 60 -i corpus -o _fuzz -- %{exe:fuzz.exe} @@))))
+1
fuzz/fuzz.ml
+1
fuzz/fuzz.ml
···
1
+
let () = Alcobar.run "hcomp" [ Fuzz_compress.suite ]
+73
fuzz/fuzz_compress.ml
+73
fuzz/fuzz_compress.ml
···
1
+
(** Fuzz tests for CCSDS 123 compression. *)
2
+
3
+
open Alcobar
4
+
5
+
(** Roundtrip - compress then decompress must recover original data. *)
6
+
let test_roundtrip buf =
7
+
(* Fixed small image: 2 bands, 4x4 *)
8
+
let bands = 2 and width = 4 and height = 4 and bits_per_sample = 8 in
9
+
let needed = bands * width * height * 2 in
10
+
let data =
11
+
if String.length buf >= needed then
12
+
Bytes.of_string (String.sub buf 0 needed)
13
+
else
14
+
let b = Bytes.make needed '\x00' in
15
+
Bytes.blit_string buf 0 b 0 (String.length buf);
16
+
b
17
+
in
18
+
(* Clamp values to valid range *)
19
+
let max_val = (1 lsl bits_per_sample) - 1 in
20
+
for i = 0 to (needed / 2) - 1 do
21
+
let hi = Char.code (Bytes.get data (i * 2)) in
22
+
let lo = Char.code (Bytes.get data ((i * 2) + 1)) in
23
+
let v = min max_val ((hi lsl 8) lor lo) in
24
+
Bytes.set data (i * 2) (Char.chr ((v lsr 8) land 0xFF));
25
+
Bytes.set data ((i * 2) + 1) (Char.chr (v land 0xFF))
26
+
done;
27
+
let compressed = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
28
+
let decompressed =
29
+
Hcomp.decompress ~bands ~width ~height ~bits_per_sample compressed
30
+
in
31
+
if data <> decompressed then fail "roundtrip mismatch"
32
+
33
+
(** Decompress must not crash on arbitrary input. *)
34
+
let test_decompress_crash_safety buf =
35
+
let data = Bytes.of_string buf in
36
+
try
37
+
let _ =
38
+
Hcomp.decompress ~bands:1 ~width:4 ~height:4 ~bits_per_sample:8 data
39
+
in
40
+
()
41
+
with _ -> ()
42
+
43
+
(** Compressed output must be deterministic. *)
44
+
let test_deterministic buf =
45
+
let bands = 2 and width = 4 and height = 4 and bits_per_sample = 8 in
46
+
let needed = bands * width * height * 2 in
47
+
let data =
48
+
if String.length buf >= needed then
49
+
Bytes.of_string (String.sub buf 0 needed)
50
+
else
51
+
let b = Bytes.make needed '\x00' in
52
+
Bytes.blit_string buf 0 b 0 (String.length buf);
53
+
b
54
+
in
55
+
let max_val = (1 lsl bits_per_sample) - 1 in
56
+
for i = 0 to (needed / 2) - 1 do
57
+
let hi = Char.code (Bytes.get data (i * 2)) in
58
+
let lo = Char.code (Bytes.get data ((i * 2) + 1)) in
59
+
let v = min max_val ((hi lsl 8) lor lo) in
60
+
Bytes.set data (i * 2) (Char.chr ((v lsr 8) land 0xFF));
61
+
Bytes.set data ((i * 2) + 1) (Char.chr (v land 0xFF))
62
+
done;
63
+
let c1 = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
64
+
let c2 = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
65
+
if c1 <> c2 then fail "compression is not deterministic"
66
+
67
+
let suite =
68
+
( "compress",
69
+
[
70
+
test_case "roundtrip" [ bytes ] test_roundtrip;
71
+
test_case "decompress crash safety" [ bytes ] test_decompress_crash_safety;
72
+
test_case "deterministic" [ bytes ] test_deterministic;
73
+
] )
+4
fuzz/fuzz_compress.mli
+4
fuzz/fuzz_compress.mli
+33
hcomp.opam
+33
hcomp.opam
···
1
+
# This file is generated by dune, edit dune-project instead
2
+
opam-version: "2.0"
3
+
synopsis:
4
+
"CCSDS 123.0-B Lossless Multispectral and Hyperspectral Image Compression"
5
+
description:
6
+
"Predictor-based lossless compression for multi-band (hyperspectral) images following the CCSDS 123.0-B standard. Uses a local difference predictor with spatial and spectral neighbors and an adaptive entropy coder."
7
+
maintainer: ["Thomas Gazagnaire"]
8
+
authors: ["Thomas Gazagnaire"]
9
+
homepage: "https://tangled.org/gazagnaire.org/ocaml-hcomp"
10
+
bug-reports: "https://tangled.org/gazagnaire.org/ocaml-hcomp/issues"
11
+
depends: [
12
+
"dune" {>= "3.21"}
13
+
"ocaml" {>= "5.1"}
14
+
"alcotest" {>= "1.0" & with-test}
15
+
"crowbar" {>= "0.2" & with-test}
16
+
"odoc" {with-doc}
17
+
]
18
+
build: [
19
+
["dune" "subst"] {dev}
20
+
[
21
+
"dune"
22
+
"build"
23
+
"-p"
24
+
name
25
+
"-j"
26
+
jobs
27
+
"@install"
28
+
"@runtest" {with-test}
29
+
"@doc" {with-doc}
30
+
]
31
+
]
32
+
dev-repo: "git+https://tangled.org/gazagnaire.org/ocaml-hcomp"
33
+
x-maintenance-intent: ["(latest)"]
+332
lib/hcomp.ml
+332
lib/hcomp.ml
···
1
+
(* Copyright (c) 2026 Thomas Gazagnaire <thomas@gazagnaire.org>
2
+
3
+
Permission to use, copy, modify, and distribute this software for any
4
+
purpose with or without fee is hereby granted, provided that the above
5
+
copyright notice and this permission notice appear in all copies.
6
+
7
+
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
8
+
WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
9
+
MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
10
+
ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
11
+
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
12
+
ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
13
+
OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. *)
14
+
15
+
(** CCSDS 123.0-B Lossless Multispectral and Hyperspectral Image Compression.
16
+
17
+
Implements: 1. Local difference predictor using spatial and spectral
18
+
neighbors 2. Sample-adaptive entropy coder (Rice-like variable-length
19
+
coding) *)
20
+
21
+
(* ---- Bitstream I/O ---- *)
22
+
23
+
module Bitstream : sig
24
+
type writer
25
+
26
+
val create : int -> writer
27
+
val write_bits : writer -> int -> int -> unit
28
+
val write_unary : writer -> int -> unit
29
+
val contents : writer -> bytes
30
+
31
+
type reader
32
+
33
+
val of_bytes : bytes -> reader
34
+
val read_bits : reader -> int -> int
35
+
val read_unary : reader -> int
36
+
end = struct
37
+
type writer = { mutable buf : bytes; mutable pos : int }
38
+
39
+
let create cap = { buf = Bytes.make (max 16 cap) '\x00'; pos = 0 }
40
+
41
+
let ensure_capacity w nbits =
42
+
let needed = (w.pos + nbits + 7) / 8 in
43
+
if needed > Bytes.length w.buf then begin
44
+
let new_cap = max needed (Bytes.length w.buf * 2) in
45
+
let new_buf = Bytes.make new_cap '\x00' in
46
+
Bytes.blit w.buf 0 new_buf 0 (Bytes.length w.buf);
47
+
w.buf <- new_buf
48
+
end
49
+
50
+
let write_bits w nbits value =
51
+
ensure_capacity w nbits;
52
+
for i = nbits - 1 downto 0 do
53
+
let bit = (value lsr i) land 1 in
54
+
let byte_idx = w.pos / 8 in
55
+
let bit_idx = 7 - (w.pos mod 8) in
56
+
let cur = Char.code (Bytes.get w.buf byte_idx) in
57
+
Bytes.set w.buf byte_idx (Char.chr (cur lor (bit lsl bit_idx)));
58
+
w.pos <- w.pos + 1
59
+
done
60
+
61
+
let write_unary w n =
62
+
ensure_capacity w (n + 1);
63
+
for _ = 1 to n do
64
+
(* Write n zeros *)
65
+
w.pos <- w.pos + 1
66
+
done;
67
+
(* Write terminating 1 *)
68
+
let byte_idx = w.pos / 8 in
69
+
let bit_idx = 7 - (w.pos mod 8) in
70
+
let cur = Char.code (Bytes.get w.buf byte_idx) in
71
+
Bytes.set w.buf byte_idx (Char.chr (cur lor (1 lsl bit_idx)));
72
+
w.pos <- w.pos + 1
73
+
74
+
let contents w =
75
+
let len = (w.pos + 7) / 8 in
76
+
Bytes.sub w.buf 0 len
77
+
78
+
type reader = { data : bytes; mutable rpos : int; max_bits : int }
79
+
80
+
let of_bytes data = { data; rpos = 0; max_bits = Bytes.length data * 8 }
81
+
82
+
let read_bits r nbits =
83
+
let value = ref 0 in
84
+
for _ = 1 to nbits do
85
+
if r.rpos < r.max_bits then begin
86
+
let byte_idx = r.rpos / 8 in
87
+
let bit_idx = 7 - (r.rpos mod 8) in
88
+
let bit = (Char.code (Bytes.get r.data byte_idx) lsr bit_idx) land 1 in
89
+
value := (!value lsl 1) lor bit;
90
+
r.rpos <- r.rpos + 1
91
+
end
92
+
else value := !value lsl 1
93
+
done;
94
+
!value
95
+
96
+
let read_unary r =
97
+
let n = ref 0 in
98
+
while
99
+
r.rpos < r.max_bits
100
+
&& begin
101
+
let byte_idx = r.rpos / 8 in
102
+
let bit_idx = 7 - (r.rpos mod 8) in
103
+
let bit = (Char.code (Bytes.get r.data byte_idx) lsr bit_idx) land 1 in
104
+
r.rpos <- r.rpos + 1;
105
+
bit = 0
106
+
end
107
+
do
108
+
incr n
109
+
done;
110
+
!n
111
+
end
112
+
113
+
(* ---- Sample access helpers ---- *)
114
+
115
+
(** Image data is stored in band-interleaved-by-pixel (BIP) order with 16-bit
116
+
big-endian samples. *)
117
+
118
+
let get_sample (data : bytes) ~bands ~width ~band ~x ~y =
119
+
let idx = ((((y * width) + x) * bands) + band) * 2 in
120
+
let hi = Char.code (Bytes.get data idx) in
121
+
let lo = Char.code (Bytes.get data (idx + 1)) in
122
+
(hi lsl 8) lor lo
123
+
124
+
let set_sample (data : bytes) ~bands ~width ~band ~x ~y value =
125
+
let idx = ((((y * width) + x) * bands) + band) * 2 in
126
+
Bytes.set data idx (Char.chr ((value lsr 8) land 0xFF));
127
+
Bytes.set data (idx + 1) (Char.chr (value land 0xFF))
128
+
129
+
(* ---- Local difference predictor ---- *)
130
+
131
+
(** Predict a sample using spatial and spectral neighbors.
132
+
133
+
The predictor uses:
134
+
- North (x, y-1): spatial neighbor above
135
+
- West (x-1, y): spatial neighbor to the left
136
+
- Northwest (x-1, y-1): spatial diagonal neighbor
137
+
- Previous band (band-1) at same position: spectral neighbor
138
+
139
+
Prediction formula (simplified from CCSDS 123.0-B):
140
+
- If first band, first pixel: mid-range value
141
+
- If first band: spatial prediction (median of N, W, N+W-NW)
142
+
- Otherwise: weighted combination of spatial prediction and spectral
143
+
neighbor *)
144
+
145
+
let predict ~(data : bytes) ~bands ~width ~height:(_height : int) ~band ~x ~y
146
+
~bits_per_sample =
147
+
ignore _height;
148
+
let mid = 1 lsl (bits_per_sample - 1) in
149
+
if band = 0 then begin
150
+
(* First band: spatial-only prediction *)
151
+
if x = 0 && y = 0 then mid
152
+
else if y = 0 then get_sample data ~bands ~width ~band ~x:(x - 1) ~y
153
+
else if x = 0 then get_sample data ~bands ~width ~band ~x ~y:(y - 1)
154
+
else begin
155
+
let n = get_sample data ~bands ~width ~band ~x ~y:(y - 1) in
156
+
let w = get_sample data ~bands ~width ~band ~x:(x - 1) ~y in
157
+
let nw = get_sample data ~bands ~width ~band ~x:(x - 1) ~y:(y - 1) in
158
+
(* Median predictor (same as PNG "Paeth" without abs distance) *)
159
+
let p = n + w - nw in
160
+
let max_val = (1 lsl bits_per_sample) - 1 in
161
+
max 0 (min max_val p)
162
+
end
163
+
end
164
+
else begin
165
+
(* Subsequent bands: blend spatial prediction with spectral neighbor *)
166
+
let spectral = get_sample data ~bands ~width ~band:(band - 1) ~x ~y in
167
+
if x = 0 && y = 0 then spectral
168
+
else begin
169
+
let spatial =
170
+
if y = 0 then get_sample data ~bands ~width ~band ~x:(x - 1) ~y
171
+
else if x = 0 then get_sample data ~bands ~width ~band ~x ~y:(y - 1)
172
+
else begin
173
+
let n = get_sample data ~bands ~width ~band ~x ~y:(y - 1) in
174
+
let w = get_sample data ~bands ~width ~band ~x:(x - 1) ~y in
175
+
let nw = get_sample data ~bands ~width ~band ~x:(x - 1) ~y:(y - 1) in
176
+
let p = n + w - nw in
177
+
let max_val = (1 lsl bits_per_sample) - 1 in
178
+
max 0 (min max_val p)
179
+
end
180
+
in
181
+
(* Weight: 3/4 spectral + 1/4 spatial *)
182
+
((spectral * 3) + spatial) / 4
183
+
end
184
+
end
185
+
186
+
(* ---- Mapped residual (for entropy coding) ---- *)
187
+
188
+
(** Map a signed prediction residual to an unsigned value. Uses the standard
189
+
interleaving: 0, -1, 1, -2, 2, ... -> 0, 1, 2, 3, 4, ... This is sometimes
190
+
called "zigzag" encoding. *)
191
+
let map_residual r = if r >= 0 then 2 * r else (-2 * r) - 1
192
+
193
+
let unmap_residual m = if m mod 2 = 0 then m / 2 else -(m / 2) - 1
194
+
195
+
(* ---- Adaptive entropy coder (Rice-like) ---- *)
196
+
197
+
(** The entropy coder maintains a running estimate of the optimal Rice parameter
198
+
k per band, adapting it based on recent residual magnitudes.
199
+
200
+
Encoding of mapped residual m with parameter k:
201
+
- quotient q = m >> k (written in unary)
202
+
- remainder r = m & ((1 << k) - 1) (written in k bits)
203
+
204
+
The parameter k is adapted: accumulate |residuals| and count samples, then k
205
+
= floor(log2(accumulator / count)). *)
206
+
207
+
type entropy_state = { mutable accum : int; mutable count : int }
208
+
209
+
let create_entropy_state bits_per_sample =
210
+
(* Initialize accumulator to bias toward middle k values *)
211
+
let initial_k = max 1 (bits_per_sample / 2) in
212
+
{ accum = 1 lsl initial_k; count = 1 }
213
+
214
+
let compute_k st =
215
+
if st.count = 0 then 0
216
+
else begin
217
+
let ratio = st.accum / st.count in
218
+
if ratio <= 0 then 0
219
+
else
220
+
let rec log2 n acc = if n <= 1 then acc else log2 (n lsr 1) (acc + 1) in
221
+
log2 ratio 0
222
+
end
223
+
224
+
let update_entropy_state st mapped_val =
225
+
st.accum <- st.accum + mapped_val;
226
+
st.count <- st.count + 1;
227
+
(* Halve periodically to keep adaptation responsive *)
228
+
if st.count >= 64 then begin
229
+
st.accum <- st.accum / 2;
230
+
st.count <- st.count / 2
231
+
end
232
+
233
+
let encode_sample bw st mapped_val =
234
+
let k = compute_k st in
235
+
let q = mapped_val lsr k in
236
+
let r = mapped_val land ((1 lsl k) - 1) in
237
+
(* Cap unary length to avoid pathological cases *)
238
+
let max_unary = 32 in
239
+
if q < max_unary then begin
240
+
Bitstream.write_unary bw q;
241
+
if k > 0 then Bitstream.write_bits bw k r
242
+
end
243
+
else begin
244
+
(* Escape: write max_unary zeros + 1, then full value *)
245
+
Bitstream.write_unary bw max_unary;
246
+
Bitstream.write_bits bw 16 mapped_val
247
+
end;
248
+
update_entropy_state st mapped_val
249
+
250
+
let decode_sample br st =
251
+
let k = compute_k st in
252
+
let max_unary = 32 in
253
+
let q = Bitstream.read_unary br in
254
+
let mapped_val =
255
+
if q < max_unary then begin
256
+
let r = if k > 0 then Bitstream.read_bits br k else 0 in
257
+
(q lsl k) lor r
258
+
end
259
+
else Bitstream.read_bits br 16
260
+
in
261
+
update_entropy_state st mapped_val;
262
+
mapped_val
263
+
264
+
(* ---- Public API ---- *)
265
+
266
+
let compress ~bands ~width ~height ~bits_per_sample data =
267
+
let expected_len = bands * width * height * 2 in
268
+
if Bytes.length data <> expected_len then
269
+
invalid_arg
270
+
(Printf.sprintf "Hcomp.compress: expected %d bytes, got %d" expected_len
271
+
(Bytes.length data));
272
+
if bits_per_sample < 1 || bits_per_sample > 16 then
273
+
invalid_arg
274
+
(Printf.sprintf "Hcomp.compress: bits_per_sample must be 1..16, got %d"
275
+
bits_per_sample);
276
+
let bw = Bitstream.create (expected_len / 2) in
277
+
(* Header *)
278
+
Bitstream.write_bits bw 16 bands;
279
+
Bitstream.write_bits bw 16 width;
280
+
Bitstream.write_bits bw 16 height;
281
+
Bitstream.write_bits bw 8 bits_per_sample;
282
+
(* One entropy state per band for spectral adaptation *)
283
+
let states =
284
+
Array.init bands (fun _ -> create_entropy_state bits_per_sample)
285
+
in
286
+
(* Encode in BIP order: for each pixel, for each band *)
287
+
for y = 0 to height - 1 do
288
+
for x = 0 to width - 1 do
289
+
for band = 0 to bands - 1 do
290
+
let sample = get_sample data ~bands ~width ~band ~x ~y in
291
+
let predicted =
292
+
predict ~data ~bands ~width ~height ~band ~x ~y ~bits_per_sample
293
+
in
294
+
let residual = sample - predicted in
295
+
let mapped = map_residual residual in
296
+
encode_sample bw states.(band) mapped
297
+
done
298
+
done
299
+
done;
300
+
Bitstream.contents bw
301
+
302
+
let decompress ~bands ~width ~height ~bits_per_sample data =
303
+
let br = Bitstream.of_bytes data in
304
+
(* Read and verify header *)
305
+
let h_bands = Bitstream.read_bits br 16 in
306
+
let h_width = Bitstream.read_bits br 16 in
307
+
let h_height = Bitstream.read_bits br 16 in
308
+
let h_bps = Bitstream.read_bits br 8 in
309
+
if
310
+
h_bands <> bands || h_width <> width || h_height <> height
311
+
|| h_bps <> bits_per_sample
312
+
then invalid_arg "Hcomp.decompress: header mismatch";
313
+
let output = Bytes.make (bands * width * height * 2) '\x00' in
314
+
let states =
315
+
Array.init bands (fun _ -> create_entropy_state bits_per_sample)
316
+
in
317
+
let max_val = (1 lsl bits_per_sample) - 1 in
318
+
for y = 0 to height - 1 do
319
+
for x = 0 to width - 1 do
320
+
for band = 0 to bands - 1 do
321
+
let mapped = decode_sample br states.(band) in
322
+
let residual = unmap_residual mapped in
323
+
let predicted =
324
+
predict ~data:output ~bands ~width ~height ~band ~x ~y
325
+
~bits_per_sample
326
+
in
327
+
let sample = max 0 (min max_val (predicted + residual)) in
328
+
set_sample output ~bands ~width ~band ~x ~y sample
329
+
done
330
+
done
331
+
done;
332
+
output
+33
lib/hcomp.mli
+33
lib/hcomp.mli
···
1
+
(** CCSDS 123.0-B Lossless and Near-Lossless Multispectral and Hyperspectral
2
+
Image Compression.
3
+
4
+
Predictor-based compression for multi-band (hyperspectral) images. The
5
+
algorithm uses a local difference predictor with spatial and spectral
6
+
neighbors followed by an adaptive entropy coder (sample-adaptive, similar to
7
+
Rice coding). *)
8
+
9
+
val compress :
10
+
bands:int -> width:int -> height:int -> bits_per_sample:int -> bytes -> bytes
11
+
(** [compress ~bands ~width ~height ~bits_per_sample data] compresses multi-band
12
+
image data.
13
+
14
+
@param bands Number of spectral bands.
15
+
@param width Image width in pixels (samples per line).
16
+
@param height Image height in pixels (number of lines).
17
+
@param bits_per_sample Bits per sample (1--16).
18
+
@param data
19
+
Raw image data in band-interleaved-by-pixel (BIP) order: for each pixel
20
+
(y, x), all [bands] samples are stored consecutively as 16-bit big-endian
21
+
unsigned integers. Total length must be [bands * width * height * 2]
22
+
bytes. *)
23
+
24
+
val decompress :
25
+
bands:int -> width:int -> height:int -> bits_per_sample:int -> bytes -> bytes
26
+
(** [decompress ~bands ~width ~height ~bits_per_sample compressed] decompresses
27
+
data produced by {!compress}.
28
+
29
+
@param bands Number of spectral bands.
30
+
@param width Image width in pixels.
31
+
@param height Image height in pixels.
32
+
@param bits_per_sample Bits per sample (1--16).
33
+
@param compressed Compressed bitstream. *)
+167
test/test.ml
+167
test/test.ml
···
1
+
(* Copyright (c) 2026 Thomas Gazagnaire <thomas@gazagnaire.org>
2
+
3
+
Permission to use, copy, modify, and distribute this software for any
4
+
purpose with or without fee is hereby granted, provided that the above
5
+
copyright notice and this permission notice appear in all copies.
6
+
7
+
THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
8
+
WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
9
+
MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
10
+
ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
11
+
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
12
+
ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
13
+
OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. *)
14
+
15
+
(** Helper: create multi-band image data in BIP order with 16-bit big-endian
16
+
samples. *)
17
+
let make_image ~bands ~width ~height f =
18
+
let data = Bytes.make (bands * width * height * 2) '\x00' in
19
+
for y = 0 to height - 1 do
20
+
for x = 0 to width - 1 do
21
+
for band = 0 to bands - 1 do
22
+
let v = f ~band ~x ~y in
23
+
let idx = ((((y * width) + x) * bands) + band) * 2 in
24
+
Bytes.set data idx (Char.chr ((v lsr 8) land 0xFF));
25
+
Bytes.set data (idx + 1) (Char.chr (v land 0xFF))
26
+
done
27
+
done
28
+
done;
29
+
data
30
+
31
+
let check_roundtrip ~bands ~width ~height ~bits_per_sample data =
32
+
let compressed = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
33
+
let decompressed =
34
+
Hcomp.decompress ~bands ~width ~height ~bits_per_sample compressed
35
+
in
36
+
if data <> decompressed then
37
+
Alcotest.failf
38
+
"roundtrip failed: input and output differ (bands=%d width=%d height=%d \
39
+
bps=%d)"
40
+
bands width height bits_per_sample
41
+
42
+
(* Single band, constant image. *)
43
+
let test_single_band_constant () =
44
+
let bands = 1 and width = 8 and height = 8 and bits_per_sample = 8 in
45
+
let data = make_image ~bands ~width ~height (fun ~band:_ ~x:_ ~y:_ -> 128) in
46
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data
47
+
48
+
(* Single band, gradient. *)
49
+
let test_single_band_gradient () =
50
+
let bands = 1 and width = 8 and height = 8 and bits_per_sample = 8 in
51
+
let data = make_image ~bands ~width ~height (fun ~band:_ ~x ~y:_ -> x * 32) in
52
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data
53
+
54
+
(* Multi-band constant. *)
55
+
let test_multi_band_constant () =
56
+
let bands = 4 and width = 8 and height = 8 and bits_per_sample = 8 in
57
+
let data =
58
+
make_image ~bands ~width ~height (fun ~band ~x:_ ~y:_ -> 50 * (band + 1))
59
+
in
60
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data
61
+
62
+
(* Multi-band with spatial gradient. *)
63
+
let test_multi_band_gradient () =
64
+
let bands = 3 and width = 16 and height = 16 and bits_per_sample = 8 in
65
+
let data =
66
+
make_image ~bands ~width ~height (fun ~band ~x ~y ->
67
+
(((x + y) * 8) + (band * 20)) mod 256)
68
+
in
69
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data
70
+
71
+
(* 16-bit samples. *)
72
+
let test_16bit_samples () =
73
+
let bands = 2 and width = 8 and height = 8 and bits_per_sample = 16 in
74
+
let data =
75
+
make_image ~bands ~width ~height (fun ~band ~x ~y ->
76
+
((x * 1000) + (y * 500) + (band * 10000)) mod 65536)
77
+
in
78
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data
79
+
80
+
(* All zeros. *)
81
+
let test_all_zeros () =
82
+
let bands = 3 and width = 4 and height = 4 and bits_per_sample = 8 in
83
+
let data = make_image ~bands ~width ~height (fun ~band:_ ~x:_ ~y:_ -> 0) in
84
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data
85
+
86
+
(* Single pixel image. *)
87
+
let test_single_pixel () =
88
+
let bands = 5 and width = 1 and height = 1 and bits_per_sample = 12 in
89
+
let data =
90
+
make_image ~bands ~width ~height (fun ~band ~x:_ ~y:_ -> band * 100)
91
+
in
92
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data
93
+
94
+
(* Spectral correlation: bands are similar with small offsets. *)
95
+
let test_spectral_correlation () =
96
+
let bands = 8 and width = 8 and height = 8 and bits_per_sample = 10 in
97
+
let data =
98
+
make_image ~bands ~width ~height (fun ~band ~x ~y ->
99
+
let base = ((x * 30) + (y * 30)) mod 1024 in
100
+
(base + (band * 5)) mod 1024)
101
+
in
102
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data
103
+
104
+
(* Checkerboard pattern per band. *)
105
+
let test_checkerboard_multi_band () =
106
+
let bands = 2 and width = 8 and height = 8 and bits_per_sample = 8 in
107
+
let data =
108
+
make_image ~bands ~width ~height (fun ~band ~x ~y ->
109
+
if (x + y + band) mod 2 = 0 then 200 else 50)
110
+
in
111
+
check_roundtrip ~bands ~width ~height ~bits_per_sample data
112
+
113
+
(* Verify compression ratio for correlated multi-band data. *)
114
+
let test_compression_ratio () =
115
+
let bands = 4 and width = 16 and height = 16 and bits_per_sample = 8 in
116
+
let data =
117
+
make_image ~bands ~width ~height (fun ~band ~x ~y ->
118
+
let base = (x + y) * 255 / 30 in
119
+
(base + (band * 2)) mod 256)
120
+
in
121
+
let compressed = Hcomp.compress ~bands ~width ~height ~bits_per_sample data in
122
+
let ratio =
123
+
Float.of_int (Bytes.length compressed) /. Float.of_int (Bytes.length data)
124
+
in
125
+
if ratio >= 1.0 then
126
+
Alcotest.failf
127
+
"compression ratio %.3f >= 1.0 for correlated multi-band data \
128
+
(compressed=%d, original=%d)"
129
+
ratio (Bytes.length compressed) (Bytes.length data)
130
+
131
+
(* Invalid input size should raise. *)
132
+
let test_invalid_size () =
133
+
let bands = 2 and width = 4 and height = 4 and bits_per_sample = 8 in
134
+
let data = Bytes.make 10 '\x00' in
135
+
match Hcomp.compress ~bands ~width ~height ~bits_per_sample data with
136
+
| _ -> Alcotest.fail "expected Invalid_argument for wrong data size"
137
+
| exception Invalid_argument _ -> ()
138
+
139
+
(* Invalid bits_per_sample should raise. *)
140
+
let test_invalid_bps () =
141
+
let bands = 1 and width = 2 and height = 2 in
142
+
let data = Bytes.make (1 * 2 * 2 * 2) '\x00' in
143
+
(match Hcomp.compress ~bands ~width ~height ~bits_per_sample:0 data with
144
+
| _ -> Alcotest.fail "expected Invalid_argument for bps=0"
145
+
| exception Invalid_argument _ -> ());
146
+
match Hcomp.compress ~bands ~width ~height ~bits_per_sample:17 data with
147
+
| _ -> Alcotest.fail "expected Invalid_argument for bps=17"
148
+
| exception Invalid_argument _ -> ()
149
+
150
+
let suite =
151
+
( "hcomp",
152
+
[
153
+
("single band constant roundtrip", `Quick, test_single_band_constant);
154
+
("single band gradient roundtrip", `Quick, test_single_band_gradient);
155
+
("multi-band constant roundtrip", `Quick, test_multi_band_constant);
156
+
("multi-band gradient roundtrip", `Quick, test_multi_band_gradient);
157
+
("16-bit samples roundtrip", `Quick, test_16bit_samples);
158
+
("all zeros roundtrip", `Quick, test_all_zeros);
159
+
("single pixel roundtrip", `Quick, test_single_pixel);
160
+
("spectral correlation roundtrip", `Quick, test_spectral_correlation);
161
+
("checkerboard multi-band roundtrip", `Quick, test_checkerboard_multi_band);
162
+
("compression ratio", `Quick, test_compression_ratio);
163
+
("invalid input size", `Quick, test_invalid_size);
164
+
("invalid bits_per_sample", `Quick, test_invalid_bps);
165
+
] )
166
+
167
+
let () = Alcotest.run "hcomp" [ suite ]