M-Labs/nalgebra
3
0
Fork 1
Code Issues Pull Requests Packages Projects Releases Wiki Activity

convolution for static and dynamic

Browse Source
This commit is contained in:
Marc Haubenstock 2019-07-28 13:39:43 +02:00
parent 8b85ba1081
commit a62b58b529
2 changed files with 117 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

106
src/linalg/convolution.rs
View File

@ -2,9 +2,9 @@ use std::cmp;
use crate::base::allocator::Allocator; use crate::base::allocator::Allocator;
use crate::base::default_allocator::DefaultAllocator; use crate::base::default_allocator::DefaultAllocator;
use crate::base::dimension::{Dim, DimAdd, DimDiff, DimSub, DimSum}; use crate::base::dimension::{Dim, DimAdd, DimDiff, DimSub, DimSum, DimName, DimMul};
use crate::storage::Storage; use crate::storage::Storage;
use crate::{zero, RealField, Vector, VectorN, U1, Scalar, Matrix, DMatrix}; use crate::{zero, RealField, Vector, VectorN, U1, Scalar, Matrix, MatrixMN, DMatrix};
use crate::alga::general::Field; use crate::alga::general::Field;
impl<N: RealField, D1: Dim, S1: Storage<N, D1>> Vector<N, D1, S1> { impl<N: RealField, D1: Dim, S1: Storage<N, D1>> Vector<N, D1, S1> {
@ -130,7 +130,7 @@ impl<N: RealField, D1: Dim, S1: Storage<N, D1>> Vector<N, D1, S1> {
} }
} }
// TODO: @Investigate -> Only implemented for DMatrix for now as images are usually DMatrix
impl<N: RealField> DMatrix<N> { impl<N: RealField> DMatrix<N> {
/// Returns the convolution of the target vector and a kernel. /// Returns the convolution of the target vector and a kernel.
/// ///
@ -141,7 +141,7 @@ impl<N: RealField> DMatrix<N> {
/// # Errors /// # Errors
/// Inputs must satisfy `self.shape() >= kernel.shape() > 0`. /// Inputs must satisfy `self.shape() >= kernel.shape() > 0`.
/// ///
pub fn mat_convolve_full<R1, C1, S1>( pub fn dmat_convolve_full<R1, C1, S1>(
&self, &self,
kernel: Matrix<N, R1, C1, S1>, //TODO: Would be nice to have an IsOdd trait. As kernels could be of even size atm kernel: Matrix<N, R1, C1, S1>, //TODO: Would be nice to have an IsOdd trait. As kernels could be of even size atm
) -> DMatrix<N> ) -> DMatrix<N>
@ -150,10 +150,10 @@ impl<N: RealField> DMatrix<N> {
C1: Dim, C1: Dim,
S1: Storage<N, R1, C1> S1: Storage<N, R1, C1>
{ {
let mat_rows = self.nrows(); let mat_rows = self.nrows() as i32;
let mat_cols = self.ncols(); let mat_cols = self.ncols() as i32;
let ker_rows = kernel.data.shape().0.value(); let ker_rows = kernel.data.shape().0.value() as i32;
let ker_cols = kernel.data.shape().1.value(); let ker_cols = kernel.data.shape().1.value() as i32;
if ker_rows == 0 || ker_rows > mat_rows || ker_cols == 0|| ker_cols > mat_cols { if ker_rows == 0 || ker_rows > mat_rows || ker_cols == 0|| ker_cols > mat_cols {
panic!( panic!(
@ -166,29 +166,30 @@ impl<N: RealField> DMatrix<N> {
let kernel_size = ker_rows; let kernel_size = ker_rows;
let kernel_min = kernel_size/2; let kernel_min = kernel_size/2;
let zero = zero::<N>(); let zero = zero::<N>();
let mut conv = DMatrix::from_element(mat_cols, mat_rows, zero); let mut conv = DMatrix::<N>::zeros(mat_cols as usize, mat_rows as usize);
for i in 0..mat_rows { for i in 0..mat_rows {
for j in 0..mat_cols { for j in 0..mat_cols {
for k_i in 0..kernel_size { for k_i in 0..kernel_size {
for k_j in 0..kernel_size { for k_j in 0..kernel_size {
let i_matrix = (i + k_i - kernel_min) as i32; let i_matrix = i + k_i - kernel_min;
let j_matrix = (j + k_j - kernel_min) as i32; let j_matrix = j + k_j - kernel_min;
let is_i_in_range = i_matrix >=0 && i_matrix < mat_rows as i32; let is_i_in_range = i_matrix >=0 && i_matrix < mat_rows;
let is_j_in_range = j_matrix >=0 && j_matrix < mat_cols as i32; let is_j_in_range = j_matrix >=0 && j_matrix < mat_cols;
let convolved_value = let convolved_value =
match is_i_in_range && is_j_in_range { match is_i_in_range && is_j_in_range {
true => { true => {
let pixel_value = *self.index((i_matrix as usize, j_matrix as usize)); let pixel_value = *self.index((i_matrix as usize, j_matrix as usize));
let kernel_value = *kernel.index((k_i,k_j)); let kernel_value = *kernel.index((k_i as usize,k_j as usize));
kernel_value*pixel_value kernel_value*pixel_value
} }
//TODO: More behaviour on borders
false => zero false => zero
}; };
*conv.index_mut((i,j)) = convolved_value; *conv.index_mut((i as usize,j as usize)) += convolved_value;
} }
} }
@ -198,7 +199,80 @@ impl<N: RealField> DMatrix<N> {
conv conv
} }
//TODO: rest //TODO: rest ?
}
impl<N: RealField, R1: Dim +DimName, C1: Dim +DimName> MatrixMN<N, R1, C1> where DefaultAllocator: Allocator<N, R1, C1> {
/// Returns the convolution of the target vector and a kernel.
///
/// # Arguments
///
/// * `kernel` - A Matrix with rows > 0 and cols > 0
///
/// # Errors
/// Inputs must satisfy `self.shape() >= kernel.shape() > 0`.
///
pub fn smat_convolve_full<R2, C2, S1>(
&self,
kernel: Matrix<N, R2, C2, S1>, //TODO: Would be nice to have an IsOdd trait. As kernels could be of even size atm
) -> MatrixMN<N, R1, C1>
where
R2: Dim,
C2: Dim,
S1: Storage<N, R2, C2>
{
let mat_rows = self.nrows() as i32;
let mat_cols = self.ncols() as i32;
let ker_rows = kernel.data.shape().0.value() as i32;
let ker_cols = kernel.data.shape().1.value() as i32;
if ker_rows == 0 || ker_rows > mat_rows || ker_cols == 0|| ker_cols > mat_cols {
panic!(
"convolve_full expects `self.nrows() >= kernel.nrows() > 0 and self.ncols() >= kernel.ncols() > 0 `, \
rows received {} and {} respectively. \
cols received {} and {} respectively.",
mat_rows, ker_rows, mat_cols, ker_cols);
}
let kernel_size = ker_rows;
let kernel_min = kernel_size/2;
let zero = zero::<N>();
let mut conv = MatrixMN::<N,R1,C1>::zeros();
for i in 0..mat_rows {
for j in 0..mat_cols {
for k_i in 0..kernel_size {
for k_j in 0..kernel_size {
let i_matrix = i + k_i - kernel_min;
let j_matrix = j + k_j - kernel_min;
let is_i_in_range = i_matrix >=0 && i_matrix < mat_rows;
let is_j_in_range = j_matrix >=0 && j_matrix < mat_cols;
let convolved_value =
match is_i_in_range && is_j_in_range {
true => {
let pixel_value = *self.index((i_matrix as usize, j_matrix as usize));
let kernel_value = *kernel.index((k_i as usize,k_j as usize));
kernel_value*pixel_value
}
//TODO: More behaviour on borders
false => zero
};
*conv.index_mut((i as usize,j as usize)) += convolved_value;
}
}
}
}
conv
}
//TODO: rest ?
} }

49
tests/linalg/convolution.rs
View File

@ -1,4 +1,4 @@
use na::{Vector2,Vector3,Vector4,Vector5,DVector, DMatrix}; use na::{Vector2,Vector3,Vector4,Vector5,DVector, DMatrix, Matrix5, Matrix3};
use std::panic; use std::panic;
// //
@ -163,29 +163,34 @@ fn convolve_valid_check(){
// >>> convolve([1,2,3,4],[1,2],"same") // >>> convolve([1,2,3,4],[1,2],"same")
// array([ 1, 4, 7, 10]) // array([ 1, 4, 7, 10])
#[test] #[test]
fn convolve_same_dmat_check(){ fn convolve_same_mat_check(){
let actual_s = Matrix5::from_vec( vec![3.0,4.0,4.0,4.0,3.0,4.0,5.0,5.0,5.0,4.0,4.0,5.0,5.0,5.0,4.0,4.0,5.0,5.0,5.0,4.0,3.0,4.0,4.0,4.0,3.0]);
let expected_s = Matrix5::from_element(1.0).smat_convolve_full(Matrix3::from_vec(vec![0.0,1.0,0.0,1.0,1.0,1.0,0.0,1.0,0.0]));
assert!(relative_eq!(actual_s, expected_s, epsilon = 1.0e-7));
let actual_d = DMatrix::from_vec(5,5, vec![3.0,4.0,4.0,4.0,3.0,4.0,5.0,5.0,5.0,4.0,4.0,5.0,5.0,5.0,4.0,4.0,5.0,5.0,5.0,4.0,3.0,4.0,4.0,4.0,3.0]); let actual_d = DMatrix::from_vec(5,5, vec![3.0,4.0,4.0,4.0,3.0,4.0,5.0,5.0,5.0,4.0,4.0,5.0,5.0,5.0,4.0,4.0,5.0,5.0,5.0,4.0,3.0,4.0,4.0,4.0,3.0]);
let expected_d = DMatrix::from_element(5,5,1.0).mat_convolve_full(DMatrix::from_vec(3,3,vec![0.0,1.0,0.0,1.0,1.0,1.0,0.0,1.0,0.0])); let expected_d = DMatrix::from_element(5,5,1.0).dmat_convolve_full(DMatrix::from_vec(3,3,vec![0.0,1.0,0.0,1.0,1.0,1.0,0.0,1.0,0.0]));
assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7)); assert!(relative_eq!(actual_d, expected_d, epsilon = 1.0e-7));
// // Panic Tests // Panic Tests
// // These really only apply to dynamic sized vectors // These really only apply to dynamic sized vectors
// assert!( assert!(
// panic::catch_unwind(|| { panic::catch_unwind(|| {
// DVector::from_vec(vec![1.0, 2.0]).convolve_same(DVector::from_vec(vec![1.0, 2.0, 3.0, 4.0])); DMatrix::from_element(2,2,1.0).dmat_convolve_full(DMatrix::from_vec(3,3,vec![0.0,1.0,0.0,1.0,1.0,1.0,0.0,1.0,0.0]));
// }).is_err() }).is_err()
// ); );
//
// assert!( assert!(
// panic::catch_unwind(|| { panic::catch_unwind(|| {
// DVector::<f32>::from_vec(vec![]).convolve_same(DVector::from_vec(vec![1.0, 2.0, 3.0, 4.0])); DMatrix::from_element(0,0,1.0).dmat_convolve_full(DMatrix::from_vec(3,3,vec![0.0,1.0,0.0,1.0,1.0,1.0,0.0,1.0,0.0]));
// }).is_err() }).is_err()
// ); );
//
// assert!( assert!(
// panic::catch_unwind(|| { panic::catch_unwind(|| {
// DVector::from_vec(vec![1.0, 2.0, 3.0, 4.0]).convolve_same(DVector::<f32>::from_vec(vec![])); DMatrix::from_element(5,5,1.0).dmat_convolve_full(DMatrix::from_vec(0,0,vec![]));
// }).is_err() }).is_err()
// ); );
} }