Re: Error in `optimize!` #50

mdhe1248 · 2017-04-26T17:00:35Z

Hi,
I wanted to update deformation. I ran the script below in flash.wustl.edu:/mnt/donghoon_030/20170412

using JLD
using Images, Unitful, FixedSizeArrays, JLD, AxisArrays
using BlockRegistration

fnmm = "/mnt/donghoon_030/20170412/exp2_20170412" #mm `base` name without "_i"`

#### Load data
Es, cs, Qs, knots, mmis = jldopen(string(fnmm, "_", 1, ".mm"), mmaparrays=true) do file
    read(file, "Es"), read(file, "cs"), read(file, "Qs"), read(file, "knots"), read(file, "mmis")
end

#### Calculate deformation (It may take several hours)
λ_old = 1e-6
λt_old = 1e-5
ap = AffinePenalty(knots, λ)
ϕs_old, penalty = fixed_λ(cs, Qs, knots, ap, λt, mmis; show_trace=true, iterations=5)

λ = 1e-7
λt = 1e-3
ϕs, penalty = fixed_λ(cs, Qs, knots, ap, λt, mmis; show_trace=true, iterations=5)

#### Optimization
ϕ_1, fval_1, fval0_1 = optimize!(ϕs, ϕs_old, ap, mmis);
ϕ_2, fval_2, fval0_2 = optimize!(ϕs, ϕs_old, ap, λt,  mmis);

Is this correct implementation? Should it update ϕs based on ϕs_old, ap, mmis?

I obtained this error:

julia> ϕ, fval, fval0 = optimize!(ϕ, ϕs_old, ap, mmis);
ERROR: MethodError: no method matching maxshift(::Float32)
Closest candidates are:
  maxshift(::CenterIndexedArrays.CenterIndexedArray{ND<:RegisterCore.NumDenom,N,A}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterCore.jl:226
  maxshift{T,N}(::CachedInterpolations.CachedInterpolation{T,N,M,O}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:1044
 in #_optimize!#5(::Array{Any,1}, ::Function, ::RegisterOptimize.DeformOpt{Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},RegisterPenalty.AffinePenalty{Float64,3},Array{Float32,8}}, ::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{Float32,8}, ::Float64, ::Int64) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:432
 in #optimize!#4(::Float64, ::Int64, ::Array{Any,1}, ::Function, ::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{Float32,8}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:426
 in optimize!(::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{Float32,8}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:425

The text was updated successfully, but these errors were encountered:

timholy · 2017-04-26T19:01:13Z

The only reason to update the deformation is if you warp the moving image and recalculate all of the mismatch-related parameters for the warped image (Es, cs, Qs, and mmis).

timholy · 2017-04-26T19:03:04Z

And if you're doing whole-experiment registration, that would basically require that you write out a disk file with the intermediate warped image.

mdhe1248 · 2017-04-26T19:25:50Z

Ah, I just wanted to briefly test it if optimize! works.
I will try from the beginning with small data set.
The procedure would be then (I only listed some key parameters),

Calculate mm using the original image: obtain mmis_old and other parameters.
Obtain ϕs_old using fixed_λ. And warp the original image.
Calculate mm using warped image: obtain mmis, and other parameters.
Get ap. Get ϕs using fixed_λ.
Run optimized_ϕs, fval, fval0 = optimize!(ϕs, ϕs_old, ap, mmis)
Warp the original image using an optimized_ϕs.

timholy · 2017-04-26T19:51:59Z

Yes, that's right. No guarantees that it will be any better than ϕs_old (one can always hope...)

timholy · 2017-04-26T19:56:53Z

I should also say that this isn't heavily tested so you're likely to get the same issue. Is mmis an array of CenterIndexedArrays (or something similar)? If it's an array of real values, check out the various versions of fixed_lambda (in RegisterOptimize) that convert from an Array{Tf} to the mmisc version.

mdhe1248 · 2017-04-26T20:14:47Z

That was my concern.

julia> typeof(mmis)  
Array{Float32,8}

I think this is it:
The 3rd and 4th lines from the bottom convert mmis -> mmisr -> mmisc.

function fixed_λ{Tf<:Number,T,N}(cs::Array{Tf}, Qs::Array{Tf}, knots::NTuple{N}, ap::AffinePenalty{T,N}, λt, mmis::Array{Tf}; kwargs...)
    csr = reinterpret(Vec{N,Tf}, cs, tail(size(cs)))
    Qsr = reinterpret(Mat{N,N,Tf}, Qs, tail(tail(size(Qs))))
    if length(mmis) > 10^7
        L = length(mmis)*sizeof(Tf)/1024^3
        @printf "The mismatch data are %0.2f GB in size.\n  During optimization, the initial function evaluations may be limited by I/O and\n  could be very slow. Later evaluations should be faster.\n" L
    end
    ND = NumDenom{Tf}
    mmisr = reinterpret(ND, mmis, tail(size(mmis)))
    mmisc = cachedinterpolators(mmisr, N, ntuple(d->(size(mmisr,d)+1)>>1, N))
    fixed_λ(csr, Qsr, knots, ap, λt, mmisc; kwargs...)
end

and maxshift function is being used here line 4.

function _optimize!(objective, ϕ, dp, mmis, tol, print_level; kwargs...)
    uvec = u_as_vec(ϕ)
    T = eltype(uvec)
    mxs = maxshift(first(mmis))

    solver = IpoptSolver(;hessian_approximation="limited-memory",
                         print_level=print_level,
                         sb="yes",
                         tol=tol, kwargs...)
    m = MathProgBase.NonlinearModel(solver)
    ub1 = T[mxs...] - T(RegisterFit.register_half)
    ub = repeat(ub1, outer=[div(length(uvec), length(ub1))])
    MathProgBase.loadproblem!(m, length(uvec), 0, -ub, ub, T[], T[], :Min, objective)
    MathProgBase.setwarmstart!(m, uvec)
    fval0 = MathProgBase.eval_f(objective, uvec)
    isfinite(fval0) || error("Initial value must be finite")
    MathProgBase.optimize!(m)

    stat = MathProgBase.status(m)
    stat == :Optimal || warn("Solution was not optimal")
    uopt = MathProgBase.getsolution(m)
    fval = MathProgBase.getobjval(m)
    _copy!(ϕ, uopt)
    ϕ, fval, fval0
end

mdhe1248 · 2017-04-26T20:17:23Z

mmisc seems to be passed into this function:

"""
`ϕs, penalty = fixed_λ(cs, Qs, knots, affinepenalty, λt, mmis)`
computes an optimal vector-of-deformations `ϕs` for an image sequence,
using an temporal penalty coefficient `λt`.
"""
function fixed_λ{T,N,_}(cs::AbstractArray{Vec{N,T}}, Qs::AbstractArray{Mat{N,N,T}}, knots::NTuple{N}, ap::AffinePenalty{_,N}, λt, mmis; mu_init=0.1, kwargs...)
    λtT = T(λt)
    apT = convert(AffinePenalty{T,N}, ap)
    maxshift = map(x->(x-1)>>1, size(first(mmis)))
    print("Calculating initial guess (this may take a while)...")
    u0, isconverged = initial_deformation(apT, λtT, cs, Qs)
    println("done")
    if !isconverged
        Base.warn_once("initial_deformation failed to converge with λ = ", ap.λ, ", λt = ", λt)
    end
    uclamp!(u0, maxshift)
    colons = ntuple(ColonFun, Val{N})
    ϕs = [GridDeformation(u0[colons..., i], knots) for i = 1:size(u0)[end]]
    local mismatch
    println("Starting optimization.") 
    optimize!(ϕs, identity, apT, λtT, mmis; kwargs...)
end

timholy · 2017-04-27T00:51:24Z

Yes, what I'm basically saying is that you might have to imitate the creation of mmisc and pass that. The issue is that saving it to disk destroys some of the "julia" structure and that's why, when you load it, it comes back as a plain array.

mdhe1248 · 2017-04-27T05:01:36Z

Okay,
I changed added one more optimize! function.

function optimize!{Tf<:Number, T, N}(ϕ, ϕ_old, dp::AffinePenalty{T,N}, mmis::Array{Tf})
    ND = NumDenom{Tf}
    mmisr = reinterpret(ND, mmis, tail(size(mmis)));
    mmisc = cachedinterpolators(mmisr, N, ntuple(d->(size(mmisr,d)+1)>>1, N))
    optimize!(ϕ, ϕ_old, dp::DeformationPenalty, mmisc) 
end

Would it be a suitable solution? hm. Is dp::DeformationPenalty the same as dp::AffinePenalty? I may delete ::DeformationPenalty.

After I run this modification, I got another error. I haven't looked at it carefully yet.

ERROR: MethodError: no method matching createProblem(::Int64, ::Array{Float32,1}, ::Array{Float32,1}, ::Int64, ::Array{Float32,1}, ::Array{Float32,1}, ::Int64, ::Int64, ::Ipopt.#eval_f_cb#4{RegisterOptimize.DeformOpt{Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},Arr
ay{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},RegisterPenalty.AffinePenalty{Float64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}}}, ::Ipopt.#eval_g_cb#6{RegisterOptimize.DeformOpt{Array{RegisterDeformation.GridDefor
mation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},RegisterPenalty.AffinePenalty{Float64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}}}, :
:Ipopt.#eval_grad_f_cb#5{RegisterOptimize.DeformOpt{Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},RegisterPenalty.AffinePenalty{Float64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}}}, ::Ipopt.#eval_jac_g_cb#7{RegisterOptimize.DeformOpt{Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Fl
oat32},3},LinSpace{Float64}},1},RegisterPenalty.AffinePenalty{Float64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}}}, ::Void)
Closest candidates are:
  createProblem(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Int64, ::Any, ::Any, ::Any, ::Any, ::Any) 
at /home/donghoon/.julia/v0.5/Ipopt/src/Ipopt.jl:171
  createProblem(::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Array{Float64,1}, ::Array{Float64,1}, ::Int64, ::Int64, ::Any, ::Any, ::Any, ::Any) at /hom
e/donghoon/.julia/v0.5/Ipopt/src/Ipopt.jl:171
 in loadproblem!(::Ipopt.IpoptMathProgModel, ::Int64, ::Int64, ::Array{Float32,1}, ::Array{Float32,1}, ::Array{Float32,1}, ::Array{Float32,1}, ::Symbol, ::RegisterOpti
mize.DeformOpt{Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},Array{RegisterDeformation.GridDeformat
ion{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},RegisterPenalty.AffinePenalty{Float64,3},Array{CachedInterpolations.CachedInterpolation{Reg
isterCore.NumDenom{Float32},3,7,(21,21,4)},4}}) at /home/donghoon/.julia/v0.5/Ipopt/src/IpoptSolverInterface.jl:92
 in #_optimize!#5(::Array{Any,1}, ::Function, ::RegisterOptimize.DeformOpt{Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},
LinSpace{Float64}},1},Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1},RegisterPenalty.AffinePenalty{F
loat64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}}, ::Array{RegisterDeformation.GridDeformation{Float32,3,Array
{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.Num
Denom{Float32},3,7,(21,21,4)},4}, ::Float64, ::Int64) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:441
 in #optimize!#4(::Float64, ::Int64, ::Array{Any,1}, ::Function, ::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace
{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Flo
at64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/Registe
rOptimize.jl:426
 in #__optimize!#20(::Float64, ::Int64, ::Array{Any,1}, ::Function, ::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSp
ace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{
Float64,3}, ::Array{Float32,8}) at ./REPL[91]:6
 in __optimize!(::Array{RegisterDeformation.GridDeformation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDe
formation{Float32,3,Array{FixedSizeArrays.Vec{3,Float32},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{Float32,8}) at ./REPL[91]:2

timholy · 2017-04-27T08:55:13Z

Would it be a suitable solution?

Yes!

Is dp::DeformationPenalty the same as dp::AffinePenalty? I may delete ::DeformationPenalty

You don't need the ::DeformationPenalty annotation on the last of your function (it doesn't do anything). DeformationPenalty is the abstract type:

julia> supertype(AffinePenalty)
RegisterPenalty.DeformationPenalty{T,N}

We used to use a different penalty but switched to the affine penalty. (See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3745275/ for a sense of just how many choices there are...)

After I run this modification, I got another error.

Looks like Ipopt (which is a C++ optimizer called by Julia) only supports Float64.

mdhe1248 · 2017-04-27T16:56:29Z

Converting Es, cs, Qs, mmis into Float64 resolved the error above.
for example, I ran the below:

cs = convert(Array{Float64, length(size(cs))}, cs)

Now I get another error:

julia> optimize!(ϕs, ϕs_old, ap, mmis)
ERROR: MethodError: no method matching zero(::Type{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}})
Closest candidates are:
  zero(::Type{Base.LibGit2.Oid}) at libgit2/oid.jl:88
  zero(::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuildItem}) at pkg/resolve/versionweight.jl:80
  zero(::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuild}) at pkg/resolve/versionweight.jl:120
  ...
 in penalty!(::Void, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.Gr
idDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedIn
terpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}, ::BitArray{4}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterPenalty.jl:80
 in penalty!(::Void, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.Gr
idDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedIn
terpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterPenalty.jl:79
 in #_optimize!#7(::Array{Any,1}, ::Function, ::RegisterOptimize.DeformOpt{Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},L
inSpace{Float64}},1},Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1},RegisterPenalty.AffinePenalty{Flo
at64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{Fi
xedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDeno
m{Float32},3,7,(21,21,4)},4}, ::Float64, ::Int64) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:451
 in #optimize!#4(::Float64, ::Int64, ::Array{Any,1}, ::Function, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{
Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float
64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOp
timize.jl:426
 in optimize!(::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDefor
mation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpola
tion{RegisterCore.NumDenom{Float32},3,7,(21,21,4)},4}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:425
 in optimize!(::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDefor
mation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{Float32,8}) at /home/donghoon/.julia
/v0.5/BlockRegistration/src/RegisterOptimize.jl:434

It seems zero(eltype(ϕs)) gives the error message.

timholy · 2017-04-27T17:38:43Z

Maybe replace that line with local val.

timholy · 2017-04-27T17:41:09Z

Oh, wait, you'll hit an error on the last line. Perhaps you need something like this:

numbertype{G<:GridDeformation}(ϕs::AbstractVector{G}) = numbertype(first(ϕs))
numbertype(ϕ) = eltype(ϕ)

and use T = numbertype(ϕ) instead of T = eltype(ϕ) on line 79.

mdhe1248 · 2017-04-27T19:22:52Z

I got another error in compose.
This seems like a type error like the one above.

julia> ϕ_1, fval_1, fval0_1 = optimize!(ϕs, ϕs_old, ap, mmis);
ERROR: MethodError: no method matching compose(::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1})
 in penalty!(::Void, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,7,(21,21,4)},4}, ::BitArray{4}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterPenalty.jl:86
 in penalty!(::Void, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,7,(21,21,4)},4}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterPenalty.jl:79
 in #_optimize!#7(::Array{Any,1}, ::Function, ::RegisterOptimize.DeformOpt{Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1},Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1},RegisterPenalty.AffinePenalty{Float64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,7,(21,21,4)},4}}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,7,(21,21,4)},4}, ::Float64, ::Int64) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:451
 in #optimize!#4(::Float64, ::Int64, ::Array{Any,1}, ::Function, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,7,(21,21,4)},4}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:426
 in optimize!(::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,7,(21,21,4)},4}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:425
 in optimize!(::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{Float64,8}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:434

timholy · 2017-04-27T19:47:04Z

It's basically the same thing: compose has been defined for a single deformation but not an array of deformations. You could add a compose method (to GridDeformations) like this:

compose{G1<:GridDeformation, G2<:GridDeformation}(ϕs_old::AbstractVector{G1}, ϕs_new::AbstractVector{G2}) = map(compose, ϕs_old, ϕs_new)

mdhe1248 · 2017-04-27T20:14:54Z

I added new method to GridDeformations. map(f, old, new) seems an elegant way.

I got another error.
This is also type error. I am slowly getting familiar with types. Still, I am not sure what to do here.

julia> compose(ϕs_old, ϕs)
ERROR: MethodError: no method matching compose(::RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}, ::RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}})
Closest candidates are:
  compose{T1,T2,N,A<:Interpolations.AbstractInterpolation{T,N,IT<:Union{Interpolations.InterpolationType,Tuple{Vararg{Interpolations.InterpolationType,N}}},GT<:Union{Interpolations.GridType,Interpolations.NoInterp,Tuple{Vararg{Union{Interpolations.GridType,Interpolations.NoInterp},N}}}}}(::RegisterDeformation.GridDeformation{T1,N,A<:Interpolations.AbstractInterpolation,L}, ::RegisterDeformation.GridDeformation{T2,N,A<:AbstractArray,L}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterDeformation.jl:295
  compose{T,N}(::Function, ::RegisterDeformation.GridDeformation{T,N,A<:AbstractArray,L}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterDeformation.jl:320
 in collect(::Base.Generator{Base.Zip2{Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1},Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}},Base.##3#4{RegisterDeformation.#compose}}) at ./array.jl:307
 in compose(::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}, ::Array{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},1}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterDeformation.jl:317

Also, if I look at compose, this looks for phi_old.u and phi_old.knots.
I am not sure my phi_old and phi_new has .u and .knots?

function compose{T1,T2,N,A<:AbstractInterpolation}(ϕ_old::GridDeformation{T1,N,A}, ϕ_new::GridDeformation{T2,N})
    u, knots = ϕ_old.u, ϕ_old.knots
    ϕ_new.knots == knots || error("Not yet implemented for incommensurate knots")
    unew = ϕ_new.u
    sz = map(length, knots)
    x = knot(knots, 1)
    out = _compose(u, unew, x, 1)
    ucomp = similar(u, typeof(out))
    TG = Mat{N,N,eltype(out)}
    g = Array{TG}(size(u))
    gtmp = Vector{typeof(out)}(N)
    eyeN = eye(TG)
    for I in CartesianRange(sz)
        x = knot(knots, I)
        dx = lookup(unew, x, I)
        y = x + dx
        ucomp[I] = dx + vecindex(u, y)
        vecgradient!(gtmp, u, y)
        g[I] = convert(TG, gtmp) + eyeN
    end
    GridDeformation(ucomp, knots), g
end

timholy · 2017-04-27T20:22:57Z

What is eltype(ϕs_old) and eltype(ϕs_new)?

One possibility is that you have to call interpolate!(ϕ_old) where ϕ_old is a single deformation. Beware that this overwrite the information stored in ϕ_old, so use copy if you don't want that to happen.

mdhe1248 · 2017-04-27T22:32:04Z

s_old and s_new have the same type.

julia> eltype(ϕs_old) 
RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}

Hm..
There is Interpolations.interpolate!. Is this the right function to use?
It does not change the type of ϕs[1]

julia> eltype(ϕs[1])
Float64
julia> eltype(Interpolations.interpolate!(ϕs[1]))
Float64

timholy · 2017-04-28T02:19:40Z

You have to have RegisterDeformation loaded (it's defined there, as an extension of the one in Interpolations).

mdhe1248 · 2017-04-28T03:37:48Z

Ah, this seems work.

julia> compose((Interpolations.interpolate(ϕs_old[1])), ϕs[1])

This is because of the way compose is defined:

function compose{T1,T2,N,A<:AbstractInterpolation}(ϕ_old::GridDeformation{T1,N,A}, ϕ_new::GridDeformation{T2,N})
...

I needed to interpolate ϕ_old but not ϕ_new.

I guess I also need to change

compose{G1<:GridDeformation, G2<:GridDeformation}(ϕs_old::AbstractVector{G1}, ϕs_new::AbstractVector{G2}) = map(compose, ϕs_old, ϕs_new)

to

compose{G1<:GridDeformation, G2<:GridDeformation}(ϕs_old::AbstractVector{G1}, ϕs_new::AbstractVector{G2}) = map(compose, map(Interpolations.interpolate!, ϕs_old), ϕs_new)

Do you think this change is reasonable? I will try it and see if I get error anyway.

mdhe1248 · 2017-04-28T05:01:43Z

After the modification above, I got another error.
I kept having errors after some revisions. So I tried this (inspired by interpolate!).

julia> optimize!(ϕs[1], Interpolations.interpolate!(ϕs_old[1]), ap, mmis) 

ERROR: u should have length 63480, but length(u) = 6348
 in penalty!(::Void, ::Array{FixedSizeArrays.Vec{3,Float64},3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}, ::Bi
tArray{3}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterPenalty.jl:192
 in penalty!(::Void, ::RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}, ::RegisterDeformation.GridDeformation{
Float64,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{FixedSizeArrays.Vec{3,Float64},3,Array{FixedSizeArray
s.Vec{3,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},Interpolations.BSpline{Interpolations.Quadratic{I
nterpolations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{Float64}}, ::RegisterPenalty.AffinePenalty{Float64
,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}, ::BitArray{3}) at /home/donghoon/.julia/v0.5/BlockRegistration/
src/RegisterPenalty.jl:93
 in eval_f(::RegisterOptimize.DeformOpt{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}},RegisterDeformation.Gr
idDeformation{Float64,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{FixedSizeArrays.Vec{3,Float64},3,Array{
FixedSizeArrays.Vec{3,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},Interpolations.BSpline{Interpolatio
ns.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{Float64}},RegisterPenalty.AffinePena
lty{Float64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}}, ::Array{Float64,1}) at /home/donghoon/.julia/v0.5/Bloc
kRegistration/src/RegisterOptimize.jl:501
 in #_optimize!#7(::Array{Any,1}, ::Function, ::RegisterOptimize.DeformOpt{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpa
ce{Float64}},RegisterDeformation.GridDeformation{Float64,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{Fixe
dSizeArrays.Vec{3,Float64},3,Array{FixedSizeArrays.Vec{3,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},
Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{F
loat64}},RegisterPenalty.AffinePenalty{Float64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}}, ::RegisterDeformati
on.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.Cache
dInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}, ::Float64, ::Int64) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:451
 in #optimize!#4(::Float64, ::Int64, ::Array{Any,1}, ::Function, ::RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float
64}}, ::RegisterDeformation.GridDeformation{Float64,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{FixedSize
Arrays.Vec{3,Float64},3,Array{FixedSizeArrays.Vec{3,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},Inter
polations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{Float6
4}}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}) at /home/donghoon/.
julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:426
 in optimize!(::RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}, ::RegisterDeformation.GridDeformation{Float64
,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{FixedSizeArrays.Vec{3,Float64},3,Array{FixedSizeArrays.Vec{3
,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},Interpolations.BSpline{Interpolations.Quadratic{Interpol
ations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{Float64}}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::
Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.j
l:425
 in optimize!(::RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}, ::RegisterDeformation.GridDeformation{Float64
,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{FixedSizeArrays.Vec{3,Float64},3,Array{FixedSizeArrays.Vec{3
,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},Interpolations.BSpline{Interpolations.Quadratic{Interpol
ations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{Float64}}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::
Array{Float64,7}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:434

I got another error. I think but it is more promising. I think the problem comes from mismatch between parameter size; while I am using ϕs[1] and ϕs_old[1], I use full ap and mmis. I am not sure though. mmis[:,:,:,:,:,:,:,1] also failed:

julia> optimize!(ϕs[1], Interpolations.interpolate!(ϕs_old[1]), ap, mmis[:,:,:,:,:,:,:,1])
ERROR: Initial value must be finite
 in #_optimize!#7(::Array{Any,1}, ::Function, ::RegisterOptimize.DeformOpt{RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpa
ce{Float64}},RegisterDeformation.GridDeformation{Float64,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{Fixe
dSizeArrays.Vec{3,Float64},3,Array{FixedSizeArrays.Vec{3,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},
Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{F
loat64}},RegisterPenalty.AffinePenalty{Float64,3},Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}}, ::RegisterDeformati
on.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.Cache
dInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}, ::Float64, ::Int64) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:452
 in #optimize!#4(::Float64, ::Int64, ::Array{Any,1}, ::Function, ::RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float
64}}, ::RegisterDeformation.GridDeformation{Float64,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{FixedSize
Arrays.Vec{3,Float64},3,Array{FixedSizeArrays.Vec{3,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},Inter
polations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{Float6
4}}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}) at /home/donghoon/.
julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:426
 in optimize!(::RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}, ::RegisterDeformation.GridDeformation{Float64
,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{FixedSizeArrays.Vec{3,Float64},3,Array{FixedSizeArrays.Vec{3
,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},Interpolations.BSpline{Interpolations.Quadratic{Interpol
ations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{Float64}}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::
Array{CachedInterpolations.CachedInterpolation{RegisterCore.NumDenom{Float64},3,6,(21,21,4)},3}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.j
l:425
 in optimize!(::RegisterDeformation.GridDeformation{Float64,3,Array{FixedSizeArrays.Vec{3,Float64},3},LinSpace{Float64}}, ::RegisterDeformation.GridDeformation{Float64
,3,Interpolations.ScaledInterpolation{FixedSizeArrays.Vec{3,Float64},3,Interpolations.BSplineInterpolation{FixedSizeArrays.Vec{3,Float64},3,Array{FixedSizeArrays.Vec{3
,Float64},3},Interpolations.BSpline{Interpolations.Quadratic{Interpolations.InPlace}},Interpolations.OnCell,0},Interpolations.BSpline{Interpolations.Quadratic{Interpol
ations.InPlace}},Interpolations.OnCell,Tuple{LinSpace{Float64},LinSpace{Float64},LinSpace{Float64}}},LinSpace{Float64}}, ::RegisterPenalty.AffinePenalty{Float64,3}, ::
Array{Float64,7}) at /home/donghoon/.julia/v0.5/BlockRegistration/src/RegisterOptimize.jl:434

timholy · 2017-04-28T07:55:11Z

Hard to tell without a little more information. Let's check together? (Or post your complete script/data file?)

mdhe1248 · 2017-04-28T12:52:32Z

The error comes from _optimize! in RegisterOptimize.jl

    fval0 = MathProgBase.eval_f(objective, uvec)
    isfinite(fval0) || error("Initial value must be finite")

I am uploading my script. I will talk to you later today.
My data (mm file) is in flash.wustl.edu:/mnt/donghoon_030/20170412. It will take ~5 minutes to finish the below.

using JLD
using Images, Unitful, FixedSizeArrays, JLD, AxisArrays
using BlockRegistration

fnmm = "/mnt/donghoon_030/20170412/exp2_20170412" #mm `base` name without "_i"`

#### Load data
Es, cs, Qs, knots, mmis = jldopen(string(fnmm, "_", 1, ".mm"), mmaparrays=true) do file
    read(file, "Es", "cs", "Qs", "knots", "mmis")
end;

#knots = map(x -> convert(LinSpace{Float32}, x), knots)
Es = convert(Array{Float64, length(size(Es))}, Es);
cs = convert(Array{Float64, length(size(cs))}, cs);
Qs = convert(Array{Float64, length(size(Qs))}, Qs);
mmis = convert(Array{Float64, length(size(mmis))}, mmis);

### Calculate deformation: Use same data for test.
λ_old = 1e-6
λt_old = 1e-5
ap = AffinePenalty(knots, λ_old)
ϕs_old, penalty_old = fixed_λ(cs, Qs, knots, ap, λt_old, mmis; show_trace=true, iterations=3)

λ = 1e-7
λt = 1e-3
ϕs, penalty = fixed_λ(cs, Qs, knots, ap, λt, mmis; show_trace=true, iterations=3)

#### Optimization
ϕ_1, fval_1, fval0_1 = optimize!(ϕs[1], Interpolations.interpolate!(ϕs_old), ap, mmis);
#ϕ_2, fval_2, fval0_2 = optimize!(ϕs, ϕs_old, ap, λt,  mmis);

mdhe1248 · 2017-04-28T23:16:38Z

Hi Tim,

Could you teach me how to interpolate deformation?
For example, I have 3X3X2 deformation. I want to interpolate and have 5X6X3 deformation.

I thought I could do something like this:

knots = (linspace(1.0,300.0,3),linspace(1.0,300.0,3),linspace(1.0,20.0,2));
itp = interpolate(knots, ϕs[1].u)

#new_u = [itp[1,1,1] itp[1.5, 1,1] itp[2, 1, 1] ...

timholy · 2017-04-29T09:49:41Z

You've got the right idea, but it's even easier. Something like this:

using Interpolations
ϕi = interpolate(ϕ)   # yes, this does work!
# Now, all we have to do is evaluate ϕi at the new knot locations:
knots_new = (linspace(1,300,5), linspace(1,300,6), linspace(1,20,3))
sz = map(length, knots_new)
unew = similar(ϕ.u, sz)
for I in CartesianRange(sz)
    x = map((k,i)->k[i], knots_new, I.I)
    unew[I] = ϕi.u[x...]
end
ϕnew = GridDeformation(unew, knots_new)

timholy · 2017-04-29T10:46:08Z

See #52.

mdhe1248 · 2017-04-30T16:38:56Z

Thank you Tim, the script above works well. map and CartesianRange for looping seems clever to me.
I warped an image using interpolated ϕs. The result was not as quite good as the warped image from non-interpolated ϕs. Is it expected? Otherwise, I may have made some mistake.

timholy · 2017-04-30T19:30:06Z

I didn't test the script above, but #52 is well-tested so I'd recommend using it. There can be very subtle differences, e.g., using ϕi and ϕnew from #52's new tests (check the code), we have this:

julia> ϕni = interpolate(ϕnew);

julia> ϕi[3.2,4.77]
Vec(7.177265532251683,8.416)

julia> ϕni[3.2,4.77]
Vec(7.188572268098131,8.395451023838106)

So there's not identical, but they are extremely close, so close I doubt you could see it. Maybe only around the edges? To say anything more, I'd need you to be more specific about "not quite as good."

mdhe1248 · 2017-05-01T00:01:50Z

Ah, sorry Tim, I think I made mistake. I re-ran my test script and found little difference.
First, I looked at the warped images by eye. Later, I calculated errors (sum of squares).

I will try optimization and see if I don't get any errors.

mdhe1248 · 2017-05-01T13:01:23Z

Everything worked without error. So I can close this report.
Here is what I did.

calculated initial mismatch.
got initial_ϕs (later, I call it ϕs_old) it and an intermediate warped image.
calculated mismatch from the intermediate image.
obtained new ϕs

ϕs, penalty = fixed_λ(cs_new, Qs_new, knots_new, ap, λt, mmis_new; ϕs_old = ϕs_old, show_trace=true, iterations=10000);

ran ϕs_c = compose(ϕs, new_ϕs)
warped the original image using ϕs_c.
compare the result with the image warped twice.

The error seems lower when image was warped twice. Does it mean it is better to warp multiple times rather than warp one time with optimized ϕ? Is it possible to have a bug?
I will briefly test several more things and close this issue.

mdhe1248 · 2017-05-01T13:27:36Z

I wanted to make sure compose worked well.
I compared warped images:
Warping was done with ϕs_c or just ϕs without compose. For each warped image, I calculated error between a moving image and a fixed image.
Yes, compose seems working. ϕs_c reduced errors overall. I don't understand why this is not better than "multiple warping".
I will close this thread.

timholy · 2017-05-01T14:27:46Z

I'd like to understand better the exact error computation you're showing here. Is the error the penalty output of fixed_λ, or is it mismatch0? If the latter, are you using the same normalization (:intensity or :pixels) both for optimizing and for measuring? The penalties could be different, and optimizing for one might not optimize for the other.

Finally, what are your grid size and lambdas? One concern might be that composition is just not sufficiently accurate if there are substantive discontinuities in the deformation. We might have to go to very large grid sizes, unfortunately, to make this work properly.

timholy · 2017-05-01T14:30:23Z

Also, in that second plot, in "only new ϕ" what are you comparing against? Here's my guess:

combined ϕ is for comparing moving against fixed
only new ϕ is for comparing moving_warped_once against fixed

If that's not right, please let me know.

mdhe1248 · 2017-05-01T15:21:37Z

Here, I used calculated error using "mean of squares". This is because the number of NaN is different in different images.
Original -> non- warped image
warp1 -> warped image using coarse grid (3,3,2)
combined ϕ -> warped using composed ϕ. (grid size was 3,3,2 first, 6,6,4 second.)
multi-warp -> warp twice.

multi warp still seems better.

timholy mentioned this issue Apr 29, 2017

Implement regridding for GridDeformations #52

Merged

mdhe1248 closed this as completed May 1, 2017

timholy reopened this May 1, 2017

timholy mentioned this issue May 1, 2017

Split original image creation in warpgrid out into knotgrid #53

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Re: Error in `optimize!` #50

Re: Error in `optimize!` #50

mdhe1248 commented Apr 26, 2017

timholy commented Apr 26, 2017

timholy commented Apr 26, 2017

mdhe1248 commented Apr 26, 2017

timholy commented Apr 26, 2017

timholy commented Apr 26, 2017

mdhe1248 commented Apr 26, 2017

mdhe1248 commented Apr 26, 2017

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017 •

edited

Loading

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017

timholy commented Apr 27, 2017

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017

timholy commented Apr 28, 2017

mdhe1248 commented Apr 28, 2017 •

edited

Loading

mdhe1248 commented Apr 28, 2017 •

edited

Loading

timholy commented Apr 28, 2017

mdhe1248 commented Apr 28, 2017 •

edited

Loading

mdhe1248 commented Apr 28, 2017

timholy commented Apr 29, 2017 •

edited

Loading

timholy commented Apr 29, 2017

mdhe1248 commented Apr 30, 2017

timholy commented Apr 30, 2017

mdhe1248 commented May 1, 2017

mdhe1248 commented May 1, 2017 •

edited

Loading

mdhe1248 commented May 1, 2017

timholy commented May 1, 2017

timholy commented May 1, 2017

mdhe1248 commented May 1, 2017

Re: Error in optimize! #50

Re: Error in optimize! #50

Comments

mdhe1248 commented Apr 26, 2017

timholy commented Apr 26, 2017

timholy commented Apr 26, 2017

mdhe1248 commented Apr 26, 2017

timholy commented Apr 26, 2017

timholy commented Apr 26, 2017

mdhe1248 commented Apr 26, 2017

mdhe1248 commented Apr 26, 2017

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017 • edited Loading

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017

timholy commented Apr 27, 2017

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017

timholy commented Apr 27, 2017

mdhe1248 commented Apr 27, 2017

timholy commented Apr 28, 2017

mdhe1248 commented Apr 28, 2017 • edited Loading

mdhe1248 commented Apr 28, 2017 • edited Loading

timholy commented Apr 28, 2017

mdhe1248 commented Apr 28, 2017 • edited Loading

mdhe1248 commented Apr 28, 2017

timholy commented Apr 29, 2017 • edited Loading

timholy commented Apr 29, 2017

mdhe1248 commented Apr 30, 2017

timholy commented Apr 30, 2017

mdhe1248 commented May 1, 2017

mdhe1248 commented May 1, 2017 • edited Loading

mdhe1248 commented May 1, 2017

timholy commented May 1, 2017

timholy commented May 1, 2017

mdhe1248 commented May 1, 2017

Re: Error in `optimize!` #50

Re: Error in `optimize!` #50

mdhe1248 commented Apr 27, 2017 •

edited

Loading

mdhe1248 commented Apr 28, 2017 •

edited

Loading

mdhe1248 commented Apr 28, 2017 •

edited

Loading

mdhe1248 commented Apr 28, 2017 •

edited

Loading

timholy commented Apr 29, 2017 •

edited

Loading

mdhe1248 commented May 1, 2017 •

edited

Loading