Merge pull request #32 from JuliaHCI/ml/geometry

Project.toml

@@ -1,7 +1,7 @@
 name = "HCIToolbox"
 uuid = "b6cd55e5-4d02-4e12-b82c-005f67e784bf"
 authors = ["Miles Lucas <[email protected]>"]
-version = "0.4.3"
+version = "0.4.4"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"

docs/Project.toml

@@ -1,9 +1,7 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-LossFunctions = "30fc2ffe-d236-52d8-8643-a9d8f7c094a7"
-Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-Documenter = "0.25"
+Documenter = "0.26"

docs/make.jl

@@ -22,7 +22,12 @@ makedocs(;
     ),
     pages = [
         "Home" => "index.md",
-        "API/Reference" => "api.md"
+        "Processing" => "processing.md",
+        "Spectral Processing" => "sdi.md",
+        "Geometries" => "geometry.md",
+        "Signal Injection" => "inject.md",
+        "Utilites" => "utils.md",
+        "Index" => "api.md"
     ]
 )
 

docs/src/api.md

@@ -1,14 +1,6 @@
-# API/Reference
+# Index
 
-For now, here is a dump of all documented functions and types.
-
-## Index
+Here is a dump of all documented functions and types
 
 ```@index
 ```
-
-## API/Reference
-
-```@autodocs
-Modules = [HCIToolbox, HCIToolbox.Kernels]
-```

docs/src/geometry.md

@@ -0,0 +1,50 @@
+# Geometries
+
+A common way of analyzing HCI data is to look at specific spatial extents, such as annuli. These geometries can be thought of as *spatially filtering* the input data.
+
+For example, to create a geometry that is a concentric circle with radius `r`, we could filter a single frame like this
+
+```julia
+frame = ones(101, 101)
+idxs = CartesianIndices(frame)
+radius = 10
+center = (51, 51)
+# only include indices that are within circle
+idxs_inside_circle = filter(idxs) do idx
+    translated = idx.I .- center
+    dist = sqrt(sum(abs2, translated))
+    return dist ≤ radius
+end
+```
+
+using these filtered indices we can *mask* the data with something like
+
+```julia
+masked = zero(frame)
+masked[idxs_inside_circle] = frames[idxs_inside_circle]
+```
+
+more useful, though, is filtering the data. If we think of the frame as a sample of pixels unrolled into a vector, we can filter that vector and only use the pixels that are within the mask.
+
+```julia
+filtered = frame[idxs_inside_circle]
+```
+
+This is very convenient for statistical algorithms wince we are *filtering* the data instead of just *masking* it, which greatly reduces the number of pixels. For example, the circle defined above only uses 4% of the data, so why waste time processing the rest?
+
+## Index
+
+```@index
+Pages = ["geometry.md"]
+```
+
+## API/Reference
+
+```@docs
+AnnulusView
+MultiAnnulusView
+eachannulus
+inverse
+inverse!
+copyto!
+```

docs/src/index.md

@@ -6,6 +6,7 @@ CurrentModule = HCIToolbox
 
 [![GitHub](https://img.shields.io/badge/Code-GitHub-black.svg)](https://github.com/juliahci/HCIToolbox.jl)
 [![Build Status](https://github.com/juliahci/HCIToolbox.jl/workflows/CI/badge.svg?branch=master)](https://github.com/juliahci/HCIToolbox.jl/actions)
+[![PkgEval](https://juliaci.github.io/NanosoldierReports/pkgeval_badges/H/HCIToolbox.svg)](https://juliaci.github.io/NanosoldierReports/pkgeval_badges/report.html)
 [![Coverage](https://codecov.io/gh/juliahci/HCIToolbox.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/juliahci/HCIToolbox.jl)
 [![License](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
 

docs/src/inject.md

@@ -0,0 +1,21 @@
+# Signal Injection
+
+The following functions are used for injecting fake signal into HCI data.
+
+## Index
+
+```@index
+Pages = ["inject.md"]
+```
+
+## API/Reference
+
+```@docs
+inject
+inject!
+construct
+Kernels
+Kernels.Gaussian
+Kernels.AiryDisk
+Kernels.Moffat
+```

docs/src/processing.md

@@ -0,0 +1,27 @@
+# Processing
+
+There are a variety of functions defined by HCIToolbox to process HCI data. In particular, routines like derotating cubes, collapsing data cubes, scaling SDI data, etc. are available.
+
+!!! tip "Multi-threading"
+    Many of the methods that work on cubes of data try to multi-thread the operations along the time axis. Make sure you set the environment variable `JULIA_NUM_THREADS` before staring your runtime to take advantage of this.
+
+## Index
+
+```@index
+Pages = ["processing.md"]
+```
+
+## API/Reference
+
+```@docs
+collapse
+collapse!
+crop
+cropview
+derotate
+derotate!
+expand
+flatten
+shift_frame
+shift_frame!
+```

docs/src/sdi.md

@@ -0,0 +1,19 @@
+# Spectral Processing
+
+The following functions are used specifically for processing spectral differential imaging (SDI) tensors/cubes.
+
+## Index
+
+```@index
+Pages = ["sdi.md"]
+```
+
+## API/Reference
+
+```@docs
+invscale
+invscale_and_collapse
+scale_list
+scale
+scale_and_stack
+```

docs/src/utils.md

@@ -0,0 +1,20 @@
+# Utilities
+
+The following routines are general utility or otherwise don't fit into a category
+
+## Index
+
+```@index
+Pages = ["utils.md"]
+```
+
+## API/Reference
+
+```@docs
+mask_circle
+mask_circle!
+mask_annulus
+mask_annulus!
+normalize_par_angles
+normalize_par_angles!
+```

src/HCIToolbox.jl

@@ -12,7 +12,7 @@ export derotate,
        mask_circle,
        mask_annulus!,
        mask_annulus,
-       normalize_par_angles, 
+       normalize_par_angles,
        normalize_par_angles!,
        inject,
        inject!,
@@ -22,7 +22,13 @@ export derotate,
        shift_frame!,
        scale_and_stack,
        invscale_and_collapse,
-       scale_list
+       scale_list,
+       # geometry
+       AnnulusView,
+       MultiAnnulusView,
+       eachannulus,
+       inverse!,
+       inverse
 
 # Utilities for dealing with cubes like derotating and median combining
 include("morphology.jl")
@@ -34,4 +40,6 @@ include("angles.jl")
 
 include("scaling.jl")
 
+include("geometry/geometry.jl")
+
 end

src/geometry/annulus.jl

@@ -0,0 +1,125 @@
+
+struct AnnulusView{T,N,M<:AbstractArray{T,N},IT} <: AbstractArray{T,N}
+    parent::M
+    rmin::Float64
+    rmax::Float64
+    indices::IT
+    fill::T
+end
+
+"""
+    AnnulusView(cube::AbstractArray{T,3};
+                inner=0, outer=last(size(parent))/2 + 0.5,
+                fill=0)
+
+Cut out an annulus with inner radius `inner` and outer radius `outer`. Values that fall outside of this region will be replaced with `fill`. This does not copy any data, it is merely a view into the data.
+"""
+function AnnulusView(parent::AbstractArray{T,3};
+                     inner=0,
+                     outer=(size(parent, 3) + 1) / 2,
+                     fill=zero(T)) where T
+    # check inputs
+    0 ≤ inner < outer || error("Invalid annulus region [$inner, $outer]")
+    time_axis = axes(parent, 1)
+    space_indices = CartesianIndices((axes(parent, 2), axes(parent, 3)))
+    space_axis = filter(idx -> inside_annulus(inner, outer, center(parent)[2:3], idx), space_indices)
+    return AnnulusView(parent, Float64(inner), Float64(outer), (time_axis, space_axis), T(fill))
+end
+
+Base.parent(view::AnnulusView) = view.parent
+Base.size(view::AnnulusView) = size(parent(view))
+Base.copy(view::AnnulusView) = AnnulusView(copy(parent(view)), view.rmin, view.rmax, view.indices, view.fill)
+
+inside_annulus(view::AnnulusView, args...) = inside_annulus(view.rmin, view.rmax, center(view)[2:3], args...)
+
+@propagate_inbounds function Base.getindex(view::AnnulusView{T,N}, idx::Vararg{<:Integer,N}) where {T,N}
+    @boundscheck checkbounds(parent(view), idx...)
+    ifelse(inside_annulus(view, idx...), convert(T, parent(view)[idx...]), view.fill)
+end
+
+@propagate_inbounds function Base.setindex!(view::AnnulusView{T,N}, val, idx::Vararg{<:Integer,N}) where {T,N}
+    @boundscheck checkbounds(parent(view), idx...)
+    if inside_annulus(view, idx...)
+        parent(view)[idx...] = val
+    else
+        view.fill
+    end
+end
+
+"""
+    (::AnnulusView)(asview=false)
+
+Return the pixels that fall within the annulus as a matrix. This matrix is equivalent to unrolling each frame and then spatially filtering the pixels outside the annulus. If `asview` is true, the returned values will be a view of the parent array instead of a copy.
+
+# Examples
+```jldoctest
+julia> ann = AnnulusView(ones(10, 101, 101); inner=5, outer=20);
+
+julia> X = ann();
+
+julia> size(X)
+(10, 1188)
+```
+"""
+function (view::AnnulusView)(asview=false)
+    if asview
+        @view parent(view)[view.indices...]
+    else
+        parent(view)[view.indices...]
+    end
+end
+
+"""
+    copyto!(::AnnulusView, mat::AbstractMatrix)
+
+Copy the pixels from `mat` into the pixels in the annulus. `mat` should have the same size as the matrix output from [`AnnulusView`](@ref)
+
+# Examples
+```jldoctest
+julia> ann = AnnulusView(ones(10, 101, 101); inner=5, outer=20);
+
+julia> X = ann();
+
+julia> new_ann = copyto!(ann, -X);
+
+julia> new_ann() ≈ -X
+true
+```
+"""
+function Base.copyto!(view::AnnulusView, mat::AbstractMatrix)
+    inverse!(view, view.parent, mat)
+    return view
+end
+
+"""
+    inverse!(::AnnulusView, out, mat)
+
+In-place version of [`inverse`](@ref) that fills `out` in-place.
+"""
+function inverse!(view::AnnulusView, out, mat)
+    out[view.indices...] = mat
+    return out
+end
+
+
+"""
+    inverse(::AnnulusView, mat::AbstractMatrix)
+
+Generate a cube similar to the view with the pixels from `mat`. `mat` should have the same size as the matrix output from [`AnnulusView`](@ref)
+
+# Examples
+```jldoctest
+julia> ann = AnnulusView(ones(10, 101, 101); inner=5, outer=20);
+
+julia> X = ann();
+
+julia> out = inverse(ann, -X);
+
+julia> out ≈ -ann
+true
+```
+"""
+function inverse(view::AnnulusView, mat)
+    out = fill!(similar(parent(view)), view.fill)
+    inverse!(view, out, mat)
+end

src/geometry/geometry.jl

@@ -0,0 +1,14 @@
+using CoordinateTransformations
+using StaticArrays
+using Base: @propagate_inbounds
+
+inside_annulus(rmin, rmax, center, idx...) = inside_annulus(rmin, rmax, center, SVector(idx[end-1:end]))
+inside_annulus(rmin, rmax, center, idx::CartesianIndex{2}) = inside_annulus(rmin, rmax, center, SVector(idx.I))
+function inside_annulus(rmin, rmax, center::AbstractVector, position::AbstractVector)
+    Δ = center - position
+    r = sqrt(sum(abs2, Δ))
+    return rmin ≤ r ≤ rmax
+end
+
+include("annulus.jl")
+include("multiannulus.jl")