Histogram equalization for color images

Cargo.toml

@@ -27,10 +27,10 @@ directories = "2.0.2"
 filetime = "0.2.10"
 gcd = "2.0.1"
 hound = "3.4.0"
+lab = "0.8.1"
 line_drawing = "0.8.0"
 log = "0.4.8"
 image = { version = "0.23.4", features = ["png"] }
-imageproc = "0.21.0"
 reqwest = { version = "0.10.4", features = ["blocking"] }
 rustfft = "3.0.1"
 satellite = { git = "/~https://github.com/richinfante/satellite-rs", rev = "1f95726" }

src/gui/work.rs

@@ -459,7 +459,6 @@ pub fn process() {
             );
             callback(noaa_apt::process(
                 &mut context,
-                &settings,
                 &signal,
                 contrast_adjustment,
                 rotate,

src/imageext.rs

@@ -0,0 +1,139 @@
+/// Some extra utilities for working with images, that use or complement
+/// available functions from `image` crate
+
+use image::{RgbaImage, SubImage, Pixel, GenericImage, GenericImageView};
+use lab::Lab;
+use std::convert::TryInto;
+
+/// A set of per-channel histograms from an image with 8 bits per channel.
+pub struct ChannelHistogram {
+    /// Per-channel histograms.
+    pub channels: Vec<[u32; 256]>,
+}
+
+/// A set of per-channel cumulative histograms from an image with 8 bits per channel.
+pub struct CumulativeChannelHistogram {
+    /// Per-channel cumulative histograms.
+    pub channels: Vec<[u32; 256]>,
+}
+
+/// Equalize the histogram of the grayscale (but still Rgba image with
+/// R = G = B, A = 255), by equalizing the histogram of one of channels (R),
+/// and using that for all the other (G, B). Alpha channel is not modified.
+pub fn equalize_histogram_grayscale(sub_image: &mut SubImage<&mut RgbaImage>) {
+    // since it's a grayscale image (R = G = B, A = 255), use R channel histogram:
+    let hist = cumulative_histogram_rgba(sub_image).channels[0];
+    let total = hist[255] as f32;
+
+    for y in 0..sub_image.height() {
+        for x in 0..sub_image.width() {
+            let p = sub_image.get_pixel_mut(x, y);
+
+            // Each histogram has length 256 and RgbaImage has 8 bits per pixel
+            let fraction = hist[p.channels()[0] as usize] as f32 / total;
+
+            // apply f to channels r, g, b and apply g to alpha channel
+            p.apply_with_alpha(
+                // for R, G, B, use equalized values:
+                |_| (255. * fraction) as u8,
+                // for A, leave unmodified
+                |alpha| alpha
+            );
+        }
+    }
+}
+
+/// Equalize the histogram of the color subimage by converting Rgb -> Lab,
+/// equalizing the L (lightness) histogram, and converting back Lab -> Rgb.
+pub fn equalize_histogram_color(sub_image: &mut SubImage<&mut RgbaImage>) {
+    let mut lab_pixels: Vec<Lab> = rgb_to_lab(sub_image);
+
+    let lab_hist = cumulative_histogram_lab(&lab_pixels);
+    let total = lab_hist[100] as f32;
+
+    lab_pixels.iter_mut().for_each(|p: &mut Lab| {
+        // casting p.l from f32 to usize rounds towards 0
+        // l is in range [0..100] inclusive, lab_hist has lenght 101
+        let fraction = lab_hist[p.l as usize] as f32 / total;
+        p.l = 100. * fraction;
+    });
+    lab_to_rgb_mut(&lab_pixels, sub_image);
+}
+
+/// Returns a vector of Lab pixel values, alpha channel value is not used.
+fn rgb_to_lab(sub_image: &mut SubImage<&mut RgbaImage>) -> Vec<Lab> {
+    sub_image.pixels().map(|(_x, _y, p)| {
+        let rgb: [u8; 3] = p.channels()[..3].try_into().unwrap();
+        Lab::from_rgb(&rgb)
+    }).collect()
+}
+
+/// Converts Lab to Rgb and modifies the R, B, G values of pixels
+/// in the original subimage. The value of the alpha channel is unmodified.
+fn lab_to_rgb_mut(lab_pixels: &Vec<Lab>, sub_image: &mut SubImage<&mut RgbaImage>) {
+    let rgb_pixels: Vec<[u8; 3]> = lab_pixels.iter()
+        .map(|x: &Lab| x.to_rgb()).collect();
+
+    let height = sub_image.height();
+    let width = sub_image.width();
+
+    for y in 0..height {
+        for x in 0..width {
+            let p = sub_image.get_pixel_mut(x, y);
+            let [r, g, b] = rgb_pixels[(y * width + x) as usize];
+            let a = p.channels()[3]; // get original alpha channel
+            *p = Pixel::from_channels(r, g, b, a);
+        }
+    }
+}
+
+/// Calculates the cumulative histograms for each channel of the subimage.
+fn cumulative_histogram_rgba(
+    sub_image: &mut SubImage<&mut RgbaImage>
+) -> CumulativeChannelHistogram {
+    let mut hist = histogram_rgba(sub_image);
+    for c in 0..hist.channels.len() {
+        for i in 1..hist.channels[c].len() {
+            hist.channels[c][i] += hist.channels[c][i - 1];
+        }
+    }
+    CumulativeChannelHistogram {
+        channels: hist.channels,
+    }
+}
+
+/// Calculates the histograms for each channel of the subimage.
+fn histogram_rgba(sub_image: &mut SubImage<&mut RgbaImage>) -> ChannelHistogram {
+    let mut hist = vec![[0u32; 256]; 4];
+
+    sub_image.pixels().for_each(|(_x, _y, p)| {
+        for (i, c) in p.channels().iter().enumerate() {
+            hist[i][*c as usize] += 1;
+        }
+    });
+    ChannelHistogram { channels: hist }
+}
+
+/// Calculates the cumulative histogram using the L (lightness) channel.
+/// L values are in range [0..100] inclusive, so the resulting array
+/// has 101 elements: `[u32; 101]`
+fn cumulative_histogram_lab(lab_pixels: &Vec<Lab>) -> [u32; 101] {
+    let mut hist = histogram_lab(lab_pixels);
+    for i in 1..hist.len() {
+        hist[i] += hist[i - 1];
+    }
+    hist
+}
+
+/// Calculates the histogram using the L (lightness) channel.
+/// L values are in range [0..100] inclusive, so the resulting array
+/// has 101 elements: `[u32; 101]`.
+/// If the histogram for the other channels is needed in the future,
+/// consider defining a struct similar to `ChannelHistogram`.
+fn histogram_lab(lab_pixels: &Vec<Lab>) -> [u32; 101] {
+    let mut hist = [0u32; 101];
+    for p in lab_pixels {
+        hist[p.l as usize] += 1; // use L (lightness) channel
+    }
+    hist
+}

src/main.rs

@@ -19,6 +19,7 @@ mod context;
 mod decode;
 mod dsp;
 mod err;
+mod imageext;
 mod filters;
 mod frequency;
 mod geo;
@@ -88,7 +89,7 @@ fn inner_main() -> err::Result<()> {
 
             if !sync {
                 match contrast_adjustment {
-                    noaa_apt::Contrast::Telemetry | noaa_apt::Contrast::Histogram => 
+                    noaa_apt::Contrast::Telemetry | noaa_apt::Contrast::Histogram =>
                         warn!("Adjusting contrast without syncing, expect horrible results!"),
                     _ => ()
                 }
@@ -114,7 +115,6 @@ fn inner_main() -> err::Result<()> {
 
             let img = noaa_apt::process(
                 &mut context,
-                &settings,
                 &raw_data,
                 contrast_adjustment,
                 rotate,

src/noaa_apt.rs

@@ -17,7 +17,8 @@ use crate::map;
 use crate::misc;
 use crate::processing;
 use crate::telemetry;
-use crate::{config, wav};
+use crate::wav;
+use image::GrayImage;
 
 pub type Image = image::RgbaImage;
 
@@ -116,14 +117,13 @@ pub fn load(input_filename: &Path) -> err::Result<(Signal, Rate)> {
 
 pub fn process(
     context: &mut Context,
-    settings: &config::Settings,
     signal: &Signal,
     contrast_adjustment: Contrast,
     rotate: Rotate,
     color: Option<ColorSettings>,
     orbit: Option<OrbitSettings>,
 ) -> err::Result<Image> {
-    let (low, high) = match contrast_adjustment {
+    let (mut low, mut high) = match contrast_adjustment {
         Contrast::Telemetry => {
             context.status(0.1, "Adjusting contrast from telemetry".to_string());
 
@@ -150,28 +150,41 @@ pub fn process(
         }
     };
 
+    // for colorization & histogram equalization,
+    // always do 98% contrast adjust first, then colorize,
+    // then equalize histogram of color image if needed
+    if color.is_some() {
+        if let Contrast::Histogram = contrast_adjustment {
+            let (l, h) = misc::percent(&signal, 0.98)?;
+            low = l;
+            high = h;
+        }
+    }
+
     // --------------------
 
     context.status(0.3, "Generating image".to_string());
 
     let height = signal.len() as u32 / PX_PER_ROW;
 
-    use image::GrayImage;
-
-    let mut img: GrayImage = GrayImage::from_vec(PX_PER_ROW, height, map(&signal, low, high))
-        .ok_or_else(|| {
-            err::Error::Internal("Could not create image, wrong buffer length".to_string())
-        })?;
-
-    if let Contrast::Histogram = contrast_adjustment {
-        img = processing::histogram_equalization(&img)?;
-    }
+    // grayscale image obtained by mapping signal values to 0..255
+    // based on the selected contrast adjustment
+    let img: GrayImage = GrayImage::from_vec(
+        PX_PER_ROW, height, map_signal_u8(&signal, low, high)
+    ).ok_or_else(|| {
+        err::Error::Internal("Could not create image, wrong buffer length".to_string())
+    })?;
 
     let mut img: Image = image::DynamicImage::ImageLuma8(img).into_rgba(); // convert to RGBA
 
-    if let Some(color_settings) = color {
+    if let Some(color_settings) = &color {
         processing::false_color(&mut img, color_settings);
     }
+
+    if let Contrast::Histogram = contrast_adjustment {
+        processing::histogram_equalization(&mut img, color.is_some());
+    }
+
     // --------------------
 
     if let Some(orbit_settings) = orbit.clone() {
@@ -220,7 +233,7 @@ pub fn process(
 ///
 /// `low` becomes 0 and `high` becomes 255. Values are clamped to prevent `u8`
 /// overflow.
-fn map(signal: &Signal, low: f32, high: f32) -> Vec<u8> {
+fn map_signal_u8(signal: &Signal, low: f32, high: f32) -> Vec<u8> {
     let range = high - low;
     let raw_data: Vec<u8> = signal
         .iter()
@@ -252,6 +265,6 @@ mod tests {
         let low = 0. * 123.123 - 234.234;
         let high = 255. * 123.123 - 234.234;
 
-        assert_eq!(expected, map(&shifted_values, low, high));
+        assert_eq!(expected, map_signal_u8(&shifted_values, low, high));
     }
 }

src/processing.rs

@@ -1,13 +1,14 @@
 //! Image processing functions.
 
-use image::{GenericImageView, GenericImage, GrayImage, ImageBuffer, Rgba};
+use image::{GenericImage, Rgba, RgbaImage};
 use log::{warn, info};
 
 use crate::decode::{PX_PER_CHANNEL, PX_SYNC_FRAME, PX_SPACE_DATA, PX_CHANNEL_IMAGE_DATA};
 use crate::err;
+use crate::imageext;
 use crate::geo;
 use crate::misc;
-use crate::noaa_apt::{Image, ColorSettings, OrbitSettings, RefTime};
+use crate::noaa_apt::{ColorSettings, OrbitSettings, RefTime};
 
 
 /// Rotates the channels in place, keeping the sync bands and telemetry intact.
@@ -18,24 +19,24 @@ use crate::noaa_apt::{Image, ColorSettings, OrbitSettings, RefTime};
 /// Care is taken to leave lines from the A channel at the same height as the B
 /// channel. Otherwise there can be a vertical offset of one pixel between each
 /// channel.
-pub fn rotate(img: &mut Image) {
+pub fn rotate(img: &mut RgbaImage) {
     info!("Rotating image");
 
     // where the actual image data starts, past the sync frames and deep space band
-    let x_offset = PX_SYNC_FRAME + PX_SPACE_DATA; // !
+    let x_offset = PX_SYNC_FRAME + PX_SPACE_DATA;
 
     let mut channel_a = img.sub_image(
-        x_offset, 0, PX_CHANNEL_IMAGE_DATA, img.height() // !
+        x_offset, 0, PX_CHANNEL_IMAGE_DATA, img.height()
     );
     image::imageops::rotate180_in_place(&mut channel_a);
 
     let mut channel_b = img.sub_image(
-        x_offset + PX_PER_CHANNEL, 0, PX_CHANNEL_IMAGE_DATA, img.height() // !
+        x_offset + PX_PER_CHANNEL, 0, PX_CHANNEL_IMAGE_DATA, img.height()
     );
     image::imageops::rotate180_in_place(&mut channel_b);
 }
 
-/// Rotate image if the pass was south to north.
+/// Returns true if this was a south to north pass, and the image needs to be rotated.
 pub fn south_to_north_pass(orbit_settings: &OrbitSettings) -> err::Result<bool> {
 
     let tle = match &orbit_settings.custom_tle {
@@ -78,32 +79,31 @@ pub fn south_to_north_pass(orbit_settings: &OrbitSettings) -> err::Result<bool>
     return Ok(azimuth < PI / 4. || azimuth > 3. * PI / 4.);
 }
 
-/// Histogram equalization, for each channel separately.
-/// Works only on the grayscale image,
-/// needs to be done before the RGBA conversion.
-pub fn histogram_equalization(img: &GrayImage) -> err::Result<GrayImage> {
-    info!("Performing histogram equalization");
-
-    let mut output = GrayImage::new(img.width(), img.height());
-    let mut channel_a = img.view(0, 0, PX_PER_CHANNEL, img.height()).to_image();
-    let mut channel_b = img
-        .view(PX_PER_CHANNEL, 0, PX_PER_CHANNEL, img.height())
-        .to_image();
-
-    imageproc::contrast::equalize_histogram_mut(&mut channel_a);
-    imageproc::contrast::equalize_histogram_mut(&mut channel_b);
-
-    output.copy_from(&channel_a, 0, 0)?;
-    output.copy_from(&channel_b, PX_PER_CHANNEL, 0)?;
+/// Histogram equalization, in place, for each channel (A, B) separately.
+/// If `has_color=false`, it will treat the image as grayscale (R = G = B, A = 255).
+/// If `has_color=true`, it will convert image from Rgba to Lab, equalize the histogram
+/// for L (lightness) channel, convert back to Rgb and adjust image values accordingly.
+pub fn histogram_equalization(img: &mut RgbaImage, has_color: bool) {
+    info!("Performing histogram equalization, has color: {}", has_color);
+    let height = img.height();
+
+    let mut channel_a = img.sub_image(0, 0, PX_PER_CHANNEL, height);
+    if has_color {
+        imageext::equalize_histogram_color(&mut channel_a);
+    } else {
+        imageext::equalize_histogram_grayscale(&mut channel_a);
+    }
 
-    Ok(output)
+    let mut channel_b = img.sub_image(PX_PER_CHANNEL, 0, PX_PER_CHANNEL, height);
+    imageext::equalize_histogram_grayscale(&mut channel_b);
 }
 
+
 /// Attempts to produce a colored image from grayscale channel and IR data.
-/// Works best when contrast is set to "telemetry".
+/// Works best when contrast is set to "telemetry" or "98 percent".
 /// Needs a way to allow tweaking hardcoded values for water, land, ice
 /// and dirt detection, from the UI or command line.
-pub fn false_color(img: &mut ImageBuffer<Rgba<u8>, Vec<u8>>, color_settings: ColorSettings) {
+pub fn false_color(img: &mut RgbaImage, color_settings: &ColorSettings) {
     let water = color_settings.water_threshold;
     let vegetation = color_settings.vegetation_threshold;
     let clouds = color_settings.clouds_threshold;