Mismatch increases in this case

[Cody-G]

For some time I've been trying to get good registrations when there is inconsistency in sampling locations along the z-axis. This is common with high-speed OCPI recordings. We're doing a lot to handle this at the level of the microscope hardware and acquisition software, but it would be good if our registration algorithm could clean up the last bit of it. I wasn't satisfied with the results so i made a synthetic image volume and simulated the situation (script below). I'm seeing that the loss (aka mismatch) actually increases when registering these volumes. I didn't yet track down why that is, and this can probably be triggered by a simpler test case, but here's a demo that I know demonstrates the issue (note that it requires a utility package BinaryVolumes, another HolyLab repo.)

using BlockRegistration, BlockRegistrationScheduler, RegisterDeformation, RegisterMismatch
using AffineTransforms
using ImageView2, Images
using FixedSizeArrays
using JLD
using Interpolations

using BinaryVolumes

#utility function for running Apertures registration on a pair of images
function _reg_def(fixed, moving, knots, mxshift; lambda = 0.01)
    λ = lambda
    algorithm = Apertures[Apertures(fixed, knots, mxshift, λ, identity; pid=myid(), correctbias=false) for i = 1:1]
    mon = monitor(algorithm, (), Dict{Symbol,Any}(:u=>ArrayDecl(Array{Vec{3,Float64},3}, ([convert(Int,x.len) for x in knots]...))))
    mon[1][:mismatch] = 0.0
    fileout = tempname()
    @time driver(fileout, algorithm, moving, mon)
    u = JLD.load(fileout, "u")
    mm = JLD.load(fileout, "mismatch")
    return u, mm
end

#add jitter in z sampling location, simulating inconsistencies in piezo position when using OCPI under certain conditions
function z_jitter{T}(img::Array{T,3}, npix::Float64)
    @assert npix < 0.5 && npix >= 0.0
    etp = extrapolate(interpolate(img, BSpline(Linear()),OnGrid()), Flat())
    out = zeros(eltype(img), size(img))
    z_def = Float64[]
    for i in 1:size(img,3)
        r = (rand()*2*npix)-npix
        push!(z_def, r)
        for y = 1:size(img,2)
            for x = 1:size(img,1)
                out[x,y,i] = etp[x,y,i+r] #Interpolations doesn't currently support slicing
            end
        end
    end
    return out, z_def
end
    

#create template ellipsoid volume
radii = [10;20;15]*2
num_pad = 20
line_spacings = [6;6;6];

vmask = ellipsoid_volume(radii; fill_val=true, num_pad = num_pad);
vbound = find_boundaries(vmask; val=true);

#draw grid lines. Make density higher on one side so that there's no
#ambiguity in orientation.  Only need to do this for n-1 dimensions
sz = [size(vmask)...];
vgrid = deepcopy(vbound);
for i = 1:length(line_spacings)
    line_idxs_a = floor(Int,linspace(1,sz[i]/2, floor(Int,sz[i]/line_spacings[i])));
    line_idxs_b = floor(Int,linspace(sz[i]/2,sz[i], floor(Int, sz[i]/(1.5*line_spacings[i]))))[2:end];
    cond = x->x==true
    vgrid = draw_lines!(vgrid, vmask, i, line_idxs_a, cond; fill_val=true);
end
vgrid = map(Float64, vgrid);
fixed = vgrid

moving, z_def = z_jitter(vgrid, 0.49); #plus or minus a half frame of jitter

mxshift = (3,3,3)
gridsize = (2,2,2) 
knots = map(d->linspace(1,size(fixed,d),gridsize[d]), (1:ndims(fixed)...))

u, mm = _reg_def(fixed,moving,knots,mxshift;lambda=0.0);
d = griddeformations(u, knots)
warped = warp(moving, d[1])

#mismatch is worse after registration
imshow(colorview(RGB, moving, zeroarray, fixed))
imshow(colorview(RGB, warped, zeroarray, fixed))
@show ratio(mismatch0(moving,fixed), NaN) 
@show mm

[Cody-G]

Note that the mismatch does decrease when the grid size is made large in the z dimension (as it should be if one wants to correct this kind of jitter). I just discovered this by mistake, and it's concerning to me that it's choosing a solution that increases the mismatch.

Even with a full-sized grid spacing in z it has trouble finding the right deformation, but it's not too far off. I'm looking into improving that now.

[Cody-G]

It's pretty complicated to understand what's happening in that test case, so I've been looking for simpler ways to trigger strange behavior. I've found some surprising cases with 1-D registration. Again I'm creating the moving image by jittering the sampling locations in the one dimension. I find that

• The monitor'ed mismatch is negative in this case
• The mismatch of the warped image is actually worse than my own educated guess, so we're not finding a very good solution.

Here's an example script. If you run it a few times you'll notice that my guess is often better than the algorithm's, and the monitored mismatch is almost always negative. Could this be a problem with the way we're interpolating the mismatch?

using BlockRegistration, BlockRegistrationScheduler, RegisterDeformation, RegisterMismatch
using ImageView2, Images
using FixedSizeArrays
using JLD
using Interpolations


#utility function for running Apertures registration on a pair of images
function _reg_def(fixed, moving, knots, mxshift; lambda = 0.01)
    λ = lambda
    algorithm = Apertures[Apertures(fixed, knots, mxshift, λ, identity; pid=myid(), correctbias=false) for i = 1:1]
    mon = monitor(algorithm, (), Dict{Symbol,Any}(:u=>ArrayDecl(Array{Vec{1,Float64},1}, ([convert(Int,x.len) for x in knots]...))))
    mon[1][:mismatch] = 0.0
    fileout = tempname()
    @time driver(fileout, algorithm, moving, mon)
    u = JLD.load(fileout, "u")
    mm = JLD.load(fileout, "mismatch")
    return u, mm
end

#add jitter in sampling location, simulating inconsistencies in piezo position when using OCPI under certain conditions
function jitter{T}(img::Array{T,1}, npix::Float64)
    @assert npix < 0.5 && npix >= 0.0 #don't tear space
    etp = extrapolate(interpolate(img, BSpline(Linear()),OnGrid()), Flat())
    out = zeros(eltype(img), size(img))
    z_def = Float64[]
    for i in 1:length(img)
        r = (rand()*2*npix)-npix
        push!(z_def, r)
        out[i] = etp[i+r]
    end
    return out, z_def
end
    
fixed = zeros(6)
fixed[3:4] = 1.0

moving, z_def = jitter(fixed, 0.45);

mxshift = (1)
gridsize = (6) 
knots = map(d->linspace(1,size(fixed,d),gridsize[d]), (1:ndims(fixed)...))
u, mm = _reg_def(fixed,moving,knots,mxshift;lambda=0.0);
@show mm #the monitored mismatch is negatve.  red flag?
d = griddeformations(u, knots)
warped = warp(moving, d[1])

ideal_u = fill(0.0, size(u))
ideal_u[1,2,1] = -0.9
ideal_u[1,5,1] = 0.9

ideal_d = griddeformations(ideal_u, knots)
ideal_warped = warp(moving, ideal_d[1]);
@show ratio(mismatch0(moving,fixed), NaN)
@show ratio(mismatch0(warped,fixed), NaN) #this time it's not negative
@show ratio(mismatch0(ideal_warped,fixed), NaN) #our guessed answer is usually better

[timholy]

Thanks for the very helpful issue reports. I've started with the last report, the 1d example. I can trace the negative mismatch to a simple source: In your example, the mismatch at certain grid points looks like

z = NumDenom(0.0,1)
o = NumDenom(1.0,1)
mm = CenterIndexedArray([z, z, o])

and we can reproduce the negative penalty with a simple analog:

julia> using Interpolations
                                                                                                                                                                            
julia> penalty = [0.0, 0.0, 1.0]
3-element Array{Float64,1}:                                                                                                                                                 
 0.0                                                                                                                                                                                         
 0.0                                                                                                                                                                                         
 1.0                                                                                                                                                                                         
                                                                                                                                                                                             
julia> itp = interpolate!(penalty, BSpline(Quadratic(InPlace())), OnCell())
3-element Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0}:                  
  0.0                                                                                                                                                                                        
 -2.77556e-17
  1.0        

julia> itp[1.51]
-0.08792000000000001

You get the same thing from

itp = CenterIndexedArray(interpolate!([z, z, o], BSpline(Quadratic(InPlace())), OnCell()))
RegisterPenalty.vecindex(itp, Vec(-0.49))

In a certain sense, the returned value is correct: if you interpolate a quadratic through those three points, the minimum will occur between the first and second point, and the value will be negative.

However, it's also not hard to understand how this could be problematic. Linear interpolation wouldn't suffer from the same problem, but the issue there is that (for efficiency) we really need the gradient of the penalty during optimization, and you don't have a smooth gradient unless you use at least quadratic interpolation.

Bottom line, I'm not quite sure what to do about this. I almost wonder if we should switch to a completely different optimization algorithm and focus only on integer-valued shifts.

[timholy]

One option might be to do linear interpolation and use a subgradient method. I might be tempted to move so that the maximum step is 0.1 pixel for any grid point, and just stop when the resulting penalty increases. I'll try cooking something up along those lines.

[Cody-G]

That would make sense to me! I also thought we would have to keep the step size small if we went with something like that.

There's one other thing I should bring up before we proceed. I think there could be a fundamental problem with interpolating the mismatch. Consider fixed = [0.0; 0.0; 1.0; 1.0; 0.0], and a moving image slid half a pixel to the left moving = [0.0; 0.5; 1.0; 0.5; 0.0]. So the ideal shift of the moving image would be 0.5. If we evaluate the mismatch at integer offsets we get equal values at offsets of 0 and 1. Moreover both of the values we get are kind of mediocre compared to the perfect 0 we could get if we evaluated the mismatch directly at 0.5 (assuming linear interpolation). Instead we get a range of equally-scored solutions from 1 to 0. Edit: Well actually if we continue to use quadratic interpolation I guess they aren't equal, but I don't see a guarantee that the quadratic would help and not hurt.

So my question is, how hard would it be to move the interpolation upward into the mismatch calculation itself? I think I don't understand Fourier methods well enough to know if that makes sense. Alternatively I think it may make sense to upsample the images before running registration, though that may get expensive...

[timholy]

Keep in mind that a linear interpolation of [0.0, 0.5 0.1, 0.5, 0.0] at half-grid points is [0.25, 0.75, 0.75, 0.25]. There is no way to reconstruct the original image by interpolation.

[Cody-G]

Oh yes you are right! So we wouldn't get a perfect answer like I said. But we would still get a lower mismatch value than at the integer locations.

[timholy]

True. Interestingly, for that example, quadratic interpolation would suggest the right location for the minimum mismatch---the mismatch is symmetric around a half-grid point, so quadratic interpolation would place the minimum at 0.5.

[Cody-G]

But it seems to me like the quadratic may be hurting us as often as it helps. I ran registration just to verify that it finds the right minimum in that case, and surprisingly it doesn't:

using BlockRegistration, BlockRegistrationScheduler, RegisterDeformation, RegisterMismatch
using FixedSizeArrays
using JLD 
using Interpolations

#utility function for running Apertures registration on a pair of images
function _reg_def(fixed, moving, knots, mxshift; lambda = 0.01)
    �� = lambda
    algorithm = Apertures[Apertures(fixed, knots, mxshift, ��, identity; pid=myid(), correctbias=false) for i = 1:1]
    mon = monitor(algorithm, (), Dict{Symbol,Any}(:u=>ArrayDecl(Array{Vec{1,Float64},1}, (1,))))
    mon[1][:mismatch] = 0.0 
    fileout = tempname()
    @time driver(fileout, algorithm, moving, mon)
    u = JLD.load(fileout, "u")
    mm = JLD.load(fileout, "mismatch")
    return u, mm
end
fixed = [0.0;0.0;1.0;1.0;0.0]
moving = [0.0;0.5;1.0;0.5;0.0]
mxshift = (1) 
gridsize = (1) 
knots = ([3],)
u, mm = _reg_def(fixed,moving,knots,mxshift;lambda=0.0);
@show u   # == -0.336435

[Cody-G]

Sorry the last line should have been @show u, I updated the code to reflect that.

[timholy]

One issue I noted is that mxshift = (1,) will in fact confine the maximum deformation to 0.5 or smaller (it computes mismatch for shifts up to 1, but to evaluate the penalty you have to interpolate and for quadratic interpolation that requires that it be at least a half-step interior to the boundary). So your ideal_u is not quite a fair comparison.

But I have one or more fixes coming (still testing...).

[timholy]

Worth pointing out something else in this context. Setting lambda = 0 can be useful for debugging, and indeed this detected a bug in initial_deformation (fix coming). However, weird stuff can happen. For example, in the code I'm testing right now I got this:

julia> ϕ
RegisterDeformation.GridDeformation{Float64,1,Array{FixedSizeArrays.Vec{1,Float64},1},LinSpace{Float64}}(6-elementArray{FixedSizeArrays.Vec{1,Float64},1}:
 Vec(0.0)
 Vec(-1.0)
 Vec(0.0)
 Vec(-1.0)
 Vec(1.0)
 Vec(0.0)
,(linspace(1.0,6.0,6),))

julia> warped
6-element Array{Float64,1}:
   0.0     
 NaN       
   0.898912
   0.898912
   0.0     
   0.0     

Notice the NaN in the middle of the image. That's because the second pixel in the warped image didn't "receive" any intensity from the original image. Using a small lambda fixes that problem.

[Cody-G]

Oops accidentally closed the issue via tab+enter

Using a small lambda fixes that problem.

Ah interesting, that is good to know. One thing I haven't quite figured out is how to decide what is a small lambda for a given grid size. It seems to me that the two parameters are linked in a complex way. i.e. a lambda = 0.01 may allow little deformation on a coarse grid but significant deformation on a finer grid. Do you have any tips about how to scale lambda with the grid size? Maybe there is a normalization I could do so that I get an affine penalty on a per-cell basis rather than for an entire image? Could it be as simple as dividing lambda by prod(gridsize)? (assuming all cells have the same volume).