HolyLab / BlockRegistrationScheduler.jl

Multi-core image registration scheduler
1 stars 0 forks source link

adapting registration for 2D time series #38

Closed kellyhill90 closed 7 years ago

kellyhill90 commented 7 years ago

I am trying to adapt the registration for a 2d time series. I started from the BlockRegistrationScheduler stack-by-stack example and tried to identify all the places that 3d data needs to be converted to 2d. However, I still end up with an error that seems to be related to size mismatch. Any suggestions? Error and code both included below. (Probably not relevant, but note that I skipped steps related to preprocessing, bias adjustment, and roi selection.)

ERROR: ssize and gridsize must have the same length
aperture_grid(::Tuple{Int64,Int64,Int64}, ::Tuple{Int64,Int64}) at /home/kelly/.julia/v0.6/BlockRegistration/src/RegisterMismatchCommon.jl:153
worker(::RegisterWorkerApertures.Apertures{Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},2},Float64,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},2}, ::Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3}, ::Int64, ::Dict{Symbol,Any}) at /home/kelly/.julia/v0.6/BlockRegistrationScheduler/src/RegisterWorkerApertures.jl:163
macro expansion at /home/kelly/.julia/v0.6/BlockRegistrationScheduler/src/RegisterDriver.jl:90 [inlined]
(::RegisterDriver.##3#9{JLD.JldFile,Dict{Symbol,Any},SharedArray{Bool,1},SharedArray{Bool,1},RegisterDriver.#getnextidx#8,RemoteChannel{Channel{Any}},Array{RegisterWorkerApertures.Apertures,1},Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3},Array{Dict{Symbol,Any},1}})() at ./task.jl:335
Stacktrace:
 [1] sync_end() at ./task.jl:287
 [2] macro expansion at ./task.jl:303 [inlined]
 [3] (::RegisterDriver.##2#7{Array{RegisterWorkerApertures.Apertures,1},Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3},Array{Dict{Symbol,Any},1},Int64})(::JLD.JldFile) at /home/kelly/.julia/v0.6/BlockRegistrationScheduler/src/RegisterDriver.jl:78
 [4] #jldopen#11(::Array{Any,1}, ::Function, ::RegisterDriver.##2#7{Array{RegisterWorkerApertures.Apertures,1},Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3},Array{Dict{Symbol,Any},1},Int64}, ::String, ::Vararg{String,N} where N) at /home/kelly/.julia/v0.6/JLD/src/JLD.jl:256
 [5] driver(::String, ::Array{RegisterWorkerApertures.Apertures,1}, ::Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3}, ::Array{Dict{Symbol,Any},1}) at /home/kelly/.julia/v0.6/BlockRegistrationScheduler/src/RegisterDriver.jl:68

For reference, this is how I adapted the code. The ssize/gridsize error comes at the last line (driver).

using BlockRegistration, BlockRegistrationScheduler
using RegisterWorkerApertures

const μm = u"μm"
const s=u"s"

img0=load("testM.tif");
fixedidx = (size(img0,3)+1) >> 1
fixed0=img0[:,:,fixedidx];
ps=[1.15μm,1.15μm]
σ = 25μm
sigmahp = Float64[σ/x for x in ps]
sigmalp = [0,0]  # lowpass filtering is not currently recommended
fixed=fixed0;

gridsize=(round(Int,Float64[10μm/x for x in ps])[1],round(Int,Float64[10μm/x for x in ps])[2])
knots = map(d->linspace(1,size(fixed,d),gridsize[d]), (1:ndims(fixed)...))

mxshift=(10,20);
λ=0.01

algorithm = Apertures[Apertures(fixed, knots, mxshift, λ, correctbias=false)]

mon = monitor(algorithm, (), Dict{Symbol,Any}(:u=>ArrayDecl(Array{Vec{2,Float64},2}, gridsize)))

basename = splitext(splitdir(fn)[2])[1]
@show basename
fileout = string(basename, ".register")

@time driver(fileout, algorithm, img0, mon)
mdhe1248 commented 7 years ago

Kelly, Could you add Axis on your image? Something like this:

fn = "testM.tif"
img0=load(fn);
img = AxisArray(img0, (:x, :y, :time), (1.15μm, 1.15μm, 2s))

I think may help solving the error above. However, you may get another error.

mdhe1248 commented 7 years ago

This is the error I get:

julia> @time driver(fileout, algorithm, img, mon)
Working on 1
Planning FFTs (maximum Inf seconds)...done (0.99 seconds)
ERROR: MethodError: Cannot `convert` an object of type SVector{2,Float64} to an object of type Tuple{Float64,Float64}
This may have arisen from a call to the constructor Tuple{Float64,Float64}(...),
since type constructors fall back to convert methods.
copy!(::IndexLinear, ::Array{FixedSizeArrays.Vec{2,Float64},2}, ::IndexLinear, ::Array{SVector{2,Float64},2}) at ./abstractarray.jl:655
copy!(::SharedArray{FixedSizeArrays.Vec{2,Float64},2}, ::Array{SVector{2,Float64},2}) at ./sharedarray.jl:558
monitor!(::Dict{Symbol,Any}, ::Symbol, ::Array{SVector{2,Float64},2}) at /home/donghoon/.julia/v0.6/BlockRegistrationScheduler/src/RegisterWo
rkerShell.jl:105
worker(::RegisterWorkerApertures.Apertures{AxisArrays.AxisArray{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},2,SubArray{ColorTypes.Gr
ay{FixedPointNumbers.Normed{UInt16,16}},2,Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3},Tuple{Base.Slice{Base.OneTo{Int64}},B
ase.Slice{Base.OneTo{Int64}},Int64},true},Tuple{AxisArrays.Axis{:x,StepRangeLen{Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float6
4, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}},AxisArrays.Axis{:y,StepRangeLen{Quantity{Float64, Dimensions:{
𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}}}},Float64,StepRangeLen{Float64,
Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},2}, ::AxisArrays.AxisArray{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3,
Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3},Tuple{AxisArrays.Axis{:x,StepRangeLen{Quantity{Float64, Dimensions:{𝐋}, Units:{
μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}},AxisArrays.Axis{:y,StepRangeLen{Quantity{F
loat64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}},AxisArrays.
Axis{:time,StepRange{Quantity{Int64, Dimensions:{𝐓}, Units:{s}},Quantity{Int64, Dimensions:{𝐓},Units:{s}}}}}}, ::Int64, ::Dict{Symbol,Any}) a
t /home/donghoon/.julia/v0.6/BlockRegistrationScheduler/src/RegisterWorkerApertures.jl:190
macro expansion at /home/donghoon/.julia/v0.6/BlockRegistrationScheduler/src/RegisterDriver.jl:90 [inlined]
(::RegisterDriver.##3#9{JLD.JldFile,Dict{Symbol,Any},SharedArray{Bool,1},SharedArray{Bool,1},RegisterDriver.#getnextidx#8,RemoteChannel{Chann
el{Any}},Array{RegisterWorkerApertures.Apertures,1},AxisArrays.AxisArray{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3,Array{ColorTy
pes.Gray{FixedPointNumbers.Normed{UInt16,16}},3},Tuple{AxisArrays.Axis{:x,StepRangeLen{Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity
{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}},AxisArrays.Axis{:y,StepRangeLen{Quantity{Float64, Dimen
sions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}},AxisArrays.Axis{:time,St
epRange{Quantity{Int64, Dimensions:{𝐓}, Units:{s}},Quantity{Int64, Dimensions:{𝐓}, Units:{s}}}}}},Array{Dict{Symbol,Any},1}})() at ./task.jl:
335
Stacktrace:
 [1] sync_end() at ./task.jl:287
 [2] macro expansion at ./task.jl:303 [inlined]
 [3] (::RegisterDriver.##2#7{Array{RegisterWorkerApertures.Apertures,1},AxisArrays.AxisArray{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,
16}},3,Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3},Tuple{AxisArrays.Axis{:x,StepRangeLen{Quantity{Float64, Dimensions:{𝐋}, 
Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}},AxisArrays.Axis{:y,StepRangeLen{Qua
ntity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}},Axis
Arrays.Axis{:time,StepRange{Quantity{Int64, Dimensions:{𝐓}, Units:{s}},Quantity{Int64, Dimensions:{𝐓}, Units:{s}}}}}},Array{Dict{Symbol,Any},
1},Int64})(::JLD.JldFile) at /home/donghoon/.julia/v0.6/BlockRegistrationScheduler/src/RegisterDriver.jl:78
 [4] #jldopen#11(::Array{Any,1}, ::Function, ::RegisterDriver.##2#7{Array{RegisterWorkerApertures.Apertures,1},AxisArrays.AxisArray{ColorType
s.Gray{FixedPointNumbers.Normed{UInt16,16}},3,Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3},Tuple{AxisArrays.Axis{:x,StepRang
eLen{Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}
}}},AxisArrays.Axis{:y,StepRangeLen{Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Floa
t64, Dimensions:{𝐋}, Units:{μm}}}},AxisArrays.Axis{:time,StepRange{Quantity{Int64, Dimensions:{𝐓}, Units:{s}},Quantity{Int64, Dimensions:{𝐓},
 Units:{s}}}}}},Array{Dict{Symbol,Any},1},Int64}, ::String, ::Vararg{String,N} where N) at /home/donghoon/.julia/v0.6/JLD/src/JLD.jl:256
 [5] driver(::String, ::Array{RegisterWorkerApertures.Apertures,1}, ::AxisArrays.AxisArray{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16
}},3,Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt16,16}},3},Tuple{AxisArrays.Axis{:x,StepRangeLen{Quantity{Float64, Dimensions:{𝐋}, Un
its:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}},AxisArrays.Axis{:y,StepRangeLen{Quant
ity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}},Quantity{Float64, Dimensions:{𝐋}, Units:{μm}}}},AxisAr
rays.Axis{:time,StepRange{Quantity{Int64, Dimensions:{𝐓}, Units:{s}},Quantity{Int64, Dimensions:{𝐓}, Units:{s}}}}}}, ::Array{Dict{Symbol,Any}
,1}) at /home/donghoon/.julia/v0.6/BlockRegistrationScheduler/src/RegisterDriver.jl:68

Does anyone have an idea?

kellyhill90 commented 7 years ago

Jerry thanks for the idea. I get the same new error that you posted and haven't been able to work through a solution to it either. We'll see if anyone has another thought.

timholy commented 7 years ago

Fixes coming!

kellyhill90 commented 7 years ago

I'm trying to write a header file to open the warped output from a 2D registration. My original data was a tif file, not Imagine format. My best guess with some help from Jerry and Cody is:

props=imgoriginal.properties
props =Dict("headeronly" => true, "datafile" => imgwarped)
save("test.nhdr", imgoriginal; props = props)  #imgoriginal is just used a placeholder here; imgwarped will be loaded based on the headerfile

However the first line is unsuccessful; either ERROR: "type Array has no field properties" or ERROR: "type AxisArray has no field properties" (depending whether I loaded the tif file directly or converted it to AxisArray). If I skip the first line completely and just set the two properties (headeronly, datafile) explicitly, then the resulting header file loads as a time series of random noise.

Is there another solution to reach the end goal of opening the warped image? How to load properties of an AxisArray if not by img.properties?

mdhe1248 commented 7 years ago

Okay. Tim taught me how to rewrite a header before. I found how I did it. I did not fully test with Julia06 yet. I am working on it. It will be quick.

datafn = "warp_pre.cam"  ### Warped image file name.
save_headerfn = "warp_pre.nhdr"    #### new header name.

####if you have img with AxisArray, `ax = img.axes` may work instead of the
steps below.
psx = (1.2:1.2:512*1.2)*μm
psy = (1.2:1.2:512*1.2)*μm
ti = (2:2:16*2)*s
ax = (Axis{:x}(psx), Axis{:y}(psy), Axis{:time}(ti))

header = NRRD.headerinfo(Float16, ax)  ## Maybe in your case,` Float32`
header["datafile"] = datafn

#### Save header
open(save_headerfn, "w") do file
    write(file, magic(format"NRRD")) # Because header does not know NRRD
format.
    NRRD.write_header(file, "0004", header)  #"0004" is newer version.
"0002" is older version.
end

UPDATE:: It works on my Cannon environment. (with subtle modifications e.g. filenames, using modules)

2017-07-06 10:07 GMT-05:00 kellyhill90 notifications@github.com:

I'm trying to write a header file to open the warped output from a 2D registration. My original data was a tif file, not Imagine format. My best guess with some help from Jerry and Cody is:

props=imgoriginal.properties props =Dict("headeronly" => true, "datafile" => imgwarped) save("test.nhdr", imgoriginal; props = props) #imgoriginal is just used a placeholder here; imgwarped will be loaded based on the headerfile

However the first line is unsuccessful; either ERROR: "type Array has no field properties" or ERROR: "type AxisArray has no field properties" (depending whether I loaded the tif file directly or converted it to AxisArray). If I skip the first line completely and just set the two properties (headeronly, datafile) explicitly, then the resulting header file loads as a time series of random noise.

Is there another solution to reach the end goal of opening the warped image? How to load properties of an AxisArray if not by img.properties?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/HolyLab/BlockRegistrationScheduler/issues/38#issuecomment-313424410, or mute the thread https://github.com/notifications/unsubscribe-auth/AMvfFWtJEHINA_gUBENKFzf3xDj82W_gks5sLPgPgaJpZM4OJ9dw .

-- Donghoon (Jerry) Lee Neuroscience PhD Candidate | Lab of Timothy E. Holy Department of Anatomy and Neurobiology Graduate Program in Neuroscience Washington University in St. Louis

kellyhill90 commented 7 years ago

What packages do you use for this to work? Do I need to clone this NRRD package directly? https://github.com/JuliaIO/NRRD.jl

using Images, FileIO --> I get UndefVarError: NRRD not defined

Maybe there is another package that I need to add, clone, or update before using NRRD?

mdhe1248 commented 7 years ago

Pkg.add("NRRD") using NRRD

May solve the issue.

On Jul 6, 2017 1:36 PM, "kellyhill90" notifications@github.com wrote:

Do you think we need to clone this NRRD package directly? https://github.com/JuliaIO/NRRD.jl

Even after Using Images, FileIO, I still get UndefVarError: NRRD not defined

— You are receiving this because you commented.

Reply to this email directly, view it on GitHub https://github.com/HolyLab/BlockRegistrationScheduler/issues/38#issuecomment-313482416, or mute the thread https://github.com/notifications/unsubscribe-auth/AMvfFQB-T-xe9uGedAqx8C4MYjiAgrFDks5sLSkpgaJpZM4OJ9dw .

kellyhill90 commented 7 years ago

Fantastic. Got it working. Now I can play around with optimization. Thanks for your help, Jerry!

timholy commented 7 years ago

@mdhe1248, can I copy your example and paste to the NRRD README? Seems worth documenting somewhere (and I've been lazy about documenting it).

mdhe1248 commented 7 years ago

Sure!

On Jul 7, 2017 5:26 PM, "Tim Holy" notifications@github.com wrote:

@mdhe1248 https://github.com/mdhe1248, can I copy your example and paste to the NRRD README? Seems worth documenting somewhere (and I've been lazy about documenting it).

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/HolyLab/BlockRegistrationScheduler/issues/38#issuecomment-313807359, or mute the thread https://github.com/notifications/unsubscribe-auth/AMvfFfbrwtAE0P0yebEk4xmqScU4YCOvks5sLrCKgaJpZM4OJ9dw .