forked from agiovann/CalBlitz
-
Notifications
You must be signed in to change notification settings - Fork 0
/
DemoMotionCorrection.py
138 lines (115 loc) · 4.07 KB
/
DemoMotionCorrection.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
# -*- coding: utf-8 -*-
"""
Created on Tue Apr 26 12:13:18 2016
@author: agiovann
"""
# -*- coding: utf-8 -*-
"""
Created on Tue Jun 30 20:56:07 2015
@author: agiovann
"""
import sys
#% add required packages
import h5py
import calblitz as cb
import time
import pylab as pl
import numpy as np
#% set basic ipython functionalities
try:
pl.ion()
%load_ext autoreload
%autoreload 2
except:
print "Probably not a Ipython interactive environment"
#%% define movie
filename='movies/demoMovie_PC.tif'
frameRate=15.62;
start_time=0;
#%%
filename_py=filename[:-4]+'.npz'
filename_hdf5=filename[:-4]+'.hdf5'
filename_mc=filename[:-4]+'_mc.npz'
#%% load movie
m=cb.load(filename, fr=frameRate,start_time=start_time);
#%% load portion of the movie or only some channels
#m=cb.load(filename, fr=frameRate,start_time=start_time,subindices=range(0,1500,10));
#%% save movie hdf5 format. fastest
m.save(filename_hdf5)
#%% backend='opencv' is much faster
m.play(fr=100,gain=15.0,magnification=1,backend='pylab')
#%%
low_SNR=True
if low_SNR:
N=1000000
mn1=m.copy().bilateral_blur_2D(diameter=5,sigmaColor=10000,sigmaSpace=0)
mn1,shifts,xcorrs, template=mn1.motion_correct()
mn2=mn1.apply_shifts(shifts)
#mn1=cb.movie(np.transpose(np.array(Y_n),[2,0,1]),fr=30)
mn=cb.concatenate([mn1,mn2],axis=1)
#%% automatic parameters motion correction
max_shift_h=10; # maximum allowed shifts in y
max_shift_w=10; # maximum allowed shifts in x
m=cb.load(filename_hdf5);
m,shifts,xcorrs,template=m.motion_correct(max_shift_w=max_shift_w,max_shift_h=max_shift_h, num_frames_template=None, template = None,method='opencv')
# removie the borders
max_h,max_w= np.max(shifts,axis=0)
min_h,min_w= np.min(shifts,axis=0)
m=m.crop(crop_top=max_h,crop_bottom=-min_h+1,crop_left=max_w,crop_right=-min_w,crop_begin=0,crop_end=0)
np.savez(filename_mc,template=template,shifts=shifts,xcorrs=xcorrs,max_shift_h=max_shift_h,max_shift_w=max_shift_w)
pl.subplot(2,1,1)
pl.plot(shifts)
pl.ylabel('pixels')
pl.subplot(2,1,2)
pl.plot(xcorrs)
pl.ylabel('cross correlation')
pl.xlabel('Frames')
#%%
m.play(fr=50,gain=5.0,magnification=1,backend='opencv')
#%% IF YOU WANT MORE CONTROL USE THE FOLLOWING
# motion correct for template purpose. Just use a subset to compute the template
m=cb.load(filename_hdf5);
min_val_add=np.min(np.mean(m,axis=0)) # movies needs to be not too negative
m=m-min_val_add
every_x_frames=1 # for memory/efficiency reason can be increased
submov=m[::every_x_frames,:]
max_shift_h=10;
max_shift_w=10;
templ=np.nanmedian(submov,axis=0); # create template with portion of movie
shifts,xcorrs=submov.extract_shifts(max_shift_w=max_shift_w, max_shift_h=max_shift_h, template=templ, method='opencv') #
submov.apply_shifts(shifts,interpolation='cubic')
template=(np.nanmedian(submov,axis=0))
shifts,xcorrs=m.extract_shifts(max_shift_w=max_shift_w, max_shift_h=max_shift_h, template=template, method='opencv') #
m=m.apply_shifts(shifts,interpolation='cubic')
template=(np.median(m,axis=0))
pl.imshow(template,cmap=pl.cm.gray)
#%% now use the good template to correct
max_shift_h=10;
max_shift_w=10;
m=cb.load(filename_hdf5);
min_val_add=np.float32(np.percentile(m,.0001))
m=m-min_val_add
shifts,xcorrs=m.extract_shifts(max_shift_w=max_shift_w, max_shift_h=max_shift_h, template=template, method='opencv') #
if m.meta_data[0] is None:
m.meta_data[0]=dict()
m.meta_data[0]['shifts']=shifts
m.meta_data[0]['xcorrs']=xcorrs
m.meta_data[0]['template']=template
m.meta_data[0]['min_val_add']=min_val_add
m.save(filename_hdf5)
#%% visualize shifts
pl.subplot(2,1,1)
pl.plot(shifts)
pl.subplot(2,1,2)
pl.plot(xcorrs)
#%% reload and apply shifts
m=cb.load(filename_hdf5)
meta_data=m.meta_data[0];
shifts,min_val_add,xcorrs=[meta_data[x] for x in ['shifts','min_val_add','xcorrs']]
m=m.apply_shifts(shifts,interpolation='cubic')
#% crop borders created by motion correction
max_h,max_w= np.max(shifts,axis=0)
min_h,min_w= np.min(shifts,axis=0)
m=m.crop(crop_top=max_h,crop_bottom=-min_h+1,crop_left=max_w,crop_right=-min_w,crop_begin=0,crop_end=0)
#%% play movie
m.play(fr=50,gain=3.0,magnification=1)