Question: Does the image data type matter?

[mdhe1248]

Hi Tim and Cody, I have two questions:

1. I am wondering if data type affects registration speed? For example, share array results in faster registration than array?
2. Another question is how can I save new image into shared array? Can I use something like this line: f = open(fn, "w+") output = Mmap.mmap(f, Array{Float32,4}, (size(img)[1:3]..., length(tindex)))

Best, Jerry

[timholy]

Hi Jerry,

A SharedArray means an Array that can be used efficiently by multiple julia processes (e.g., julia -p 4 or addprocs). If you're only using one process, there is no difference in usage speed. (SharedArrays are slower to construct in the first place, but for a long-running job like registration you won't notice the few milliseconds it takes.)

You can create a SharedArray from a file or by specifying its element type and size, see ?SharedArray.

[mdhe1248]

Okay. Thanks Tim. Registration speed increases a lot when I use SharedArray and multiple julia processes.

Here might be a problem and potential solution: Problem: When I use subim from an Array, registration freezes. Solution: If I use subim from SharedArray, then the registration doesn't freeze. I always used multiple threads. Thus, I am not sure if single core registration can also solve this problem.

[timholy]

I don't understand what you mean by "registration freezes." Does it just fail to make progress? If you hit Ctrl-C, where does it break?

[mdhe1248]

Hi Tim, Ctrl-C didn't work immediately, so I killed the process outside Julia. Cody had the same issue when he used subim, and we end up with not using subim.

The process usually froze at the mismatch step: For example, I could see the messages below, but did not proceed further. From worker 2: Worker 2 is working on 1 From worker 3: Worker 3 is working on 2 From worker 4: Worker 4 is working on 3 From worker 6: Worker 6 is working on 5 From worker 5: Worker 5 is working on 4

[timholy]

How are you calling it when this happens? Can you paste your script?

[mdhe1248]

Oh, Tim, could you see this script below instead of the script above? I will delete the script above. It has many errors.

wpids = addprocs(4)  # use 8 worker processes (no CUDA)

using Images, SIUnits.ShortUnits, FixedSizeArrays, JLD, MAT
using BlockRegistration, BlockRegistrationScheduler
using RegisterWorkerApertures, RegisterMismatch
using ImagePlayer
####### mismatch
fn = "/mnt/donghoon_G001/sampleimg/TSeries-03192016-1142-029.tif" # filename ****
savedir = "/mnt/donghoon_G001/sampleimg"
fileoutext = ".register"
img0 = load(fn);
img = Image(SharedArray(fn, Float32, size(img0); mode = "r"), properties(img0));
tindex0 = 1:size(img,3) #time
fixedidx0 = Int16(ceil(length(tindex0)*0.50)); # Choose fixedindx from the tindex ****

# Select moving stacks
tindex = tindex0;
tindex = tindex0[[1; round(Int16, fixedidx0/2); fixedidx0; round(Int16, fixedidx0/2)+fixedidx0; end-1]]   #for test

# Snip out a region of interest (discard empty portions of the image stack)
roi = (1:250,1:80)

##### Select img: different 'img' (even shared array) other than img0 will give errors.
  img = img0;   # Array
#  img = subim(img0, roi..., tindex)    # Array with roi
#  img = img   # Shared array: entire image
#  img = subim(img, roi..., tindex)  # Shared array: you could alternatively select a subset of times

# Select "fixed" image 
fixedidx = find(tindex .== fixedidx0)
fixed0 = getindexim(img, "t", fixedidx[1])

# Pixel space
ps = [0.11e-6m 0.11e-6m]

# Define preprocessing
sigmahp = Float64[25e-6m/x for x in ps]
sigmalp = [2, 2]
pp = PreprocessSNF(0, sigmalp, sigmahp)
fixed = copy(pp(fixed0))  # use copy because of julia #14625

# grid size
gridspacing =5e-6m
gridsize = (ceil(Int64,[size(img)...][1:2].*ps'./gridspacing)...)
knots = map(d->linspace(1,size(fixed,d),gridsize[d]), (1:ndims(fixed)...))

# max shift
mxshift = (5, 5)

# Choose volume regularization penalty. See the README for
# BlockRegistration and `fixed_λ`.
λ = 1e-4

# Set aperture overlap.
apertureoverlap = 0.0;  #Aperture overlap percentage (between 0 and 1)
aperture_width = default_aperture_width(fixed, gridsize)
overlap = map(x->round(Int64,x*apertureoverlap), aperture_width)

# Create the worker algorithm structures. We assign one per worker process.
algorithm = Apertures[Apertures(fixed, knots, mxshift, λ, pp; overlap=overlap, pid=wpids[i], correctbias=false) for i = 1:length(wpids)]
mon = monitor(algorithm, (), Dict{Symbol,Any}(:mismatch => Float64, :u=>ArrayDecl(Array{Vec{2,Float64},2}, gridsize)))

# Define the output file and run the job
basename = splitext(splitdir(fn)[2])[1]
@show basename
fileout = string(savedir,"/", basename,fileoutext)
@time driver(fileout, algorithm, img, mon)

# Append important extra information to the file
jldopen(fileout, "r+") do io
    write(io, "imagefile", fn)
    write(io, "roi", roi)
    write(io, "fixedidx", fixedidx)
    write(io, "knots", knots)
    write(io, "sigmalp", sigmalp)
    write(io, "sigmahp", sigmahp)
    write(io, "tindex", tindex)
    write(io, "gridsize", gridsize)
    write(io, "λ", λ)
    write(io, "apertureoverlap", apertureoverlap)
    write(io, "mxshift", mxshift)
end

[timholy]

For me that script seems to be running fine: it doesn't illustrate the freezing you mentioned.

A couple of notes:

• Creating a SharedArray from a file only works if you know that the data are written to disk in exactly the same way the same bytes are stored in memory. So if there's any file formatting whatsoever (there is with a TIF) file, you cannot use SharedArray(filename, ...). To learn more, read about memory mapping. The Wikipedia page isn't really written in an introductory style, but perhaps you can find one (http://www.mathworks.com/help/matlab/import_export/overview-of-memory-mapping.html?refresh=true looks better). So I used img0 for my tests.
• One small detail is to create integer ranges rather than floating-point ranges. Observe:

julia> 1:8/2
1.0:1.0:4.0

julia> 1:8÷2
1:4

You can get that "integer division" operator, ÷, by typing "\div" and then hitting TAB.

[timholy]

Fixed on vivid by backporting https://github.com/JuliaLang/julia/pull/14625.