How can I register/align shifted image stacks?

Good evening.

I have the problem that in many microscopy experiments the sample drifts over time (mainly because a buffer solution is flowing continuously over the sample).

I am looking for a way to register/align the image stack to remove or correct for this drift (@phlavenk mentioned this issue here).

I currently do this registration in Fiji/ImageJ using the StackReg-plugin with either the translation or Rigid Body algorithms.

My questions are:

  1. Is this possible with Images.jl? Probably not (the Registration section is empty here).
  2. Has anybody ever performed image/stack registration in Julia and can tell me how to do it?

Thank you very much in advance.

Moses

might be useful. Let us know how it works out.

1 Like

So perhaps the Registration section should be updated?

I’d also give a try to SubpixelRegistration.jl. I’m very happy with this package for drift correction in electron microscopy. In my case I’ve also used StackReg in ImageJ and i used to externally call it from my julia code, but the native SubpixelRegistration just works great.

I don’t remember the reason why this package was not transfered to JuliaImages (licencing?), but maybe just having a link on the Images.jl webpage would be great for better discoverability.

Thank you for pointing me to this one. However, I have trouble installing it. The package is obviously not officially registered, so I was not able to add it directly.
When I tried add "https://github.com/HolyLab/RegisterQD.jl.git" I got the following error:

Cloning git-repo `https://github.com/HolyLab/RegisterQD.jl.git`
  Updating git-repo `https://github.com/HolyLab/RegisterQD.jl.git`
 Resolving package versions...
ERROR: Unsatisfiable requirements detected for package RegisterMismatch [3c0dd727]:
 RegisterMismatch [3c0dd727] log:
 ├─RegisterMismatch [3c0dd727] has no known versions!
 └─restricted to versions * by RegisterQD [ac24ea0c] — no versions left
   └─RegisterQD [ac24ea0c] log:
     ├─possible versions are: 0.2.0 or uninstalled
     └─RegisterQD [ac24ea0c] is fixed to version 0.2.0

Could you please tell me, how to install it? Am I missing some dependencies?
I am using Julia 1.0.4 with Atom/Juno on Ubuntu 18.04 LTS.

Thank you, too. I could at least install it, but am at a loss of how to use it. The readme is not really clear to me. When I copy the example from the readme to my julia (version 1.0.4), I get

ERROR: MethodError: no method matching setindex_shape_check(::Int64, ::Int64, ::Int64, ::Int64, ::Int64)
Closest candidates are:
  setindex_shape_check(::AbstractArray, ::Integer...) at indices.jl:179
  setindex_shape_check(::AbstractArray{#s57,1} where #s57, ::Integer, ::Integer) at indices.jl:221
  setindex_shape_check(::AbstractArray{#s57,2} where #s57, ::Integer, ::Integer) at indices.jl:225
  ...
Stacktrace:
 [1] macro expansion at ./multidimensional.jl:649 [inlined]
 [2] _unsafe_setindex!(::IndexLinear, ::Array{Float64,4}, ::Int64, ::UnitRange{Int64}, ::UnitRange{Int64}, ::UnitRange{Int64}, ::Int64) at ./multidimensional.jl:644
 [3] _setindex!(::IndexLinear, ::Array{Float64,4}, ::Int64, ::UnitRange{Int64}, ::UnitRange{Int64}, ::Vararg{Union{Real, AbstractArray},N} where N) at ./multidimensional.jl:639
 [4] setindex!(::Array{Float64,4}, ::Int64, ::UnitRange{Int64}, ::UnitRange{Int64}, ::UnitRange{Int64}, ::Vararg{Any,N} where N) at ./abstractarray.jl:998
 [5] top-level scope at none:0

Would you mind to share some dummy-code for the registration of an image (tiff)?

The following worked for me from a “clean” Julia environment:

pkg> up    # this is necessary if you haven't cloned the General registry yet
pkg> registry add https://github.com/HolyLab/HolyLabRegistry.git
pkg> add RegisterQD

If all you ever need to do is translation, then I imagine SubpixelRegistration should be fine. RegisterQD does translation, but its novelty (if you can call it that) is pretty reliable global rigid (rotation + translation) registration, and somewhat reliable affine registration.

Thank you very much for the hint of adding the HolyLabRegistry.

For anybody else being on a version before 1.1 (like me):
You need to clone the registry to your local drive first. Otherwise, the name of the Package is not recognised. It is described here.
If I got it right registry add is only available from version 1.1 onwards. (In my case it registry throws an error for an unknown command.)

@tim.holy Could you possibly point me to a description of the syntax or a tutorial of RegisterQD.jl?

Using Atom/Juno I found a documentation of the commands, like this one:
tform, mm = qd_rigid(fixed, moving, mxshift, mxrot; SD=I, minwidth_rot=default_minwidth_rot(fixed, SD), thresh=thresh, initial_tfm=IdentityTransformation(), kwargs...)

However, this is quite confusing to me.

  1. Which of the options or variables is the input image stack and which format is expected?
    (EDIT: O.K., I think fixed is the reference image and moving the one I want to align with the reference.)
  2. Can I directly asign an image (tiff) file to fixed and moving?
  3. tform is possibly the transformed result, but what is mm?
  4. I need to run through the stack slice by slice, correct?

Do you by any chance have some sample/demo code you could share?

For you and others: My current workflow in Fiji/ImageJ looks like this:

  1. I have two timelapse stacks, representing two channels.
  2. I align the brighter stack with MultiStackReg (a plugin that basically calls StackReg, but is able to save the transformation matrix) using RigidBody.
  3. I apply the transformation matrix saved in step 2 two the second channel stack. This is necessary, because the lower signal/noise ratio of this channel usually results in completely different alignments, if I call the StackReg plugin for this stack.

I would like to transfer this to Julia for some additional analysis.

Did you do ?qd_rigid? There is quite a lot of help in the functions of the package. Though it looks like some things are still not documented, e.g. mm is the mismatch, essentially the mean-square difference.

Here’s an example using the “mri” image in TestImages.jl. (You don’t expect alignment of these images, this is just for demonstration purposes.)

julia> using TestImages, RegisterQD

julia> mri = testimage("mri");

julia> fixed = mri[:,:,1];

julia> for i = 2:size(mri,3)
           moving = mri[:,:,i]
           tfm, mm = qd_rigid(fixed, moving, (5,5), 0.1, print_interval=typemax(Int))
           @show i tfm mm
       end
i = 2
tfm = AffineMap([0.9999999989634472 -4.5531370764914994e-5; 4.5531370764914994e-5 0.9999999989634472], [0.5040017264556756, 0.2990075519981027])
mm = 0.09511831219416511
i = 3
tfm = AffineMap([0.9999994381696169 -0.0010600285140643243; 0.0010600285140643243 0.9999994381696169], [2.5874336682326287, 0.4090526341672306])
mm = 0.2142679963090935
i = 4
tfm = AffineMap([0.9999981403042801 -0.0019285714871368261; 0.0019285714871368261 0.9999981403042801], [4.830450516702656, -0.7747954743088441])
mm = 0.2715531457442805

and so on. Use warp to apply a transformation to the image:

julia> using Images, ImageView

julia> moving = mri[:,:,5];

julia> tfm, mm = qd_rigid(fixed, moving, (5,5), 0.1, print_interval=typemax(Int))
(AffineMap([0.9996029561869612 0.02817676316201408; -0.02817676316201408 0.9996029561869612], [5.002439810122244, 5.099462334279371]), 0.2706807127431129)

julia> movw = warp(moving, tfm, axes(fixed));

julia> imshow(colorview(RGB, fixed, moving, zeroarray); name="original")

julia> imshow(colorview(RGB, fixed, movw, zeroarray); name="rigid");

original rigid
(Left is original, right is rigid-warping.)

1 Like

I did not exactly call ?qd_rigid, but used the Documentation-pane provided by Atom/Juno. The content is the same as the REPL help command.

Thank you for the sample. Things are much clearer now. I hope I can test it tomorrow.

So, I finally managed to test it properly on one stack I have. The result is amazing, at least in comparison to Fiji’s StackReg-Plugin. However, it currently takes much more time (Fiji ~ 1.5 minutes, Julia ~ 70 minutes). But this can be due to the @show macro inculded in the sample code.

@tim.holy: Is there any included command to actually save the tfm-values for each slice to a matrix/file? Otherwise, I will just write them to a file individually, so that I can load them later again.

1 Like

Glad you are pleased with the results. The amount of time it takes is subject to quite a few factors, and you may be able to reduce it. See, for example, https://github.com/HolyLab/RegisterQD.jl/issues/21#issuecomment-537295913.

You can use JLD, JLD2, or Julia’s own serializer to save the transformation arrays. Or you can convert them to affine transformations and just store the matrices using CSV, HDF5, MAT, or some other format.

I currently have Julia running in a VM. That’s probably the main reason for the slow performance. I will check the issues you pointed out. However, I don’t mind if the script needs more time, as long as the result is fine.

Two more questions:

  1. Using qd_rigid I get a result that look like this:
    tfm = AffineMap([1.0 -6.76146e-5; 6.76146e-5 1.0], [0.441512, -0.42513])
    Am I right to assume that the first 4 values define the rotation, while the second set of two values defines the translation? According to the help, the center coordinate would usually be 0,0
  2. Would CoordinateTransformations.jl be the correct choice for your suggested conversion? (Sorry, I/we have never done something like this before.)

I currently have Julia running in a VM. That’s probably the main reason for the slow performance.

Might be, might not. RegisterQD uses QuadDIRECT to do global optimization. That may account for the better results. But it’s also much more expensive than local optimization. Moreover, there is no good criterion for termination, short of being able to prove there isn’t a better minimum (which can require quite exhaustive exploration). If you need to speed it up, termination would be the first thing I’d look at. But it sounds like it’s good enough as-is for your needs.

Regarding your questions, yes to both.

All of this is super cool, I somehow missed these two packages. I can’t help thinking that we truly have all the tools we might ever need to build some kind of camera calibration package… All we really need is time.

@tim.holy do you have any metrics of how this performs on brains compared to different software or is that a forthcoming paper? Ideally we would have a comparison similar to what’s done with MindBoggle.

@tim.holy: Do you have any hint for an image registration package (in 3D) for non-rigid image registration?

Here are some options:

2 Likes

Does blockRegistration offer an option to access the displacement field or is that fundamentally different?