These are chat archives for thunder-project/thunder

2nd
Jun 2016
Nikita Vladimirov
@nvladimus
Jun 02 2016 18:22
@d-v-b , have you tried registration with new Thunder? I am running into mismatch between dimensions because squeeze() keeps one extra dim:
reference = datCrop[refR[0]:refR[1],:,:,:].mean() # taking 20 stacks average as a ref stack
reference.shape
(1, 500, 500, 20)
datCrop.shape
(1200, 500, 500, 20)
register = CrossCorr()
model = register.fit(datCrop, reference=reference)

Exception: Image shape (500, 500, 20) and reference shape (1, 500, 500, 20) must match
the thunder squeeze() keeps 1-st dimention in ref. image no matter what :(
Jeremy Freeman
@freeman-lab
Jun 02 2016 18:26
@nvladimus try calling toarray() on the reference
so do reference = datCrop[refR[0]:refR[1],:,:,:].mean().toarray()
i realize this squeeze() thing has become confusing, and we might need to change it
Nikita Vladimirov
@nvladimus
Jun 02 2016 18:28
it works, thanks, @freeman-lab
Nikita Vladimirov
@nvladimus
Jun 02 2016 18:45

sorry, I am stuck again. In older Thunder version, registration model transformations were lists, and I used to extract the xy shift from them like

xy_shifts = np.array([model.transformations[x].delta for x in range(len(model.transformations))])

This would give me an array of xy displacements for each stack, ordered by Z. But now it seems to be a dictionary or some key-value pairs. How do I access them now?

Nikita Vladimirov
@nvladimus
Jun 02 2016 18:50

My full code would be, in old Thunder, like

zStart = 10 #starting plane of registration 
zEnd = 30 #ending plane
window = 500 #width and height of the registration ROI
datCrop = dat.crop((dims[0]/3,dims[1]/2,zStart), \
                       (dims[0]/3 + window,dims[1]/2 + window,zEnd)) # take sub-stack, (length, width, z)
regCrop = Registration('planarcrosscorr').prepare(datCrop,startIdx=refR[0], stopIdx=refR[1])
regParamsCrop = regCrop.fit(datCrop.medianFilter(3))

    # extrapolate all missing planes from known registration shifts
xysCrop = np.array([regParamsCrop.transformations[x].delta for x in range(len(regParamsCrop.transformations))])
xysNew = np.zeros((numFrames,dims[2],2))
for i in range(zStart): #extrapolate
        xysNew[:,i,:] = np.median(xysCrop,axis = 1)
for i in range(zStart,zEnd): # populate with real values
        xysNew[:,i,:] = xysCrop[:,i-zStart,:]
for i in range(zEnd,dims[2]): #extrapolate
        xysNew[:,i,:] = np.median(xysCrop,axis = 1)  

# apply registration to full data    
for x in range(xysNew.shape[0]):
        regParamsCrop.transformations[x].delta = list(xysNew[x,:,:])
dat = regParamsCrop.transform(dat)

xys = np.array([regParamsCrop.transformations[x].delta for x in range(len(regParamsCrop.transformations))])

distsX, distsY = xys[:,:,0], xys[:,:,1]

and then I would plot the distsX, distsY to know how much the sample shifted during experiment. Is there any way to do the same thing in new Thunder?

Nikita Vladimirov
@nvladimus
Jun 02 2016 19:13
I know this is a lot of code, sorry.
I just need an example how to extract displacements from the model.transformations, if possible.
Kyle
@kr-hansen
Jun 02 2016 20:07

@nvladimus You should be able to use the toarray() method mentioned above. Do something like:
model.toarray()[:,0] for x and
model.toarray()[:,1] for y values.

.toarray() will make it into a numpy array (or bolt array possibly with distributed data, I'm not sure off the top of my head) so you don't need to call np.array() on it.

At least for the registration method with Thunder 1.0, the examples have been updated on the docs page (http://docs.thunder-project.org/) under the registration tutorial. I was able to work out most of the changes by looking through those tutorials.

Nikita Vladimirov
@nvladimus
Jun 02 2016 20:15
cool, thanks! I will try now.
Nikita Vladimirov
@nvladimus
Jun 02 2016 20:49
sorry, the new registration method does not make much sense to me. It used to register stacks plane-by-plane, and for each stack I would get 2 * N shifts (N being number of planes in a stack, shifts were in x, y). The new registration method does it by 3D volume, and for every stack returns 3 shifts, now in x,y,z, but they do not agree with previously calculated (correct) plane-by-plane shifts.
e.g. planes were shifted by 15 pixels in x, but the new method reports shifting by 25 pixels
Davis Bennett
@d-v-b
Jun 02 2016 20:52
what is the new method? I haven't used registration in thunder 1.0
Nikita Vladimirov
@nvladimus
Jun 02 2016 20:53
# run registration on cropped data
zStart = 10 #starting plane of registration 
zEnd = 30 #ending plane
window = 500 #width and height of the registration ROI
datCrop = dat[:,dims[0]/3:dims[0]/3 + window,dims[1]/2:dims[1]/2 + window,zStart:zEnd] # take sub-stack, (length, width, z)
reference = datCrop[refR[0]:refR[1],:,:,:].mean().toarray()

register = CrossCorr()
model = register.fit(datCrop.median_filter(3), reference=reference)

deltas = model.toarray()
plt.plot(deltas[:,1]) #plot x shifts for every time point
Davis Bennett
@d-v-b
Jun 02 2016 20:56
so it looks like unless you specify otherwise the registration is 3D
Jeremy Freeman
@freeman-lab
Jun 02 2016 20:57
@d-v-b yeah CrossCorr now has an option axis which allows you to restrict the registration to a subset of axes
Nikita Vladimirov
@nvladimus
Jun 02 2016 20:57
ah, ok, I was looking for this option in register.fit() function. Will try now.
Jeremy Freeman
@freeman-lab
Jun 02 2016 20:58
yeah i think CrossCorr(axis=2) should do the "planar" registration from before
yeah in general arguments that control how the algorithm behaviors are part of the constructor
probably need to better document this
Davis Bennett
@d-v-b
Jun 02 2016 20:58
and it wouldn't surprise me if 3D cross-correlation performs worse for a single plane than 2D; cross-correlation is good for aligning photographs, but not great for fluorescence imaging (in my experience)
This message was deleted
Jeremy Freeman
@freeman-lab
Jun 02 2016 21:01
oh yeah you definitely want to use the planar stuff
the goal of the new argument was to make it more general
whereas before we just sort of defaulted to xy but not z as the default
Nikita Vladimirov
@nvladimus
Jun 02 2016 21:39
no, it still prefers +25 or -25 values for the shifts. Can it be because of wrong type casting, from int16?
Nikita Vladimirov
@nvladimus
Jun 02 2016 21:49
casting to float64 does not help either
Davis Bennett
@d-v-b
Jun 02 2016 21:52
are you saying that it only returns shift values that are +- 25 pixels?
Nikita Vladimirov
@nvladimus
Jun 02 2016 21:52
I think calculating the ref stack is messed up, I am looking into it.
Jeremy Freeman
@freeman-lab
Jun 02 2016 21:55
ok cool
Nikita Vladimirov
@nvladimus
Jun 02 2016 21:56

ok, there is something wrong in how I load the data: dat = td.images.frombinary(rawDir, shape = dims, ext='stack', npartitions=numFrames, engine=sc)
then I try to plot one plane:

plt.figure(figsize = (10,20))
plt.imshow(dat[0,:,:,20].toarray())

and I get some mosaic image instead of fish brain.

Jeremy Freeman
@freeman-lab
Jun 02 2016 21:57
well that definitely sounds like the issue!
Nikita Vladimirov
@nvladimus
Jun 02 2016 21:57
the dimensions look right, but the data gets jumbled for some reason.
Davis Bennett
@d-v-b
Jun 02 2016 21:58
what order are you using for dims?
Jeremy Freeman
@freeman-lab
Jun 02 2016 21:58
yeah it could be you need to reverse them
Davis Bennett
@d-v-b
Jun 02 2016 21:58
could be [x y z] vs [z y x]
Nikita Vladimirov
@nvladimus
Jun 02 2016 21:58
[2048, 1024, 41]
[x,y,z]
ok, will try
aha! Now it looks right.
Nikita Vladimirov
@nvladimus
Jun 02 2016 22:04
can I swap axes in the dat?
like in numpy array?