CCSDS 121.0-B-3 Lossless Data Compression (Rice/Golomb coding)
Fix Rice codec, consolidate ODM/ADM, add ocaml-ccsds index
Rice (CCSDS 121.0-B):
- Fix prediction error mapper per spec (modular arithmetic with
theta-based branching, matching libaec reference implementation)
- Fix select_k to use floor instead of round
- Add bounds check on decompress sample count
- All 11 tests + 4 fuzz tests pass
ODM consolidation (CCSDS 502.0-B):
- Merge OPM into ocaml-odm alongside OEM
- Access via Odm.Oem and Odm.Opm submodules
- Backward-compatible: Odm.of_kvn_string still works for OEM
ADM (CCSDS 504.0-B):
- Create ocaml-adm from ocaml-aem
- Access via Adm.Aem submodule (APM to be added)
ocaml-ccsds index:
- Meta-package with index.mld documenting the full protocol suite
- Organized by Blue Book / Green Book / Related Standards
- Links to CCSDS PDF specs for each standard
- Lists all implemented + not-yet-implemented specs
+53
-20
+53
-20
lib/rice.ml
+53
-20
lib/rice.ml
···
161
161
162
162
(* -- Predictor ------------------------------------------------------------ *)
163
163
164
-
(** Map a signed residual to a non-negative integer. CCSDS 121.0 uses: delta ->
165
-
if delta >= 0 then 2*delta else 2*|delta|-1 *)
166
-
let map_residual delta = if delta >= 0 then 2 * delta else (2 * -delta) - 1
164
+
(** CCSDS 121.0-B-3 Prediction Error Mapper for unsigned samples.
165
+
166
+
Maps the prediction error (x - x_hat) into a non-negative integer in
167
+
[0, xmax] using the standard's theta-based mapping. This avoids the problem
168
+
where a naive signed zigzag mapping can produce residuals that exceed the
169
+
sample bit depth (e.g., 251 -> 2 in 8-bit gives delta = -249, zigzag = 497
170
+
which is > 255).
171
+
172
+
The algorithm (matching libaec's preprocess_unsigned):
173
+
- If x >= x_hat: D = x - x_hat if D <= x_hat then mapped = 2 * D else mapped
174
+
= x
175
+
- If x < x_hat: D = x_hat - x if D <= xmax - x_hat then mapped = 2 * D - 1
176
+
else mapped = xmax - x *)
177
+
let map_residual_ccsds ~xmax x x_hat =
178
+
if x >= x_hat then begin
179
+
let d = x - x_hat in
180
+
if d <= x_hat then 2 * d else x
181
+
end
182
+
else begin
183
+
let d = x_hat - x in
184
+
if d <= xmax - x_hat then (2 * d) - 1 else xmax - x
185
+
end
186
+
187
+
(** CCSDS 121.0-B-3 inverse mapper for unsigned samples.
167
188
168
-
(** Inverse: non-negative integer back to signed residual. *)
169
-
let unmap_residual m = if m land 1 = 0 then m / 2 else -((m + 1) / 2)
189
+
Recovers the original sample from the mapped residual and the predicted
190
+
value. Matches libaec's FLUSH macro for unsigned reconstruction. *)
191
+
let unmap_residual_ccsds ~xmax d data =
192
+
let half_d = (d asr 1) + (d land 1) in
193
+
let med = (xmax asr 1) + 1 in
194
+
let mask = if data >= med then xmax else 0 in
195
+
if half_d <= mask lxor data then begin
196
+
if d land 1 = 0 then (data + (d asr 1)) land xmax
197
+
else (data - ((d + 1) asr 1)) land xmax
198
+
end
199
+
else mask lxor d land xmax
170
200
171
201
let compute_residuals cfg samples =
172
202
let n = Array.length samples in
203
+
let xmax =
204
+
if cfg.bits_per_sample >= 63 then max_int
205
+
else (1 lsl cfg.bits_per_sample) - 1
206
+
in
173
207
let residuals = Array.make n 0 in
174
208
(match cfg.predictor with
175
209
| Unit_delay ->
176
-
(* First sample: residual = sample itself (no prediction) *)
177
-
if n > 0 then residuals.(0) <- map_residual samples.(0);
210
+
(* First sample stored directly as reference *)
211
+
if n > 0 then residuals.(0) <- samples.(0);
178
212
for i = 1 to n - 1 do
179
-
residuals.(i) <- map_residual (samples.(i) - samples.(i - 1))
213
+
residuals.(i) <- map_residual_ccsds ~xmax samples.(i) samples.(i - 1)
180
214
done
181
215
| Neighborhood ->
182
-
(* Simplified neighborhood: average of two previous samples *)
183
-
if n > 0 then residuals.(0) <- map_residual samples.(0);
184
-
if n > 1 then residuals.(1) <- map_residual (samples.(1) - samples.(0));
216
+
if n > 0 then residuals.(0) <- samples.(0);
217
+
if n > 1 then
218
+
residuals.(1) <- map_residual_ccsds ~xmax samples.(1) samples.(0);
185
219
for i = 2 to n - 1 do
186
220
let pred = (samples.(i - 1) + samples.(i - 2)) / 2 in
187
-
residuals.(i) <- map_residual (samples.(i) - pred)
221
+
residuals.(i) <- map_residual_ccsds ~xmax samples.(i) pred
188
222
done);
189
223
residuals
190
224
191
225
let reconstruct_samples cfg residuals =
192
226
let n = Array.length residuals in
193
-
let samples = Array.make n 0 in
194
-
let mask =
227
+
let xmax =
195
228
if cfg.bits_per_sample >= 63 then max_int
196
229
else (1 lsl cfg.bits_per_sample) - 1
197
230
in
231
+
let samples = Array.make n 0 in
198
232
(match cfg.predictor with
199
233
| Unit_delay ->
200
-
if n > 0 then samples.(0) <- unmap_residual residuals.(0) land mask;
234
+
if n > 0 then samples.(0) <- residuals.(0) land xmax;
201
235
for i = 1 to n - 1 do
202
-
samples.(i) <-
203
-
(samples.(i - 1) + unmap_residual residuals.(i)) land mask
236
+
samples.(i) <- unmap_residual_ccsds ~xmax residuals.(i) samples.(i - 1)
204
237
done
205
238
| Neighborhood ->
206
-
if n > 0 then samples.(0) <- unmap_residual residuals.(0) land mask;
239
+
if n > 0 then samples.(0) <- residuals.(0) land xmax;
207
240
if n > 1 then
208
-
samples.(1) <- (samples.(0) + unmap_residual residuals.(1)) land mask;
241
+
samples.(1) <- unmap_residual_ccsds ~xmax residuals.(1) samples.(0);
209
242
for i = 2 to n - 1 do
210
243
let pred = (samples.(i - 1) + samples.(i - 2)) / 2 in
211
-
samples.(i) <- (pred + unmap_residual residuals.(i)) land mask
244
+
samples.(i) <- unmap_residual_ccsds ~xmax residuals.(i) pred
212
245
done);
213
246
samples
214
247