diff --git a/.claude/freshen-package-status b/.claude/freshen-package-status new file mode 100644 index 0000000..b2bf205 --- /dev/null +++ b/.claude/freshen-package-status @@ -0,0 +1,11 @@ +DONE: design review +DONE: API review +DONE: update .gitignore +DONE: format with runic +DONE: add Aqua.jl +DONE: remove deprecations +DONE: add ExplicitImports.jl +DONE: limit struct mutability +DONE: improve test coverage +DONE: add and improve docstrings +TODO: add or improve documentation diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 65a8e0e..1dae340 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -13,17 +13,15 @@ jobs: fail-fast: false matrix: version: - - '1.6' - - '1.1' + - '1.10' - '1' -# - 'nightly' os: - ubuntu-latest arch: - x64 steps: - - uses: actions/checkout@v2 - - uses: julia-actions/setup-julia@v1 + - uses: actions/checkout@v4 + - uses: julia-actions/setup-julia@v2 with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} @@ -42,6 +40,6 @@ jobs: - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v4 with: file: lcov.info diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml index f549621..14c6cea 100644 --- a/.github/workflows/Documenter.yml +++ b/.github/workflows/Documenter.yml @@ -10,11 +10,14 @@ jobs: name: Documentation runs-on: ubuntu-latest steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - name: Install HolyLab Registry run: julia -e 'using Pkg; Pkg.Registry.add([RegistrySpec(name="General"),RegistrySpec(url="https://github.com/HolyLab/HolyLabRegistry")]);' - - uses: julia-actions/julia-buildpkg@latest - - uses: julia-actions/julia-docdeploy@latest + - uses: julia-actions/setup-julia@v2 + with: + version: '1' + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-docdeploy@v1 env: GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} diff --git a/.gitignore b/.gitignore index 49541d3..605dc84 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,7 @@ *.jl.mem deps/deps.jl Manifest.toml +Manifest-v*.toml +docs/Manifest.toml +docs/Manifest-v*.toml +docs/build/ diff --git a/Project.toml b/Project.toml index 3c7fa97..7d765f1 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "RegisterCore" uuid = "67712758-55e7-5c3c-8e85-dda1d7758434" +version = "1.0.0" authors = ["Tim Holy "] -version = "0.2.5" [deps] CenterIndexedArrays = "46a7138f-0d70-54e1-8ada-fb8296f91f24" @@ -10,16 +10,23 @@ ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" Requires = "ae029012-a4dd-5104-9daa-d747884805df" [compat] +Aqua = "0.8" CenterIndexedArrays = "0.2" +ExplicitImports = "1" ImageCore = "0.8.1, 0.9, 0.10" ImageFiltering = "0.5, 0.6, 0.7" +ImageMetadata = "0.9" +Interpolations = "0.13, 0.14, 0.15, 0.16" Requires = "1" +Test = "1" julia = "^1.1" [extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7" ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["ImageMetadata", "Interpolations", "Test"] +test = ["Aqua", "ExplicitImports", "ImageMetadata", "Interpolations", "Test"] diff --git a/README.md b/README.md index 86c7c8c..727942c 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,10 @@ # RegisterCore +[![CI](https://github.com/HolyLab/RegisterCore.jl/actions/workflows/CI.yml/badge.svg)](https://github.com/HolyLab/RegisterCore.jl/actions/workflows/CI.yml) [![codecov](https://codecov.io/gh/HolyLab/RegisterCore.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/HolyLab/RegisterCore.jl) +[![Aqua QA](https://juliatesting.github.io/Aqua.jl/dev/assets/badge.svg)](https://github.com/JuliaTesting/Aqua.jl) +[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://HolyLab.github.io/RegisterCore.jl/stable) +[![](https://img.shields.io/badge/docs-dev-blue.svg)](https://HolyLab.github.io/RegisterCore.jl/dev) This package defines low-level types and operations used for a larger image registration infrastructure. -The documentation gives something of a big-picture overview: - -[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://HolyLab.github.io/RegisterCore.jl/dev) +The documentation gives a big-picture overview of the key concepts. diff --git a/docs/Project.toml b/docs/Project.toml index 3fcbfb7..1814eb3 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,4 +2,4 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" [compat] -Documenter = "~0.23" +Documenter = "1" diff --git a/docs/make.jl b/docs/make.jl index 89c69bb..5ef95cf 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -7,7 +7,11 @@ makedocs( prettyurls = get(ENV, "CI", nothing) == "true" ), modules = [RegisterCore], - pages = ["index.md", "api.md"] + checkdocs = :exports, + pages = [ + "Overview" => "index.md", + "API Reference" => "api.md", + ] ) deploydocs( diff --git a/docs/src/api.md b/docs/src/api.md index 59971f0..6bb51d7 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -1,15 +1,22 @@ # API +```@docs +RegisterCore +``` + ## Types ```@docs NumDenom +MismatchArray +ColonFun +PreprocessSNF ``` -## Functions related to NumDenom +## MismatchArray functions ```@docs -indmin_mismatch +argmin_mismatch maxshift mismatcharrays ratio @@ -20,7 +27,7 @@ separate ```@docs highpass -PreprocessSNF +highpass! paddedview trimmedview ``` diff --git a/docs/src/index.md b/docs/src/index.md index 2be9039..470daf9 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -74,7 +74,7 @@ corresponds to an identical shift in the opposite direction. ## NumDenom -Mismatch computations actually return two numbers, conventionaly +Mismatch computations actually return two numbers, conventionally called `num` and `denom` packed into a type called `NumDenom`. `num` represents the "numerator" of the mismatch, and always holds the sum-of-squared-differences. In our example above, @@ -88,6 +88,31 @@ julia> sum((img[1:end-1,:] - img[2:end,:]).^2) 24 ``` +`NumDenom` supports vector-space arithmetic, which enables interpolation +of mismatch arrays. A brief example without `RegisterMismatch`: + +```jldoctest +julia> using RegisterCore + +julia> a = NumDenom(3.0, 10.0) +NumDenom(3.0,10.0) + +julia> b = NumDenom(1.0, 6.0) +NumDenom(1.0,6.0) + +julia> a + b +NumDenom(4.0,16.0) + +julia> 0.5 * a +NumDenom(1.5,5.0) + +julia> ratio(a, 0.0) +0.3 + +julia> ratio(NumDenom(0.0, 1e-6), 0.0) +0.0 +``` + One key thing to note is that this sum-of-squared-differences includes *only overlapping pixels*. In the example above, `D[2,0]` is `NumDenom(16, 22)`, and the 16 diff --git a/src/RegisterCore.jl b/src/RegisterCore.jl index a8affe7..099d7ff 100644 --- a/src/RegisterCore.jl +++ b/src/RegisterCore.jl @@ -1,8 +1,9 @@ module RegisterCore -using CenterIndexedArrays -using ImageCore, ImageFiltering -using Requires +using CenterIndexedArrays: CenterIndexedArrays, CenterIndexedArray +using ImageCore: ImageCore, Gray, gray +using ImageFiltering: ImageFiltering, KernelFactors, NA, imfilter +using Requires: Requires, @require import Base: +, -, *, / @@ -13,8 +14,9 @@ export ColonFun, PreprocessSNF, # functions + argmin_mismatch, highpass, - indmin_mismatch, + highpass!, maxshift, mismatcharrays, ratio, @@ -23,40 +25,51 @@ export trimmedview """ -The major functions/types exported by this module are: - -- `NumDenom` and `MismatchArray`: packed pair representation of - `(num,denom)` mismatch data -- `separate`: splits a `NumDenom` array into its component `num,denom` arrays -- `indmin_mismatch`: find the location of the minimum mismatch -- `highpass`: highpass filter an image before performing registration -- `PreprocessSNF`: shot-noise standardization and filtering - +Core types and utilities for image registration mismatch computations. + +**Types** +- [`NumDenom`](@ref): a packed `(numerator, denominator)` pair; supports vector-space + arithmetic for use with `Interpolations.jl` +- [`MismatchArray`](@ref): a `CenterIndexedArray` of `NumDenom` elements representing + mismatch over a range of shifts +- [`PreprocessSNF`](@ref): shot-noise–filtering image preprocessor (bias subtract, + square-root transform, band-pass filter) +- [`ColonFun`](@ref): callable singleton returning `Colon()`, for type-stable slicing + +**Functions** +- [`separate`](@ref): unpack a `NumDenom` or `MismatchArray` into `(num, denom)` arrays +- [`ratio`](@ref): convert `NumDenom` to a scalar ratio, with threshold masking +- [`maxshift`](@ref): return the maximum-shift half-size of a `MismatchArray` +- [`mismatcharrays`](@ref): pack array-of-arrays pairs into an array of `MismatchArray`s +- [`argmin_mismatch`](@ref): find the shift index of the minimum mismatch +- [`highpass`](@ref) / [`highpass!`](@ref): high-pass filter an image (Gaussian-based, NaN-safe) +- [`paddedview`](@ref) / [`trimmedview`](@ref): extend/trim a `SubArray` to/from its parent """ RegisterCore """ x = NumDenom(num, denom) -`NumDenom{T}` is an object containing a `(num,denom)` pair. -`x.num` is `num` and `x.denom` is `denom`. +A packed `(num, denom)` pair with element type `T`. Access fields as `x.num` and +`x.denom`. Arguments of type `Gray` are automatically unwrapped; mismatched numeric +types are promoted. -Algebraically, `NumDenom` objects act like 2-vectors, and can be added and multiplied -by scalars: +`NumDenom` objects act as 2-vectors under arithmetic — addition and scalar +multiplication operate component-wise on both fields: - nd1 + nd2 = NumDenom(nd1.num + nd2.num, nd1.denom + nd2.denom) - 2*nd = NumDenom(2*nd.num, 2*nd.denom) + nd1 + nd2 → NumDenom(nd1.num + nd2.num, nd1.denom + nd2.denom) + nd1 - nd2 → NumDenom(nd1.num - nd2.num, nd1.denom - nd2.denom) + c * nd → NumDenom(c*nd.num, c*nd.denom) + nd / c → NumDenom(nd.num/c, nd.denom/c) -Note that this is *not* what you'd get from normal arithmetic with ratios, where, e.g., -`2*nd` would be expected to produce `NumDenom(2*nd.num, nd.denom)`. -The reason for calling these `*` and `+` is for use in `Interpolations.jl`, because it -allows interpolation to be performed on "both arrays" at once without -recomputing the interpolation coefficients. See the documentation for information -about how this is used for performing aperturered mismatch computations. +This vector-space algebra (rather than ratio algebra) is intentional: it allows +`Interpolations.jl` to interpolate numerator and denominator jointly without +recomputing interpolation weights, enabling efficient aperture-windowed mismatch +computations. -As a consequence, there is no `convert(Float64, nd::NumDenom)` method, because the algebra -above breaks any pretense that `NumDenom` numbers are somehow equivalent to ratios. -If you want to convert to a ratio, see [`ratio`](@ref). +As a consequence, `convert(Float64, ::NumDenom)` is deliberately not defined, because +the component-wise arithmetic breaks the ratio interpretation. Use [`ratio`](@ref) to +convert to a scalar ratio. """ struct NumDenom{T <: Number} num::T @@ -95,20 +108,38 @@ function Base.show(io::IO, p::NumDenom) return end +""" + MismatchArray{ND,N,A} + +A `CenterIndexedArray` whose elements are [`NumDenom`](@ref) pairs, representing +packed numerator/denominator mismatch data over a range of shifts. + +The axes are centered at zero, so `axes(D, k)` runs from `-maxshift(D)[k]` to +`+maxshift(D)[k]`; an index value of `(i, j, ...)` corresponds to a shift of +`(i, j, ...)` pixels. + + D = MismatchArray(num::AbstractArray, denom::AbstractArray) + +Pack same-size arrays `num` and `denom` into a `MismatchArray` with centered axes of +half-size `size(num) .÷ 2`. Element type is `NumDenom{promote_type(eltype(num), eltype(denom))}`. + + D = MismatchArray(T, dims::Dims) + D = MismatchArray(T, dims::Integer...) + +Allocate an uninitialized `MismatchArray` with element type `NumDenom{T}` and the given +dimensions. Useful for pre-allocating output before filling with `copyto!`. +""" const MismatchArray{ND <: NumDenom, N, A} = CenterIndexedArray{ND, N, A} """ - mxs = maxshift(D) + mxs = maxshift(D::MismatchArray) -Return the `maxshift` value used to compute the mismatch array `D`. +Return the half-size of the `MismatchArray` `D` as a tuple of integers — i.e., the +maximum shift (in pixels, per dimension) that was searched when computing the mismatch. +Equivalent to `D.halfsize`. """ maxshift(A::MismatchArray) = A.halfsize -""" -`numdenom = MismatchArray(num, denom)` packs the array-pair -`(num,denom)` into a single `MismatchArray`. This is useful -preparation for interpolation. -""" function (::Type{M})(num::AbstractArray, denom::AbstractArray) where {M <: MismatchArray} size(num) == size(denom) || throw(DimensionMismatch("num and denom must have the same size")) T = promote_type(eltype(num), eltype(denom)) @@ -143,9 +174,16 @@ end # The next are mostly used just for testing """ -`mms = mismatcharrays(nums, denoms)` packs array-of-arrays num/denom pairs as an array-of-MismatchArrays. + mms = mismatcharrays(nums, denom::AbstractArray{<:Number}) + mms = mismatcharrays(nums, denoms::AbstractArray{<:AbstractArray}) + +Pack array-of-arrays `(num, denom)` pairs into an `Array{MismatchArray}`. -`mms = mismatcharrays(nums, denom)`, for `denom` a single array, uses the same `denom` array for all `nums`. +The first form uses the same `denom` array for every entry in `nums`. The second form +pairs each `nums[i]` with the corresponding `denoms[i]`; `nums` and `denoms` must have +the same size. + +Returns an `Array{MismatchArray}` with the same shape as `nums`. """ function mismatcharrays(nums::AbstractArray{A}, denom::AbstractArray{T}) where {A <: AbstractArray, T <: Number} first = true @@ -178,8 +216,17 @@ function mismatcharrays(nums::AbstractArray{A1}, denoms::AbstractArray{A2}) wher end """ -`num, denom = separate(mm)` splits an `AbstractArray{NumDenom}` into separate -numerator and denominator arrays. + num, denom = separate(data::AbstractArray{NumDenom{T}}) + num, denom = separate(mm::MismatchArray) + nums, denoms = separate(mma::AbstractArray{<:MismatchArray}) + +Split packed [`NumDenom`](@ref) data into separate numerator and denominator arrays. + +- For a plain `AbstractArray{NumDenom{T}}`, returns a pair of `Array{T}` with the same + size and linear indexing as `data`. +- For a [`MismatchArray`](@ref), returns a pair of `CenterIndexedArray{T}`, preserving + the centered axes (so index ranges correspond to shift values). +- For an array of `MismatchArray`s, returns a pair of `Array{CenterIndexedArray{T}}`. """ function separate(data::AbstractArray{NumDenom{T}}) where {T} num = Array{T}(undef, size(data)) @@ -208,17 +255,25 @@ function separate(mma::AbstractArray{M}) where {M <: MismatchArray} end """ - r = ratio(nd::NumDenom, thresh, fillval=NaN) + r = ratio(nd::NumDenom, thresh; fillval=NaN) + r = ratio(r::Real, thresh; fillval=NaN) -Return `nd.num/nd.denom`, unless `nd.denom < thresh`, in which case return `fillval` converted -to the same type as the ratio. -Choosing a `thresh` of zero will always return the ratio. +Return the ratio `nd.num/nd.denom`, or `fillval` (converted to the ratio type) when +`nd.denom < thresh`. Setting `thresh = 0` always returns the ratio. + +The second form accepts a plain `Real` and returns it unchanged, regardless of `thresh` +and `fillval` — `thresh` and `fillval` are silently ignored. This allows callers to +handle both [`NumDenom`](@ref) and pre-computed ratio arrays uniformly without +branching on the element type. """ -@inline function ratio(nd::NumDenom{T}, thresh, fillval = convert(T, NaN)) where {T} +@inline function ratio(nd::NumDenom{T}, thresh; fillval = convert(T, NaN)) where {T} r = nd.num / nd.denom return nd.denom < thresh ? oftype(r, fillval) : r end -ratio(r::Real, thresh, fillval = NaN) = r +ratio(r::Real, thresh; fillval = NaN) = r + +@deprecate(ratio(nd::NumDenom{T}, thresh, fillval) where {T}, ratio(nd, thresh; fillval=fillval)) +@deprecate(ratio(r::Real, thresh, fillval), ratio(r, thresh)) (::Type{M})(::Type{T}, dims::Dims) where {M <: MismatchArray, T} = CenterIndexedArray{NumDenom{T}}(undef, dims) (::Type{M})(::Type{T}, dims::Integer...) where {M <: MismatchArray, T} = CenterIndexedArray{NumDenom{T}}(undef, dims) @@ -236,13 +291,23 @@ end #### Utility functions #### """ -`index = indmin_mismatch(numdenom, thresh)` returns the location of -the minimum value of what is effectively `num./denom`. However, it -considers only those points for which `denom .> thresh`; moreover, it -will never choose an edge point. `index` is a CartesianIndex into the -arrays. + index = argmin_mismatch(numdenom::MismatchArray, thresh::Real) + index = argmin_mismatch(r::CenterIndexedArray{<:Number}) + +Return the `CartesianIndex` of the minimum mismatch, excluding edge points. + +The first form operates on a [`MismatchArray`](@ref): it computes `num/denom` at each +interior point and finds the minimum among points where `denom > thresh`. The second +form operates on a pre-computed ratio array `r` of plain numbers and finds the interior +minimum unconditionally. + +"Edge points" are the outermost index value on each side of every dimension; they are +always excluded from consideration. + +If no valid point is found (e.g., all `denom ≤ thresh`), returns a zero +`CartesianIndex` — check for this case before using the result as an array index. """ -function indmin_mismatch(numdenom::MismatchArray{NumDenom{T}, N}, thresh::Real) where {T, N} +function argmin_mismatch(numdenom::MismatchArray{NumDenom{T}, N}, thresh::Real) where {T, N} imin = CartesianIndex(ntuple(d -> 0, Val(N))) rmin = typemax(T) threshT = convert(T, thresh) @@ -259,7 +324,7 @@ function indmin_mismatch(numdenom::MismatchArray{NumDenom{T}, N}, thresh::Real) return imin end -function indmin_mismatch(r::CenterIndexedArray{T, N}) where {T <: Number, N} +function argmin_mismatch(r::CenterIndexedArray{T, N}) where {T <: Number, N} imin = CartesianIndex(ntuple(d -> 0, Val(N))) rmin = typemax(T) @inbounds for I in CartesianIndices(map(trimedges, axes(r))) @@ -277,17 +342,30 @@ trimedges(r::AbstractUnitRange) = (first(r) + 1):(last(r) - 1) ### Miscellaneous """ -`datahp = highpass([T], data, sigma)` returns a highpass-filtered -version of `data`, with all negative values truncated at 0. The -highpass is computed by subtracting a lowpass-filtered version of -data, using Gaussian filtering of width `sigma`. As it is based on -`Image.jl`'s Gaussian filter, it gracefully handles `NaN` values. - -If you do not wish to highpass-filter along a particular axis, put -`Inf` into the corresponding slot in `sigma`. - -You may optionally specify the element type of the result, which for -`Integer` or `FixedPoint` inputs defaults to `Float32`. + datahp = highpass([T], data, sigma) + highpass!(out, data, sigma) + highpass!(data, sigma) + +Return (or write in-place) a high-pass–filtered version of `data` with negative values +clamped to zero. The high-pass is computed by subtracting a Gaussian-smoothed copy of +`data` (via `ImageFiltering.jl`), which gracefully handles `NaN` values. + +`sigma` must be an iterable (e.g., a tuple or vector) with one width per dimension of +`data`. To skip filtering along a particular axis, set the corresponding entry to `Inf` +(the input is then returned as-is, converted to `Array{T}`, with no subtraction or +clamping applied). + +For `highpass`, the optional first argument `T` sets the element type of the output: +- `highpass(T, data, sigma)` — allocates an output of element type `T` +- `highpass(data, sigma)` — `T` defaults to `eltype(data)` for `AbstractFloat` inputs, + or `Float32` for `Integer`/`FixedPoint` inputs + +For `highpass!`, the output element type is `eltype(out)`: +- `highpass!(out, data, sigma)` — writes result into `out`; `out` and `data` may be + distinct arrays (useful for pre-allocated buffers in hot loops) +- `highpass!(data, sigma)` — filters `data` in-place + +See also [`PreprocessSNF`](@ref) for a combined shot-noise–filtering preprocessor. """ function highpass(::Type{T}, data::AbstractArray, sigma) where {T} if any(isinf, sigma) @@ -302,26 +380,55 @@ highpass(data::AbstractArray{T}, sigma) where {T <: AbstractFloat} = highpass(T, highpass(data::AbstractArray, sigma) = highpass(Float32, data, sigma) """ -`pp = PreprocessSNF(bias, sigmalp, sigmahp)` constructs an object that -can be used to pre-process an image as `pp(img)`. The "SNF" part of -the name means "shot-noise filtered," meaning that this preprocessor -is specifically designed for situations in which you are dominated by -shot noise (i.e., from photon-counting statistics). - -The processing is of the form -``` - imgout = bandpass(√max(0,img-bias)) -``` -i.e., the image is bias-subtracted, square-root transformed (to turn -shot noise into constant variance), and then band-pass filtered using -Gaussian filters of width `sigmalp` (for the low-pass) and `sigmahp` -(for the high-pass). You can pass `sigmalp=zeros(n)` to skip low-pass -filtering, and `sigmahp=fill(Inf, n)` to skip high-pass filtering. + highpass!(out, data, sigma) + highpass!(data, sigma) + +In-place variant of [`highpass`](@ref). See that function for full documentation. +""" +function highpass!(out::AbstractArray, data::AbstractArray, sigma) + T = eltype(out) + if any(isinf, sigma) + out .= data + else + out .= data .- imfilter(T, data, KernelFactors.IIRGaussian(T, (sigma...,)), NA()) + end + out[out .< 0] .= 0 + return out +end +highpass!(data::AbstractArray, sigma) = highpass!(data, data, sigma) + +""" + pp = PreprocessSNF(bias, sigmalp, sigmahp) + +Construct a shot-noise–filtered preprocessor. Call it as `pp(img)` to preprocess an +image. All arguments are stored as `Float32` (or `Vector{Float32}`), so `Float64` +inputs are silently converted. + +The "SNF" name stands for "shot-noise filtered": this preprocessor is designed for +photon-counting (Poisson) data, where the variance equals the mean. + +Processing steps: +1. Subtract `bias` and clamp to zero: `A′ = max(0, img - bias)` +2. Square-root transform to stabilize variance: `A″ = √A′` +3. High-pass filter with Gaussian width `sigmahp` (subtracts a low-pass Gaussian) +4. Low-pass filter with Gaussian width `sigmalp` + +`sigmalp` and `sigmahp` must each be a vector of length `ndims(img)`. Pass +`sigmalp = zeros(ndims(img))` to skip low-pass filtering; pass +`sigmahp = fill(Inf, ndims(img))` to skip high-pass filtering. + +When called on a `SubArray`, `pp` extends the view to the full parent (via +[`paddedview`](@ref)) before processing, then trims the result back to the original +index range (via [`trimmedview`](@ref)). This preserves any padding context around +the sub-region. + +If `ImageMetadata` is loaded, `pp` also accepts `ImageMeta` arrays and propagates +image properties to the output. """ mutable struct PreprocessSNF # Shot-noise filtered - bias::Float32 - sigmalp::Vector{Float32} - sigmahp::Vector{Float32} + const bias::Float32 + const sigmalp::Vector{Float32} + const sigmahp::Vector{Float32} end # PreprocessSNF(bias::T, sigmalp, sigmahp) = PreprocessSNF{T}(bias, T[sigmalp...], T[sigmahp...]) @@ -350,9 +457,17 @@ end """ -`Apad = paddedview(A)`, for a SubArray `A`, returns a SubArray that -extends to the full parent along any non-sliced dimensions of the -parent. See also [`trimmedview`](@ref). + Apad = paddedview(A::SubArray) + +Return a `SubArray` of `A`'s parent that extends to the full parent size along every +dimension of `A` that was indexed by a `UnitRange`. Dimensions indexed by a scalar are +still dropped (as usual for `SubArray`), and dimensions indexed by a `Slice` are +unchanged. + +Only `Slice`, `Real`, and `UnitRange` index types are supported; any other index type +throws an error. + +See also [`trimmedview`](@ref). """ paddedview(A::SubArray) = _paddedview(A, (), (), A.indices...) _paddedview(A::SubArray{T, N, P, I}, newindexes, newsize) where {T, N, P, I} = @@ -371,10 +486,18 @@ pdsize(A, newsize, d, i::Real) = newsize pdsize(A, newsize, d, i::UnitRange) = tuple(newsize..., size(A, d)) """ -`B = trimmedview(Bpad, A::SubArray)` returns a SubArray `B` with -`axes(B) = axes(A)`. `Bpad` must have the same size as `paddedview(A)`. + B = trimmedview(Bpad::AbstractArray, A::SubArray) + +Return a view `B` of `Bpad` with `axes(B) == axes(A)`. + +`Bpad` must be an `AbstractArray` with the same size as `paddedview(A)` — typically +it is the result of an operation applied to `paddedview(A)`. Dimensions of `A` that +were indexed by a scalar are skipped (they are dropped in `A`); all other dimensions +are sliced with the original index range from `A`. + +See also [`paddedview`](@ref). """ -function trimmedview(Bpad, A::SubArray) +function trimmedview(Bpad::AbstractArray, A::SubArray) ndims(Bpad) == ndims(A) || throw(DimensionMismatch("dimensions $(ndims(Bpad)) and $(ndims(A)) of Bpad and A must match")) return _trimmedview(Bpad, A.parent, 1, (), A.indices...) end @@ -389,7 +512,13 @@ _trimmedview(Bpad, P, d, newindexes) = view(Bpad, newindexes...) return _trimmedview(Bpad, P, d + 1, (newindexes..., index), indexes...) end -# For faster and type-stable slicing +""" + ColonFun()(i::Int) -> Colon() + +A callable singleton that returns `Colon()` for any integer argument. +Used in place of `Colon()` directly where type-stable, composable slicing +is needed (e.g., constructing index tuples programmatically). +""" struct ColonFun end ColonFun(::Int) = Colon() @@ -399,6 +528,4 @@ function __init__() end end -include("deprecations.jl") - end diff --git a/src/deprecations.jl b/src/deprecations.jl deleted file mode 100644 index 1e6ab13..0000000 --- a/src/deprecations.jl +++ /dev/null @@ -1,5 +0,0 @@ -import Base: one -@deprecate one(::Type{NumDenom{T}}) where {T} oneunit(NumDenom{T}) -@deprecate one(p::NumDenom) oneunit(p) - -@deprecate ratio(mm::AbstractArray, args...) ratio.(mm, args...) diff --git a/test/runtests.jl b/test/runtests.jl index 61c9bac..8d79909 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,8 @@ using RegisterCore using CenterIndexedArrays, ImageCore, ImageMetadata, Interpolations using Test +using Aqua +using ExplicitImports @testset "NumDenom and arrays" begin nd = NumDenom(3.5, 10) @@ -58,17 +60,30 @@ using Test @test round(NumDenom{Int}, NumDenom(1.2, 4.8)) === NumDenom{Int}(1, 5) + # promote_rule for NumDenom + plain Number (line 80) + @test promote_type(NumDenom{Int}, Float64) == NumDenom{Float64} + + # identity convert when T matches (line 83) + @test convert(NumDenom{Float32}, NumDenom(1.5f0, 2.0f0)) === NumDenom(1.5f0, 2.0f0) + + # vararg integer dims constructor (line 225) + mm2 = MismatchArray(Float32, 3, 3) + @test axes(mm2) == Base.IdentityUnitRange.((-1:1, -1:1)) + + # ColonFun constructor (line 395) + @test ColonFun(5) === Colon() + # Finding the location of the minimum numer = [5, 4, 3, 4.5, 7] .* [2, 1, 1.5, 2, 3]' numer[1, 5] = -1 # on the edge, so it shouldn't be selected denom = ones(5, 5) mma = MismatchArray(numer, denom) - @test indmin_mismatch(mma, 0) == CartesianIndex((0, -1)) + @test argmin_mismatch(mma, 0) == CartesianIndex((0, -1)) denom = reshape(float(1:25), 5, 5) mma = MismatchArray(numer, denom) - @test indmin_mismatch(mma, 0) == CartesianIndex((0, 1)) + @test argmin_mismatch(mma, 0) == CartesianIndex((0, 1)) rat = ratio.(mma, 0.5) - @test indmin_mismatch(rat) == CartesianIndex((0, 1)) + @test argmin_mismatch(rat) == CartesianIndex((0, 1)) end @testset "mismatcharrays and separate" begin @@ -97,6 +112,20 @@ end af = float.(a) ahp = highpass(af, (2,)) @test all(x -> abs(x) < 1.0e-12, ahp) + + # in-place variants + out = zeros(Float64, 4) + highpass!(out, af, (2,)) + @test all(x -> abs(x) < 1.0e-12, out) + + af2 = copy(af) + highpass!(af2, (2,)) + @test all(x -> abs(x) < 1.0e-12, af2) + + # inf sigma (skip filtering) + out2 = zeros(Float32, 4) + highpass!(out2, a, (Inf,)) + @test out2 == Float32.(a) end @testset "PreprocessSNF" begin @@ -118,6 +147,21 @@ end Bmeta = ImageMeta(B, date = "today") @test isa(pp(Bmeta), ImageMeta) + + # PreprocessSNF on a SubArray (lines 336–338) + Bbig = fill(1000, 21, 17) + pp2 = PreprocessSNF(100, [0, 0], [Inf, Inf]) + S = view(Bbig, 2:20, 2:16) + ppS = pp2(S) + @test size(ppS) == size(S) +end + +@testset "Aqua" begin + Aqua.test_all(RegisterCore) +end + +@testset "ExplicitImports" begin + @test check_no_implicit_imports(RegisterCore) === nothing end @testset "Padding and trimming" begin @@ -140,4 +184,9 @@ end S = view(A, :, 1:3, 2) Spad = paddedview(S) @test Spad == A[:, :, 2] + + # pdindex error for unsupported index type (line 368) + A3 = reshape(1:8, 2, 2, 2) + S_bad = view(A3, [1, 2], :, 1) + @test_throws ErrorException paddedview(S_bad) end