diff --git a/src/RegisterWorkerRigid.jl b/src/RegisterWorkerRigid.jl index 405acd4..01f2627 100644 --- a/src/RegisterWorkerRigid.jl +++ b/src/RegisterWorkerRigid.jl @@ -8,7 +8,7 @@ using RegisterWorkerShell import RegisterWorkerShell: worker -export Rigid +export Rigid, RigidGridStart type Rigid{A<:AbstractArray,T} <: AbstractWorker fixed::A @@ -65,4 +65,55 @@ function worker(algorithm::Rigid, img, tindex, mon) mon end + +type RigidGridStart{A<:AbstractArray,T} <: AbstractWorker + fixed::A + thresh::T + maxradians::Vector{Float64} + rgridsz::Vector{Int} + mxshift::Vector{Int} + SD::Matrix{T} + params + workerpid::Int +end + +""" +`RigidGridStart(fixed, maxradians, rgridsz, mxshift; thresh_fac=(0.5)^ndims(fixed), thresh=nothing, SD = eye(ndims(fixed)), pid=1, kwargs...)` +initializes the `RigidGridStart` registration algorithm. This algorithm works much like `Rigid` except that instead of initializing with a +principle axis transformation it intializes with a grid search of rotations+shifts. Possible shifts are specified by `mxshift` while rotations +are specified by `maxradians` and `rgridsz`. See docs for `RegisterOptimize.rotation_gridsearch` for more details. +""" +function RigidGridStart(fixed, maxradians, rgridsz, mxshift; thresh_fac=(0.5)^ndims(fixed), thresh=0, SD = eye(ndims(fixed)), pid=1, kwargs...) + nimages(fixed) == 1 || error("Register to a single image") + SDm = convert(Matrix, SD) + T = eltype(SDm) + T = T <: AbstractFloat ? T : Float64 + cthresh = thresh == nothing ? thresh_fac*sum(abs2, fixed) : thresh + params = Dict{Symbol,Any}(kwargs) + RigidGridStart{typeof(fixed),T}(fixed, convert(T, cthresh), [maxradians...], [rgridsz...], [mxshift...], SDm, params, pid) +end + +function init!(algorithm::RigidGridStart) + eval(:(using RegisterMismatch)) +end + +function worker(algorithm::RigidGridStart, img, tindex, mon) + moving = getindex_t(img, tindex) + if any(x->x>1, algorithm.rgridsz) + tfm, best_mm = rotation_gridsearch(algorithm.fixed, moving, algorithm.mxshift, algorithm.maxradians, algorithm.rgridsz, algorithm.SD) + else + tfm = tformeye(ndims(moving)) + end + print("Beginning iterative optimization with this transform, returned by grid search:\n $tfm") + tfm, mismatch = optimize_rigid(algorithm.fixed, moving, tfm, [size(algorithm.fixed)...]/2, algorithm.SD, thresh=algorithm.thresh; print_level=get(algorithm.params, :print_level, 0)) + # There are no Rigid parameters that are expected as outputs, + # so no need to call monitor!(mon, algorithm) + monitor!(mon, :tform, tfm) + monitor!(mon, :mismatch, mismatch) + if haskey(mon, :warped) + monitor!(mon, :warped, transform(moving, tfm)) + end + mon +end + end # module diff --git a/test/rigid.jl b/test/rigid.jl index 30d1ebb..0514f6d 100644 --- a/test/rigid.jl +++ b/test/rigid.jl @@ -3,6 +3,7 @@ using BlockRegistration using BlockRegistrationScheduler, RegisterDriver, RegisterWorkerShell, RegisterWorkerRigid using Base.Test +#Rigid fixed = testimage("cameraman") tfm = tformrotate(pi/12) moving = transform(fixed, tfm) @@ -16,3 +17,23 @@ ptfm = mon[:tform]*tfm @test ≈(ptfm.scalefwd, eye(2), atol=0.01) @test all(abs.(ptfm.offset) .< 0.2) @test mon[:mismatch] < 0.001 + +#RigidGridStart +fixed = testimage("cameraman") +tfm = tformrotate(pi/10) +moving = transform(fixed, tfm) +moving[find(x->isnan(x), moving)] = zero(eltype(moving)) + +maxradians = pi/10 +rgridsz = 3 +mxshift = (3,3) +alg = RigidGridStart(fixed, maxradians, rgridsz, mxshift; print_level=5) +#alg = RigidGridStart(fixed, maxradians, 1, mxshift; print_level=5) +mon = monitor(alg, ()) +mon[:tform] = nothing +mon[:mismatch] = 0.0 +mon = driver(alg, moving, mon) +ptfm = mon[:tform]*tfm +@test ≈(ptfm.scalefwd, eye(2), atol=0.01) +@test all(abs.(ptfm.offset) .< 0.2) +@test mon[:mismatch] < 0.001