GitHub - HolyLab/RegisterDeformation.jl: Representation of spatial deformations (it’s linked from here)
1 Like
Sorry to bug you yet another time. I had a look at the deformation package and it seems that that focuses on the part where one applies the transformation and not the registration part of it. Or is that the part where BlockRegistration plugs in and solves the inverse problem?
Sorry, just tried it our and I seem to got it. The ϕs in the cookbook is actually the deformed grid and your link just was the documentation on the GridDeformation type.
It’s a whole suite of tools. You asked about the displacement field, which is implemented and documented in RegisterDeformation. The other packages are built on top of that, and optimize the displacement field to approximately minimize the mean-square mismatch. The order in this list is roughly the “stacking order”: RegisterOptimize is the top-level package which minimizes the mismatch as a function of the displacement field.
Just step your way through the demo in the BlockRegistration documentation, I think it will become far clearer then.
2 Likes
yes, thanks. Yesterday I had issues installing the HolyLabRegistry and was not able to step through. Now it works fine!
1 Like
Great. Keep in mind that most of BlockRegistration was designed to handle the multi-terabtye datasets of light sheet microscopy, and it compromises on quality for the sake of speed. The key tradeoff is best documented in Guaranteed improvement? · Issue #8 · HolyLab/RegisterOptimize.jl · GitHub. If you want high quality, expect to have to polish its solutions a bit (and there’s a whole documentation page on what you can do to polish them, though it doesn’t mention options like “use a different package” e.g. ANTSRegistration.jl). On the other hand, it has some fairly innovative features that mimic some of the benefits of global optimization, so it might be quite useful for getting into a good configuration for polishing.
1 Like
I am not sure I fully understand what you mean by quality. If I understand correctly you are using a low resolution deformation grid, which inherently leads to regularization since the deformation cannot be arbitrarily fine. But isn’t it possible to change the blocksize and in that way control the resolution of the deformation?
As stated in that issue link, the modeling of the effect of the deformation is “prospective” rather than “retrospective”: it tries to predict the effect of applying the deformation using the fourier shift theorem without actually calculating it precisely. That introduces some additional inaccuracy but makes it blazingly fast (for a voxelwise mean square error algorithm).
While the reasons are not fully described in the Documenter documentation, the issue of accuracy is extensively discussed in Improving alignment · BlockRegistration.jl
ok, so since time is no issue in my application it seems I want to go with the RegisterHindsight approach. I will give this a spin and show how far I get. After going through the cookbook I am pretty impressed what I have seen so far! Well done.
Just out of interest: Is there a mathematical/algorithmical paper out there on the approach you are using here? Or is that new stuff?
1 Like
It’s new. I always seem to prefer to write software over papers.
1 Like
The community is appreciating that! 
My colleagues in my group are always joking that there have to be multiple Tim Holys since you have created such a multitude of awesome Julia packages.
3 Likes
@tim.holy, thanks for the great work! Some day I’d like to have a really thorough comparison of brain registration using a Julia based pipeline compared to an ANTs pipeline. I don’t know if I’ll ever get the time to take the lead on something like that, but maybe I can convince someone in the biophysics department to do it at some point.
@tim.holy We are now frequently using the QuadDirect algorithm in our lab. However, I got a strange issue on one computer, after I ran the update-command for the packages. As far as I can remember, RegisterQD was not among the updated packages. Neither did I update my julia installation (still version 1.0.5).
Minimal (non-)working sample
julia> using RegisterQD
julia> pseudofixed=rand(10,10)
julia> tfm, mm = qd_rigid(pseudofixed, pseudofixed, (5,5), 0.1, print_interval=typemax(Int))
On a machine that was not updated, the last command generates the expected result
(AffineMap([1.0 -0.0; 0.0 1.0], [0.0, 0.0]), 0.0)
However, the updated machine first throws a deprecation warning and then an error
julia> tfm, mm = qd_rigid(pseudofixed, pseudofixed, (5,5), 0.1, print_interval=typemax(Int))
┌ Warning: `CenterIndexedArray(::Type{T}, dims) where T` is deprecated, use `CenterIndexedArray{T}(undef, dims...)` instead.
│ caller = Type at RegisterCore.jl:360 [inlined]
└ @ Core ~/.julia/packages/RegisterCore/0P3yV/src/RegisterCore.jl:360
ERROR: UndefVarError: shiftrange not defined
Stacktrace:
[1] #mismatch#5(::Symbol, ::Function, ::Type{Float64}, ::SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::SubArray{Float64,2,OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::Tuple{Int64,Int64}) at /home/moses/.julia/packages/RegisterMismatch/kKJng/src/RegisterMismatch.jl:120
[2] #mismatch at ./subarray.jl:0 [inlined]
[3] #mismatch#1 at /home/moses/.julia/packages/RegisterMismatchCommon/Q3oNz/src/RegisterMismatchCommon.jl:15 [inlined]
[4] #mismatch at ./none:0 [inlined]
[5] #best_shift#1(::Symbol, ::CoordinateTransformations.IdentityTransformation, ::Function, ::SubArray{Float64,2,Array{Float64,2},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::SubArray{Float64,2,OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}},Tuple{UnitRange{Int64},UnitRange{Int64}},false}, ::Tuple{Int64,Int64}, ::Float64) at /home/moses/.julia/packages/RegisterQD/ZBH17/src/util.jl:29
[6] #best_shift at ./none:0 [inlined]
[7] #rigid_mm_fast#14(::CoordinateTransformations.IdentityTransformation, ::Function, ::Array{Float64,1}, ::Tuple{Int64,Int64}, ::Array{Float64,2}, ::Array{Float64,2}, ::Float64, ::LinearAlgebra.UniformScaling{Bool}) at /home/moses/.julia/packages/RegisterQD/ZBH17/src/rigid.jl:74
[8] #rigid_mm_fast at ./none:0 [inlined]
[9] f at /home/moses/.julia/packages/RegisterQD/ZBH17/src/rigid.jl:92 [inlined]
[10] split!(::QuadDIRECT.Box{Float64,1}, ::getfield(RegisterQD, Symbol("#f#17")){LinearAlgebra.UniformScaling{Bool},CoordinateTransformations.IdentityTransformation,Float64,Array{Float64,2},Array{Float64,2},Tuple{Int64,Int64}}, ::Array{Float64,1}, ::Int64, ::Array{Float64,1}, ::Float64, ::Float64, ::Float64, ::Float64) at /home/moses/.julia/packages/QuadDIRECT/uCBvD/src/algorithm.jl:48
[11] _init(::Function, ::Array{Float64,1}, ::Tuple{Array{Float64,1}}, ::Array{Float64,1}, ::Array{Float64,1}) at /home/moses/.julia/packages/QuadDIRECT/uCBvD/src/algorithm.jl:31
[12] init at /home/moses/.julia/packages/QuadDIRECT/uCBvD/src/algorithm.jl:20 [inlined]
[13] #analyze#29(::Base.Iterators.Pairs{Symbol,Any,NTuple{5,Symbol},NamedTuple{(:minwidth, :print_interval, :maxevals, :atol, :rtol),Tuple{Array{Float64,1},Int64,Float64,Int64,Float64}}}, ::Function, ::Function, ::Tuple{Array{Float64,1}}, ::Array{Float64,1}, ::Array{Float64,1}) at /home/moses/.julia/packages/QuadDIRECT/uCBvD/src/algorithm.jl:515
[14] (::getfield(QuadDIRECT, Symbol("#kw##analyze")))(::NamedTuple{(:minwidth, :print_interval, :maxevals, :atol, :rtol),Tuple{Array{Float64,1},Int64,Float64,Int64,Float64}}, ::typeof(QuadDIRECT.analyze), ::Function, ::Tuple{Array{Float64,1}}, ::Array{Float64,1}, ::Array{Float64,1}) at ./none:0
[15] #_analyze#7(::Base.Iterators.Pairs{Symbol,Any,NTuple{5,Symbol},NamedTuple{(:minwidth, :print_interval, :maxevals, :atol, :rtol),Tuple{Array{Float64,1},Int64,Float64,Int64,Float64}}}, ::Function, ::Function, ::Array{Float64,1}, ::Array{Float64,1}) at /home/moses/.julia/packages/RegisterQD/ZBH17/src/util.jl:207
[16] (::getfield(RegisterQD, Symbol("#kw##_analyze")))(::NamedTuple{(:minwidth, :print_interval, :maxevals, :atol, :rtol),Tuple{Array{Float64,1},Int64,Float64,Int64,Float64}}, ::typeof(RegisterQD._analyze), ::Function, ::Array{Float64,1}, ::Array{Float64,1}) at ./none:0
[17] #qd_rigid_coarse#16(::LinearAlgebra.UniformScaling{Bool}, ::CoordinateTransformations.IdentityTransformation, ::Float64, ::Base.Iterators.Pairs{Symbol,Int64,Tuple{Symbol},NamedTuple{(:print_interval,),Tuple{Int64}}}, ::Function, ::Array{Float64,2}, ::Array{Float64,2}, ::Tuple{Int64,Int64}, ::Array{Float64,1}, ::Array{Float64,1}) at /home/moses/.julia/packages/RegisterQD/ZBH17/src/rigid.jl:95
[18] #qd_rigid_coarse at ./none:0 [inlined]
[19] #qd_rigid#20(::LinearAlgebra.UniformScaling{Bool}, ::Array{Float64,1}, ::Float64, ::CoordinateTransformations.IdentityTransformation, ::Int64, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Array{Float64,2}, ::Array{Float64,2}, ::Tuple{Int64,Int64}, ::Float64) at /home/moses/.julia/packages/RegisterQD/ZBH17/src/rigid.jl:189
[20] (::getfield(RegisterQD, Symbol("#kw##qd_rigid")))(::NamedTuple{(:print_interval,),Tuple{Int64}}, ::typeof(qd_rigid), ::Array{Float64,2}, ::Array{Float64,2}, ::Tuple{Int64,Int64}, ::Float64) at ./none:0
[21] top-level scope at none:0
As I did not update julia itself, I do not understand that I get a deprecation. As far as I can see, RegisterQD has not seen any updates recently, either.
Thanks for the report, this turned out to be quite the mess. As you presumably know, RegisterQD uses HolyLabRegistry in addition to General, as it requires a number of other lab packages. This turned out to be a problem stemming from the bad old days before we got @GunnarFarneback’s wonderful LocalRegistry, when our releases were crafted by editing HolyLabRegistry manually. When I switched to LocalRegistry, a mistake I’d made months prior led me to inadvertently overwriting version 0.2.1 of RegisterMismatchCommon.
I believe I’ve untangled this successfully, both for older releases and the latest. You said you’re using Julia 1.0.5; the more recent dependencies won’t work on Julia 1.0, but I believe I sorted out all the retrospective bounds. I have tested installation from a “clean” package depot on both Julia 1.4 and Julia 1.0, and it all seems to work. Please try updating your packages and let me know if it works. (If you get Pkg errors try removing and then re-adding packages.)
2 Likes
Thank you for the fast solution. It works again.
Yes, I am still on version 1.0.5, because this is currently marked as LTS.
@tim.holy thanks for this really cool package! I’d like to use RegisterQD to rigidly transform my data to a common brain atlas, however the pixel units are modestly different (0.88um vs 1.00um). Right now I use CMTK using commands described in Figure S2 to perform a rigid followed by non-rigid step, but I’m attracted to the global optimality that your algorithm strives for as well as easy integration with Julia code.
I am using AxisArrays with properly specified units but I think RegisterQD assumes units are consistent between fixed and “moving”? What’s the recommended approach here–resample the data so axes have same units? I’m planning on just using imresize unless you would recommend a different way.
Thank you!
Edit: I think things are working properly with qd_affine although need to test on higher quality data ;). One thing I’m confused by is what the second value in the returned tuple is, so for example for this return value:
(AffineMap([1.2079968439617563 -0.0601021283777167; -0.062482595377332216 1.2539712136955676], [-160.24406067251462, -428.64275616465136]), 0.8871258679764823), I understand that the first 2x2 matrix is rotation / stretch, the second 2 vector is transaltion, but what is the third number (0.88)?
You shouldn’t need to resample, the package is already set up to handle that. Just set the SD argument to qd_rigid. You can probably use ps = pixelspacing(img); SD = ps./maximum(ps) to get a numerical value free of physical units.
That’s the mm (mismatch) that is described at some places and not others (oops) in the docs. It’s essentially the average square difference, considering only the “inbounds” voxels (i.e., only those voxels for which the fixed and aligned moving images overlap).
Thanks! For SD, should I use the pixel spacing from “fixed” or “moving”? Ie the atlas is (0.78um, 0.78um, 2um) and my data is (1um, 1um, 1um)
Oh, sorry, I misunderstood; I thought they were the same but not isotropic. (The definition of a “rigid transformation” needs to change if space is not sampled isotropically.) I’m afraid you’ll have to resample to make them agree with one another, but they don’t have to be isotropic.
@tim.holy ANTsRegistration.jl is a great package, thanks for your contribution. I would have the following question: while imgw = register(fixed, moving, pipeline; kwargs...) allows to obtain the registered image as output, is there a simple way to get the corresponding geometrical-transformation information (which would be akin to the imregtform function in Matlab)? In case this is relevant to a potential solution, I am specifically interested in the simplest case of rigid registration (more specifically translation). Any help would be very much appreciated.
1 Like