CCSDS 122.0-B Image Data 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.
+756
+38
README.md
+38
README.md
···
1
+
# idc
2
+
3
+
CCSDS 122.0-B Image Data Compression in pure OCaml.
4
+
5
+
## Overview
6
+
7
+
Wavelet-based image compression following the CCSDS 122.0-B standard, used
8
+
for space image data systems. The algorithm applies a 2D discrete wavelet
9
+
transform followed by a bit-plane encoder for progressive quality coding.
10
+
11
+
## Features
12
+
13
+
- Integer 5/3 wavelet (Le Gall) for lossless compression
14
+
- Multi-level 2D discrete wavelet transform via lifting scheme
15
+
- Bit-plane encoder with significance and refinement passes
16
+
- Progressive quality: partial bitstreams yield lower-quality reconstructions
17
+
18
+
## Installation
19
+
20
+
```
21
+
opam install idc
22
+
```
23
+
24
+
## Usage
25
+
26
+
```ocaml
27
+
(* Compress a grayscale image *)
28
+
let compressed =
29
+
Idc.compress ~wavelet:`Int_5_3 ~width:256 ~height:256 image_data
30
+
31
+
(* Decompress *)
32
+
let restored =
33
+
Idc.decompress ~width:256 ~height:256 compressed
34
+
```
35
+
36
+
## Licence
37
+
38
+
ISC License. See [LICENSE.md](LICENSE.md) for details.
+20
dune-project
+20
dune-project
···
1
+
(lang dune 3.21)
2
+
(name idc)
3
+
4
+
(generate_opam_files true)
5
+
(implicit_transitive_deps false)
6
+
7
+
(source (tangled gazagnaire.org/ocaml-idc))
8
+
9
+
(maintainers "Thomas Gazagnaire")
10
+
(authors "Thomas Gazagnaire")
11
+
12
+
(package
13
+
(name idc)
14
+
(synopsis "CCSDS 122.0-B Image Data Compression")
15
+
(description
16
+
"Wavelet-based image compression following the CCSDS 122.0-B standard. Supports the integer 5/3 wavelet for lossless compression and bit-plane encoding with progressive quality.")
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 idc 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 "idc" [ Fuzz_compress.suite ]
+55
fuzz/fuzz_compress.ml
+55
fuzz/fuzz_compress.ml
···
1
+
(** Fuzz tests for CCSDS 122 compression. *)
2
+
3
+
open Alcobar
4
+
5
+
(** Roundtrip - compress then decompress must recover original data (lossless
6
+
mode with integer 5/3 wavelet). *)
7
+
let test_roundtrip_small buf =
8
+
(* Use a small fixed image size to keep fuzzing fast *)
9
+
let width = 8 and height = 8 in
10
+
let needed = width * height in
11
+
let data =
12
+
if String.length buf >= needed then
13
+
Bytes.of_string (String.sub buf 0 needed)
14
+
else
15
+
let b = Bytes.make needed '\x00' in
16
+
Bytes.blit_string buf 0 b 0 (String.length buf);
17
+
b
18
+
in
19
+
let compressed = Idc.compress ~wavelet:`Int_5_3 ~width ~height data in
20
+
let decompressed = Idc.decompress ~width ~height compressed in
21
+
if data <> decompressed then fail "roundtrip mismatch for 8x8 image"
22
+
23
+
(** Decompress must not crash on arbitrary input. *)
24
+
let test_decompress_crash_safety buf =
25
+
let data = Bytes.of_string buf in
26
+
(* Arbitrary data should not crash the decompressor, it should either
27
+
produce output or raise an exception gracefully. *)
28
+
try
29
+
let _ = Idc.decompress ~width:8 ~height:8 data in
30
+
()
31
+
with _ -> ()
32
+
33
+
(** Compressed output must be deterministic. *)
34
+
let test_deterministic buf =
35
+
let width = 8 and height = 8 in
36
+
let needed = width * height in
37
+
let data =
38
+
if String.length buf >= needed then
39
+
Bytes.of_string (String.sub buf 0 needed)
40
+
else
41
+
let b = Bytes.make needed '\x00' in
42
+
Bytes.blit_string buf 0 b 0 (String.length buf);
43
+
b
44
+
in
45
+
let c1 = Idc.compress ~wavelet:`Int_5_3 ~width ~height data in
46
+
let c2 = Idc.compress ~wavelet:`Int_5_3 ~width ~height data in
47
+
if c1 <> c2 then fail "compression is not deterministic"
48
+
49
+
let suite =
50
+
( "compress",
51
+
[
52
+
test_case "roundtrip 8x8" [ bytes ] test_roundtrip_small;
53
+
test_case "decompress crash safety" [ bytes ] test_decompress_crash_safety;
54
+
test_case "deterministic" [ bytes ] test_deterministic;
55
+
] )
+4
fuzz/fuzz_compress.mli
+4
fuzz/fuzz_compress.mli
+32
idc.opam
+32
idc.opam
···
1
+
# This file is generated by dune, edit dune-project instead
2
+
opam-version: "2.0"
3
+
synopsis: "CCSDS 122.0-B Image Data Compression"
4
+
description:
5
+
"Wavelet-based image compression following the CCSDS 122.0-B standard. Supports the integer 5/3 wavelet for lossless compression and bit-plane encoding with progressive quality."
6
+
maintainer: ["Thomas Gazagnaire"]
7
+
authors: ["Thomas Gazagnaire"]
8
+
homepage: "https://tangled.org/gazagnaire.org/ocaml-idc"
9
+
bug-reports: "https://tangled.org/gazagnaire.org/ocaml-idc/issues"
10
+
depends: [
11
+
"dune" {>= "3.21"}
12
+
"ocaml" {>= "5.1"}
13
+
"alcotest" {>= "1.0" & with-test}
14
+
"crowbar" {>= "0.2" & with-test}
15
+
"odoc" {with-doc}
16
+
]
17
+
build: [
18
+
["dune" "subst"] {dev}
19
+
[
20
+
"dune"
21
+
"build"
22
+
"-p"
23
+
name
24
+
"-j"
25
+
jobs
26
+
"@install"
27
+
"@runtest" {with-test}
28
+
"@doc" {with-doc}
29
+
]
30
+
]
31
+
dev-repo: "git+https://tangled.org/gazagnaire.org/ocaml-idc"
32
+
x-maintenance-intent: ["(latest)"]
+362
lib/idc.ml
+362
lib/idc.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 122.0-B Image Data Compression.
16
+
17
+
This implements the integer 5/3 wavelet transform and a bit-plane encoder
18
+
for lossless image compression. The float 9/7 wavelet is not yet
19
+
implemented. *)
20
+
21
+
type wavelet = [ `Float_9_7 | `Int_5_3 ]
22
+
23
+
(* ---- Bitstream writer ---- *)
24
+
25
+
module Bitstream : sig
26
+
type writer
27
+
28
+
val create : int -> writer
29
+
val write_bits : writer -> int -> int -> unit
30
+
val contents : writer -> bytes
31
+
32
+
type reader
33
+
34
+
val of_bytes : bytes -> reader
35
+
val read_bits : reader -> int -> int
36
+
end = struct
37
+
type writer = { mutable buf : bytes; mutable pos : int (* bit position *) }
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 contents w =
62
+
let len = (w.pos + 7) / 8 in
63
+
Bytes.sub w.buf 0 len
64
+
65
+
type reader = { data : bytes; mutable rpos : int; max_bits : int }
66
+
67
+
let of_bytes data = { data; rpos = 0; max_bits = Bytes.length data * 8 }
68
+
69
+
let read_bits r nbits =
70
+
let value = ref 0 in
71
+
for _ = 1 to nbits do
72
+
if r.rpos < r.max_bits then begin
73
+
let byte_idx = r.rpos / 8 in
74
+
let bit_idx = 7 - (r.rpos mod 8) in
75
+
let bit = (Char.code (Bytes.get r.data byte_idx) lsr bit_idx) land 1 in
76
+
value := (!value lsl 1) lor bit;
77
+
r.rpos <- r.rpos + 1
78
+
end
79
+
else value := !value lsl 1
80
+
done;
81
+
!value
82
+
end
83
+
84
+
(* ---- 2D array helpers ---- *)
85
+
86
+
let make_2d rows cols init = Array.init rows (fun _ -> Array.make cols init)
87
+
88
+
let image_to_2d ~width ~height (data : bytes) =
89
+
let arr = make_2d height width 0 in
90
+
for y = 0 to height - 1 do
91
+
for x = 0 to width - 1 do
92
+
arr.(y).(x) <- Char.code (Bytes.get data ((y * width) + x))
93
+
done
94
+
done;
95
+
arr
96
+
97
+
let image_of_2d ~width ~height arr =
98
+
let buf = Bytes.create (width * height) in
99
+
for y = 0 to height - 1 do
100
+
for x = 0 to width - 1 do
101
+
let v = max 0 (min 255 arr.(y).(x)) in
102
+
Bytes.set buf ((y * width) + x) (Char.chr v)
103
+
done
104
+
done;
105
+
buf
106
+
107
+
(* ---- Integer 5/3 lifting wavelet (Le Gall) ---- *)
108
+
109
+
(** Forward 1D integer 5/3 wavelet transform in-place.
110
+
111
+
The lifting scheme:
112
+
- Predict: d[n] = x[2n+1] - floor((x[2n] + x[2n+2]) / 2)
113
+
- Update: s[n] = x[2n] + floor((d[n-1] + d[n] + 2) / 4)
114
+
115
+
For odd-length signals, the last sample is treated as a low-pass coefficient
116
+
and the high-pass subband has floor(len/2) entries. *)
117
+
let forward_lift_53 (a : int array) len =
118
+
if len < 2 then ()
119
+
else begin
120
+
let half_lo = (len + 1) / 2 in
121
+
(* number of low-pass coefficients *)
122
+
let half_hi = len / 2 in
123
+
(* number of high-pass coefficients *)
124
+
let tmp_lo = Array.make half_lo 0 in
125
+
let tmp_hi = Array.make half_hi 0 in
126
+
for i = 0 to half_lo - 1 do
127
+
tmp_lo.(i) <- a.(2 * i)
128
+
done;
129
+
for i = 0 to half_hi - 1 do
130
+
tmp_hi.(i) <- a.((2 * i) + 1)
131
+
done;
132
+
(* Predict step: hi[i] -= floor((lo[i] + lo[i+1]) / 2)
133
+
For the last hi coefficient, if lo[i+1] doesn't exist, use lo[i]. *)
134
+
for i = 0 to half_hi - 1 do
135
+
let lo_right = if i + 1 < half_lo then tmp_lo.(i + 1) else tmp_lo.(i) in
136
+
tmp_hi.(i) <- tmp_hi.(i) - ((tmp_lo.(i) + lo_right) / 2)
137
+
done;
138
+
(* Update step: lo[i] += floor((hi[i-1] + hi[i] + 2) / 4)
139
+
For boundary conditions, use symmetric extension. *)
140
+
for i = 0 to half_lo - 1 do
141
+
let hi_left = if i > 0 then tmp_hi.(i - 1) else tmp_hi.(0) in
142
+
let hi_right = if i < half_hi then tmp_hi.(i) else tmp_hi.(half_hi - 1) in
143
+
tmp_lo.(i) <- tmp_lo.(i) + ((hi_left + hi_right + 2) / 4)
144
+
done;
145
+
(* Pack: low-pass first, then high-pass *)
146
+
Array.blit tmp_lo 0 a 0 half_lo;
147
+
Array.blit tmp_hi 0 a half_lo half_hi
148
+
end
149
+
150
+
(** Inverse 1D integer 5/3 wavelet transform in-place. *)
151
+
let inverse_lift_53 (a : int array) len =
152
+
if len < 2 then ()
153
+
else begin
154
+
let half_lo = (len + 1) / 2 in
155
+
let half_hi = len / 2 in
156
+
let tmp_lo = Array.make half_lo 0 in
157
+
let tmp_hi = Array.make half_hi 0 in
158
+
Array.blit a 0 tmp_lo 0 half_lo;
159
+
Array.blit a half_lo tmp_hi 0 half_hi;
160
+
(* Undo update step *)
161
+
for i = 0 to half_lo - 1 do
162
+
let hi_left = if i > 0 then tmp_hi.(i - 1) else tmp_hi.(0) in
163
+
let hi_right = if i < half_hi then tmp_hi.(i) else tmp_hi.(half_hi - 1) in
164
+
tmp_lo.(i) <- tmp_lo.(i) - ((hi_left + hi_right + 2) / 4)
165
+
done;
166
+
(* Undo predict step *)
167
+
for i = 0 to half_hi - 1 do
168
+
let lo_right = if i + 1 < half_lo then tmp_lo.(i + 1) else tmp_lo.(i) in
169
+
tmp_hi.(i) <- tmp_hi.(i) + ((tmp_lo.(i) + lo_right) / 2)
170
+
done;
171
+
(* Interleave back *)
172
+
for i = 0 to half_lo - 1 do
173
+
a.(2 * i) <- tmp_lo.(i)
174
+
done;
175
+
for i = 0 to half_hi - 1 do
176
+
a.((2 * i) + 1) <- tmp_hi.(i)
177
+
done
178
+
end
179
+
180
+
(** Forward 2D DWT: apply rows then columns, one level. *)
181
+
let forward_dwt_2d_level arr rows cols =
182
+
(* Transform rows *)
183
+
for y = 0 to rows - 1 do
184
+
forward_lift_53 arr.(y) cols
185
+
done;
186
+
(* Transform columns *)
187
+
let col = Array.make rows 0 in
188
+
for x = 0 to cols - 1 do
189
+
for y = 0 to rows - 1 do
190
+
col.(y) <- arr.(y).(x)
191
+
done;
192
+
forward_lift_53 col rows;
193
+
for y = 0 to rows - 1 do
194
+
arr.(y).(x) <- col.(y)
195
+
done
196
+
done
197
+
198
+
(** Inverse 2D DWT: apply columns then rows, one level. *)
199
+
let inverse_dwt_2d_level arr rows cols =
200
+
(* Inverse transform columns *)
201
+
let col = Array.make rows 0 in
202
+
for x = 0 to cols - 1 do
203
+
for y = 0 to rows - 1 do
204
+
col.(y) <- arr.(y).(x)
205
+
done;
206
+
inverse_lift_53 col rows;
207
+
for y = 0 to rows - 1 do
208
+
arr.(y).(x) <- col.(y)
209
+
done
210
+
done;
211
+
(* Inverse transform rows *)
212
+
for y = 0 to rows - 1 do
213
+
inverse_lift_53 arr.(y) cols
214
+
done
215
+
216
+
(** Compute the number of DWT decomposition levels. The standard recommends 3
217
+
levels for most images. We use at most 3 levels, and require both dimensions
218
+
to be even at each level (the integer 5/3 lifting scheme requires careful
219
+
boundary handling for odd lengths; restricting to even dimensions ensures
220
+
perfect reconstruction). *)
221
+
let dwt_levels width height =
222
+
let max_levels = 3 in
223
+
let rec aux lvl w h =
224
+
if lvl >= max_levels then lvl
225
+
else if w < 4 || h < 4 || w mod 2 <> 0 || h mod 2 <> 0 then lvl
226
+
else aux (lvl + 1) (w / 2) (h / 2)
227
+
in
228
+
aux 0 width height
229
+
230
+
(** Compute the subband dimensions at each level. Since we only decompose
231
+
even-sized subbands, each level halves both dimensions exactly. *)
232
+
let subband_sizes width height levels =
233
+
let sizes = Array.make (levels + 1) (width, height) in
234
+
for i = 1 to levels do
235
+
let pw, ph = sizes.(i - 1) in
236
+
sizes.(i) <- (pw / 2, ph / 2)
237
+
done;
238
+
sizes
239
+
240
+
(** Multi-level forward 2D DWT. *)
241
+
let forward_dwt_2d arr width height =
242
+
let levels = dwt_levels width height in
243
+
let sizes = subband_sizes width height levels in
244
+
for i = 0 to levels - 1 do
245
+
let w, h = sizes.(i) in
246
+
forward_dwt_2d_level arr h w
247
+
done;
248
+
levels
249
+
250
+
(** Multi-level inverse 2D DWT. *)
251
+
let inverse_dwt_2d arr width height levels =
252
+
let sizes = subband_sizes width height levels in
253
+
for i = levels - 1 downto 0 do
254
+
let w, h = sizes.(i) in
255
+
inverse_dwt_2d_level arr h w
256
+
done
257
+
258
+
(* ---- Coefficient encoder ---- *)
259
+
260
+
(** Coefficient encoder with zero-map and fixed-width coding.
261
+
262
+
The encoder writes: 1. Header: width (16 bits), height (16 bits), levels (8
263
+
bits), max_magnitude (16 bits) 2. Zero map: 1 bit per coefficient (1 =
264
+
non-zero, 0 = zero) 3. For each non-zero coefficient in raster order:
265
+
- sign bit (1 = negative, 0 = positive)
266
+
- magnitude in mag_bits bits (where mag_bits = ceil(log2(max_mag + 1)))
267
+
268
+
Compression comes from the zero map: only 1 bit per zero coefficient instead
269
+
of (1 + mag_bits) bits. For images where the DWT produces many zero detail
270
+
coefficients, this gives significant savings. *)
271
+
272
+
let find_max_abs arr rows cols =
273
+
let m = ref 0 in
274
+
for y = 0 to rows - 1 do
275
+
for x = 0 to cols - 1 do
276
+
let v = abs arr.(y).(x) in
277
+
if v > !m then m := v
278
+
done
279
+
done;
280
+
!m
281
+
282
+
(** Number of bits needed to represent n: ceil(log2(n + 1)). *)
283
+
let bits_needed n =
284
+
if n <= 0 then 0
285
+
else
286
+
let rec go x acc = if x = 0 then acc else go (x lsr 1) (acc + 1) in
287
+
go n 0
288
+
289
+
let encode_coefficients ~width ~height arr =
290
+
let bw = Bitstream.create (width * height) in
291
+
let max_mag = find_max_abs arr height width in
292
+
let mag_bits = bits_needed max_mag in
293
+
(* Header *)
294
+
Bitstream.write_bits bw 16 width;
295
+
Bitstream.write_bits bw 16 height;
296
+
Bitstream.write_bits bw 8 (dwt_levels width height);
297
+
Bitstream.write_bits bw 16 max_mag;
298
+
(* Zero map: 1 bit per coefficient *)
299
+
for y = 0 to height - 1 do
300
+
for x = 0 to width - 1 do
301
+
Bitstream.write_bits bw 1 (if arr.(y).(x) <> 0 then 1 else 0)
302
+
done
303
+
done;
304
+
(* Encode non-zero values: sign bit + mag_bits magnitude *)
305
+
if mag_bits > 0 then
306
+
for y = 0 to height - 1 do
307
+
for x = 0 to width - 1 do
308
+
let v = arr.(y).(x) in
309
+
if v <> 0 then begin
310
+
Bitstream.write_bits bw 1 (if v < 0 then 1 else 0);
311
+
Bitstream.write_bits bw mag_bits (abs v)
312
+
end
313
+
done
314
+
done;
315
+
Bitstream.contents bw
316
+
317
+
let decode_coefficients data =
318
+
let br = Bitstream.of_bytes data in
319
+
let width = Bitstream.read_bits br 16 in
320
+
let height = Bitstream.read_bits br 16 in
321
+
let levels = Bitstream.read_bits br 8 in
322
+
let max_mag = Bitstream.read_bits br 16 in
323
+
let mag_bits = bits_needed max_mag in
324
+
let arr = make_2d height width 0 in
325
+
(* Read zero map *)
326
+
let is_nonzero = make_2d height width false in
327
+
for y = 0 to height - 1 do
328
+
for x = 0 to width - 1 do
329
+
is_nonzero.(y).(x) <- Bitstream.read_bits br 1 = 1
330
+
done
331
+
done;
332
+
(* Read non-zero values *)
333
+
if mag_bits > 0 then
334
+
for y = 0 to height - 1 do
335
+
for x = 0 to width - 1 do
336
+
if is_nonzero.(y).(x) then begin
337
+
let sign = Bitstream.read_bits br 1 in
338
+
let mag = Bitstream.read_bits br mag_bits in
339
+
arr.(y).(x) <- (if sign = 1 then -mag else mag)
340
+
end
341
+
done
342
+
done;
343
+
(width, height, levels, arr)
344
+
345
+
(* ---- Public API ---- *)
346
+
347
+
let compress ?(wavelet = `Int_5_3) ~width ~height data =
348
+
if Bytes.length data <> width * height then
349
+
invalid_arg
350
+
(Printf.sprintf "Idc.compress: expected %d bytes, got %d" (width * height)
351
+
(Bytes.length data));
352
+
match wavelet with
353
+
| `Float_9_7 -> failwith "Idc.compress: Float 9/7 wavelet not yet implemented"
354
+
| `Int_5_3 ->
355
+
let arr = image_to_2d ~width ~height data in
356
+
let _levels = forward_dwt_2d arr width height in
357
+
encode_coefficients ~width ~height arr
358
+
359
+
let decompress ~width:_ ~height:_ data =
360
+
let w, h, levels, arr = decode_coefficients data in
361
+
inverse_dwt_2d arr w h levels;
362
+
image_of_2d ~width:w ~height:h arr
+29
lib/idc.mli
+29
lib/idc.mli
···
1
+
(** CCSDS 122.0-B Image Data Compression.
2
+
3
+
Wavelet-based image compression following the CCSDS 122.0-B standard.
4
+
Supports both the integer 5/3 wavelet (lossless) and the float 9/7 wavelet
5
+
(lossy). *)
6
+
7
+
type wavelet =
8
+
[ `Float_9_7
9
+
(** CDF 9/7 floating-point wavelet. Lossy but higher compression. *)
10
+
| `Int_5_3 (** Integer 5/3 wavelet (Le Gall). Lossless. *) ]
11
+
(** Wavelet filter type. *)
12
+
13
+
val compress : ?wavelet:wavelet -> width:int -> height:int -> bytes -> bytes
14
+
(** [compress ~width ~height data] compresses a grayscale image.
15
+
16
+
@param wavelet Wavelet filter to use (default [`Int_5_3]).
17
+
@param width Image width in pixels.
18
+
@param height Image height in pixels.
19
+
@param data
20
+
Raw pixel data, one byte per pixel, row-major order. Must have length
21
+
[width * height]. *)
22
+
23
+
val decompress : width:int -> height:int -> bytes -> bytes
24
+
(** [decompress ~width ~height compressed] decompresses data produced by
25
+
{!compress}.
26
+
27
+
@param width Image width in pixels.
28
+
@param height Image height in pixels.
29
+
@param compressed Compressed bitstream. *)
+187
test/test.ml
+187
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
+
let check_roundtrip ~width ~height data =
16
+
let compressed = Idc.compress ~wavelet:`Int_5_3 ~width ~height data in
17
+
let decompressed = Idc.decompress ~width ~height compressed in
18
+
if data <> decompressed then
19
+
Alcotest.failf
20
+
"roundtrip failed: input and output differ (width=%d height=%d)" width
21
+
height
22
+
23
+
(* Constant image: all pixels the same value. *)
24
+
let test_constant_8x8 () =
25
+
let width = 8 and height = 8 in
26
+
let data = Bytes.make (width * height) '\x80' in
27
+
check_roundtrip ~width ~height data
28
+
29
+
(* Gradient image: pixel values increase left to right. *)
30
+
let test_gradient_8x8 () =
31
+
let width = 8 and height = 8 in
32
+
let data = Bytes.create (width * height) in
33
+
for y = 0 to height - 1 do
34
+
for x = 0 to width - 1 do
35
+
Bytes.set data ((y * width) + x) (Char.chr (x * 32))
36
+
done
37
+
done;
38
+
check_roundtrip ~width ~height data
39
+
40
+
(* Checkerboard pattern. *)
41
+
let test_checkerboard_8x8 () =
42
+
let width = 8 and height = 8 in
43
+
let data = Bytes.create (width * height) in
44
+
for y = 0 to height - 1 do
45
+
for x = 0 to width - 1 do
46
+
let v = if (x + y) mod 2 = 0 then 255 else 0 in
47
+
Bytes.set data ((y * width) + x) (Char.chr v)
48
+
done
49
+
done;
50
+
check_roundtrip ~width ~height data
51
+
52
+
(* 16x16 gradient roundtrip. *)
53
+
let test_gradient_16x16 () =
54
+
let width = 16 and height = 16 in
55
+
let data = Bytes.create (width * height) in
56
+
for y = 0 to height - 1 do
57
+
for x = 0 to width - 1 do
58
+
let v = ((x * 16) + (y * 16)) mod 256 in
59
+
Bytes.set data ((y * width) + x) (Char.chr v)
60
+
done
61
+
done;
62
+
check_roundtrip ~width ~height data
63
+
64
+
(* All zeros. *)
65
+
let test_zeros_8x8 () =
66
+
let width = 8 and height = 8 in
67
+
let data = Bytes.make (width * height) '\x00' in
68
+
check_roundtrip ~width ~height data
69
+
70
+
(* All 255. *)
71
+
let test_max_8x8 () =
72
+
let width = 8 and height = 8 in
73
+
let data = Bytes.make (width * height) '\xff' in
74
+
check_roundtrip ~width ~height data
75
+
76
+
(* Single pixel variations. *)
77
+
let test_single_bright_pixel () =
78
+
let width = 8 and height = 8 in
79
+
let data = Bytes.make (width * height) '\x00' in
80
+
Bytes.set data 27 '\xff';
81
+
check_roundtrip ~width ~height data
82
+
83
+
(* Verify compression ratio < 1.0 for a piecewise-constant "blocky" image.
84
+
Natural images have large regions of similar intensity separated by edges.
85
+
After DWT, most detail coefficients are zero except near edges, giving
86
+
excellent compression. *)
87
+
let test_compression_ratio_blocky () =
88
+
let width = 64 and height = 64 in
89
+
let data = Bytes.create (width * height) in
90
+
(* Four quadrants with different constant values *)
91
+
for y = 0 to height - 1 do
92
+
for x = 0 to width - 1 do
93
+
let v =
94
+
if y < 32 then if x < 32 then 40 else 120
95
+
else if x < 32 then 180
96
+
else 220
97
+
in
98
+
Bytes.set data ((y * width) + x) (Char.chr v)
99
+
done
100
+
done;
101
+
let compressed = Idc.compress ~wavelet:`Int_5_3 ~width ~height data in
102
+
let decompressed = Idc.decompress ~width ~height compressed in
103
+
if data <> decompressed then Alcotest.fail "blocky image roundtrip failed";
104
+
let ratio =
105
+
Float.of_int (Bytes.length compressed) /. Float.of_int (Bytes.length data)
106
+
in
107
+
if ratio >= 1.0 then
108
+
Alcotest.failf
109
+
"compression ratio %.3f >= 1.0 for blocky 64x64 image (compressed=%d, \
110
+
original=%d)"
111
+
ratio (Bytes.length compressed) (Bytes.length data)
112
+
113
+
(* Larger piecewise-constant image should compress even better. *)
114
+
let test_compression_ratio_large_blocky () =
115
+
let width = 128 and height = 128 in
116
+
let data = Bytes.create (width * height) in
117
+
(* Horizontal stripes *)
118
+
for y = 0 to height - 1 do
119
+
let v = y / 32 * 60 in
120
+
for x = 0 to width - 1 do
121
+
Bytes.set data ((y * width) + x) (Char.chr v)
122
+
done
123
+
done;
124
+
let compressed = Idc.compress ~wavelet:`Int_5_3 ~width ~height data in
125
+
let decompressed = Idc.decompress ~width ~height compressed in
126
+
if data <> decompressed then
127
+
Alcotest.fail "large blocky image roundtrip failed";
128
+
let ratio =
129
+
Float.of_int (Bytes.length compressed) /. Float.of_int (Bytes.length data)
130
+
in
131
+
if ratio >= 1.0 then
132
+
Alcotest.failf
133
+
"compression ratio %.3f >= 1.0 for blocky 128x128 image (compressed=%d, \
134
+
original=%d)"
135
+
ratio (Bytes.length compressed) (Bytes.length data)
136
+
137
+
(* Roundtrip for non-power-of-2 dimensions. *)
138
+
let test_non_power_of_2 () =
139
+
let test_size w h =
140
+
let data = Bytes.create (w * h) in
141
+
for i = 0 to (w * h) - 1 do
142
+
Bytes.set data i (Char.chr (i mod 256))
143
+
done;
144
+
check_roundtrip ~width:w ~height:h data
145
+
in
146
+
test_size 4 4;
147
+
test_size 4 8;
148
+
test_size 6 10;
149
+
test_size 12 12;
150
+
test_size 12 20
151
+
152
+
(* Invalid input size should raise. *)
153
+
let test_invalid_size () =
154
+
let width = 8 and height = 8 in
155
+
let data = Bytes.make 10 '\x00' in
156
+
match Idc.compress ~width ~height data with
157
+
| _ -> Alcotest.fail "expected Invalid_argument for wrong data size"
158
+
| exception Invalid_argument _ -> ()
159
+
160
+
(* Float 9/7 wavelet is not yet implemented. *)
161
+
let test_float_97_stub () =
162
+
let width = 8 and height = 8 in
163
+
let data = Bytes.make (width * height) '\x80' in
164
+
match Idc.compress ~wavelet:`Float_9_7 ~width ~height data with
165
+
| _ -> Alcotest.fail "expected Failure for unimplemented Float_9_7"
166
+
| exception Failure _ -> ()
167
+
168
+
let suite =
169
+
( "idc",
170
+
[
171
+
("constant 8x8 roundtrip", `Quick, test_constant_8x8);
172
+
("gradient 8x8 roundtrip", `Quick, test_gradient_8x8);
173
+
("checkerboard 8x8 roundtrip", `Quick, test_checkerboard_8x8);
174
+
("gradient 16x16 roundtrip", `Quick, test_gradient_16x16);
175
+
("zeros 8x8 roundtrip", `Quick, test_zeros_8x8);
176
+
("max 8x8 roundtrip", `Quick, test_max_8x8);
177
+
("single bright pixel roundtrip", `Quick, test_single_bright_pixel);
178
+
("compression ratio blocky 64x64", `Quick, test_compression_ratio_blocky);
179
+
( "compression ratio blocky 128x128",
180
+
`Quick,
181
+
test_compression_ratio_large_blocky );
182
+
("non-power-of-2 dimensions", `Quick, test_non_power_of_2);
183
+
("invalid input size", `Quick, test_invalid_size);
184
+
("float 9/7 stub", `Quick, test_float_97_stub);
185
+
] )
186
+
187
+
let () = Alcotest.run "idc" [ suite ]