Direct Convolution Package

Table of Contents

1 Introduction

This package goal is to compute convolution products

\begin{equation} \label{eq:main} \gamma[k]=\sum\limits_{i\in\Omega^\alpha}\alpha[i]\beta[k+\lambda i],\text{ with }\lambda\in\mathbb{Z}^* \end{equation}

a_offset.png

using direct (no FFT) methods.

2 API documentation

2.1 Linear Filter

Index: [L] LinearFilter [f] fcoef [l] length [o] offset [r] range

  • LinearFilter 
abstract type LinearFilter{T<:Number}

Abstract type defining a linear filter. A linear filter is defined by its coefficients and by its domain

linearFilter.jl:7, back to index 

function LinearFilter(c::AbstractArray{T,1},offset::Int)::LinearFilter  where {T}

Creates a linear filter from its coefficients and an offset

The offset is the position of the filter coefficient to be aligned with zero, see range(…).

Example: 

f=LinearFilter([0:5;],4);
hcat([range(f);],fcoef(f))

6×2 Array{Int64,2}:
 -4  0
 -3  1
 -2  2
 -1  3
  0  4
  1  5

linearFilter.jl:104, back to index 

function LinearFilter(c::AbstractArray{T,1})::LinearFilter  where {T}

Creates a centered linear filter, it must have an odd number of coefficients, \(2n+1\) and is centered by construction (offset=n)

Example: 

f=LinearFilter([0:4;]);
hcat([range(f);],fcoef(f))

5×2 Array{Int64,2}:
 -2  0
 -1  1
  0  2
  1  3
  2  4

linearFilter.jl:118, back to index

  • fcoef 
fcoef(c::LinearFilter)

Returns filter coefficients

linearFilter.jl:13, back to index

  • length 
length(c::LinearFilter)::Int

Returns filter length

linearFilter.jl:17, back to index

  • offset 
offset(c::LinearFilter)::Int

Returns filter offset

Caveat: the first position is 0 (and not 1)

See: range(…)

linearFilter.jl:20, back to index

  • range 
range(c::LinearFilter)::UnitRange

Returns filter range \(\Omega\)

Filter support is defined by \[ \Omega_\alpha = [ -\text{offset}(\alpha) , \text{size}(\alpha) -\text{offset}(\alpha) - 1 ] \]

linearFilter.jl:36, back to index

2.2 Convolution functions

These are the main functions of the package, allowing to compute Eq. \ref{eq:main}.

2.2.1 Boundary extensions

Index: [B] BoundaryExtension [C] ConstantBE [M] MirrorBE [P] PeriodicBE [Z] ZeroPaddingBE

  • BoundaryExtension 
 abstract type BoundaryExtension

Used for tag dispatching, parent of available boundary extensions 

 subtypes(BoundaryExtension)

 4-element Array{Union{DataType, UnionAll},1}:
  DirectConvolution.ConstantBE   
  DirectConvolution.MirrorBE     
  DirectConvolution.PeriodicBE   
  DirectConvolution.ZeroPaddingBE

directConvolution.jl:11, back to index

  • ConstantBE 
 struct ConstantBE <: BoundaryExtension

directConvolution.jl:21, back to index

  • MirrorBE 
 struct MirrorBE <: BoundaryExtension

directConvolution.jl:25, back to index

  • PeriodicBE 
 struct PeriodicBE <: BoundaryExtension

directConvolution.jl:23, back to index

  • ZeroPaddingBE 
 struct ZeroPaddingBE <: BoundaryExtension

directConvolution.jl:19, back to index

2.2.2 Convolution computation

Index: [d] directConv, directConv!, directConv2D!, directCrossCorrelation, directCrossCorrelation2D

  • directConv 
 function directConv(α::LinearFilter{T},

                     λ::Int64,

                     β::AbstractArray{T,1},

                     ::Type{LeftBE}=ZeroPaddingBE,
                     ::Type{RightBE}=ZeroPaddingBE) where {T <: Number,
                                                           LeftBE <: BoundaryExtension,
                                                           RightBE <: BoundaryExtension}

Computes a convolution.

Convenience function that allocate \(\gamma\) and compute all its component using directConv!(…)

Returns: \(\gamma\) a created vector of length identical to the \(\beta\) one.

Example: 

 β=[1:15;];
 γ=ones(Int,15);
 α=LinearFilter([0,0,1],0);
 γ=directConv(α,1,β);
 hcat([1:length(γ);],γ)'

 2×15 Array{Int64,2}:
  1  2  3  4  5  6  7   8   9  10  11  12  13  14  15
  3  4  5  6  7  8  9  10  11  12  13  14  15   0   0

directConvolution.jl:296, back to index 

 function directConv(α::LinearFilter{T},
                     β::AbstractArray{T,1},

                     ::Type{LeftBE}=ZeroPaddingBE,
                     ::Type{RightBE}=ZeroPaddingBE) where {T <: Number,
                                                           LeftBE <: BoundaryExtension,
                                                           RightBE <: BoundaryExtension}

Computes a convolution.

This is a convenience function where \(\lambda=-1\)

Returns: \(\gamma\) a created vector of length identical to the \(\beta\) one.

directConvolution.jl:342, back to index

  • directConv! 
 function directConv!(α::LinearFilter{T},
                      λ::Int,

                      β::AbstractArray{T,1},

                      γ::AbstractArray{T,1},
                      Ωγ::UnitRange{Int},

                      ::Type{LeftBE}=ZeroPaddingBE,
                      ::Type{RightBE}=ZeroPaddingBE;

                      accumulate::Bool=false)::Void where {T <: Number,
                                                           LeftBE <: BoundaryExtension,
                                                           RightBE <: BoundaryExtension}

Computes a convolution.

Inplace modification of \(\gamma[k], k\in\Omega_\gamma\). \[ \gamma[k]=\sum\limits_{i\in\Omega^\alpha}\alpha[i]\beta[k+\lambda i],\text{ with }\lambda\in\mathbb{Z}^* \] If \(k\notin \Omega_\gamma\), \(\gamma[k]\) is unmodified.

If accumulate=false then an erasing step \(\gamma[k]=0, k\in\Omega_\gamma\) is performed before computation.

If \(\lambda=-1\) you compute a convolution, if \(\lambda=+1\) you compute a cross-correlation.

Example: 

 β=[1:15;];
 γ=ones(Int,15);
 α=LinearFilter([0,0,1],0);
 directConv!(α,1,β,γ,5:10);
 hcat([1:length(γ);],γ)'

 2×15 Array{Int64,2}:
  1  2  3  4  5  6  7   8   9  10  11  12  13  14  15
  1  1  1  1  7  8  9  10  11  12   1   1   1   1   1

directConvolution.jl:239, back to index

  • directConv2D! 
 function directConv2D!(α_I::LinearFilter{T},
                        λ_I::Int,
                        α_J::LinearFilter{T},
                        λ_J::Int,

                        β::AbstractArray{T,2},

                        min_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                        max_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                        min_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                        max_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE)::Void where {T<:Number}

Computes a 2D (separable) convolution.

For general information about parameters, see directConv!(…)

α_I must be interpreted as filter for running index I

CAVEAT: the result overwrites β

TODO: @parallel

directConvolution.jl:384, back to index

  • directCrossCorrelation 
 function directCrossCorrelation(α::LinearFilter{T},
                                 β::AbstractArray{T,1},

                                 ::Type{LeftBE}=ZeroPaddingBE,
                                 ::Type{RightBE}=ZeroPaddingBE) where {T <: Number,
                                                                       LeftBE <: BoundaryExtension,
                                                                       RightBE <: BoundaryExtension}

Computes a cross-correlation

This is a convenience function where \(\lambda=+1\)

Returns: \(\gamma\) a created vector of length identical to the \(\beta\) one.

directConvolution.jl:361, back to index

  • directCrossCorrelation2D 
 function directCrossCorrelation2D(α_I::LinearFilter{T},
                                   α_J::LinearFilter{T},

                                   β::AbstractArray{T,2},

                                   min_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                                   max_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                                   min_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                                   max_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE)::Array{T,2} where {T<:Number}

Computes a 2D cross-correlation

This is a wrapper that calls directConv2D!(…)

Note: β is not modified, instead the function returns the result.

directConvolution.jl:455, back to index

2.3 Savitzky-Golay Filters

Index: [S] SG_Filter [a] apply_SG_filter, apply_SG_filter2D [f] filter [l] length [m] maxDerivativeOrder [p] polynomialOrder

  • SG_Filter 
struct SG_Filter{T<:AbstractFloat,N}

A structure to store Savitzky-Golay filters.

SG_Filter.jl:30, back to index 

function SG_Filter(T::DataType=Float64;halfWidth::Int=5,degree::Int=2)::SG_Filter

Creates a set of Savitzky-Golay filters

  • filter length is 2*halfWidth+1
  • polynomial degree is degree

SG_Filter.jl:58, back to index

  • apply_SG_filter 
function apply_SG_filter(signal::Array{T,1},
                         sg::SG_Filter{T};
                         derivativeOrder::Int=0,
                         left_BE::Type{<:BoundaryExtension}=ConstantBE,
                         right_BE::Type{<:BoundaryExtension}=ConstantBE) where {T<:AbstractFloat}

Applies SG filter to 1D signal

Returns:

  • the smoothed signal

SG_Filter.jl:87, back to index

  • apply_SG_filter2D 
function apply_SG_filter2D(signal::Array{T,2},
                           sg_I::SG_Filter{T},
                           sg_J::SG_Filter{T};
                           derivativeOrder_I::Int=0,
                           derivativeOrder_J::Int=0,
                           min_I_BE::Type{<:BoundaryExtension}=ConstantBE,
                           max_I_BE::Type{<:BoundaryExtension}=ConstantBE,
                           min_J_BE::Type{<:BoundaryExtension}=ConstantBE,
                           max_J_BE::Type{<:BoundaryExtension}=ConstantBE) where {T<:AbstractFloat}

Applies SG filter to 2D signal

Returns:

  • the smoothed signal

SG_Filter.jl:106, back to index

  • filter 
function filter(sg::SG_Filter{T,N};derivativeOrder::Int=0) where {T<:AbstractFloat,N}

Returns the filter to be used to compute the smoothed derivatives of order derivativeOrder.

SG_Filter.jl:37, back to index

  • length 
length(sg::SG_Filter{T,N}) where {T<:AbstractFloat,N}

Returns filter length, this is an odd number, see SG_Filter(…)

SG_Filter.jl:45, back to index

  • maxDerivativeOrder 
maxDerivativeOrder(sg::SG_Filter{T,N}) where {T<:AbstractFloat,N}

Maximum order of the smoothed derivatives we can compute with sg

SG_Filter.jl:49, back to index

  • polynomialOrder 
polynomialOrder(sg::SG_Filter{T,N}) where {T<:AbstractFloat,N}

Returns the degree of the polynomial used to construct the Savitzky-Golay filters, see SG_Filter(…).

SG_Filter.jl:54, back to index

2.4 Undecimated Wavelet Transform

2.4.1 UDWT Filters

Index: [U] UDWT_Filter, UDWT_Filter_Biorthogonal [t] tildeψ_filter, tildeϕ_filter [ψ] ψ_filter [ϕ] ϕ_filter

  • UDWT_Filter 
 abstract type UDWT_Filter{T<:Number} <: UDWT_Filter_Biorthogonal{T}

A specialization of UDWT_Filter_Biorthogonal for orthogonal filters.

For orthogonal filters we have: \(\phi=\tilde{\phi}\) and \(\psi=\tilde{\psi}\)

udwt.jl:27, back to index

  • UDWT_Filter_Biorthogonal 
 abstract type UDWT_Filter_Biorthogonal{T<:Number}

Abstract type defining the \(\phi\), \(\psi\), \(\tilde{\phi}\) and \(\tilde{\psi}\) filters associated to an undecimated biorthogonal wavelet transform

udwt.jl:9, back to index

  • tildeψ_filter 
 tildeψ_filter(c::UDWT_Filter_Biorthogonal)::LinearFilter

udwt.jl:22, back to index 

 tildeψ_filter(c::UDWT_Filter)::LinearFilter

udwt.jl:40, back to index

  • tildeϕ_filter 
 tildeϕ_filter(c::UDWT_Filter_Biorthogonal)::LinearFilter

udwt.jl:20, back to index 

 tildeϕ_filter(c::UDWT_Filter)::LinearFilter

udwt.jl:37, back to index

  • ψ_filter 
 ψ_filter(c::UDWT_Filter_Biorthogonal)::LinearFilter

udwt.jl:18, back to index

  • ϕ_filter 
 ϕ_filter(c::UDWT_Filter_Biorthogonal)::LinearFilter

udwt.jl:16, back to index

2.4.2 UDWT Computational subroutines

Index: [U] UDWT [i] inverse_udwt, inverse_udwt! [l] length [s] scale [u] udwt

  • UDWT 
 struct UDWT{T<:Number}

A structure to store 1D UDWT

udwt.jl:80, back to index 

 UDWT{T}(filter::UDWT_Filter_Biorthogonal{T};
             n::Int=0,
             scale::Int=0) where {T<:Number}

Creates an instance

Parameters:

  • filter: used filter
  • scale : max scale
  • n: signal length

udwt.jl:91, back to index

  • inverse_udwt 
 function inverse_udwt(udwt_domain::UDWT{T})::Array{T,1} where {T<:Number}

Performs an 1D inverse undecimated wavelet transform

Returns: a vector containing the reconstructed signal.

udwt.jl:234, back to index

  • inverse_udwt! 
 function inverse_udwt!(udwt_domain::UDWT{T},
                        reconstructed_signal::AbstractArray{T,1}) where {T<:Number}

Performs an 1D inverse undecimated wavelet transform

Caveat: uses a pre-allocated vector reconstructed_signal

udwt.jl:176, back to index

  • length 
 length(udwt::UDWT)::Int

Returns expected signal length

udwt.jl:110, back to index

  • scale 
 scale(udwt::UDWT)::Int

Returns max scale

udwt.jl:107, back to index

  • udwt 
 function udwt(signal::AbstractArray{T,1},
               filter::UDWT_Filter_Biorthogonal{T};
               scale::Int=3) where {T<:Number}

Performs an 1D undecimated wavelet transform

\[(\mathcal{W}_{j+1}f)[u]=(\bar{g}_j*\mathcal{V}_{j}f)[u]\] \[(\mathcal{V}_{j+1}f)[u]=(\bar{h}_j*\mathcal{V}_{j}f)[u]\]

udwt.jl:114, back to index

3 Unit tests

Test Summary:     | Pass  Total
DirectConvolution |   38     38

Author: picaud

Created: 2018-06-25 Mon 18:51

Validate