Mercurial > repos > gregor.m > spyboat

comparison spyboat_cli.py @ 3:4d7f30a7e2f6 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit d070f1ba04a5141a65487f45b29c22767639e44b"
author gregor.m
date Tue, 24 Nov 2020 13:06:26 +0000
parents
children a4c6fcf2c456
comparison
equal deleted inserted replaced
2:c59d1373230e 3:4d7f30a7e2f6
1 #!/usr/bin/env python
2
3 ## Gets interfaced by Galaxy or can be used for bash scripting
4 import argparse
5 import sys, os
6 import logging
7
8 from skimage import io
9 from numpy import float32
10
11 import spyboat
12 import output_report
13
14 logging.basicConfig(level=logging.INFO, stream=sys.stdout, force=True)
15 logger = logging.getLogger('wrapper')
16
17 # ----------command line parameters ---------------
18
19 parser = argparse.ArgumentParser(description='Process some arguments.')
20
21 # I/O
22 parser.add_argument('--input_path', help="Input movie location", required=True)
23 parser.add_argument('--phase_out', help='Phase output file name', required=True)
24 parser.add_argument('--period_out', help='Period output file name', required=True)
25 parser.add_argument('--power_out', help='Power output file name', required=True)
26 parser.add_argument('--amplitude_out', help='Amplitude output file name', required=True)
27 parser.add_argument('--preprocessed_out', help="Preprocessed-input output file name, 'None'", required=False)
28
29
30 # (Optional) Multiprocessing
31
32 parser.add_argument('--ncpu', help='Number of processors to use',
33 required=False, type=int, default=1)
34
35 # Optional spatial downsampling
36 parser.add_argument('--rescale', help='Rescale the image by a factor given in %%, None means no rescaling',
37 required=False, type=int)
38 # Optional Gaussian smoothing
39 parser.add_argument('--gauss_sigma', help='Gaussian smoothing parameter, None means no smoothing', required=False, type=float)
40
41 # Wavelet Analysis Parameters
42 parser.add_argument('--dt', help='Sampling interval', required=True, type=float)
43 parser.add_argument('--Tmin', help='Smallest period', required=True, type=float)
44 parser.add_argument('--Tmax', help='Biggest period', required=True, type=float)
45 parser.add_argument('--nT', help='Number of periods to scan for', required=True, type=int)
46
47 parser.add_argument('--Tcutoff', help='Sinc cut-off period, disables detrending if not set', required=False, type=float)
48 parser.add_argument('--win_size', help='Sliding window size for amplitude normalization, None means no normalization',
49 required=False, type=float)
50
51 # Optional masking
52 parser.add_argument('--masking', help="Set to either 'dynamic', 'static' or 'None' which is the default", default='None', required=False, type=str)
53
54 parser.add_argument('--mask_frame',
55 help="The frame of the input movie to create a static mask from, needs masking set to 'static'",
56 required=False, type=int)
57
58
59 parser.add_argument('--mask_thresh', help='The threshold of the mask, all pixels with less than this value get masked (if masking enabled).',
60 required=False, type=float,
61 default=0)
62
63 # output html report/snapshots
64 parser.add_argument('--html_fname', help="Name of the html report.",
65 default='OutputReport.html', required=False, type=str)
66
67 parser.add_argument('--report_img_path', help="For the html report, to be set in Galaxy. Without galaxy leave at cwd!", default='.', required=False, type=str)
68
69 parser.add_argument('--version', action='version', version='0.1.0')
70
71 arguments = parser.parse_args()
72
73 logger.info("Received following arguments:")
74 for arg in vars(arguments):
75 logger.info(f'{arg} -> {getattr(arguments, arg)}')
76
77 # ------------Read the input----------------------------------------
78 try:
79 movie = spyboat.open_tif(arguments.input_path)
80 except FileNotFoundError:
81 logger.critical(f"Couldn't open {arguments.input_path}, check movie storage directory!")
82
83 sys.exit(1)
84
85 # -------- Do (optional) spatial downsampling ---------------------------
86
87 scale_factor = arguments.rescale
88
89 # defaults to None
90 if not scale_factor:
91 logger.info('No downsampling requested..')
92
93 elif 0 < scale_factor < 100:
94 logger.info(f'Downsampling the movie to {scale_factor:d}% of its original size..')
95 movie = spyboat.down_sample(movie, scale_factor / 100)
96 else:
97 raise ValueError('Scale factor must be between 0 and 100!')
98
99 # -------- Do (optional) pre-smoothing -------------------------
100 # note that downsampling already is a smoothing operation..
101
102 # check if pre-smoothing requested
103 if not arguments.gauss_sigma:
104 logger.info('No pre-smoothing requested..')
105 else:
106 logger.info(f'Pre-smoothing the movie with Gaussians, sigma = {arguments.gauss_sigma:.2f}..')
107
108 movie = spyboat.gaussian_blur(movie, arguments.gauss_sigma)
109
110 # ----- Set up Masking before processing ----
111
112 mask = None
113 if arguments.masking == 'static':
114 if not arguments.mask_frame:
115 logger.critical("Frame number for static masking is missing!")
116 sys.exit(1)
117
118 if (arguments.mask_frame > movie.shape[0]) or (arguments.mask_frame < 0):
119 logger.critical(f'Requested frame does not exist, input only has {movie.shape[0]} frames.. exiting')
120 sys.exit(1)
121
122 else:
123 logger.info(f'Creating static mask from frame {arguments.mask_frame} with threshold {arguments.mask_thresh}')
124 mask = spyboat.create_static_mask(movie, arguments.mask_frame,
125 arguments.mask_thresh)
126 elif arguments.masking == 'dynamic':
127 logger.info(f'Creating dynamic mask with threshold {arguments.mask_thresh}')
128 mask = spyboat.create_dynamic_mask(movie, arguments.mask_thresh)
129
130 else:
131 logger.info('No masking requested..')
132
133 # ------ Retrieve wavelet parameters ---------------------------
134
135 Wkwargs = {'dt': arguments.dt,
136 'Tmin': arguments.Tmin,
137 'Tmax': arguments.Tmax,
138 'nT': arguments.nT,
139 'T_c' : arguments.Tcutoff, # defaults to None
140 'win_size' : arguments.win_size # defaults to None
141 }
142
143 # --- start parallel processing ---
144
145 results = spyboat.run_parallel(movie, arguments.ncpu, **Wkwargs)
146
147 # --- masking? ---
148
149 if mask is not None:
150 # mask all output movies (in place!)
151 for key in results:
152 logger.info(f'Masking {key}')
153 spyboat.apply_mask(results[key], mask, fill_value=-1)
154
155 # --- Produce Output HTML Report Figures/png's ---
156
157 # create the directory, yes we have to do that ourselves :)
158 # galaxy then magically renders the html from that directory
159 try:
160
161 if arguments.report_img_path != '.':
162 logger.info(f'Creating report directory {arguments.report_img_path}')
163 os.mkdir(arguments.report_img_path)
164
165 # make 6 times 4 snapshots each
166 Nsnap = 6
167 NFrames = movie.shape[0]
168 # show only frames at least one Tmin
169 # away from the edge (-effects)
170 start_frame = int(Wkwargs['Tmin']/Wkwargs['dt'])
171
172 if (start_frame > NFrames//2):
173 logger.warning("Smallest period already is larger than half the observation time!")
174 # set to 0 in this case
175 start_frame = 0
176
177 frame_increment = int( (NFrames - 2*start_frame) / Nsnap)
178 snapshot_frames = range(start_frame, NFrames - start_frame, frame_increment)
179
180 for snapshot_frame in snapshot_frames:
181 output_report.produce_snapshots(movie, results, snapshot_frame, Wkwargs, img_path=arguments.report_img_path)
182
183 output_report.produce_distr_plots(results, Wkwargs, img_path=arguments.report_img_path)
184
185 output_report.create_html(snapshot_frames, arguments.html_fname)
186
187 except FileExistsError as e:
188 logger.critical(f"Could not create html report directory: {repr(e)}")
189
190
191 # --- save out result movies ---
192
193 # save phase movie
194 io.imsave(arguments.phase_out, results['phase'], plugin="tifffile")
195 logger.info(f'Written {arguments.phase_out}')
196 # save period movie
197 io.imsave(arguments.period_out, results['period'], plugin="tifffile")
198 logger.info(f'Written {arguments.period_out}')
199 # save power movie
200 io.imsave(arguments.power_out, results['power'], plugin="tifffile")
201 logger.info(f'Written {arguments.power_out}')
202 # save amplitude movie
203 io.imsave(arguments.amplitude_out, results['amplitude'], plugin="tifffile")
204 logger.info(f'Written {arguments.amplitude_out}')
205
206 # save out the probably pre-processed (scaled and blurred) input movie for
207 # direct comparison to results and coordinate mapping etc.
208 if arguments.preprocessed_out:
209 io.imsave(arguments.preprocessed_out, movie, plugin='tifffile')
210 logger.info(f'Written {arguments.preprocessed_out}')