[feat](trx-rs): remove NOAA APT decoder
Co-Authored-By: Claude Opus 4.6 <noreply@anthropic.com> Signed-off-by: Stan Grams <sjg@haxx.space>
This commit is contained in:
@@ -4,8 +4,7 @@
|
||||
|
||||
//! Shared PNG image encoding for weather satellite decoders.
|
||||
//!
|
||||
//! Both NOAA APT and Meteor-M LRPT decoders produce PNG output through
|
||||
//! this common module.
|
||||
//! The Meteor-M LRPT decoder produces PNG output through this module.
|
||||
|
||||
use std::io::Cursor;
|
||||
|
||||
|
||||
@@ -4,17 +4,12 @@
|
||||
|
||||
//! Weather satellite image decoders.
|
||||
//!
|
||||
//! This crate provides decoders for two weather satellite transmission formats:
|
||||
//!
|
||||
//! - **NOAA APT** ([`noaa`]): Automatic Picture Transmission from NOAA-15/18/19
|
||||
//! on 137 MHz using FM/AM subcarrier modulation at 4160 samples/sec.
|
||||
//!
|
||||
//! - **Meteor-M LRPT** ([`lrpt`]): Low Rate Picture Transmission from
|
||||
//! Meteor-M N2-3/N2-4 using QPSK modulation at 72 kbps with CCSDS framing.
|
||||
//! This crate provides a decoder for Meteor-M LRPT (Low Rate Picture
|
||||
//! Transmission) from Meteor-M N2-3/N2-4 using QPSK modulation at 72 kbps
|
||||
//! with CCSDS framing.
|
||||
|
||||
pub mod image_enc;
|
||||
pub mod lrpt;
|
||||
pub mod noaa;
|
||||
|
||||
/// Current time in milliseconds since UNIX epoch.
|
||||
pub(crate) fn now_ms() -> i64 {
|
||||
|
||||
@@ -1,353 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2026 Stan Grams <sjg@haxx.space>
|
||||
//
|
||||
// SPDX-License-Identifier: BSD-2-Clause
|
||||
|
||||
//! APT (Automatic Picture Transmission) demodulator and line decoder.
|
||||
//!
|
||||
//! Weather satellite APT signal chain:
|
||||
//! FM-demodulated audio → 2400 Hz AM subcarrier → envelope → 4160 Hz image
|
||||
//!
|
||||
//! Frame layout at 4160 Hz (2080 samples = 0.5 s per line, 2 lines/sec):
|
||||
//! [SyncA 39][SpaceA 47][ImageA 909][TelA 45][SyncB 39][SpaceB 47][ImageB 909][TelB 45]
|
||||
|
||||
use num_complex::Complex;
|
||||
use rustfft::FftPlanner;
|
||||
use std::sync::Arc;
|
||||
|
||||
pub const APT_RATE: u32 = 4160;
|
||||
pub const LINE_SAMPLES: usize = 2080;
|
||||
|
||||
// Line layout offsets (samples into a line at APT_RATE Hz)
|
||||
pub const SYNC_A_LEN: usize = 39;
|
||||
const SPACE_A_LEN: usize = 47;
|
||||
pub const IMAGE_A_LEN: usize = 909;
|
||||
const TEL_A_LEN: usize = 45;
|
||||
const SYNC_B_LEN: usize = 39;
|
||||
const SPACE_B_LEN: usize = 47;
|
||||
pub const IMAGE_B_LEN: usize = 909;
|
||||
|
||||
pub const IMAGE_A_OFFSET: usize = SYNC_A_LEN + SPACE_A_LEN; // 86
|
||||
pub const IMAGE_B_OFFSET: usize =
|
||||
IMAGE_A_OFFSET + IMAGE_A_LEN + TEL_A_LEN + SYNC_B_LEN + SPACE_B_LEN; // 1126
|
||||
|
||||
// FFT block size for Hilbert-based AM demodulation
|
||||
const BLOCK_SIZE: usize = 4096;
|
||||
|
||||
// Sync detection parameters
|
||||
const SYNC_THRESHOLD: f32 = 0.15;
|
||||
const SYNC_SEARCH_LOCKED: usize = 12; // ±samples around expected sync position when locked
|
||||
const MAX_BAD_SYNC_LINES: u32 = 8; // unlock after this many low-confidence lines
|
||||
|
||||
/// Telemetry block length (samples per channel).
|
||||
pub const TEL_LEN: usize = 45;
|
||||
|
||||
/// A decoded APT line: raw pixel arrays for both image channels plus telemetry.
|
||||
#[derive(Clone)]
|
||||
pub struct RawLine {
|
||||
pub pixels_a: Box<[u8; IMAGE_A_LEN]>,
|
||||
pub pixels_b: Box<[u8; IMAGE_B_LEN]>,
|
||||
/// Telemetry block A (45 samples, normalised to 0-255).
|
||||
pub tel_a: Box<[u8; TEL_LEN]>,
|
||||
/// Telemetry block B (45 samples, normalised to 0-255).
|
||||
pub tel_b: Box<[u8; TEL_LEN]>,
|
||||
pub line_no: u32,
|
||||
}
|
||||
|
||||
/// Sync A reference pattern at APT_RATE Hz.
|
||||
///
|
||||
/// 1040 Hz square wave (period = 4 samples): alternating pairs hi/lo.
|
||||
fn sync_a_ref() -> [f32; SYNC_A_LEN] {
|
||||
let mut p = [0.0f32; SYNC_A_LEN];
|
||||
for (i, v) in p.iter_mut().enumerate() {
|
||||
// 7 cycles of alternating pairs, rest is zero (end-of-sync blank)
|
||||
*v = if i < 28 && (i % 4) < 2 { 1.0 } else { -1.0 };
|
||||
}
|
||||
p
|
||||
}
|
||||
|
||||
/// Compute normalised cross-correlation of `buf[offset..]` with the sync A
|
||||
/// reference pattern. Returns a value approximately in `[-1.0, 1.0]`.
|
||||
fn sync_score(buf: &[f32], offset: usize) -> f32 {
|
||||
if offset + SYNC_A_LEN > buf.len() {
|
||||
return 0.0;
|
||||
}
|
||||
let ref_pat = sync_a_ref();
|
||||
let window = &buf[offset..offset + SYNC_A_LEN];
|
||||
let mean = window.iter().sum::<f32>() / SYNC_A_LEN as f32;
|
||||
let rms =
|
||||
(window.iter().map(|&x| (x - mean) * (x - mean)).sum::<f32>() / SYNC_A_LEN as f32).sqrt();
|
||||
if rms < 1e-7 {
|
||||
return 0.0;
|
||||
}
|
||||
window
|
||||
.iter()
|
||||
.zip(ref_pat.iter())
|
||||
.map(|(&s, &r)| (s - mean) * r)
|
||||
.sum::<f32>()
|
||||
/ (SYNC_A_LEN as f32 * rms)
|
||||
}
|
||||
|
||||
/// Find the offset in `buf[0..search_len]` with the highest sync A score.
|
||||
/// Returns `(offset, score)`.
|
||||
fn find_best_sync(buf: &[f32], search_len: usize) -> (usize, f32) {
|
||||
let limit = search_len.min(buf.len().saturating_sub(SYNC_A_LEN));
|
||||
let mut best_off = 0usize;
|
||||
let mut best_score = f32::NEG_INFINITY;
|
||||
for off in 0..=limit {
|
||||
let s = sync_score(buf, off);
|
||||
if s > best_score {
|
||||
best_score = s;
|
||||
best_off = off;
|
||||
}
|
||||
}
|
||||
(best_off, best_score)
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// AM demodulator (Hilbert-based via rustfft)
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
/// Converts PCM at `sample_rate` to APT envelope samples at `APT_RATE` Hz.
|
||||
///
|
||||
/// Uses an FFT-based Hilbert transform to obtain the analytic signal, then
|
||||
/// extracts the AM envelope. Processes input in non-overlapping blocks of
|
||||
/// `BLOCK_SIZE` samples.
|
||||
pub struct AptDemod {
|
||||
fft_fwd: Arc<dyn rustfft::Fft<f32>>,
|
||||
fft_inv: Arc<dyn rustfft::Fft<f32>>,
|
||||
/// Lower bin index of the 2400 Hz ± 1040 Hz bandpass filter.
|
||||
k_lo: usize,
|
||||
/// Upper bin index of the 2400 Hz ± 1040 Hz bandpass filter.
|
||||
k_hi: usize,
|
||||
/// Input sample accumulation buffer.
|
||||
in_buf: Vec<f32>,
|
||||
/// Fractional position into the next input block for the resampler.
|
||||
resamp_phase: f64,
|
||||
/// Input samples consumed per APT_RATE output sample.
|
||||
resamp_step: f64,
|
||||
/// Output envelope buffer at APT_RATE Hz.
|
||||
pub out: Vec<f32>,
|
||||
}
|
||||
|
||||
impl AptDemod {
|
||||
pub fn new(sample_rate: u32) -> Self {
|
||||
let mut planner = FftPlanner::new();
|
||||
let fft_fwd = planner.plan_fft_forward(BLOCK_SIZE);
|
||||
let fft_inv = planner.plan_fft_inverse(BLOCK_SIZE);
|
||||
|
||||
let fs = sample_rate as f64;
|
||||
// Bandpass around 2400 Hz carrier, ±1040 Hz (APT image bandwidth)
|
||||
let k_lo = ((1360.0 * BLOCK_SIZE as f64 / fs).floor() as usize).max(1);
|
||||
let k_hi = ((3440.0 * BLOCK_SIZE as f64 / fs).ceil() as usize).min(BLOCK_SIZE / 2);
|
||||
|
||||
Self {
|
||||
fft_fwd,
|
||||
fft_inv,
|
||||
k_lo,
|
||||
k_hi,
|
||||
in_buf: Vec::new(),
|
||||
resamp_phase: 0.0,
|
||||
resamp_step: sample_rate as f64 / APT_RATE as f64,
|
||||
out: Vec::new(),
|
||||
}
|
||||
}
|
||||
|
||||
/// Push raw PCM samples; envelope output accumulates in `self.out`.
|
||||
pub fn push(&mut self, samples: &[f32]) {
|
||||
self.in_buf.extend_from_slice(samples);
|
||||
while self.in_buf.len() >= BLOCK_SIZE {
|
||||
// Drain exactly BLOCK_SIZE samples into a stack array
|
||||
let block: Vec<f32> = self.in_buf.drain(..BLOCK_SIZE).collect();
|
||||
self.process_block(&block);
|
||||
}
|
||||
}
|
||||
|
||||
fn process_block(&mut self, block: &[f32]) {
|
||||
// Forward FFT
|
||||
let mut spectrum: Vec<Complex<f32>> = block.iter().map(|&s| Complex::new(s, 0.0)).collect();
|
||||
self.fft_fwd.process(&mut spectrum);
|
||||
|
||||
// Bandpass + analytic signal:
|
||||
// - Keep positive-freq band [k_lo, k_hi], doubled (for single-sideband power).
|
||||
// - Zero all other bins (negative freqs and out-of-band positives).
|
||||
// IFFT of the resulting one-sided spectrum ≈ analytic signal of the
|
||||
// bandpass-filtered input; its magnitude is the AM envelope.
|
||||
for (k, bin) in spectrum.iter_mut().enumerate() {
|
||||
if k >= self.k_lo && k <= self.k_hi {
|
||||
*bin *= 2.0;
|
||||
} else {
|
||||
*bin = Complex::ZERO;
|
||||
}
|
||||
}
|
||||
|
||||
// Inverse FFT → complex analytic signal
|
||||
self.fft_inv.process(&mut spectrum);
|
||||
let scale = 1.0_f32 / BLOCK_SIZE as f32;
|
||||
|
||||
// Magnitude = AM envelope; resample to APT_RATE via linear interpolation
|
||||
let n = BLOCK_SIZE as f64;
|
||||
while self.resamp_phase + 1.0 < n {
|
||||
let i = self.resamp_phase as usize;
|
||||
let frac = (self.resamp_phase - i as f64) as f32;
|
||||
let s0 = spectrum[i].norm() * scale;
|
||||
let s1 = spectrum[i + 1].norm() * scale;
|
||||
self.out.push(s0 + frac * (s1 - s0));
|
||||
self.resamp_phase += self.resamp_step;
|
||||
}
|
||||
self.resamp_phase -= n;
|
||||
if self.resamp_phase < 0.0 {
|
||||
self.resamp_phase = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
pub fn reset(&mut self) {
|
||||
self.in_buf.clear();
|
||||
self.out.clear();
|
||||
self.resamp_phase = 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
// ---------------------------------------------------------------------------
|
||||
// Sync tracker and line assembler
|
||||
// ---------------------------------------------------------------------------
|
||||
|
||||
/// Consumes resampled APT envelope samples (at `APT_RATE` Hz), detects
|
||||
/// sync A markers, and assembles decoded image lines.
|
||||
pub struct SyncTracker {
|
||||
buf: Vec<f32>,
|
||||
locked: bool,
|
||||
bad_sync_count: u32,
|
||||
line_no: u32,
|
||||
pub lines: Vec<RawLine>,
|
||||
}
|
||||
|
||||
impl Default for SyncTracker {
|
||||
fn default() -> Self {
|
||||
Self::new()
|
||||
}
|
||||
}
|
||||
|
||||
impl SyncTracker {
|
||||
pub fn new() -> Self {
|
||||
Self {
|
||||
buf: Vec::new(),
|
||||
locked: false,
|
||||
bad_sync_count: 0,
|
||||
line_no: 0,
|
||||
lines: Vec::new(),
|
||||
}
|
||||
}
|
||||
|
||||
/// Push envelope samples; fully assembled lines accumulate in `self.lines`.
|
||||
pub fn push(&mut self, samples: &[f32]) {
|
||||
self.buf.extend_from_slice(samples);
|
||||
self.drain();
|
||||
}
|
||||
|
||||
fn drain(&mut self) {
|
||||
loop {
|
||||
if !self.locked {
|
||||
// Need 2 × LINE_SAMPLES to have room for a full line after scan
|
||||
if self.buf.len() < 2 * LINE_SAMPLES {
|
||||
break;
|
||||
}
|
||||
let (pos, score) = find_best_sync(&self.buf, LINE_SAMPLES);
|
||||
if score >= SYNC_THRESHOLD {
|
||||
self.buf.drain(0..pos);
|
||||
self.locked = true;
|
||||
// Fall through to locked extraction below
|
||||
} else {
|
||||
// No convincing sync yet; discard half a line and retry
|
||||
self.buf.drain(0..LINE_SAMPLES / 2);
|
||||
continue;
|
||||
}
|
||||
}
|
||||
|
||||
// Locked mode: need LINE_SAMPLES + search window to be available
|
||||
if self.buf.len() < LINE_SAMPLES + SYNC_SEARCH_LOCKED {
|
||||
break;
|
||||
}
|
||||
|
||||
// Refine within a small window around the expected line start
|
||||
let (drift, score) = find_best_sync(&self.buf, SYNC_SEARCH_LOCKED);
|
||||
if score >= SYNC_THRESHOLD * 0.5 {
|
||||
// Accept refined position
|
||||
if drift > 0 {
|
||||
self.buf.drain(0..drift);
|
||||
}
|
||||
self.bad_sync_count = 0;
|
||||
} else {
|
||||
self.bad_sync_count += 1;
|
||||
if self.bad_sync_count >= MAX_BAD_SYNC_LINES {
|
||||
self.locked = false;
|
||||
self.bad_sync_count = 0;
|
||||
// Discard the current suspect region and restart search
|
||||
self.buf.drain(0..LINE_SAMPLES / 2);
|
||||
continue;
|
||||
}
|
||||
// Keep going at the expected position (don't drift on bad sync)
|
||||
}
|
||||
|
||||
if self.buf.len() < LINE_SAMPLES {
|
||||
break;
|
||||
}
|
||||
self.extract_line();
|
||||
}
|
||||
}
|
||||
|
||||
fn extract_line(&mut self) {
|
||||
let samples: Vec<f32> = self.buf.drain(0..LINE_SAMPLES).collect();
|
||||
|
||||
// Per-line normalisation: scale to [0, 255] using the 2nd–98th percentile
|
||||
// of this line's values to clip noise and hot pixels.
|
||||
let mut sorted: Vec<f32> = samples.clone();
|
||||
sorted.sort_unstable_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
|
||||
let p_lo = sorted[(sorted.len() * 2 / 100).max(0)];
|
||||
let p_hi = sorted[(sorted.len() * 98 / 100).min(sorted.len() - 1)];
|
||||
let range = (p_hi - p_lo).max(1e-6);
|
||||
|
||||
let norm = |v: f32| -> u8 { ((v - p_lo) / range * 255.0).round().clamp(0.0, 255.0) as u8 };
|
||||
|
||||
let mut pixels_a = Box::new([0u8; IMAGE_A_LEN]);
|
||||
for (i, p) in pixels_a.iter_mut().enumerate() {
|
||||
*p = norm(samples[IMAGE_A_OFFSET + i]);
|
||||
}
|
||||
let mut pixels_b = Box::new([0u8; IMAGE_B_LEN]);
|
||||
for (i, p) in pixels_b.iter_mut().enumerate() {
|
||||
*p = norm(samples[IMAGE_B_OFFSET + i]);
|
||||
}
|
||||
|
||||
// Extract telemetry blocks (adjacent to image data)
|
||||
let tel_a_offset = IMAGE_A_OFFSET + IMAGE_A_LEN; // right after image A
|
||||
let tel_b_offset = IMAGE_B_OFFSET + IMAGE_B_LEN; // right after image B
|
||||
let mut tel_a = Box::new([0u8; TEL_LEN]);
|
||||
for (i, p) in tel_a.iter_mut().enumerate() {
|
||||
if tel_a_offset + i < LINE_SAMPLES {
|
||||
*p = norm(samples[tel_a_offset + i]);
|
||||
}
|
||||
}
|
||||
let mut tel_b = Box::new([0u8; TEL_LEN]);
|
||||
for (i, p) in tel_b.iter_mut().enumerate() {
|
||||
if tel_b_offset + i < LINE_SAMPLES {
|
||||
*p = norm(samples[tel_b_offset + i]);
|
||||
}
|
||||
}
|
||||
|
||||
self.lines.push(RawLine {
|
||||
pixels_a,
|
||||
pixels_b,
|
||||
tel_a,
|
||||
tel_b,
|
||||
line_no: self.line_no,
|
||||
});
|
||||
self.line_no += 1;
|
||||
}
|
||||
|
||||
pub fn reset(&mut self) {
|
||||
self.buf.clear();
|
||||
self.locked = false;
|
||||
self.bad_sync_count = 0;
|
||||
self.line_no = 0;
|
||||
self.lines.clear();
|
||||
}
|
||||
}
|
||||
@@ -1,31 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2026 Stan Grams <sjg@haxx.space>
|
||||
//
|
||||
// SPDX-License-Identifier: BSD-2-Clause
|
||||
|
||||
//! APT image assembly.
|
||||
//!
|
||||
//! Standard output layout: channel A (visible / IR-A) on the left half and
|
||||
//! channel B (IR-B / IR) on the right half, stacked vertically by line number.
|
||||
|
||||
use super::apt::{RawLine, IMAGE_A_LEN, IMAGE_B_LEN};
|
||||
|
||||
/// Assemble decoded APT lines into a PNG image.
|
||||
///
|
||||
/// Returns the PNG bytes, or `None` if `lines` is empty or encoding fails.
|
||||
/// Width = `IMAGE_A_LEN + IMAGE_B_LEN` (1818 px), height = number of lines.
|
||||
pub fn encode_png(lines: &[RawLine]) -> Option<Vec<u8>> {
|
||||
if lines.is_empty() {
|
||||
return None;
|
||||
}
|
||||
|
||||
let width = (IMAGE_A_LEN + IMAGE_B_LEN) as u32;
|
||||
let height = lines.len() as u32;
|
||||
let mut pixels: Vec<u8> = Vec::with_capacity((width * height) as usize);
|
||||
|
||||
for line in lines {
|
||||
pixels.extend_from_slice(line.pixels_a.as_ref());
|
||||
pixels.extend_from_slice(line.pixels_b.as_ref());
|
||||
}
|
||||
|
||||
crate::image_enc::encode_grayscale_png(width, height, pixels)
|
||||
}
|
||||
@@ -1,166 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2026 Stan Grams <sjg@haxx.space>
|
||||
//
|
||||
// SPDX-License-Identifier: BSD-2-Clause
|
||||
|
||||
//! NOAA APT (Automatic Picture Transmission) weather satellite image decoder.
|
||||
//!
|
||||
//! Decodes the APT format broadcast by NOAA-15 (137.620 MHz),
|
||||
//! NOAA-18 (137.9125 MHz) and NOAA-19 (137.100 MHz).
|
||||
//!
|
||||
//! # Signal chain
|
||||
//!
|
||||
//! The input is FM-demodulated audio containing a 2400 Hz AM subcarrier.
|
||||
//! The decoder:
|
||||
//! 1. Extracts the AM envelope via a FFT-based Hilbert transform (rustfft).
|
||||
//! 2. Resamples to 4160 Hz (the APT image sample rate).
|
||||
//! 3. Detects line sync markers (1040 Hz alternating pattern).
|
||||
//! 4. Assembles image lines (2080 samples each) and extracts both channels.
|
||||
//!
|
||||
//! Call [`AptDecoder::process_samples`] with each audio batch, then
|
||||
//! [`AptDecoder::finalize`] when the pass ends to obtain JPEG bytes.
|
||||
|
||||
pub mod apt;
|
||||
mod image_enc;
|
||||
pub mod telemetry;
|
||||
|
||||
use apt::{AptDemod, SyncTracker};
|
||||
use telemetry::{Satellite, SensorChannel};
|
||||
|
||||
/// Completed APT image returned by [`AptDecoder::finalize`].
|
||||
pub struct AptImage {
|
||||
/// PNG-encoded image bytes.
|
||||
pub png: Vec<u8>,
|
||||
/// Number of decoded image lines.
|
||||
pub line_count: u32,
|
||||
/// Millisecond timestamp when the first line was decoded.
|
||||
pub first_line_ms: i64,
|
||||
/// Identified satellite, if telemetry was decodable.
|
||||
pub satellite: Satellite,
|
||||
/// Detected sensor channel for sub-channel A.
|
||||
pub sensor_a: SensorChannel,
|
||||
/// Detected sensor channel for sub-channel B.
|
||||
pub sensor_b: SensorChannel,
|
||||
}
|
||||
|
||||
/// Top-level NOAA APT decoder.
|
||||
///
|
||||
/// Feed audio samples with [`process_samples`] and call [`finalize`] at
|
||||
/// pass end to retrieve the assembled JPEG.
|
||||
pub struct AptDecoder {
|
||||
demod: AptDemod,
|
||||
sync: SyncTracker,
|
||||
first_line_ms: Option<i64>,
|
||||
}
|
||||
|
||||
impl AptDecoder {
|
||||
pub fn new(sample_rate: u32) -> Self {
|
||||
Self {
|
||||
demod: AptDemod::new(sample_rate),
|
||||
sync: SyncTracker::new(),
|
||||
first_line_ms: None,
|
||||
}
|
||||
}
|
||||
|
||||
/// Process a batch of PCM samples (float32, mono or will be treated as-is).
|
||||
///
|
||||
/// Returns the number of new lines decoded in this batch.
|
||||
pub fn process_samples(&mut self, samples: &[f32]) -> u32 {
|
||||
self.demod.push(samples);
|
||||
|
||||
let before = self.sync.lines.len() as u32;
|
||||
|
||||
// Move accumulated envelope output into the sync tracker
|
||||
if !self.demod.out.is_empty() {
|
||||
let envelope = std::mem::take(&mut self.demod.out);
|
||||
self.sync.push(&envelope);
|
||||
}
|
||||
|
||||
let after = self.sync.lines.len() as u32;
|
||||
let new_lines = after - before;
|
||||
|
||||
if new_lines > 0 && self.first_line_ms.is_none() {
|
||||
self.first_line_ms = Some(crate::now_ms());
|
||||
}
|
||||
|
||||
new_lines
|
||||
}
|
||||
|
||||
/// Total number of lines decoded so far.
|
||||
pub fn line_count(&self) -> u32 {
|
||||
self.sync.lines.len() as u32
|
||||
}
|
||||
|
||||
/// Encode all accumulated lines as a JPEG image and return the result.
|
||||
///
|
||||
/// Performs telemetry extraction, radiometric calibration (when enough
|
||||
/// lines are available for a full 128-line telemetry frame), and
|
||||
/// histogram equalisation before JPEG encoding.
|
||||
///
|
||||
/// Returns `None` if no lines have been decoded yet.
|
||||
/// Does **not** reset the decoder; call [`reset`] afterwards if needed.
|
||||
pub fn finalize(&self) -> Option<AptImage> {
|
||||
if self.sync.lines.is_empty() {
|
||||
return None;
|
||||
}
|
||||
|
||||
// Extract telemetry for calibration and satellite identification
|
||||
let tel = telemetry::extract_telemetry(&self.sync.lines);
|
||||
|
||||
// Clone lines so we can apply calibration without mutating decoder state
|
||||
let mut lines = self.sync.lines.clone();
|
||||
|
||||
let (satellite, sensor_a, sensor_b) = if let Some(ref tf) = tel {
|
||||
// Apply radiometric calibration using telemetry wedge LUTs
|
||||
for line in &mut lines {
|
||||
telemetry::calibrate_line_a(&mut line.pixels_a, &tf.cal_lut_a);
|
||||
telemetry::calibrate_line_b(&mut line.pixels_b, &tf.cal_lut_b);
|
||||
}
|
||||
(tf.satellite, tf.sensor_a, tf.sensor_b)
|
||||
} else {
|
||||
(
|
||||
Satellite::Unknown,
|
||||
SensorChannel::Unknown,
|
||||
SensorChannel::Unknown,
|
||||
)
|
||||
};
|
||||
|
||||
// Apply histogram equalisation per-channel for contrast enhancement
|
||||
let mut all_a: Vec<u8> = lines
|
||||
.iter()
|
||||
.flat_map(|l| l.pixels_a.iter().copied())
|
||||
.collect();
|
||||
let mut all_b: Vec<u8> = lines
|
||||
.iter()
|
||||
.flat_map(|l| l.pixels_b.iter().copied())
|
||||
.collect();
|
||||
telemetry::histogram_equalize(&mut all_a);
|
||||
telemetry::histogram_equalize(&mut all_b);
|
||||
|
||||
// Write equalised pixels back
|
||||
let width_a = apt::IMAGE_A_LEN;
|
||||
let width_b = apt::IMAGE_B_LEN;
|
||||
for (i, line) in lines.iter_mut().enumerate() {
|
||||
line.pixels_a
|
||||
.copy_from_slice(&all_a[i * width_a..(i + 1) * width_a]);
|
||||
line.pixels_b
|
||||
.copy_from_slice(&all_b[i * width_b..(i + 1) * width_b]);
|
||||
}
|
||||
|
||||
let png = image_enc::encode_png(&lines)?;
|
||||
Some(AptImage {
|
||||
png,
|
||||
line_count: lines.len() as u32,
|
||||
first_line_ms: self.first_line_ms.unwrap_or_else(crate::now_ms),
|
||||
satellite,
|
||||
sensor_a,
|
||||
sensor_b,
|
||||
})
|
||||
}
|
||||
|
||||
/// Clear all state; ready to decode a fresh pass.
|
||||
pub fn reset(&mut self) {
|
||||
self.demod.reset();
|
||||
self.sync.reset();
|
||||
self.first_line_ms = None;
|
||||
}
|
||||
}
|
||||
@@ -1,396 +0,0 @@
|
||||
// SPDX-FileCopyrightText: 2026 Stan Grams <sjg@haxx.space>
|
||||
//
|
||||
// SPDX-License-Identifier: BSD-2-Clause
|
||||
|
||||
//! APT telemetry frame parsing, satellite identification, and channel detection.
|
||||
//!
|
||||
//! Each APT line contains two 45-sample telemetry blocks (one per channel).
|
||||
//! The telemetry frame repeats every 128 lines and contains 16 wedges of
|
||||
//! 8 lines each. Wedges 1-8 carry calibration reference levels, wedge 9
|
||||
//! carries the channel ID, and wedges 10-15 carry thermal calibration data.
|
||||
//! Wedge 16 is the "zero modulation" reference (black body equivalent).
|
||||
|
||||
use super::apt::{RawLine, IMAGE_A_LEN, IMAGE_B_LEN};
|
||||
|
||||
/// Lines per telemetry frame (128 lines = 16 wedges x 8 lines each).
|
||||
pub const FRAME_LINES: usize = 128;
|
||||
|
||||
/// Lines per wedge.
|
||||
pub const WEDGE_LINES: usize = 8;
|
||||
|
||||
/// Number of wedges in a telemetry frame.
|
||||
pub const NUM_WEDGES: usize = 16;
|
||||
|
||||
/// The 8 calibration step values defined by the APT spec (wedges 1-8).
|
||||
/// These represent known modulation levels from 1/8 to 8/8 of full scale.
|
||||
pub const WEDGE_STEPS: [f32; 8] = [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0];
|
||||
|
||||
/// NOAA AVHRR sensor channel assignments.
|
||||
///
|
||||
/// The NOAA APT format transmits two channels simultaneously. Which sensors
|
||||
/// are mapped to channel A and B depends on the satellite and illumination.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub enum SensorChannel {
|
||||
/// Channel 1: Visible (0.58 - 0.68 um)
|
||||
Visible1,
|
||||
/// Channel 2: Near-IR (0.725 - 1.0 um)
|
||||
NearIr2,
|
||||
/// Channel 3A: Near-IR (1.58 - 1.64 um) — daytime only on NOAA-15/18/19
|
||||
NearIr3A,
|
||||
/// Channel 3B: Mid-IR thermal (3.55 - 3.93 um)
|
||||
MidIr3B,
|
||||
/// Channel 4: Thermal IR (10.30 - 11.30 um)
|
||||
ThermalIr4,
|
||||
/// Channel 5: Thermal IR (11.50 - 12.50 um) — not on NOAA-15 APT
|
||||
ThermalIr5,
|
||||
/// Unknown / could not be determined.
|
||||
Unknown,
|
||||
}
|
||||
|
||||
impl std::fmt::Display for SensorChannel {
|
||||
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
|
||||
match self {
|
||||
SensorChannel::Visible1 => write!(f, "1-VIS"),
|
||||
SensorChannel::NearIr2 => write!(f, "2-NIR"),
|
||||
SensorChannel::NearIr3A => write!(f, "3A-NIR"),
|
||||
SensorChannel::MidIr3B => write!(f, "3B-MIR"),
|
||||
SensorChannel::ThermalIr4 => write!(f, "4-TIR"),
|
||||
SensorChannel::ThermalIr5 => write!(f, "5-TIR"),
|
||||
SensorChannel::Unknown => write!(f, "unknown"),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Identified NOAA satellite.
|
||||
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
|
||||
pub enum Satellite {
|
||||
Noaa15,
|
||||
Noaa18,
|
||||
Noaa19,
|
||||
Unknown,
|
||||
}
|
||||
|
||||
impl std::fmt::Display for Satellite {
|
||||
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
|
||||
match self {
|
||||
Satellite::Noaa15 => write!(f, "NOAA-15"),
|
||||
Satellite::Noaa18 => write!(f, "NOAA-18"),
|
||||
Satellite::Noaa19 => write!(f, "NOAA-19"),
|
||||
Satellite::Unknown => write!(f, "Unknown"),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Wedge 9 channel-ID values for each satellite.
|
||||
///
|
||||
/// The channel ID wedge has a distinctive grey level that encodes which
|
||||
/// AVHRR sensor channel is being transmitted on that APT sub-channel.
|
||||
/// Values are approximate normalised levels (0.0 - 1.0).
|
||||
///
|
||||
/// Reference: NOAA KLM User's Guide, Section 4.2 (APT format).
|
||||
///
|
||||
/// Channel A mapping:
|
||||
/// Wedge 9 ≈ step 1 (1/8) → channel 1 (VIS)
|
||||
/// Wedge 9 ≈ step 2 (2/8) → channel 2 (NIR)
|
||||
/// Wedge 9 ≈ step 3 (3/8) → channel 3A (NIR, daytime)
|
||||
///
|
||||
/// Channel B mapping:
|
||||
/// Wedge 9 ≈ step 4 (4/8) → channel 3B (MIR)
|
||||
/// Wedge 9 ≈ step 5 (5/8) → channel 4 (TIR)
|
||||
/// Wedge 9 ≈ step 6 (6/8) → channel 5 (TIR)
|
||||
fn wedge9_to_sensor(normalised: f32) -> SensorChannel {
|
||||
// Map to nearest step (1/8 increments)
|
||||
let step = (normalised * 8.0).round() as u8;
|
||||
match step {
|
||||
1 => SensorChannel::Visible1,
|
||||
2 => SensorChannel::NearIr2,
|
||||
3 => SensorChannel::NearIr3A,
|
||||
4 => SensorChannel::MidIr3B,
|
||||
5 => SensorChannel::ThermalIr4,
|
||||
6 => SensorChannel::ThermalIr5,
|
||||
_ => SensorChannel::Unknown,
|
||||
}
|
||||
}
|
||||
|
||||
/// Extracted telemetry data from one complete 128-line frame.
|
||||
#[derive(Debug, Clone)]
|
||||
pub struct TelemetryFrame {
|
||||
/// Mean pixel value for each of the 16 wedges (normalised 0.0 - 1.0).
|
||||
pub wedge_means_a: [f32; NUM_WEDGES],
|
||||
pub wedge_means_b: [f32; NUM_WEDGES],
|
||||
/// Detected sensor channel for sub-channel A.
|
||||
pub sensor_a: SensorChannel,
|
||||
/// Detected sensor channel for sub-channel B.
|
||||
pub sensor_b: SensorChannel,
|
||||
/// Calibration mapping: maps raw pixel [0,255] → calibrated [0.0, 1.0]
|
||||
/// using wedges 1-8 as known reference levels.
|
||||
pub cal_lut_a: [u8; 256],
|
||||
pub cal_lut_b: [u8; 256],
|
||||
/// Identified satellite (from channel pairing heuristics).
|
||||
pub satellite: Satellite,
|
||||
}
|
||||
|
||||
/// Extract telemetry from raw lines, requiring at least one full 128-line frame.
|
||||
///
|
||||
/// Picks the best complete frame (highest overall signal quality) and parses
|
||||
/// wedge values from the telemetry blocks.
|
||||
pub fn extract_telemetry(lines: &[RawLine]) -> Option<TelemetryFrame> {
|
||||
if lines.len() < FRAME_LINES {
|
||||
return None;
|
||||
}
|
||||
|
||||
// Use the middle complete frame for best quality (avoids pass start/end noise)
|
||||
let num_frames = lines.len() / FRAME_LINES;
|
||||
let frame_idx = num_frames / 2;
|
||||
let frame_start = frame_idx * FRAME_LINES;
|
||||
let frame = &lines[frame_start..frame_start + FRAME_LINES];
|
||||
|
||||
// Extract wedge means from telemetry blocks.
|
||||
// Each wedge spans 8 lines; we average the telemetry samples across those lines.
|
||||
let mut wedge_means_a = [0.0f32; NUM_WEDGES];
|
||||
let mut wedge_means_b = [0.0f32; NUM_WEDGES];
|
||||
|
||||
for wedge_idx in 0..NUM_WEDGES {
|
||||
let line_start = wedge_idx * WEDGE_LINES;
|
||||
let mut sum_a = 0.0f32;
|
||||
let mut sum_b = 0.0f32;
|
||||
let mut count = 0u32;
|
||||
|
||||
for line_offset in 0..WEDGE_LINES {
|
||||
let line = &frame[line_start + line_offset];
|
||||
for &v in line.tel_a.as_ref() {
|
||||
sum_a += v as f32;
|
||||
count += 1;
|
||||
}
|
||||
for &v in line.tel_b.as_ref() {
|
||||
sum_b += v as f32;
|
||||
}
|
||||
}
|
||||
|
||||
if count > 0 {
|
||||
wedge_means_a[wedge_idx] = sum_a / count as f32 / 255.0;
|
||||
wedge_means_b[wedge_idx] = sum_b / count as f32 / 255.0;
|
||||
}
|
||||
}
|
||||
|
||||
// Detect sensor channels from wedge 9 (index 8)
|
||||
let sensor_a = wedge9_to_sensor(wedge_means_a[8]);
|
||||
let sensor_b = wedge9_to_sensor(wedge_means_b[8]);
|
||||
|
||||
// Build calibration LUTs from wedges 1-8
|
||||
let cal_lut_a = build_calibration_lut(&wedge_means_a);
|
||||
let cal_lut_b = build_calibration_lut(&wedge_means_b);
|
||||
|
||||
// Identify satellite from channel pairing
|
||||
let satellite = identify_satellite(sensor_a, sensor_b);
|
||||
|
||||
Some(TelemetryFrame {
|
||||
wedge_means_a,
|
||||
wedge_means_b,
|
||||
sensor_a,
|
||||
sensor_b,
|
||||
cal_lut_a,
|
||||
cal_lut_b,
|
||||
satellite,
|
||||
})
|
||||
}
|
||||
|
||||
/// Build a 256-entry calibration look-up table from wedge means.
|
||||
///
|
||||
/// Wedges 1-8 (indices 0-7) represent known reference levels at 1/8 to 8/8.
|
||||
/// We fit a piecewise linear mapping from observed pixel values to calibrated
|
||||
/// output levels, producing a corrected 0-255 output.
|
||||
fn build_calibration_lut(wedge_means: &[f32; NUM_WEDGES]) -> [u8; 256] {
|
||||
let mut lut = [0u8; 256];
|
||||
|
||||
// Collect (observed_pixel_value, target_normalised) pairs from wedges 1-8
|
||||
let mut pairs: Vec<(f32, f32)> = Vec::with_capacity(8);
|
||||
for i in 0..8 {
|
||||
let observed = wedge_means[i] * 255.0;
|
||||
let target = WEDGE_STEPS[i];
|
||||
pairs.push((observed, target));
|
||||
}
|
||||
|
||||
// Sort by observed value
|
||||
pairs.sort_by(|a, b| a.0.partial_cmp(&b.0).unwrap_or(std::cmp::Ordering::Equal));
|
||||
|
||||
// Deduplicate (if two wedges map to nearly the same observed value)
|
||||
pairs.dedup_by(|a, b| (a.0 - b.0).abs() < 0.5);
|
||||
|
||||
if pairs.len() < 2 {
|
||||
// Not enough calibration data — return identity
|
||||
for (i, v) in lut.iter_mut().enumerate() {
|
||||
*v = i as u8;
|
||||
}
|
||||
return lut;
|
||||
}
|
||||
|
||||
// Piecewise linear interpolation
|
||||
for (i, entry) in lut.iter_mut().enumerate() {
|
||||
let x = i as f32;
|
||||
let calibrated = if x <= pairs[0].0 {
|
||||
pairs[0].1
|
||||
} else if x >= pairs[pairs.len() - 1].0 {
|
||||
pairs[pairs.len() - 1].1
|
||||
} else {
|
||||
let mut cal = pairs[0].1;
|
||||
for w in pairs.windows(2) {
|
||||
if x >= w[0].0 && x <= w[1].0 {
|
||||
let t = (x - w[0].0) / (w[1].0 - w[0].0).max(1e-6);
|
||||
cal = w[0].1 + t * (w[1].1 - w[0].1);
|
||||
break;
|
||||
}
|
||||
}
|
||||
cal
|
||||
};
|
||||
*entry = (calibrated * 255.0).round().clamp(0.0, 255.0) as u8;
|
||||
}
|
||||
|
||||
lut
|
||||
}
|
||||
|
||||
/// Identify the satellite based on channel pairing heuristics.
|
||||
///
|
||||
/// Typical APT channel pairings:
|
||||
/// - NOAA-15: Ch A = 2 (NIR), Ch B = 4 (TIR) daytime;
|
||||
/// Ch A = 3A (NIR), Ch B = 4 (TIR) alternate daytime
|
||||
/// - NOAA-18: Ch A = 1 (VIS), Ch B = 4 (TIR) daytime;
|
||||
/// Ch A = 3A (NIR), Ch B = 4 (TIR) alternate
|
||||
/// - NOAA-19: Ch A = 2 (NIR), Ch B = 4 (TIR) daytime
|
||||
///
|
||||
/// Night passes typically transmit Ch 3B or Ch 4 on channel A.
|
||||
fn identify_satellite(sensor_a: SensorChannel, sensor_b: SensorChannel) -> Satellite {
|
||||
match (sensor_a, sensor_b) {
|
||||
// NOAA-18 typically sends VIS ch1 on A
|
||||
(SensorChannel::Visible1, SensorChannel::ThermalIr4) => Satellite::Noaa18,
|
||||
// NOAA-15 and NOAA-19 both send NIR ch2 on A; distinguish by B channel
|
||||
(SensorChannel::NearIr2, SensorChannel::ThermalIr4) => {
|
||||
// Both NOAA-15 and NOAA-19 use this pairing; cannot easily distinguish
|
||||
// without orbital data. Default to NOAA-19 (most common active).
|
||||
Satellite::Noaa19
|
||||
}
|
||||
(SensorChannel::NearIr3A, SensorChannel::ThermalIr4) => Satellite::Noaa15,
|
||||
(SensorChannel::NearIr2, SensorChannel::ThermalIr5) => Satellite::Noaa19,
|
||||
_ => Satellite::Unknown,
|
||||
}
|
||||
}
|
||||
|
||||
/// Apply calibration LUT to a line's pixel data (in-place).
|
||||
pub fn calibrate_line_a(pixels: &mut [u8; IMAGE_A_LEN], lut: &[u8; 256]) {
|
||||
for p in pixels.iter_mut() {
|
||||
*p = lut[*p as usize];
|
||||
}
|
||||
}
|
||||
|
||||
/// Apply calibration LUT to a line's pixel data (in-place).
|
||||
pub fn calibrate_line_b(pixels: &mut [u8; IMAGE_B_LEN], lut: &[u8; 256]) {
|
||||
for p in pixels.iter_mut() {
|
||||
*p = lut[*p as usize];
|
||||
}
|
||||
}
|
||||
|
||||
/// Apply histogram equalisation to an image channel for contrast enhancement.
|
||||
pub fn histogram_equalize(pixels: &mut [u8]) {
|
||||
if pixels.is_empty() {
|
||||
return;
|
||||
}
|
||||
|
||||
// Build histogram
|
||||
let mut hist = [0u32; 256];
|
||||
for &p in pixels.iter() {
|
||||
hist[p as usize] += 1;
|
||||
}
|
||||
|
||||
// Compute CDF
|
||||
let mut cdf = [0u32; 256];
|
||||
cdf[0] = hist[0];
|
||||
for i in 1..256 {
|
||||
cdf[i] = cdf[i - 1] + hist[i];
|
||||
}
|
||||
|
||||
// Find minimum non-zero CDF value
|
||||
let cdf_min = cdf.iter().copied().find(|&v| v > 0).unwrap_or(0);
|
||||
let total = pixels.len() as u32;
|
||||
let denom = (total - cdf_min).max(1);
|
||||
|
||||
// Build equalisation LUT
|
||||
let mut lut = [0u8; 256];
|
||||
for i in 0..256 {
|
||||
lut[i] = ((cdf[i].saturating_sub(cdf_min) as f64 / denom as f64) * 255.0).round() as u8;
|
||||
}
|
||||
|
||||
// Apply
|
||||
for p in pixels.iter_mut() {
|
||||
*p = lut[*p as usize];
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::*;
|
||||
|
||||
#[test]
|
||||
fn test_wedge9_to_sensor() {
|
||||
assert_eq!(wedge9_to_sensor(0.125), SensorChannel::Visible1);
|
||||
assert_eq!(wedge9_to_sensor(0.25), SensorChannel::NearIr2);
|
||||
assert_eq!(wedge9_to_sensor(0.375), SensorChannel::NearIr3A);
|
||||
assert_eq!(wedge9_to_sensor(0.5), SensorChannel::MidIr3B);
|
||||
assert_eq!(wedge9_to_sensor(0.625), SensorChannel::ThermalIr4);
|
||||
assert_eq!(wedge9_to_sensor(0.75), SensorChannel::ThermalIr5);
|
||||
assert_eq!(wedge9_to_sensor(0.0), SensorChannel::Unknown);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_identify_satellite() {
|
||||
assert_eq!(
|
||||
identify_satellite(SensorChannel::Visible1, SensorChannel::ThermalIr4),
|
||||
Satellite::Noaa18
|
||||
);
|
||||
assert_eq!(
|
||||
identify_satellite(SensorChannel::NearIr2, SensorChannel::ThermalIr4),
|
||||
Satellite::Noaa19
|
||||
);
|
||||
assert_eq!(
|
||||
identify_satellite(SensorChannel::NearIr3A, SensorChannel::ThermalIr4),
|
||||
Satellite::Noaa15
|
||||
);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_calibration_lut_identity_on_insufficient_data() {
|
||||
let mut means = [0.0f32; NUM_WEDGES];
|
||||
// All zeros → insufficient data → identity LUT
|
||||
let lut = build_calibration_lut(&means);
|
||||
for i in 0..256 {
|
||||
assert_eq!(lut[i], i as u8);
|
||||
}
|
||||
|
||||
// One non-zero wedge still insufficient (need ≥ 2 distinct)
|
||||
means[0] = 0.5;
|
||||
let lut = build_calibration_lut(&means);
|
||||
// Still degenerate
|
||||
assert!(lut[0] == lut[0]); // trivially true, but confirms no panic
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_histogram_equalize_uniform() {
|
||||
// Uniform distribution should remain roughly unchanged
|
||||
let mut pixels: Vec<u8> = (0..=255).collect();
|
||||
histogram_equalize(&mut pixels);
|
||||
// After equalization, values should span full range
|
||||
assert_eq!(*pixels.first().unwrap(), 0);
|
||||
assert_eq!(*pixels.last().unwrap(), 255);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_sensor_channel_display() {
|
||||
assert_eq!(format!("{}", SensorChannel::Visible1), "1-VIS");
|
||||
assert_eq!(format!("{}", SensorChannel::ThermalIr4), "4-TIR");
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn test_satellite_display() {
|
||||
assert_eq!(format!("{}", Satellite::Noaa15), "NOAA-15");
|
||||
assert_eq!(format!("{}", Satellite::Noaa19), "NOAA-19");
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user