Add my experience#42
Conversation
I didn't include any function or script. Should I? I thought my script or functions were not suitable for this page.
| sigmahp = Float64[sigma/x for x in ps] #Highpass filter | ||
| sigmalp = [3, 3, 0] #Lowpass filter | ||
| λ = 1e-4 | ||
|
|
There was a problem hiding this comment.
Maybe show the function call that these parameters are used in?
There was a problem hiding this comment.
I thought it was redundant with the example above. I will add function calls sometime soon and see if it is informative.
| 1. Choose good images for registration: | ||
| I selected images that do not show evoked calcium activity. There are still spontaneous activity. | ||
|
|
||
| 2. Median or quantile filtering across time: |
There was a problem hiding this comment.
Of the images or something else?
There was a problem hiding this comment.
I didn't use median filter in Images. This is because I didn't filter across entire time. For example, if I think stacks at around t = 30 and t = 60 are good, I median filtered stacks at t = 27:33 and t = 57:63, respectively. However, running median filter across entire time should be okay. It may take longer time.
There was a problem hiding this comment.
I didn't mean the Images package, I meant "please clarify whether you're median-filtering the raw data or the deformation that comes from registration."
There was a problem hiding this comment.
Okay I will clarify that.
I thought having images with less noise might give better results on mismatch calculation. Is this mathematically the same as medial-filtering the deformation?
There was a problem hiding this comment.
Definitely not. I'm not sure which would be better. But reducing noise to the point that pixels are meaningful is a good idea. Another approach would be to use restrict and work with smaller images, since restrict performs averaging and hence also reduces noise.
| This is to reduce noise and spontaneous activity. | ||
|
|
||
| 3. Replace too high intensity pixels with NaN: | ||
| Sometimes, I observed high-intensity objects moving around tissue surface. While these things could be biologically important features, this is disastrous for registration. I ended up with replacing high intensity object (or pixels) with NaN. |
There was a problem hiding this comment.
Manually in specific regions, or automatically based purely on intensity? Meaning, was it just img[img .> thresh] = NaN?
There was a problem hiding this comment.
Yes. It was just img[img .> thresh] = NaN
I wonder replacing many pixels with NaN will still be fine. I also wonder how block moves if entire pixels are NaN.
There was a problem hiding this comment.
I think there are checks for this, hopefully it should work as long as there are still some informative blocks.
Example script. This might be too lengthy?
| algorithm = Apertures[Apertures(fixed, knots, mxshift, λ, pp; pid=wpids[i], correctbias=false) for i = 1:length(wpids)] #Notice that correctbias = false | ||
|
|
||
| #### In warping step, | ||
| ϕ_s = medfilt(griddeformations(u, knots), 1) #Basically there is no filtering. Notice the input parameter is 1. |
There was a problem hiding this comment.
ϕ_s1 = medfilt(griddeformations(u, knots), 1)
ϕ_s2 = griddeformations(u, knots)
isequal(ϕ_s1, ϕ_s2) # this gives `false`
#However,
[isequal(ϕ_s1[i].u, ϕ_s2[i].u) for i in eachindex(ϕ_s1)] # returns trues
[isequal(ϕ_s1[i].knots, ϕ_s2[i].knots) for i in eachindex(ϕ_s1)] # returns trues
#Actually,
ϕ_s3 = griddeformations(u, knots)
isequal(ϕ_s2, ϕ_s3) # this returns false!
`mappedarray` is used for intensity cut.
| #### 3. High intensity thresholding | ||
| img0 = load("exp1_med.imagine") #Load the filtered image | ||
| img1 = mappedarray(val -> val > 140 ? NaN : val, img0); # Replace high intensity pixels with NaN. 140 is a threshold. NaN is a value for replacement. | ||
| img1 = AxisArray(img1, (:x, :l, :z, :time), (0.5770μm, 0.5770μm, 5μm, 2s)); # Assign axes again. |
There was a problem hiding this comment.
Despite img0 has axis information, axis should be re-definition for img1 for the downstream process. e.g.) timeaxis(img1) returns nothing. Is there good ways to copy axes from other image?
There was a problem hiding this comment.
I've been noticing this as a problem too. The manual approach would be
img1 = AxisArray(mappedarray(f, img0.data), axes(img0))But I've even been considering making mappedarray do this automatically. At the moment I'm holding back, though (any time you don't do the "obvious" thing it's a little bit scary), and looking for an idea for how to make the order of wrappers more easily controlled by the user.
There was a problem hiding this comment.
Okay.
I will just replace two lines of my script with your manual approach.
| tindex0 = [20; stimidx; nimages(img0)-50] #this will be used for registration | ||
|
|
||
| #### 2. Temporal median filtering (Run only one time) | ||
| tmedian_filter(Float32, "exp1_med.cam", img0, collect(-3:3), tindex0) #in Jerry_RegisterUtils; See `?tmedian_filter` |
There was a problem hiding this comment.
Is it possible to use mappedarray here?
There was a problem hiding this comment.
No, mappedarray operates element-by-element. You could write a different type if you wanted. (See the code for MappedArrays, it is almost trivial.)
There was a problem hiding this comment.
I will stick to my code for now; I will just comment StreamingContainer in my script.
To add axis information after mappedarray.
2 participants
I didn't include any function or script. Should I? I thought my script or functions were not suitable for this page.