v0.3.0: allow users to choose between CPU and GPU (#45)

Also consolidated test/Project.toml into src/Project.toml

NEWS.md

@@ -1,3 +1,9 @@
+# version 3.0
+
+## Breaking changes
+
+- RegisterQD now requires an extra step because of enhancements that support either CPU and GPU processing. For CPU processing (formerly the only option), now you must manually load the [RegisterMismatch package](https://github.com/HolyLab/RegisterMismatch.jl): `using RegisterMismatch, RegisterQD`. For GPU processing, you should instead load the [RegisterMismatchCuda package](https://github.com/HolyLab/RegisterMismatchCuda.jl): `using RegisterMismatchCuda, RegisterQD`. *Note that loading both mismatch packages in the same session will cause method conflicts.* Both mismatch packages are registered in the publicly-available [HolyLabRegistry](https://github.com/HolyLab/HolyLabRegistry), and users are advised to add that registry. 
+
 # version 0.2
 
 ## Breaking changes

Project.toml

@@ -1,6 +1,6 @@
 name = "RegisterQD"
 uuid = "ac24ea0c-1830-11e9-18d4-81f172323054"
-version = "0.2.7"
+version = "0.3.0"
 
 [deps]
 CenterIndexedArrays = "46a7138f-0d70-54e1-8ada-fb8296f91f24"
@@ -14,13 +14,17 @@ PaddedViews = "5432bcbf-9aad-5242-b902-cca2824c8663"
 QuadDIRECT = "dae52e8d-d666-5120-a592-9e15c33b8d7a"
 RegisterCore = "67712758-55e7-5c3c-8e85-dda1d7758434"
 RegisterDeformation = "c19381b7-cf49-59d7-881c-50dfbd227eaf"
-RegisterMismatch = "3c0dd727-6833-5f5d-a1e8-c0d421935c74"
+RegisterMismatchCommon = "abb2e897-52bf-5d28-a379-6ca321e3b878"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
+AxisArrays = "0.3, 0.4"
 CenterIndexedArrays = "0.2"
 CoordinateTransformations = "0.5, 0.6"
+Distributions = "0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
+ImageMagick = "0.7, 1"
+ImageMetadata = "0.9"
 Images = "0.20, 0.21, 0.22, 0.23, 0.24, 0.25"
 Interpolations = "0.12, 0.13"
 MappedArrays = "0.2, 0.3, 0.4"
@@ -32,17 +36,21 @@ RegisterDeformation = "0.3, 0.4"
 RegisterMismatch = "0.3, 0.4"
 Rotations = "0.12, 0.13, 1"
 StaticArrays = "0.11, 0.12, 1"
+TestImages = "0.5, 0.6, 1"
+Unitful = "0.17, 0.18, 1"
 julia = "1"
 
 [extras]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ImageMagick = "6218d12a-5da1-5696-b52f-db25d2ecc6d1"
 ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+RegisterMismatch = "3c0dd727-6833-5f5d-a1e8-c0d421935c74"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["Test", "ImageMagick", "TestImages", "Random", "AxisArrays", "ImageMetadata", "Unitful", "Distributions"]
+test = ["Test", "ImageMagick", "TestImages", "Random", "AxisArrays", "ImageMetadata", "Unitful", "Distributions", "LinearAlgebra", "RegisterMismatch"]

README.md

@@ -5,6 +5,9 @@
 RegisterQD performs image registration using the global optimization routine [QuadDIRECT](https://github.com/timholy/QuadDIRECT.jl).
 Unlike many other registration packages, this is not "greedy" descent based on an initial guess---it attempts to find the globally-optimal alignment of your images.
 
+To use this package, users must choose between using either the CPU or the GPU. For CPU processing, you must manually load the [RegisterMismatch package](https://github.com/HolyLab/RegisterMismatch.jl): `using RegisterMismatch, RegisterQD`. For GPU processing, you should instead load the [RegisterMismatchCuda package](https://github.com/HolyLab/RegisterMismatchCuda.jl): `using RegisterMismatchCuda, RegisterQD`. *Note that loading both mismatch packages in the same session will cause method conflicts.* Both mismatch packages are registered in the publicly-available [HolyLabRegistry](https://github.com/HolyLab/HolyLabRegistry), and users are advised to add that registry. 
+In the current absense of Github resources for GPU code, "gpu_test.jl" should be run on your personal machine as required. 
+
 This package exports the following registration functions:
 - `qd_translate`: register images by shifting one with respect to another (translations only)
 - `qd_rigid`: register images using rotations and translations

src/RegisterQD.jl

@@ -1,8 +1,8 @@
 module RegisterQD
 
 using Images, CoordinateTransformations, QuadDIRECT
-using RegisterMismatch
-using RegisterCore #just for indmin_mismatch?
+using RegisterMismatchCommon 
+using RegisterCore 
 using RegisterDeformation, PaddedViews, MappedArrays
 using Rotations
 using Interpolations, CenterIndexedArrays, StaticArrays, OffsetArrays

src/rigid.jl

@@ -187,6 +187,7 @@ function qd_rigid(fixed, moving, mxshift::VecLike, mxrot::Union{Number,VecLike};
     if initial_tfm == IdentityTransformation() || isrotation(initial_tfm.linear)
     else
         @show "WARNING: initial_tfm is not a rigid transformation"
+        # @warn "initial_tfm is not a rigid transformation"
     end
     fixed, moving = float(fixed), float(moving)
     if presmoothed

test/gpu_tests.jl

@@ -0,0 +1,161 @@
+using CUDA, RegisterMismatchCuda, RegisterQD
+using LinearAlgebra
+using ImageMagick
+using Distributions #what is this?
+using RegisterQD.StaticArrays
+using RegisterQD.Interpolations
+using RegisterQD.Images
+using RegisterQD.CoordinateTransformations
+using RegisterQD.Rotations
+using RegisterQD.OffsetArrays
+using RegisterQD.RegisterMismatchCommon
+
+using TestImages
+using Test
+
+# CUDA.GPUArrays.default_scalar_indexing[] = CUDA.GPUArrays.ScalarDisallowed
+# CUDA.GPUArrays.default_scalar_indexing[] = CUDA.GPUArrays.ScalarAllowed
+
+
+#Helper to generate test image pairs
+function fixedmov(img, tfm)
+    img = float(img)
+    img2 = warp(img,tfm)
+    inds = OffsetArrays.IdentityUnitRange.(intersect.(axes(img), axes(img2)))
+    fixed = img[inds...]
+    moving = img2[inds...]
+    return fixed, moving
+end
+
+#Helpers to convert Transformations to AffineMaps
+to_affine(tfm::Translation) = AffineMap(Matrix{Float64}(LinearAlgebra.I, length(tfm.translation), length(tfm.translation)), tfm.translation)
+to_affine(tfm::LinearMap) = AffineMap(Matrix{Float64}(LinearAlgebra.I, length(tfm.translation), length(tfm.translation)), tfm.translation)
+to_affine(tfm::AffineMap) = tfm
+
+#Helper to test that a found transform is (roughly) the inverse of the original transform. 
+function tfmtest(tfm, tfminv)
+    comp = to_affine(tfm ∘ tfminv)  #should be the identity transform
+    diagtol = 0.005
+    offdiagtol = 0.005
+    vtol = 0.1
+    @test all(x->(1-diagtol < x < 1+diagtol), diag(comp.linear))
+    @test all(x->(-offdiagtol < x < offdiagtol), comp.linear.-Matrix(Diagonal(diag(comp.linear))))
+    @test all(x-> x.<vtol, abs.(comp.translation))
+end
+
+# helper for wrapping CuArrays in Offset Arrays
+# cu_wrap(img::OffsetArray) = OffsetArray(CuArray(img.parent), img.offsets)
+# cu_wrap(img) = CuArray(img)
+
+#rigid tests 
+img = Float32.(testimage("cameraman"))
+SD = Matrix{Float64}(LinearAlgebra.I, 2, 2)
+tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.3)) #no distortion for now
+# fixed, moving = cu_wrap.(fixedmov(centered(img), tfm));
+fixed, moving = fixedmov(centered(img), tfm)
+mxshift = (100,100) #make sure this isn't too small
+mxrot = (0.5,)
+minwidth_rot = fill(0.002, 3)
+
+
+tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002);
+
+tfmtest(tfm, tform) 
+@test mm<0.001
+
+
+img3D = Float32.(testimage("mri-stack"))
+SD3D = Matrix{Float64}(LinearAlgebra.I, 3,3 )
+tfm3D = Translation(@SVector([15, 10, 2]))∘LinearMap(RotXYZ(0.05, 0.02, 0)) #no distortion for now
+fixed, moving = fixedmov(centered(img3D), tfm3D)
+mxshift = (30, 30, 5) #make sure this isn't too small
+mxrot = (0.1, 0.1, 0.1)
+minwidth_rot = fill(0.002, 3) 
+
+tform2, mm2 = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD3D, maxevals=1000, rtol=0, fvalue=0.0002);
+tfmtest(tfm3D, tform2) 
+@test mm2<0.01
+
+#Coppied over from qd_random
+@testset "QuadDIRECT tests with standard images" begin
+    img = testimage("cameraman");
+
+    #Translation (subpixel)
+    tfm = Translation(@SVector([14.3, 17.6]))
+    fixed, moving = fixedmov(img, tfm)
+    mxshift = (100,100) #make sure this isn't too small
+    tform, mm = qd_translate(fixed, moving, mxshift; maxevals=1000, rtol=0, fvalue=0.0003)
+    tfmtest(tfm, tform)
+
+    #Rigid transform
+    SD = Matrix{Float64}(LinearAlgebra.I, 2, 2)
+    tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.3)) #no distortion for now
+    fixed, moving = fixedmov(centered(img), tfm)
+    mxshift = (100,100) #make sure this isn't too small
+    mxrot = (0.5,)
+    minwidth_rot = fill(0.002, 3)
+    tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002)
+    tfmtest(tfm, tform)
+    #with anisotropic sampling
+    SD = Matrix(Diagonal([0.5; 1.0]))
+    tfm = Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.3)*SD)
+    fixed, moving = fixedmov(centered(img), tfm)
+    tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=SD, maxevals=1000, rtol=0, fvalue=0.0002)
+    tfmtest(tfm, arrayscale(tform, SD))
+
+    #Affine transform
+    tfm = Translation(@SVector([14, 17]))∘LinearMap(RotMatrix(0.01))
+    #make it harder with nonuniform scaling
+    scale = @SMatrix [1.005 0; 0 0.995]
+    SD = Matrix{Float64}(LinearAlgebra.I, 2, 2)
+    tfm = AffineMap(tfm.linear*scale, tfm.translation)
+    mxshift = (100,100) #make sure this isn't too small
+    fixed, moving = fixedmov(centered(img), tfm)
+    tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals=1000, rtol=0, fvalue=0.0002)
+    tfmtest(tfm, tform)
+
+    #with anisotropic sampling
+    SD = Matrix(Diagonal([0.5; 1.0]))
+    tfm = Translation(@SVector([14.3, 17.8]))∘LinearMap(RotMatrix(0.1)) #Translation(@SVector([14.3, 17.8]))∘LinearMap(SD\RotMatrix(0.01)*SD)
+    scale = @SMatrix [1.005 0; 0 0.995]
+    tfm = AffineMap(tfm.linear*scale, tfm.translation)
+    tfm = arrayscale(tfm, SD)
+    fixed, moving = fixedmov(centered(img), tfm)
+    tform, mm = qd_affine(fixed, moving, mxshift; SD = SD, maxevals=10000, rtol=0, fvalue=0.0002, ndmax = 0.25)
+    tform2 = arrayscale(tform, SD)
+    tfmtest(tfm, tform2)
+end #tests with standard images
+
+@testset "Quadratic interpolation (issue #7)" begin
+    samplefrom(n) = rand(Poisson(n))
+
+    img = restrict(restrict(testimage("cameraman")))[2:end-1,2:end-1]
+    # Convert to "photons" so we can mimic shot noise
+    np = 1000  # maximum number of photons per pixel
+    img = round.(Int, np.*gray.(img))
+    fixed  = samplefrom.(img)
+    moving = samplefrom.(img)
+    ff = qsmooth(fixed)
+
+    tform, mm = qd_translate(fixed, moving, (5, 5); print_interval=typemax(Int))
+    tformq, mmq = qd_translate(ff, moving, (5, 5); presmoothed=true, print_interval=typemax(Int))
+    @test all(abs.(tformq.translation) .< abs.(tform.translation))
+
+    tform, mm = qd_rigid(fixed, moving, (5, 5), 0.1; print_interval=typemax(Int))
+    tformq, mmq = qd_rigid(ff, moving, (5, 5), 0.1; presmoothed=true, print_interval=typemax(Int))
+    @test norm(tformq.linear-I) < norm(tform.linear-I)
+    @test norm(tformq.translation) < norm(tform.translation)
+
+    tform, mm = qd_affine(fixed, moving, (5, 5); print_interval=typemax(Int))
+    tformq, mmq = qd_affine(ff, moving, (5, 5); presmoothed=true, print_interval=typemax(Int))
+    @test mmq < mm
+
+    # Test that we exactly reconstruct `qsmooth` with `presmoothed=true`
+    tformq, mmq = qd_translate(ff, fixed, (5, 5); presmoothed=true, print_interval=typemax(Int))
+    @test all(iszero, tformq.translation)
+    @test mmq < 1e-10
+    tformq, mmq = qd_rigid(ff, fixed, (5, 5), 0.1; presmoothed=true, print_interval=typemax(Int))
+    @test mmq < 1e-8
+    tformq, mmq = qd_affine(ff, fixed, (5, 5); presmoothed=true, print_interval=typemax(Int))
+    @test mmq < 1e-6   # on 32-bit systems this can't be 1e-8, not quite sure why
+end

test/qd_random.jl

@@ -4,7 +4,6 @@ using RegisterQD.Interpolations
 using RegisterQD.Images
 using RegisterQD.CoordinateTransformations
 using RegisterQD.Rotations
-using RegisterQD.RegisterMismatch
 using RegisterQD: _abs2
 using Random
 

test/qd_standard.jl

@@ -8,7 +8,6 @@ using RegisterQD.Images
 using RegisterQD.CoordinateTransformations
 using RegisterQD.Rotations
 using RegisterQD.OffsetArrays
-using RegisterQD.RegisterMismatch
 
 using Test, TestImages
 using Random

test/runtests.jl

@@ -1,5 +1,5 @@
 using ImageMagick
-using RegisterQD
+using RegisterQD, RegisterMismatch
 using Test
 
 include("util.jl")