-
-
Notifications
You must be signed in to change notification settings - Fork 11
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Added working otsu and mean thresholding #15
Merged
Merged
Changes from all commits
Commits
Show all changes
7 commits
Select commit
Hold shift + click to select a range
65b4a1e
Merge branch 'develop'
xd009642 79aba53
Merge branch 'develop'
xd009642 1043954
Added working otsu and mean thresholding
jmetz 52a0557
Added working otsu and mean thresholding
jmetz 591d984
Merge branch 'otsu' of github.com:jmetz/ndarray-vision into otsu
jmetz bc36007
Addressed comments on PR
jmetz b96faf3
Added changelog entry for threshold, removed unused debugit depenency
jmetz File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,273 @@ | ||
use crate::core::{ColourModel, Image}; | ||
use crate::core::{PixelBound}; | ||
use crate::processing::*; | ||
use ndarray::prelude::*; | ||
use num_traits::cast::{FromPrimitive}; | ||
use num_traits::cast::{ToPrimitive}; | ||
use ndarray_stats::QuantileExt; | ||
use ndarray_stats::HistogramExt; | ||
use ndarray_stats::histogram::{Grid, Bins, Edges}; | ||
use num_traits::{Num, NumAssignOps}; | ||
use std::marker::PhantomData; | ||
|
||
// Development | ||
#[cfg(test)] | ||
use assert_approx_eq::assert_approx_eq; | ||
#[cfg(test)] | ||
use noisy_float::types::n64; | ||
|
||
|
||
/// Runs the Otsu Thresholding algorithm on a type T | ||
pub trait ThresholdOtsuExt<T> { | ||
/// Output type, this is different as Otsu outputs a binary image | ||
type Output; | ||
|
||
/// Run the Otsu threshold detection algorithm with the | ||
/// given parameters. Due to Otsu being specified as working | ||
/// on greyscale images all current implementations | ||
/// assume a single channel image returning an error otherwise. | ||
fn threshold_otsu(&self) -> Result<Self::Output, Error>; | ||
} | ||
|
||
/// Runs the Mean Thresholding algorithm on a type T | ||
pub trait ThresholdMeanExt<T> { | ||
/// Output type, this is different as Otsu outputs a binary image | ||
type Output; | ||
|
||
/// Run the Otsu threshold detection algorithm with the | ||
/// given parameters. Due to Otsu being specified as working | ||
/// on greyscale images all current implementations | ||
/// assume a single channel image returning an error otherwise. | ||
fn threshold_mean(&self) -> Result<Self::Output, Error>; | ||
} | ||
|
||
impl<T, C> ThresholdOtsuExt<T> for Image<T, C> | ||
where | ||
Image<T, C>: Clone, | ||
T: Copy + Clone + Ord + Num + NumAssignOps + ToPrimitive + FromPrimitive + PixelBound, | ||
C: ColourModel, | ||
{ | ||
type Output = Image<bool, C>; | ||
|
||
fn threshold_otsu(&self) -> Result<Self::Output, Error> { | ||
let data = self.data.threshold_otsu()?; | ||
Ok(Self::Output { | ||
data, | ||
model: PhantomData, | ||
}) | ||
} | ||
} | ||
|
||
impl<T> ThresholdOtsuExt<T> for Array3<T> | ||
where | ||
T: Copy + Clone + Ord + Num + NumAssignOps + ToPrimitive + FromPrimitive | ||
{ | ||
type Output = Array3<bool>; | ||
|
||
fn threshold_otsu(&self) -> Result<Self::Output, Error> { | ||
if self.shape()[2] > 1 { | ||
Err(Error::ChannelDimensionMismatch) | ||
} else { | ||
let value = calculate_threshold_otsu(&self)?; | ||
let mask = apply_threshold(self, value); | ||
Ok(mask) | ||
} | ||
} | ||
} | ||
|
||
/// | ||
/// Calculates Otsu's threshold | ||
/// Works per channel, but currently | ||
/// assumes grayscale (see the error above if number of channels is > 1 | ||
/// i.e. single channel; otherwise we need to output all 3 threshold values). | ||
/// Todo: Add optional nbins | ||
/// | ||
fn calculate_threshold_otsu<T>(mat: &Array3<T>) -> Result<f64, Error> | ||
where | ||
T: Copy + Clone + Ord + Num + NumAssignOps + ToPrimitive + FromPrimitive | ||
{ | ||
let mut threshold = 0.0; | ||
let n_bins = 255; | ||
for c in mat.axis_iter(Axis(2)) { | ||
let scale_factor = (n_bins) as f64 | ||
/(c.max().unwrap().to_f64().unwrap()); | ||
let edges_vec: Vec<u8> = (0..n_bins).collect(); | ||
let grid = Grid::from(vec![Bins::new(Edges::from(edges_vec))]); | ||
|
||
// get the histogram | ||
let flat = Array::from_iter(c.iter()).insert_axis(Axis(1)); | ||
let flat2 = flat.mapv( | ||
|x| ((*x).to_f64().unwrap() * scale_factor).to_u8().unwrap() | ||
); | ||
let hist = flat2.histogram(grid); | ||
// Straight out of wikipedia: | ||
let counts = hist.counts(); | ||
let total = counts.sum().to_f64().unwrap(); | ||
let counts = Array::from_iter(counts.iter()); | ||
// NOTE: Could use the cdf generation for skimage-esque implementation | ||
// which entails a cumulative sum of the standard histogram | ||
let mut sum_b = 0.0; | ||
let mut weight_b = 0.0; | ||
let mut maximum = 0.0; | ||
let mut level = 0.0; | ||
let mut sum_intensity = 0.0; | ||
for (index, count) in counts.indexed_iter(){ | ||
sum_intensity += (index as f64) * (*count).to_f64().unwrap(); | ||
} | ||
for (index, count) in counts.indexed_iter(){ | ||
weight_b = weight_b + count.to_f64().unwrap(); | ||
sum_b = sum_b + (index as f64)* count.to_f64().unwrap(); | ||
let weight_f = total - weight_b; | ||
if (weight_b > 0.0) && (weight_f > 0.0){ | ||
let mean_f = (sum_intensity - sum_b) / weight_f; | ||
let val = weight_b * weight_f | ||
* ((sum_b / weight_b) - mean_f) | ||
* ((sum_b / weight_b) - mean_f); | ||
if val > maximum{ | ||
level = 1.0 + (index as f64); | ||
maximum = val; | ||
} | ||
} | ||
} | ||
threshold = level as f64 / scale_factor; | ||
} | ||
Ok(threshold) | ||
} | ||
|
||
impl<T, C> ThresholdMeanExt<T> for Image<T, C> | ||
where | ||
Image<T, C>: Clone, | ||
T: Copy + Clone + Ord + Num + NumAssignOps + ToPrimitive + FromPrimitive + PixelBound, | ||
C: ColourModel, | ||
{ | ||
type Output = Image<bool, C>; | ||
|
||
fn threshold_mean(&self) -> Result<Self::Output, Error> { | ||
let data = self.data.threshold_mean()?; | ||
Ok(Self::Output { | ||
data, | ||
model: PhantomData, | ||
}) | ||
} | ||
} | ||
|
||
impl<T> ThresholdMeanExt<T> for Array3<T> | ||
where | ||
T: Copy + Clone + Ord + Num + NumAssignOps + ToPrimitive + FromPrimitive | ||
{ | ||
type Output = Array3<bool>; | ||
|
||
fn threshold_mean(&self) -> Result<Self::Output, Error> { | ||
if self.shape()[2] > 1 { | ||
Err(Error::ChannelDimensionMismatch) | ||
} else { | ||
let value = calculate_threshold_mean(&self)?; | ||
let mask = apply_threshold(self, value); | ||
Ok(mask) | ||
} | ||
} | ||
} | ||
|
||
fn calculate_threshold_mean<T>(array: &Array3<T>) -> Result<f64, Error> | ||
where | ||
T: Copy + Clone + Num + NumAssignOps + ToPrimitive + FromPrimitive | ||
{ | ||
Ok(array.sum().to_f64().unwrap() / array.len() as f64) | ||
} | ||
|
||
|
||
fn apply_threshold<T>(data: &Array3<T>, threshold: f64) -> Array3<bool> | ||
where | ||
T: Copy + Clone + Num + NumAssignOps + ToPrimitive + FromPrimitive, | ||
{ | ||
let result = data.mapv(|x| x.to_f64().unwrap() >= threshold); | ||
result | ||
} | ||
|
||
|
||
#[cfg(test)] | ||
mod tests { | ||
use super::*; | ||
use ndarray::arr3; | ||
|
||
#[test] | ||
fn threshold_apply_threshold() { | ||
let data = arr3(&[ | ||
[[0.2], [0.4], [0.0]], | ||
[[0.7], [0.5], [0.8]], | ||
[[0.1], [0.6], [0.0]], | ||
]); | ||
|
||
let expected = arr3(&[ | ||
[[false], [false], [false]], | ||
[[true], [true], [true]], | ||
[[false], [true], [false]], | ||
]); | ||
|
||
let result = apply_threshold(&data, 0.5); | ||
|
||
assert_eq!(result, expected); | ||
} | ||
|
||
#[test] | ||
fn threshold_calculate_threshold_otsu_ints() { | ||
let data = arr3(&[ | ||
[[2], [4], [0]], | ||
[[7], [5], [8]], | ||
[[1], [6], [0]], | ||
]); | ||
let result = calculate_threshold_otsu(&data).unwrap(); | ||
println!("Done {}", result); | ||
|
||
// Calculated using Python's skimage.filters.threshold_otsu | ||
// on int input array. Float array returns 2.0156... | ||
let expected = 2.0; | ||
|
||
assert_approx_eq!(result, expected, 5e-1); | ||
} | ||
|
||
#[test] | ||
fn threshold_calculate_threshold_otsu_floats() { | ||
let data = arr3(&[ | ||
[[n64(2.0)], [n64(4.0)], [n64(0.0)]], | ||
[[n64(7.0)], [n64(5.0)], [n64(8.0)]], | ||
[[n64(1.0)], [n64(6.0)], [n64(0.0)]], | ||
]); | ||
|
||
let result = calculate_threshold_otsu(&data).unwrap(); | ||
|
||
// Calculated using Python's skimage.filters.threshold_otsu | ||
// on int input array. Float array returns 2.0156... | ||
let expected = 2.0156; | ||
|
||
assert_approx_eq!(result, expected, 5e-1); | ||
} | ||
|
||
#[test] | ||
fn threshold_calculate_threshold_mean_ints() { | ||
let data = arr3(&[ | ||
[[4], [4], [4]], | ||
[[5], [5], [5]], | ||
[[6], [6], [6]], | ||
]); | ||
|
||
let result = calculate_threshold_mean(&data).unwrap(); | ||
let expected = 5.0; | ||
|
||
assert_approx_eq!(result, expected, 1e-16); | ||
} | ||
|
||
#[test] | ||
fn threshold_calculate_threshold_mean_floats() { | ||
let data = arr3(&[ | ||
[[4.0], [4.0], [4.0]], | ||
[[5.0], [5.0], [5.0]], | ||
[[6.0], [6.0], [6.0]], | ||
]); | ||
|
||
let result = calculate_threshold_mean(&data).unwrap(); | ||
let expected = 5.0; | ||
|
||
assert_approx_eq!(result, expected, 1e-16); | ||
} | ||
} |
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Why have threshold as
f64
and not typeT
?There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
In general thresholding algorithms sometimes generate thresholds of a different type to the input (consider for example the mean threshold on
ints
, which at least from a theory perspective should be a float); if casting back toT
is used, then we would need to decide on a strategy (round
orfloor
?), which would likely depend also on whether the threshold value is used for>
or>=
comparison when being applied to generate the binary result.Am certainly open to adjusting this to
T
if we can decide on a consistent strategy that is also in-line with how this is done in other frameworks.There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
That all makes sense, I'm happy with keeping it as
f64