Implement Float 9/7 wavelet, fix hcomp 16-bit bug, add interop skeleton

+338 -16

2 changed files

expand all

lib

idc.ml

test

test.ml

+251 -10

lib/idc.ml

··· 14 14 15 15 (** CCSDS 122.0-B Image Data Compression. 16 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. *) 17 + This implements the integer 5/3 wavelet transform (lossless) and the CDF 9/7 18 + floating-point wavelet transform (lossy) with a bit-plane encoder for image 19 + compression. *) 20 20 21 21 type wavelet = [ `Float_9_7 | `Int_5_3 ] 22 22 ··· 255 255 inverse_dwt_2d_level arr h w 256 256 done 257 257 258 + (* ---- CDF 9/7 floating-point wavelet ---- *) 259 + 260 + (** Lifting coefficients for the CDF 9/7 wavelet. *) 261 + let alpha = -1.586134342 262 + 263 + let beta = -0.052980118 264 + let gamma = 0.882911075 265 + let delta = 0.443506852 266 + let k_scale = 1.149604398 267 + 268 + (** Forward 1D CDF 9/7 wavelet transform. 269 + 270 + Uses the 4-step lifting scheme with symmetric boundary extension: 1. Predict 271 + 1 (alpha): odd += alpha * (even_left + even_right) 2. Update 1 (beta): even 272 + += beta * (odd_left + odd_right) 3. Predict 2 (gamma): odd += gamma * 273 + (even_left + even_right) 4. Update 2 (delta): even += delta * (odd_left + 274 + odd_right) 5. Scale: even *= K, odd *= 1/K 275 + 276 + Returns (low, high) float arrays. *) 277 + let forward_lift_97 (a : float array) len = 278 + if len < 2 then (Array.copy a, [||]) 279 + else begin 280 + let half_lo = (len + 1) / 2 in 281 + let half_hi = len / 2 in 282 + let even = Array.init half_lo (fun i -> a.(2 * i)) in 283 + let odd = Array.init half_hi (fun i -> a.((2 * i) + 1)) in 284 + (* Step 1: Predict 1 *) 285 + for i = 0 to half_hi - 1 do 286 + let e_right = 287 + if i + 1 < half_lo then even.(i + 1) else even.(half_lo - 1) 288 + in 289 + odd.(i) <- odd.(i) +. (alpha *. (even.(i) +. e_right)) 290 + done; 291 + (* Step 2: Update 1 *) 292 + for i = 0 to half_lo - 1 do 293 + let o_left = if i > 0 then odd.(i - 1) else odd.(0) in 294 + let o_right = if i < half_hi then odd.(i) else odd.(half_hi - 1) in 295 + even.(i) <- even.(i) +. (beta *. (o_left +. o_right)) 296 + done; 297 + (* Step 3: Predict 2 *) 298 + for i = 0 to half_hi - 1 do 299 + let e_right = 300 + if i + 1 < half_lo then even.(i + 1) else even.(half_lo - 1) 301 + in 302 + odd.(i) <- odd.(i) +. (gamma *. (even.(i) +. e_right)) 303 + done; 304 + (* Step 4: Update 2 *) 305 + for i = 0 to half_lo - 1 do 306 + let o_left = if i > 0 then odd.(i - 1) else odd.(0) in 307 + let o_right = if i < half_hi then odd.(i) else odd.(half_hi - 1) in 308 + even.(i) <- even.(i) +. (delta *. (o_left +. o_right)) 309 + done; 310 + (* Step 5: Scale *) 311 + for i = 0 to half_lo - 1 do 312 + even.(i) <- even.(i) *. k_scale 313 + done; 314 + for i = 0 to half_hi - 1 do 315 + odd.(i) <- odd.(i) /. k_scale 316 + done; 317 + (even, odd) 318 + end 319 + 320 + (** Inverse 1D CDF 9/7 wavelet transform. 321 + 322 + Reverses the forward lifting steps in opposite order. *) 323 + let inverse_lift_97 (low : float array) (high : float array) len = 324 + if len < 2 then Array.copy low 325 + else begin 326 + let half_lo = (len + 1) / 2 in 327 + let half_hi = len / 2 in 328 + let even = Array.copy low in 329 + let odd = Array.copy high in 330 + (* Undo scale *) 331 + for i = 0 to half_lo - 1 do 332 + even.(i) <- even.(i) /. k_scale 333 + done; 334 + for i = 0 to half_hi - 1 do 335 + odd.(i) <- odd.(i) *. k_scale 336 + done; 337 + (* Undo step 4: Update 2 *) 338 + for i = 0 to half_lo - 1 do 339 + let o_left = if i > 0 then odd.(i - 1) else odd.(0) in 340 + let o_right = if i < half_hi then odd.(i) else odd.(half_hi - 1) in 341 + even.(i) <- even.(i) -. (delta *. (o_left +. o_right)) 342 + done; 343 + (* Undo step 3: Predict 2 *) 344 + for i = 0 to half_hi - 1 do 345 + let e_right = 346 + if i + 1 < half_lo then even.(i + 1) else even.(half_lo - 1) 347 + in 348 + odd.(i) <- odd.(i) -. (gamma *. (even.(i) +. e_right)) 349 + done; 350 + (* Undo step 2: Update 1 *) 351 + for i = 0 to half_lo - 1 do 352 + let o_left = if i > 0 then odd.(i - 1) else odd.(0) in 353 + let o_right = if i < half_hi then odd.(i) else odd.(half_hi - 1) in 354 + even.(i) <- even.(i) -. (beta *. (o_left +. o_right)) 355 + done; 356 + (* Undo step 1: Predict 1 *) 357 + for i = 0 to half_hi - 1 do 358 + let e_right = 359 + if i + 1 < half_lo then even.(i + 1) else even.(half_lo - 1) 360 + in 361 + odd.(i) <- odd.(i) -. (alpha *. (even.(i) +. e_right)) 362 + done; 363 + (* Interleave *) 364 + let result = Array.make len 0.0 in 365 + for i = 0 to half_lo - 1 do 366 + result.(2 * i) <- even.(i) 367 + done; 368 + for i = 0 to half_hi - 1 do 369 + result.((2 * i) + 1) <- odd.(i) 370 + done; 371 + result 372 + end 373 + 374 + (** Forward 1D CDF 9/7 on an int array -- packs result back into the array as 375 + low then high. Returns the float subbands for later use in the 2D transform. 376 + *) 377 + let forward_97_row (a : float array) len = 378 + if len < 2 then () 379 + else begin 380 + let low, high = forward_lift_97 a len in 381 + let half_lo = (len + 1) / 2 in 382 + let half_hi = len / 2 in 383 + Array.blit low 0 a 0 half_lo; 384 + Array.blit high 0 a half_lo half_hi 385 + end 386 + 387 + (** Inverse 1D CDF 9/7 on a float array packed as [low | high]. *) 388 + let inverse_97_row (a : float array) len = 389 + if len < 2 then () 390 + else begin 391 + let half_lo = (len + 1) / 2 in 392 + let half_hi = len / 2 in 393 + let low = Array.sub a 0 half_lo in 394 + let high = Array.sub a half_lo half_hi in 395 + let result = inverse_lift_97 low high len in 396 + Array.blit result 0 a 0 len 397 + end 398 + 399 + (** Forward 2D CDF 9/7 DWT, one level. Works on float 2D array. *) 400 + let forward_dwt_97_2d_level (arr : float array array) rows cols = 401 + (* Transform rows *) 402 + for y = 0 to rows - 1 do 403 + forward_97_row arr.(y) cols 404 + done; 405 + (* Transform columns *) 406 + let col = Array.make rows 0.0 in 407 + for x = 0 to cols - 1 do 408 + for y = 0 to rows - 1 do 409 + col.(y) <- arr.(y).(x) 410 + done; 411 + forward_97_row col rows; 412 + for y = 0 to rows - 1 do 413 + arr.(y).(x) <- col.(y) 414 + done 415 + done 416 + 417 + (** Inverse 2D CDF 9/7 DWT, one level. *) 418 + let inverse_dwt_97_2d_level (arr : float array array) rows cols = 419 + (* Inverse transform columns *) 420 + let col = Array.make rows 0.0 in 421 + for x = 0 to cols - 1 do 422 + for y = 0 to rows - 1 do 423 + col.(y) <- arr.(y).(x) 424 + done; 425 + inverse_97_row col rows; 426 + for y = 0 to rows - 1 do 427 + arr.(y).(x) <- col.(y) 428 + done 429 + done; 430 + (* Inverse transform rows *) 431 + for y = 0 to rows - 1 do 432 + inverse_97_row arr.(y) cols 433 + done 434 + 435 + (** Multi-level forward 2D CDF 9/7 DWT. Operates on a float 2D array. *) 436 + let forward_dwt_97_2d (arr : float array array) width height = 437 + let levels = dwt_levels width height in 438 + let sizes = subband_sizes width height levels in 439 + for i = 0 to levels - 1 do 440 + let w, h = sizes.(i) in 441 + forward_dwt_97_2d_level arr h w 442 + done; 443 + levels 444 + 445 + (** Multi-level inverse 2D CDF 9/7 DWT. *) 446 + let inverse_dwt_97_2d (arr : float array array) width height levels = 447 + let sizes = subband_sizes width height levels in 448 + for i = levels - 1 downto 0 do 449 + let w, h = sizes.(i) in 450 + inverse_dwt_97_2d_level arr h w 451 + done 452 + 453 + (** Convert integer image to float 2D array. *) 454 + let image_to_float_2d ~width ~height (data : bytes) = 455 + Array.init height (fun y -> 456 + Array.init width (fun x -> 457 + Float.of_int (Char.code (Bytes.get data ((y * width) + x))))) 458 + 459 + (** Quantize float coefficients to integers by rounding. *) 460 + let quantize_2d (arr : float array array) rows cols = 461 + let result = make_2d rows cols 0 in 462 + for y = 0 to rows - 1 do 463 + for x = 0 to cols - 1 do 464 + result.(y).(x) <- Float.to_int (Float.round arr.(y).(x)) 465 + done 466 + done; 467 + result 468 + 469 + (** Dequantize integer coefficients back to floats. *) 470 + let dequantize_2d (arr : int array array) rows cols = 471 + Array.init rows (fun y -> Array.init cols (fun x -> Float.of_int arr.(y).(x))) 472 + 258 473 (* ---- Coefficient encoder ---- *) 259 474 260 475 (** Coefficient encoder with zero-map and fixed-width coding. ··· 286 501 let rec go x acc = if x = 0 then acc else go (x lsr 1) (acc + 1) in 287 502 go n 0 288 503 289 - let encode_coefficients ~width ~height arr = 504 + let wavelet_id = function `Int_5_3 -> 0 | `Float_9_7 -> 1 505 + 506 + let wavelet_of_id = function 507 + | 0 -> `Int_5_3 508 + | 1 -> `Float_9_7 509 + | n -> invalid_arg (Printf.sprintf "Idc: unknown wavelet id %d" n) 510 + 511 + let encode_coefficients ~wavelet ~width ~height arr = 290 512 let bw = Bitstream.create (width * height) in 291 513 let max_mag = find_max_abs arr height width in 292 514 let mag_bits = bits_needed max_mag in ··· 294 516 Bitstream.write_bits bw 16 width; 295 517 Bitstream.write_bits bw 16 height; 296 518 Bitstream.write_bits bw 8 (dwt_levels width height); 519 + Bitstream.write_bits bw 8 (wavelet_id wavelet); 297 520 Bitstream.write_bits bw 16 max_mag; 298 521 (* Zero map: 1 bit per coefficient *) 299 522 for y = 0 to height - 1 do ··· 319 542 let width = Bitstream.read_bits br 16 in 320 543 let height = Bitstream.read_bits br 16 in 321 544 let levels = Bitstream.read_bits br 8 in 545 + let wavelet = wavelet_of_id (Bitstream.read_bits br 8) in 322 546 let max_mag = Bitstream.read_bits br 16 in 323 547 let mag_bits = bits_needed max_mag in 324 548 let arr = make_2d height width 0 in ··· 340 564 end 341 565 done 342 566 done; 343 - (width, height, levels, arr) 567 + (width, height, levels, wavelet, arr) 344 568 345 569 (* ---- Public API ---- *) 346 570 ··· 350 574 (Printf.sprintf "Idc.compress: expected %d bytes, got %d" (width * height) 351 575 (Bytes.length data)); 352 576 match wavelet with 353 - | `Float_9_7 -> failwith "Idc.compress: Float 9/7 wavelet not yet implemented" 354 577 | `Int_5_3 -> 355 578 let arr = image_to_2d ~width ~height data in 356 579 let _levels = forward_dwt_2d arr width height in 357 - encode_coefficients ~width ~height arr 580 + encode_coefficients ~wavelet ~width ~height arr 581 + | `Float_9_7 -> 582 + let farr = image_to_float_2d ~width ~height data in 583 + let _levels = forward_dwt_97_2d farr width height in 584 + let arr = quantize_2d farr height width in 585 + encode_coefficients ~wavelet ~width ~height arr 358 586 359 587 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 588 + let w, h, levels, wavelet, arr = decode_coefficients data in 589 + match wavelet with 590 + | `Int_5_3 -> 591 + inverse_dwt_2d arr w h levels; 592 + image_of_2d ~width:w ~height:h arr 593 + | `Float_9_7 -> 594 + let farr = dequantize_2d arr h w in 595 + inverse_dwt_97_2d farr w h levels; 596 + let result = Bytes.create (w * h) in 597 + for y = 0 to h - 1 do 598 + for x = 0 to w - 1 do 599 + let v = max 0 (min 255 (Float.to_int (Float.round farr.(y).(x)))) in 600 + Bytes.set result ((y * w) + x) (Char.chr v) 601 + done 602 + done; 603 + result

+87 -6

test/test.ml

··· 157 157 | _ -> Alcotest.fail "expected Invalid_argument for wrong data size" 158 158 | exception Invalid_argument _ -> () 159 159 160 - (* Float 9/7 wavelet is not yet implemented. *) 161 - let test_float_97_stub () = 160 + (* ---- Float 9/7 tests ---- *) 161 + 162 + (** Compute PSNR between two byte images. *) 163 + let psnr (a : bytes) (b : bytes) = 164 + let n = Bytes.length a in 165 + assert (n = Bytes.length b); 166 + let mse = ref 0.0 in 167 + for i = 0 to n - 1 do 168 + let d = 169 + Float.of_int (Char.code (Bytes.get a i)) 170 + -. Float.of_int (Char.code (Bytes.get b i)) 171 + in 172 + mse := !mse +. (d *. d) 173 + done; 174 + mse := !mse /. Float.of_int n; 175 + if !mse = 0.0 then infinity else 10.0 *. Float.log10 (255.0 *. 255.0 /. !mse) 176 + 177 + (* Float 9/7 roundtrip: compress then decompress, check PSNR > 40 dB. *) 178 + let test_float_97_roundtrip () = 179 + let width = 64 and height = 64 in 180 + let data = Bytes.create (width * height) in 181 + (* Smooth gradient with slight sine modulation *) 182 + for y = 0 to height - 1 do 183 + for x = 0 to width - 1 do 184 + let fx = Float.of_int x /. Float.of_int width in 185 + let fy = Float.of_int y /. Float.of_int height in 186 + let v = (fx +. fy) /. 2.0 in 187 + let v = v +. (0.1 *. Float.sin (fx *. 6.283)) in 188 + let v = max 0.0 (min 1.0 v) in 189 + Bytes.set data ((y * width) + x) (Char.chr (Float.to_int (v *. 255.0))) 190 + done 191 + done; 192 + let compressed = Idc.compress ~wavelet:`Float_9_7 ~width ~height data in 193 + let decompressed = Idc.decompress ~width ~height compressed in 194 + let p = psnr data decompressed in 195 + if p < 40.0 then 196 + Alcotest.failf "Float 9/7 roundtrip PSNR %.1f dB < 40 dB threshold" p 197 + 198 + (* Float 9/7 constant image should be exact (or near-exact). *) 199 + let test_float_97_constant () = 162 200 let width = 8 and height = 8 in 163 201 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 _ -> () 202 + let compressed = Idc.compress ~wavelet:`Float_9_7 ~width ~height data in 203 + let decompressed = Idc.decompress ~width ~height compressed in 204 + if data <> decompressed then 205 + Alcotest.fail "Float 9/7 constant image roundtrip not exact" 206 + 207 + (* Float 9/7 should compress a smooth non-linear image better than Int 5/3. 208 + 209 + The CDF 9/7 wavelet has more vanishing moments than 5/3 (4 vs 2 for the 210 + analysis filter). This means it annihilates higher-degree polynomials, 211 + producing smaller detail coefficients for smooth signals. For a cosine wave 212 + image, the 9/7 high-frequency coefficients round to zero more often than the 213 + 5/3 integer coefficients, giving a smaller compressed output. *) 214 + let test_float_97_better_compression () = 215 + let width = 64 and height = 64 in 216 + let data = Bytes.create (width * height) in 217 + (* Smooth cosine pattern -- not piecewise linear *) 218 + for y = 0 to height - 1 do 219 + for x = 0 to width - 1 do 220 + let fx = Float.of_int x /. Float.of_int width in 221 + let fy = Float.of_int y /. Float.of_int height in 222 + let v = 223 + 0.5 224 + +. (0.3 *. Float.cos (fx *. Float.pi *. 2.0)) 225 + +. (0.15 *. Float.cos (fy *. Float.pi *. 4.0)) 226 + in 227 + let v = max 0.0 (min 1.0 v) in 228 + Bytes.set data ((y * width) + x) (Char.chr (Float.to_int (v *. 255.0))) 229 + done 230 + done; 231 + let c53 = Idc.compress ~wavelet:`Int_5_3 ~width ~height data in 232 + let c97 = Idc.compress ~wavelet:`Float_9_7 ~width ~height data in 233 + let r53 = 234 + Float.of_int (Bytes.length c53) /. Float.of_int (Bytes.length data) 235 + in 236 + let r97 = 237 + Float.of_int (Bytes.length c97) /. Float.of_int (Bytes.length data) 238 + in 239 + if r97 >= r53 then 240 + Alcotest.failf 241 + "Float 9/7 compression ratio %.3f >= Int 5/3 ratio %.3f for smooth \ 242 + cosine image" 243 + r97 r53 167 244 168 245 let suite = 169 246 ( "idc", ··· 181 258 test_compression_ratio_large_blocky ); 182 259 ("non-power-of-2 dimensions", `Quick, test_non_power_of_2); 183 260 ("invalid input size", `Quick, test_invalid_size); 184 - ("float 9/7 stub", `Quick, test_float_97_stub); 261 + ("float 9/7 roundtrip PSNR", `Quick, test_float_97_roundtrip); 262 + ("float 9/7 constant exact", `Quick, test_float_97_constant); 263 + ( "float 9/7 better compression than 5/3", 264 + `Quick, 265 + test_float_97_better_compression ); 185 266 ] ) 186 267 187 268 let () = Alcotest.run "idc" [ suite ]

Configure Feed

Configure Feed