diff --git a/tierpsy/analysis/compress/compressVideo.py b/tierpsy/analysis/compress/compressVideo.py index 10179056..951bea69 100755 --- a/tierpsy/analysis/compress/compressVideo.py +++ b/tierpsy/analysis/compress/compressVideo.py @@ -122,16 +122,16 @@ def normalizeImage(img): # normalise image intensities if the data type is other # than uint8 image = image.astype(np.double) - + imax = img.max() imin = img.min() factor = 255/(imax-imin) - + imgN = ne.evaluate('(img-imin)*factor') imgN = imgN.astype(np.uint8) return imgN, (imin, imax) - + def reduceBuffer(Ibuff, is_light_background): if is_light_background: return np.min(Ibuff, axis=0) @@ -175,36 +175,36 @@ def createImgGroup(fid, name, tot_frames, im_height, im_width, is_expandable=Tru return img_dataset -def initMasksGroups(fid, expected_frames, im_height, im_width, +def initMasksGroups(fid, expected_frames, im_height, im_width, attr_params, save_full_interval, is_expandable=True): # open node to store the compressed (masked) data mask_dataset = createImgGroup(fid, "/mask", expected_frames, im_height, im_width, is_expandable) - + tot_save_full = (expected_frames // save_full_interval) + 1 full_dataset = createImgGroup(fid, "/full_data", tot_save_full, im_height, im_width, is_expandable) full_dataset._v_attrs['save_interval'] = save_full_interval - + assert all(x in ['expected_fps', 'is_light_background', 'microns_per_pixel'] for x in attr_params) set_unit_conversions(mask_dataset, **attr_params) set_unit_conversions(full_dataset, **attr_params) if is_expandable: - mean_intensity = fid.create_earray('/', + mean_intensity = fid.create_earray('/', 'mean_intensity', atom=tables.Float32Atom(), shape=(0,), expectedrows=expected_frames, filters=TABLE_FILTERS) else: - mean_intensity = fid.create_carray('/', + mean_intensity = fid.create_carray('/', 'mean_intensity', atom=tables.Float32Atom(), shape=(expected_frames,), filters=TABLE_FILTERS) - + return mask_dataset, full_dataset, mean_intensity @@ -230,13 +230,13 @@ def compressVideo(video_file, masked_image_file, mask_param, expected_fps=25, ''' #get the default values if there is any bad parameter - output = compress_defaults(masked_image_file, - expected_fps, - buffer_size = buffer_size, + output = compress_defaults(masked_image_file, + expected_fps, + buffer_size = buffer_size, save_full_interval = save_full_interval) - buffer_size = output['buffer_size'] - save_full_interval = output['save_full_interval'] + buffer_size = output['buffer_size'] + save_full_interval = output['save_full_interval'] if len(bgnd_param) > 0: is_bgnd_subtraction = True @@ -253,14 +253,14 @@ def compressVideo(video_file, masked_image_file, mask_param, expected_fps=25, # processes identifier. base_name = masked_image_file.rpartition('.')[0].rpartition(os.sep)[-1] - + # select the video reader class according to the file type. vid = selectVideoReader(video_file) - + # delete any previous if it existed with tables.File(masked_image_file, "w") as mask_fid: pass - + #Extract metadata if is_extract_timestamp: # extract and store video metadata using ffprobe @@ -270,7 +270,7 @@ def compressVideo(video_file, masked_image_file, mask_param, expected_fps=25, else: expected_frames = 1 - + # Initialize background subtraction if required if is_bgnd_subtraction: @@ -282,10 +282,10 @@ def compressVideo(video_file, masked_image_file, mask_param, expected_fps=25, frame_number = 0 full_frame_number = 0 image_prev = np.zeros([]) - + # Initialise FOV splitting if needed if is_fov_tosplit: - # masked video does not exist yet so have to initialise from data + # masked video does not exist yet so have to initialise from data # use either background or first frame if is_bgnd_subtraction: img_fov = bgnd_subtractor.bgnd.astype(np.uint8) @@ -293,17 +293,17 @@ def compressVideo(video_file, masked_image_file, mask_param, expected_fps=25, ret, img_fov = vid.read() # close and reopen the video, to restart from the beginning vid.release() - vid = selectVideoReader(video_file); - # TODO: change class creator so it only needs the video name? by using + vid = selectVideoReader(video_file); + # TODO: change class creator so it only needs the video name? by using # Tierpsy's functions such as selectVideoReader it can then read the first image by itself - + camera_serial = parse_camera_serial(masked_image_file) - + fovsplitter = FOVMultiWellsSplitter(img_fov, camera_serial=camera_serial, px2um=microns_per_pixel, **fovsplitter_param) - + # initialize timers @@ -325,26 +325,26 @@ def compressVideo(video_file, masked_image_file, mask_param, expected_fps=25, microns_per_pixel = microns_per_pixel, is_light_background = int(mask_param['is_light_background']) ) - mask_dataset, full_dataset, mean_intensity = initMasksGroups(mask_fid, + mask_dataset, full_dataset, mean_intensity = initMasksGroups(mask_fid, expected_frames, vid.height, vid.width, attr_params, save_full_interval) - + if is_bgnd_subtraction: bg_dataset = createImgGroup(mask_fid, "/bgnd", 1, vid.height, vid.width, is_expandable=False) # because we only save the one background: - bg_dataset._v_attrs['save_interval'] = vid.frame_max-vid.first_frame + 1 + bg_dataset._v_attrs['save_interval'] = vid.frame_max-vid.first_frame + 1 bg_dataset[0,:,:] = img_fov - + if vid.dtype != np.uint8: # this will worm as flags to be sure that the normalization took place. - normalization_range = mask_fid.create_earray('/', + normalization_range = mask_fid.create_earray('/', 'normalization_range', atom=tables.Float32Atom(), shape=(0, 2), expectedrows=expected_frames, filters=TABLE_FILTERS ) - + while frame_number < max_frame: ret, image = vid.read() @@ -396,23 +396,23 @@ def compressVideo(video_file, masked_image_file, mask_param, expected_fps=25, Ibuff = Ibuff[:ind_buff + 1] # mask buffer and save data into the hdf5 file - if (ind_buff == buffer_size - 1 or ret == 0) and Ibuff.size > 0: + if (ind_buff == buffer_size - 1 or ret == 0) and Ibuff.size > 0: if is_bgnd_subtraction: Ibuff_b = bgnd_subtractor.apply(Ibuff, frame_number) else: Ibuff_b = Ibuff - + #calculate the max/min in the of the buffer img_reduce = reduceBuffer(Ibuff_b, mask_param['is_light_background']) mask = getROIMask(img_reduce, **mask_param) - + Ibuff *= mask # now apply the well_mask if is MWP if is_fov_tosplit: fovsplitter.apply_wells_mask(Ibuff) # Ibuff will be modified after this - + # add buffer to the hdf5 file frame_first_buff = frame_number - Ibuff.shape[0] mask_dataset.append(Ibuff) @@ -421,20 +421,17 @@ def compressVideo(video_file, masked_image_file, mask_param, expected_fps=25, # calculate the progress and put it in a string progress_str = progressTime.get_str(frame_number) print_flush(base_name + ' ' + progress_str) - + # finish process if ret == 0: break # close the video vid.release() - + # save fovsplitting data if is_fov_tosplit: fovsplitter.write_fov_wells_to_file(masked_image_file) read_and_save_timestamp(masked_image_file) print_flush(base_name + ' Compressed video done.') - - - diff --git a/tierpsy/analysis/feat_tierpsy/get_tierpsy_features.py b/tierpsy/analysis/feat_tierpsy/get_tierpsy_features.py index 9397b04b..0f31fc15 100644 --- a/tierpsy/analysis/feat_tierpsy/get_tierpsy_features.py +++ b/tierpsy/analysis/feat_tierpsy/get_tierpsy_features.py @@ -20,7 +20,7 @@ def save_timeseries_feats_table(features_file, derivate_delta_time, fovsplitter_param={}): timeseries_features = [] fps = read_fps(features_file) - + # initialise class for splitting fov if len(fovsplitter_param) > 0: is_fov_tosplit = True @@ -30,7 +30,7 @@ def save_timeseries_feats_table(features_file, derivate_delta_time, fovsplitter_ is_fov_tosplit = False print('is fov to split?',is_fov_tosplit) - + if is_fov_tosplit: # split fov in wells masked_image_file = features_file.replace('Results','MaskedVideos') @@ -41,14 +41,14 @@ def save_timeseries_feats_table(features_file, derivate_delta_time, fovsplitter_ # well_shape=fovsplitter_param['well_shape']) fovsplitter = FOVMultiWellsSplitter(masked_image_file, **fovsplitter_param) - # store wells data in the features file + # store wells data in the features file fovsplitter.write_fov_wells_to_file(features_file) - - - + + + with pd.HDFStore(features_file, 'r') as fid: trajectories_data = fid['/trajectories_data'] - + trajectories_data_g = trajectories_data.groupby('worm_index_joined') progress_timer = TimeCounter('') base_name = get_base_name(features_file) @@ -61,39 +61,39 @@ def _display_progress(n): dd + ' Total time:' + progress_timer.get_time_str()) - + _display_progress(0) with tables.File(features_file, 'r+') as fid: - + for gg in ['/timeseries_data', '/event_durations', '/timeseries_features']: if gg in fid: fid.remove_node(gg) - - + + feat_dtypes = [(x, np.float32) for x in timeseries_all_columns] - + feat_dtypes = [('worm_index', np.int32), ('timestamp', np.int32), - ('well_name', 'S3')] + feat_dtypes - + ('well_name', 'S3')] + feat_dtypes + timeseries_features = fid.create_table( '/', 'timeseries_data', obj = np.recarray(0, feat_dtypes), filters = TABLE_FILTERS) - + if '/food_cnt_coord' in fid: food_cnt = fid.get_node('/food_cnt_coord')[:] else: food_cnt = None - + #If i find the ventral side in the multiworm case this has to change ventral_side = read_ventral_side(features_file) - + for ind_n, (worm_index, worm_data) in enumerate(trajectories_data_g): skel_id = worm_data['skeleton_id'].values - + #deal with any nan in the skeletons good_id = skel_id>=0 skel_id_val = skel_id[good_id] @@ -101,7 +101,7 @@ def _display_progress(n): args = [] for p in ('skeletons', 'widths', 'dorsal_contours', 'ventral_contours'): - + node_str = '/coordinates/' + p if node_str in fid: node = fid.get_node(node_str) @@ -114,12 +114,12 @@ def _display_progress(n): dat[good_id] = dd else: dat = None - + args.append(dat) timestamp = worm_data['timestamp_raw'].values.astype(np.int32) - - feats = get_timeseries_features(*args, + + feats = get_timeseries_features(*args, timestamp = timestamp, food_cnt = food_cnt, fps = fps, @@ -133,38 +133,38 @@ def _display_progress(n): feats['well_name'] = fovsplitter.find_well_from_trajectories_data(worm_data) else: feats['well_name'] = 'n/a' - # cast well_name to the correct type + # cast well_name to the correct type # (before shuffling columns, so it remains the last entry) # needed because for some reason this does not work: # feats['well_name'] = feats['well_name'].astype('S3') feats['_well_name'] = feats['well_name'].astype('S3') feats.drop(columns='well_name', inplace=True) feats.rename(columns={'_well_name':'well_name'}, inplace=True) - + #move the last fields to the first columns cols = feats.columns.tolist() cols = cols[-2:] + cols[:-2] cols[1], cols[2] = cols[2], cols[1] - + feats = feats[cols] - + feats['worm_index'] = feats['worm_index'].astype(np.int32) feats['timestamp'] = feats['timestamp'].astype(np.int32) feats = feats.to_records(index=False) - + timeseries_features.append(feats) _display_progress(ind_n) - + def save_feats_stats(features_file, derivate_delta_time): with pd.HDFStore(features_file, 'r') as fid: fps = fid.get_storer('/trajectories_data').attrs['fps'] timeseries_data = fid['/timeseries_data'] blob_features = fid['/blob_features'] if '/blob_features' in fid else None - + # do we need split-FOV sumaries? if 'well_name' not in timeseries_data.columns: - # for some weird reason, save_feats_stats is being called on an old + # for some weird reason, save_feats_stats is being called on an old # featuresN file without calling save_timeseries_feats_table first is_fov_tosplit = False else: @@ -176,8 +176,8 @@ def save_feats_stats(features_file, derivate_delta_time): assert all(timeseries_data['well_name']=='n/a'), \ 'Something is wrong with well naming - go check save_feats_stats' is_fov_tosplit = False - - #Now I want to calculate the stats of the video + + #Now I want to calculate the stats of the video if is_fov_tosplit: # get summary stats per well and then concatenate them all well_name_list = list(set(timeseries_data['well_name']) - set(['n/a'])) @@ -187,27 +187,27 @@ def save_feats_stats(features_file, derivate_delta_time): idx = timeseries_data['well_name'] == well # calculate stats per well tmp = get_summary_stats(timeseries_data[idx].reset_index(), - fps, - blob_features[idx].reset_index(), + fps, + blob_features[idx].reset_index(), derivate_delta_time) tmp = pd.DataFrame(zip(tmp.index, tmp), columns=['name','value']) tmp['well_name'] = well exp_feats.append(tmp) - + # now concat all exp_feats = pd.concat(exp_feats, ignore_index=True) - + else: # we don't need to split the FOV - + exp_feats = get_summary_stats(timeseries_data, - fps, - blob_features, + fps, + blob_features, derivate_delta_time) - + # save on disk # now if is_fov_tosplit exp_feats is a dataframe, otherwise a series if len(exp_feats)>0: - + # different syntax according to df or series if is_fov_tosplit: tot = max(len(x) for x in exp_feats['name']) @@ -217,7 +217,7 @@ def save_feats_stats(features_file, derivate_delta_time): tot = max(len(x) for x in exp_feats.index) dtypes = [('name', 'S{}'.format(tot)), ('value', np.float32)] exp_feats_rec = np.array(list(zip(exp_feats.index, exp_feats)), dtype = dtypes) - + # write on hdf5 file with tables.File(features_file, 'r+') as fid: for gg in ['/features_stats']: @@ -227,43 +227,43 @@ def save_feats_stats(features_file, derivate_delta_time): '/', 'features_stats', obj = exp_feats_rec, - filters = TABLE_FILTERS) + filters = TABLE_FILTERS) + - def get_tierpsy_features(features_file, derivate_delta_time = 1/3, fovsplitter_param={}): #I am adding this so if I add the parameters to calculate the features i can pass it to this function save_timeseries_feats_table(features_file, derivate_delta_time, fovsplitter_param) save_feats_stats(features_file, derivate_delta_time) - + if __name__ == '__main__': - + base_file = '/Users/lferiani/Desktop/Data_FOVsplitter/Results/drugexperiment_1hrexposure_set1_20190712_131508.22436248/metadata' - + features_file = base_file + '_featuresN.hdf5' - # restore features after previous step before testing + # restore features after previous step before testing import shutil shutil.copy(features_file.replace('.hdf5','.bk'), features_file) - + fovsplitter_param = {'total_n_wells':96, 'whichsideup':'upright', 'well_shape':'square'} - get_tierpsy_features(features_file, - derivate_delta_time = 1/3, + get_tierpsy_features(features_file, + derivate_delta_time = 1/3, fovsplitter_param=fovsplitter_param) - + # #%% # from tierpsy.features.tierpsy_features.velocities import _h_get_velocity # from tierpsy.features.tierpsy_features.helper import get_delta_in_frames -# +# # fname = '/Users/avelinojaver/Desktop/small_worms/Results/20191121_featuresN.hdf5' -# +# # delta_time = 0.3 -# -# +# +# # with pd.HDFStore(fname, 'r') as fid: # fps = fid.get_storer('/trajectories_data').attrs['fps'] # blob_features = fid['/blob_features'] @@ -273,6 +273,5 @@ def get_tierpsy_features(features_file, derivate_delta_time = 1/3, fovsplitter_p # for ind_n, (worm_index, indexes) in enumerate(trajectories_data_g.items()): # coords = blob_features.loc[indexes, ['coord_x', 'coord_y']].values # velocity = _h_get_velocity(coords, derivate_delta_frames, fps) -# +# # print(velocity.shape) - diff --git a/tierpsy/analysis/food_cnt/getFoodContourMorph.py b/tierpsy/analysis/food_cnt/getFoodContourMorph.py index a25a8996..a11645c3 100644 --- a/tierpsy/analysis/food_cnt/getFoodContourMorph.py +++ b/tierpsy/analysis/food_cnt/getFoodContourMorph.py @@ -21,6 +21,7 @@ from tierpsy.helper.misc import get_base_name + def skeletonize(img): """ OpenCV function to return a skeletonized version of img, a Mat object""" @@ -50,10 +51,10 @@ def get_patch_mask(img, min_area = None, max_area = None, block_size = None): #%% if min_area is None: min_area = max(1, int(min(img.shape)/200))**2 - + if max_area is None: max_area = (img.shape[0]*img.shape[1])/4 - + if block_size is None: block_size = int(min(img.shape)/8) block_size = block_size+1 if block_size%2==0 else block_size @@ -65,14 +66,14 @@ def get_patch_mask(img, min_area = None, max_area = None, block_size = None): cv2.THRESH_BINARY_INV, blockSize=block_size, C=3) - + mask = cv2.morphologyEx(mask_s, cv2.MORPH_CLOSE, disk(1), iterations=1) mask = cv2.erode(mask, disk(1), iterations=1) #kernel = np.array([(-1,-1,-1), (-1, 1, -1), (-1, -1,-1)]) #ss = cv2.morphologyEx(mask, cv2.MORPH_HITMISS, kernel) #mask = mask-ss #mask = cv2.morphologyEx(mask, cv2.MORPH_OPEN, disk(3)) - + mask = cv2.morphologyEx(mask, cv2.MORPH_CLOSE, disk(3), iterations=3) #plt.imshow(mask) #%% @@ -84,7 +85,7 @@ def get_patch_mask(img, min_area = None, max_area = None, block_size = None): contours, hierarchy = cv2.findContours( mask.copy(), cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)[:-2] - + # typically there are more bad contours therefore it is cheaper to draw # only the valid contours mask = np.zeros(img.shape, dtype=img.dtype) @@ -97,51 +98,51 @@ def get_patch_mask(img, min_area = None, max_area = None, block_size = None): area = cv2.contourArea(contour) if (area >= min_area): ''' - If there is a blob with a very high area it is likely - the mask captured the whole food patch. Do not fill the contour of those areas, + If there is a blob with a very high area it is likely + the mask captured the whole food patch. Do not fill the contour of those areas, otherwise the skeletonization produce weird results. ''' if area >= max_area: cv2.drawContours(mask, contours, ii, 1) else: cv2.drawContours(mask, contours, ii, 1, cv2.FILLED) - - + + #%% return mask - + def get_best_circles(mask, resize_factor = 8): ''' Get the best the best fit to a circle using the hough transform. ''' #%% - #resize image to increase speed. I don't want + #resize image to increase speed. I don't want #%% min_size = min(mask.shape) resize_factor = min_size/max(128, min_size/resize_factor) dsize = tuple(int(x/resize_factor) for x in mask.shape[::-1]) - + mask_s = cv2.dilate(mask, disk(resize_factor/2)) mask_s = cv2.resize(mask_s,dsize) #%% r_min = min(mask_s.shape) r_max = max(mask_s.shape) - + #use the first peak of the circle hough transform to initialize the food shape hough_radii = np.arange(r_min/4, r_max/2, 2) hough_res = hough_circle(mask_s, hough_radii) accums, cx, cy, radii = hough_circle_peaks(hough_res, hough_radii, total_num_peaks=9) #%% - + cx, cy, radii = [np.round(x*resize_factor).astype(np.int) for x in (cx, cy, radii)] #%% - + return list(zip(accums, cx, cy, radii)) def mask_to_food_contour(mask, n_bins = 90, frac_lowess=0.05, _is_debug=False): ''' - Estimate the food contour from a binary mask. + Estimate the food contour from a binary mask. 1) Get the best the best fit to a circle using the hough transform. 2) Transform the mask into polar coordinates centered in the fitted circle. 3) Get the closest point to the circle using 2*n_bins angles. @@ -151,70 +152,70 @@ def mask_to_food_contour(mask, n_bins = 90, frac_lowess=0.05, _is_debug=False): #%% h_res = get_best_circles(mask.copy()) _, cy0, cx0, r0 = h_res[0] - - + + (px, py) = np.where(skeletonize(mask)) - - - + + + xx = px-cx0 yy = py-cy0 - + del_r = np.sqrt(xx*xx + yy*yy) - r0 theta = np.arctan2(xx,yy) - + theta_d = np.round((theta/np.pi)*n_bins).astype(np.int) - + #g = {k:[] for k in range(-n_bins, n_bins+1)} g = {k:[] for k in np.unique(theta_d)} for k,dr in zip(theta_d, del_r): g[k].append(dr) - + _get_min = lambda g : min(g, key = lambda x:abs(x)) if g else np.nan tg, min_dr = zip(*[(k, _get_min(g)) for k,g in g.items()]) - + theta_s = np.array(tg)*np.pi/n_bins - + #increase the range for the interpolation theta_i = np.hstack((theta_s-2*np.pi, theta_s, theta_s+2*np.pi)) - + r_s = np.hstack((min_dr, min_dr, min_dr)) - + out = lowess(r_s, theta_i, frac=frac_lowess) - + #win_frac = 1/3 #filt_w = round(n_bins*win_frac) #filt_w = filt_w+1 if filt_w % 2 == 0 else filt_w #r_s = medfilt(r_s, filt_w) f = interp1d(out[:, 0], out[:, 1]) - + theta_new = np.linspace(-np.pi, np.pi, 480) r_new = f(theta_new) + r0 - + circy = r_new*np.cos(theta_new) + cy0 circx = r_new*np.sin(theta_new) + cx0 #%% if _is_debug: from skimage.draw import circle_perimeter import matplotlib.pylab as plt - + plt.figure(figsize=(5,5)) for ii, (acc, cx, cy, cr) in enumerate(h_res[0:1]): #plt.subplot(3,3,ii+1) plt.imshow(mask) cpy, cpx = circle_perimeter(cy, cx, cr) plt.plot(cpx,cpy, '.r') - - + + plt.figure() plt.plot(theta_d, del_r, '.') plt.plot(tg, min_dr) - - + + plt.figure() plt.plot(theta_i, r_s,'.') plt.plot(out[:, 0], out[:, 1], '.') - + #%% return circx, circy, h_res[0] @@ -224,67 +225,67 @@ def get_dark_mask(full_data): if full_data.shape[0] < 2: #nothing to do here returning return np.zeros((full_data.shape[1], full_data.shape[2]), np.uint8) - + #this mask shoulnd't contain many worms img_h = cv2.medianBlur(np.max(full_data, axis=0), 5) #this mask is likely to contain a lot of worms img_l = cv2.medianBlur(np.min(full_data, axis=0), 5) - + #this is the difference (the tagged pixels should be mostly worms) img_del = img_h-img_l th_d = threshold_otsu(img_del) - + #this is the maximum of the minimum pixels of the worms... th = np.max(img_l[img_del>th_d]) #this is what a darkish mask should look like dark_mask = cv2.dilate((img_h 1: - bgnd = [np.max(bgnd_o[i:i+1], axis=0) for i in range(bgnd_o.shape[0]-1)] + bgnd = [np.max(bgnd_o[i:i+1], axis=0) for i in range(bgnd_o.shape[0]-1)] else: bgnd = [np.squeeze(bgnd_o)] - + min_size = min(bgnd[0].shape) resize_factor = min(resizing_size, min_size)/min_size dsize = tuple(int(x*resize_factor) for x in bgnd[0].shape[::-1]) - + bgnd_s = [cv2.resize(x, dsize) for x in bgnd] for b_img in bgnd_s: Y_pred = get_unet_prediction(b_img, model, n_flips=1) - + if _is_debug: import matplotlib.pylab as plt plt.figure() plt.subplot(1,2,1) plt.imshow(b_img, cmap='gray') - plt.subplot(1, 2,2) + plt.subplot(1, 2,2) plt.imshow(Y_pred, interpolation='none') - + original_size = bgnd[0].shape return Y_pred, original_size, bgnd_s @@ -195,12 +195,12 @@ def get_food_prob(mask_file, model, max_bgnd_images = 2, _is_debug = False, resi def get_food_contour_nn(mask_file, model_path, _is_debug=False): ''' Get the food contour using a pretrained u-net model. - This function is faster if a preloaded model is given since it is very slow + This function is faster if a preloaded model is given since it is very slow to load the model and tensorflow. ''' - + model = load_model(model_path) - + food_prob, original_size, bgnd_images = get_food_prob(mask_file, model, _is_debug=_is_debug) #bgnd_images are only used in debug mode #%% @@ -208,7 +208,7 @@ def get_food_contour_nn(mask_file, model_path, _is_debug=False): cnts, _ = cv2.findContours(patch_m, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)[-2:] - + #pick the largest contour if not cnts: @@ -220,11 +220,11 @@ def get_food_contour_nn(mask_file, model_path, _is_debug=False): patch_m = np.zeros(patch_m.shape, np.uint8) patch_m = cv2.drawContours(patch_m, cnts , ind, color=1, thickness=cv2.FILLED) patch_m = cv2.morphologyEx(patch_m, cv2.MORPH_CLOSE, disk(3), iterations=5) - + cnts, _ = cv2.findContours(patch_m, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)[-2:] - + if len(cnts) == 1: cnts = cnts[0] elif len(cnts) > 1: @@ -232,12 +232,12 @@ def get_food_contour_nn(mask_file, model_path, _is_debug=False): cnts = max(cnts, key=cv2.contourArea) else: return np.zeros([]), food_prob, 0. - - + + hull = cv2.convexHull(cnts) hull_area = cv2.contourArea(hull) cnt_solidity = cv2.contourArea(cnts)/hull_area - + food_cnt = np.squeeze(cnts).astype(np.float) # rescale contour to be the same dimension as the original images food_cnt[:,0] *= original_size[0]/food_prob.shape[0] @@ -246,24 +246,24 @@ def get_food_contour_nn(mask_file, model_path, _is_debug=False): if _is_debug: import matplotlib.pylab as plt img = bgnd_images[0] - - + + #np.squeeze(food_cnt) patch_n = np.zeros(img.shape, np.uint8) patch_n = cv2.drawContours(patch_n, [cnts], 0, color=1, thickness=cv2.FILLED) top = img.max() bot = img.min() - img_n = (img-bot)/(top-bot) + img_n = (img-bot)/(top-bot) img_rgb = np.repeat(img_n[..., None], 3, axis=2) #img_rgb = img_rgb.astype(np.uint8) img_rgb[...,0] = ((patch_n==0)*0.5 + 0.5)*img_rgb[...,0] - + plt.figure() plt.imshow(img_rgb) - + plt.plot(hull[:,:,0], hull[:,:,1], 'r') plt.title('solidity = {:.3}'.format(cnt_solidity)) - #%% + #%% return food_cnt, food_prob, cnt_solidity @@ -271,9 +271,6 @@ def get_food_contour_nn(mask_file, model_path, _is_debug=False): if __name__ == '__main__': mask_file = '/Users/ajaver/OneDrive - Imperial College London/optogenetics/Arantza/MaskedVideos/oig8/oig-8_ChR2_control_males_3_Ch1_11052017_161018.hdf5' - + #mask_file = '/Volumes/behavgenom_archive$/Avelino/Worm_Rig_Tests/short_movies_new/MaskedVideos/Double_picking_020317/trp-4_worms6_food1-3_Set4_Pos5_Ch3_02032017_153225.hdf5' food_cnt, food_prob,cnt_solidity = get_food_contour_nn(mask_file, _is_debug=True) - - - \ No newline at end of file diff --git a/tierpsy/analysis/ske_create/getSkeletonsTables.py b/tierpsy/analysis/ske_create/getSkeletonsTables.py index 1017f35f..e2c34123 100755 --- a/tierpsy/analysis/ske_create/getSkeletonsTables.py +++ b/tierpsy/analysis/ske_create/getSkeletonsTables.py @@ -79,7 +79,7 @@ def getWormMask( # compute the thresholded mask worm_mask = worm_img < threshold if is_light_background else worm_img > threshold worm_mask = (worm_mask & (worm_img != 0)).astype(np.uint8) - + # first compute a small closing to join possible fragments of the worm. worm_mask = cv2.morphologyEx(worm_mask, cv2.MORPH_CLOSE, strel_half) @@ -91,7 +91,7 @@ def getWormMask( worm_mask = np.zeros_like(worm_mask) if worm_cnt.size > 0: cv2.drawContours(worm_mask, [worm_cnt.astype(np.int32)], 0, 1, -1) - + # let's do closing with a larger structural element to close any gaps inside the worm. # It is faster to do several iterations rather than use a single larger # strel. @@ -108,7 +108,7 @@ def getWormMask( worm_mask = np.zeros_like(worm_mask) if worm_cnt.size > 0: cv2.drawContours(worm_mask, [worm_cnt.astype(np.int32)], 0, 1, -1) - + return worm_mask, worm_cnt, cnt_area @@ -198,8 +198,8 @@ def _initSkeletonsArrays(ske_file_id, tot_rows, resampling_N, worm_midbody): # this is to initialize the arrays to one row, pytables do not accept empty arrays as initializers of carrays if tot_rows == 0: - tot_rows = 1 - + tot_rows = 1 + #define dimession of each array, it is the only part of the array that varies data_dims = {} for data_str in ['skeleton', 'contour_side1', 'contour_side2']: @@ -208,57 +208,57 @@ def _initSkeletonsArrays(ske_file_id, tot_rows, resampling_N, worm_midbody): data_dims['contour_width'] = (tot_rows, resampling_N) data_dims['width_midbody'] = (tot_rows,) data_dims['contour_area'] = (tot_rows,) - + #create and reference all the arrays def _create_array(field, dims): if '/' + field in ske_file_id: ske_file_id.remove_node('/', field) - - return ske_file_id.create_carray('/', - field, - tables.Float32Atom(dflt=np.nan), - dims, + + return ske_file_id.create_carray('/', + field, + tables.Float32Atom(dflt=np.nan), + dims, filters=TABLE_FILTERS) - + skel_arrays = {field:_create_array(field, dims) for field, dims in data_dims.items()} inram_skel_arrays = {field:np.ones(dims, dtype=np.float32)*np.nan for field, dims in data_dims.items()} - + # flags to mark if a frame was skeletonized traj_dat = ske_file_id.get_node('/trajectories_data') has_skeleton = traj_dat.cols.has_skeleton has_skeleton[:] = np.zeros_like(has_skeleton) #delete previous - + # return skel_arrays, has_skeleton return skel_arrays, has_skeleton, inram_skel_arrays -def trajectories2Skeletons(skeletons_file, +def trajectories2Skeletons(skeletons_file, masked_image_file, - resampling_N=49, - min_blob_area=50, - strel_size=5, + resampling_N=49, + min_blob_area=50, + strel_size=5, worm_midbody=(0.35, 0.65), - analysis_type="WORM", - skel_args = {'num_segments' : 24, + analysis_type="WORM", + skel_args = {'num_segments' : 24, 'head_angle_thresh' : 60} ): - + #get the index number for the width limit midbody_ind = (int(np.floor( worm_midbody[0]*resampling_N)), int(np.ceil(worm_midbody[1]*resampling_N))) - + #read trajectories data with pandas with pd.HDFStore(skeletons_file, 'r') as ske_file_id: trajectories_data = ske_file_id['/trajectories_data'] - + # extract the base name from the masked_image_file. This is used in the # progress status. base_name = masked_image_file.rpartition('.')[0].rpartition(os.sep)[-1] progress_prefix = base_name + ' Calculating skeletons.' - - - + + + # open skeleton file for append and #the compressed videos as read with tables.File(skeletons_file, "r+") as ske_file_id: @@ -271,39 +271,39 @@ def trajectories2Skeletons(skeletons_file, #invert (at least if is_light_background is true) is_light_background = not is_light_background - + #get generators to get the ROI for each frame - ROIs_generator = generateMoviesROI(masked_image_file, - trajectories_data, + ROIs_generator = generateMoviesROI(masked_image_file, + trajectories_data, bgnd_param = bgnd_param, progress_prefix = progress_prefix) # add data from the experiment info (currently only for singleworm) - with tables.File(masked_image_file, "r") as mask_fid: + with tables.File(masked_image_file, "r") as mask_fid: if '/experiment_info' in ske_file_id: ske_file_id.remove_node('/', 'experiment_info') if '/experiment_info' in mask_fid: dd = mask_fid.get_node('/experiment_info').read() ske_file_id.create_array('/', 'experiment_info', obj=dd) - - + + #initialize arrays to save the skeletons data tot_rows = len(trajectories_data) -# skel_arrays, has_skeleton = _initSkeletonsArrays(ske_file_id, - skel_arrays, has_skeleton, inram_skel_arrays = _initSkeletonsArrays(ske_file_id, - tot_rows, - resampling_N, +# skel_arrays, has_skeleton = _initSkeletonsArrays(ske_file_id, + skel_arrays, has_skeleton, inram_skel_arrays = _initSkeletonsArrays(ske_file_id, + tot_rows, + resampling_N, worm_midbody) - + # dictionary to store previous skeletons prev_skeleton = {} - + for worms_in_frame in ROIs_generator: for ind, roi_dat in worms_in_frame.items(): row_data = trajectories_data.loc[ind] worm_img, roi_corner = roi_dat skeleton_id = int(row_data['skeleton_id']) - + # get the previous worm skeletons to orient them worm_index = row_data['worm_index_joined'] if worm_index not in prev_skeleton: @@ -312,30 +312,30 @@ def trajectories2Skeletons(skeletons_file, if analysis_type == "ZEBRAFISH": output = _zebra_func(worm_img, skel_args, resampling_N) else: - _, worm_cnt, _ = getWormMask(worm_img, - row_data['threshold'], + _, worm_cnt, _ = getWormMask(worm_img, + row_data['threshold'], strel_size, - min_blob_area=row_data['area'] / 2, + min_blob_area=row_data['area'] / 2, is_light_background = is_light_background) # get skeletons output = getSkeleton(worm_cnt, prev_skeleton[worm_index], resampling_N, **skel_args) - - - + + + if output is not None and output[0].size > 0: skeleton, ske_len, cnt_side1, cnt_side2, cnt_widths, cnt_area = output prev_skeleton[worm_index] = skeleton.copy() #mark row as a valid skeleton has_skeleton[skeleton_id] = True - + # save segwrom_results # skel_arrays['skeleton_length'][skeleton_id] = ske_len # skel_arrays['contour_width'][skeleton_id, :] = cnt_widths inram_skel_arrays['skeleton_length'][skeleton_id] = ske_len inram_skel_arrays['contour_width'][skeleton_id, :] = cnt_widths - + mid_width = np.median(cnt_widths[midbody_ind[0]:midbody_ind[1]+1]) # skel_arrays['width_midbody'][skeleton_id] = mid_width inram_skel_arrays['width_midbody'][skeleton_id] = mid_width @@ -351,33 +351,33 @@ def trajectories2Skeletons(skeletons_file, inram_skel_arrays['contour_area'][skeleton_id] = cnt_area # import pdb # pdb.set_trace() - + # now write on disk for key in inram_skel_arrays: skel_arrays[key][:] = inram_skel_arrays[key].astype(np.float32) - + if __name__ == '__main__': - + import shutil from tierpsy.helper.params.tracker_param import TrackerParams - + # root_dir = '/Volumes/behavgenom$/Andre/fishVideos/' root_dir = '/Users/lferiani/Desktop/Data_FOVsplitter/short' - + #ff = 'N2_N10_F1-3_Set1_Pos7_Ch1_12112016_024337.hdf5' #ff = 'unc-9_N10_F1-3_Set1_Pos1_Ch5_17112016_193814.hdf5' #ff = 'trp-4_N1_Set3_Pos6_Ch1_19102016_172113.hdf5' #ff = 'trp-4_N10_F1-1_Set1_Pos2_Ch4_02112016_201534.hdf5' #ff = 'f3_ss_uncompressed.hdf5' ff = 'drugexperiment_1hr30minexposure_set1_bluelight_20190722_173404.22436248/metadata.hdf5' - masked_image_file = os.path.join(root_dir, 'MaskedVideos', ff) + masked_image_file = os.path.join(root_dir, 'MaskedVideos', ff) skeletons_file = os.path.join(root_dir, 'Results', ff.replace('.hdf5', '_skeletons.hdf5')) # restore skeletons from backup shutil.copy(skeletons_file.replace('.hdf5','.bk'), skeletons_file) # json_file = os.path.join(root_dir, 'f3_ss_uncompressed.json') json_file = '/Users/lferiani/Desktop/Data_FOVsplitter/loopbio_rig_96WP_upright_Hydra05.json' - + # read parameters params = TrackerParams(json_file) p = params.p_dict @@ -393,6 +393,6 @@ def trajectories2Skeletons(skeletons_file, 'skel_args' : skel_args } - + trajectories2Skeletons(skeletons_file, masked_image_file, **argkws_d) diff --git a/tierpsy/analysis/ske_create/zebrafishAnalysis/zebrafishAnalysis.py b/tierpsy/analysis/ske_create/zebrafishAnalysis/zebrafishAnalysis.py index ab736c5d..88decf0a 100644 --- a/tierpsy/analysis/ske_create/zebrafishAnalysis/zebrafishAnalysis.py +++ b/tierpsy/analysis/ske_create/zebrafishAnalysis/zebrafishAnalysis.py @@ -522,9 +522,8 @@ def getZebrafishMask(frame, config): cnt_area = cv2.contourArea(worm_cnt) output = worm_mask, worm_cnt, cnt_area, \ cleaned_mask, head_point, smoothed_points - + if output is None: return [None]*6 else: return output - diff --git a/tierpsy/analysis/traj_create/getBlobTrajectories.py b/tierpsy/analysis/traj_create/getBlobTrajectories.py index c47152d3..0539dc45 100755 --- a/tierpsy/analysis/traj_create/getBlobTrajectories.py +++ b/tierpsy/analysis/traj_create/getBlobTrajectories.py @@ -22,6 +22,7 @@ from tierpsy.helper.misc import TimeCounter, print_flush, TABLE_FILTERS + def _thresh_bw(pix_valid): # calculate otsu_threshold as lower limit. Otsu understimates the threshold. try: @@ -42,7 +43,7 @@ def _thresh_bw(pix_valid): xx = np.arange(otsu_thresh, cumhist.size) try: - # the threshold is calculated as the first pixel level above the otsu threshold + # the threshold is calculated as the first pixel level above the otsu threshold # at which there would be larger increase in the object area. hist_ratio = pix_hist[xx] / cumhist[xx] thresh = np.where( @@ -70,11 +71,11 @@ def getBufferThresh(ROI_buffer, worm_bw_thresh_factor, is_light_background, anal ''' calculate threshold using the nonzero pixels. Using the buffer instead of a single image, improves the threshold calculation, since better statistics are recovered''' - + if analysis_type == "ZEBRAFISH": # Override threshold thresh = 255 - else: + else: pix_valid = ROI_buffer[ROI_buffer != 0] @@ -93,7 +94,7 @@ def getBufferThresh(ROI_buffer, worm_bw_thresh_factor, is_light_background, anal thresh *= worm_bw_thresh_factor else: thresh = np.nan - + return thresh @@ -103,8 +104,8 @@ def _remove_corner_blobs(ROI_image): # worms ROI_valid = (ROI_image != 0).astype(np.uint8) - ROI_border_ind, _ = cv2.findContours(ROI_valid, - cv2.RETR_EXTERNAL, + ROI_border_ind, _ = cv2.findContours(ROI_valid, + cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)[-2:] @@ -120,7 +121,7 @@ def _remove_corner_blobs(ROI_image): return ROI_image def _get_blob_mask(ROI_image, thresh, thresh_block_size, is_light_background, analysis_type): - # get binary image, + # get binary image, if is_light_background: ## apply a median filter to reduce rough edges / sharpen the boundary btw worm and background ROI_image_th = cv2.medianBlur(ROI_image, 3) @@ -137,7 +138,7 @@ def _get_blob_mask(ROI_image, thresh, thresh_block_size, is_light_background, an # this case applies for example to worms where the whole body is fluorecently labeled ROI_image_th = cv2.medianBlur(ROI_image, 3) ROI_mask = ROI_image_th >= thresh - + ROI_mask &= (ROI_image != 0) ROI_mask = ROI_mask.astype(np.uint8) @@ -145,17 +146,17 @@ def _get_blob_mask(ROI_image, thresh, thresh_block_size, is_light_background, an -def getBlobContours(ROI_image, - thresh, - strel_size=(5, 5), - is_light_background=True, - analysis_type="WORM", +def getBlobContours(ROI_image, + thresh, + strel_size=(5, 5), + is_light_background=True, + analysis_type="WORM", thresh_block_size=15): - + ROI_image = _remove_corner_blobs(ROI_image) ROI_mask, thresh = _get_blob_mask(ROI_image, thresh, thresh_block_size, is_light_background, analysis_type) - + # clean it using morphological closing - make this optional by setting strel_size to 0 if np.all(strel_size): strel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, strel_size) @@ -163,8 +164,8 @@ def getBlobContours(ROI_image, # get worms, assuming each contour in the ROI is a worm - ROI_worms, hierarchy = cv2.findContours(ROI_mask, - cv2.RETR_EXTERNAL, + ROI_worms, hierarchy = cv2.findContours(ROI_mask, + cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)[-2:] @@ -172,9 +173,9 @@ def getBlobContours(ROI_image, def getBlobDimesions(worm_cnt, ROI_bbox): - + area = float(cv2.contourArea(worm_cnt)) - + worm_bbox = cv2.boundingRect(worm_cnt) bounding_box_xmin = ROI_bbox[0] + worm_bbox[0] bounding_box_xmax = bounding_box_xmin + worm_bbox[2] @@ -182,7 +183,7 @@ def getBlobDimesions(worm_cnt, ROI_bbox): bounding_box_ymax = bounding_box_ymin + worm_bbox[3] # save everything into the the proper output format - blob_bbox =(bounding_box_xmin, + blob_bbox =(bounding_box_xmin, bounding_box_xmax, bounding_box_ymin, bounding_box_ymax) @@ -195,96 +196,96 @@ def getBlobDimesions(worm_cnt, ROI_bbox): if W > L: L, W = W, L # switch if width is larger than length - + blob_dims = (CMx, CMy, L, W, angle) return blob_dims, area, blob_bbox - -def generateImages(masked_image_file, - frames=[], + +def generateImages(masked_image_file, + frames=[], bgnd_param = {}, - progress_str='', + progress_str='', progress_refresh_rate_s=20): - + #loop, save data and display progress base_name = Path(masked_image_file).stem progress_str = base_name + progress_str fps = read_fps(masked_image_file, dflt=25) - + progress_refresh_rate = fps*progress_refresh_rate_s - - - + + + with tables.File(masked_image_file, 'r') as mask_fid: mask_dataset = mask_fid.get_node("/mask") - + tot_frames = mask_dataset.shape[0] - progress_time = TimeCounter(progress_str, tot_frames) - - + progress_time = TimeCounter(progress_str, tot_frames) + + if len(bgnd_param) > 0: - - if '/bgnd' in mask_fid: + + if '/bgnd' in mask_fid: bgnd_subtractor = BackgroundSubtractorPrecalculated(masked_image_file, **bgnd_param) else: bgnd_subtractor = BackgroundSubtractorMasked(masked_image_file, **bgnd_param) else: bgnd_subtractor = None - - - + + + if len(frames) == 0: frames = range(mask_dataset.shape[0]) - + for frame_number in frames: if frame_number % progress_refresh_rate == 0: print_flush(progress_time.get_str(frame_number)) - + image = mask_dataset[frame_number] - + if bgnd_subtractor is not None: image = bgnd_subtractor.apply(image, frame_number) - + yield frame_number, image - + print_flush( progress_time.get_str(frame_number)) - - + + def generateROIBuff(masked_image_file, buffer_size, **argkws): img_generator = generateImages(masked_image_file) - + with tables.File(masked_image_file, 'r') as mask_fid: tot_frames, im_h, im_w = mask_fid.get_node("/mask").shape - + for frame_number, image in img_generator: if frame_number % buffer_size == 0: if frame_number + buffer_size > tot_frames: buffer_size = tot_frames-frame_number #change this value, otherwise the buffer will not get full image_buffer = np.zeros((buffer_size, im_h, im_w), np.uint8) - ini_frame = frame_number - - + ini_frame = frame_number + + image_buffer[frame_number-ini_frame] = image - + #compress if it is the last frame in the buffer if (frame_number+1) % buffer_size == 0 or (frame_number+1 == tot_frames): # z projection and select pixels as connected regions that were selected as worms at # least once in the masks main_mask = np.any(image_buffer, axis=0) - + # change from bool to uint since same datatype is required in # opencv main_mask = main_mask.astype(np.uint8) - - #calculate the contours, only keep the external contours (no holes) and - ROI_cnts, _ = cv2.findContours(main_mask, - cv2.RETR_EXTERNAL, + #calculate the contours, only keep the external contours (no holes) and + + ROI_cnts, _ = cv2.findContours(main_mask, + cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)[-2:] - + yield ROI_cnts, image_buffer, ini_frame - + @@ -293,7 +294,7 @@ def _cnt_to_ROIs(ROI_cnt, image_buffer, min_box_width): ROI_bbox = cv2.boundingRect(ROI_cnt) # bounding box too small to be a worm - ROI_bbox[2] and [3] are width and height if ROI_bbox[2] > min_box_width and ROI_bbox[3] > min_box_width: - # select ROI for all buffer slides + # select ROI for all buffer slides ini_x = ROI_bbox[1] fin_x = ini_x + ROI_bbox[3] ini_y = ROI_bbox[0] @@ -305,7 +306,7 @@ def _cnt_to_ROIs(ROI_cnt, image_buffer, min_box_width): return ROI_buffer, ROI_bbox def _cnt_to_props(ROI_worms, current_frame, thresh, min_area, ROI_bbox=(0,0,0,0)): - + props = [] for worm_cnt in ROI_worms: # obtain features for each worm @@ -313,19 +314,19 @@ def _cnt_to_props(ROI_worms, current_frame, thresh, min_area, ROI_bbox=(0,0,0,0) if area >= min_area: # append data to pytables only if the object is larget than min_area row = (-1, -1, current_frame, *blob_dims, area, *blob_bbox, thresh) - + props.append(row) - + return props def getBlobsData(buff_data, blob_params): #I packed input data to be able top to map the function into generateROIBuff ROI_cnts, image_buffer, frame_number = buff_data - + is_light_background, min_area, min_box_width, worm_bw_thresh_factor, \ strel_size, analysis_type, thresh_block_size = blob_params - + blobs_data = [] # examinate each region of interest for ROI_cnt in ROI_cnts: @@ -334,52 +335,52 @@ def getBlobsData(buff_data, blob_params): if ROI_buffer is not None: # calculate threshold thresh_buff = getBufferThresh(ROI_buffer, worm_bw_thresh_factor, is_light_background, analysis_type) - + for buff_ind in range(image_buffer.shape[0]): curr_ROI = ROI_buffer[buff_ind, :, :] - + # get the contour of possible worms - ROI_worms, hierarchy = getBlobContours(curr_ROI, - thresh_buff, - strel_size, + ROI_worms, hierarchy = getBlobContours(curr_ROI, + thresh_buff, + strel_size, is_light_background, - analysis_type, + analysis_type, thresh_block_size) current_frame = frame_number + buff_ind - + # make sure there are no holes in the contours. This shouldn't occur with the flag RETR_EXTERNAL assert all([hierarchy[0][x][3] == -1 for x in range(len(ROI_worms))]) - + blobs_data += _cnt_to_props(ROI_worms, current_frame, thresh_buff, min_area, ROI_bbox) - + return blobs_data def getBlobsSimple(in_data, blob_params): frame_number, image = in_data min_area, worm_bw_thresh_factor, strel_size = blob_params - - + + img_m = cv2.medianBlur(image, 3) - + valid_pix = img_m[img_m>0] if len(valid_pix) == 0: return [] - + th = _thresh_bw(valid_pix)*worm_bw_thresh_factor - + _, bw = cv2.threshold(img_m, th,255,cv2.THRESH_BINARY) if np.all(strel_size): strel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, strel_size) bw = cv2.morphologyEx(bw, cv2.MORPH_CLOSE, strel) cnts, hierarchy = cv2.findContours(bw, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)[-2:] - - + + blobs_data = _cnt_to_props(cnts, frame_number, th, min_area) return blobs_data -def getBlobsTable(masked_image_file, +def getBlobsTable(masked_image_file, trajectories_file, buffer_size = None, min_area=25, @@ -388,14 +389,14 @@ def getBlobsTable(masked_image_file, strel_size=(5,5), analysis_type="WORM", thresh_block_size=15, - n_cores_used = 1, + n_cores_used = 1, bgnd_param = {}): - + #correct strel if it is not a tuple or list if not isinstance(strel_size, (tuple,list)): strel_size = (strel_size, strel_size) assert len(strel_size) == 2 - + #read properties fps_out, _, is_light_background = read_unit_conversions(masked_image_file) expected_fps = fps_out[0] @@ -404,18 +405,18 @@ def getBlobsTable(masked_image_file, if len(bgnd_param) > 0: bgnd_param['is_light_background'] = int(is_light_background) buffer_size = traj_create_defaults(masked_image_file, buffer_size) - + def _ini_plate_worms(traj_fid, masked_image_file): # intialize main table - + int_dtypes = [('worm_index_blob', np.int), ('worm_index_joined', np.int), ('frame_number', np.int)] - dd = ['coord_x', - 'coord_y', - 'box_length', - 'box_width', + dd = ['coord_x', + 'coord_y', + 'box_length', + 'box_width', 'angle', 'area', 'bounding_box_xmin', @@ -423,9 +424,9 @@ def _ini_plate_worms(traj_fid, masked_image_file): 'bounding_box_ymin', 'bounding_box_ymax', 'threshold'] - + float32_dtypes = [(x, np.float32) for x in dd] - + plate_worms_dtype = np.dtype(int_dtypes + float32_dtypes) plate_worms = traj_fid.create_table('/', "plate_worms", @@ -433,19 +434,19 @@ def _ini_plate_worms(traj_fid, masked_image_file): "Worm feature List", filters = TABLE_FILTERS) - - + + #find if it is a mask from fluorescence and save it in the new group plate_worms._v_attrs['is_light_background'] = is_light_background plate_worms._v_attrs['expected_fps'] = expected_fps #make sure it is in a "serializable" format plate_worms._v_attrs['bgnd_param'] = bytes(json.dumps(bgnd_param), 'utf-8') - + read_and_save_timestamp(masked_image_file, trajectories_file) return plate_worms - + progress_str = ' Calculating trajectories.' if len(bgnd_param) == 0: buff_generator = generateROIBuff(masked_image_file, buffer_size, progress_str = progress_str) @@ -456,9 +457,9 @@ def _ini_plate_worms(traj_fid, masked_image_file): strel_size, analysis_type, thresh_block_size) - + f_blob_data = partial(getBlobsData, blob_params = blob_params) - + else: blob_params = (min_area, worm_bw_thresh_factor, @@ -467,18 +468,16 @@ def _ini_plate_worms(traj_fid, masked_image_file): bgnd_param = bgnd_param, progress_str = progress_str) f_blob_data = partial(getBlobsSimple, blob_params = blob_params) - - + + if n_cores_used > 1: p = mp.Pool(n_cores_used) blobs_generator = p.imap(f_blob_data, buff_generator) else: blobs_generator = map(f_blob_data, buff_generator) - + with tables.open_file(trajectories_file, mode='w') as traj_fid: plate_worms = _ini_plate_worms(traj_fid, masked_image_file) for ibuf, blobs_data in enumerate(blobs_generator): if blobs_data: plate_worms.append(blobs_data) - - \ No newline at end of file diff --git a/tierpsy/gui/QtDesignerFiles/Summarizer.ui b/tierpsy/gui/QtDesignerFiles/Summarizer.ui index 4b73b04f..f82b633d 100755 --- a/tierpsy/gui/QtDesignerFiles/Summarizer.ui +++ b/tierpsy/gui/QtDesignerFiles/Summarizer.ui @@ -7,7 +7,7 @@ 0 0 632 - 514 + 537 @@ -204,13 +204,6 @@ edited trajectories? - - - - Select features by keywords : - - - @@ -218,6 +211,9 @@ edited trajectories? + + + @@ -228,12 +224,30 @@ edited trajectories? - - + + + + Select features by keywords : + + + + + + Abbreviate feature names + + + + + + + Dorsal side annotated + + + diff --git a/tierpsy/gui/Summarizer_ui.py b/tierpsy/gui/Summarizer_ui.py index e6824ee2..a188d478 100644 --- a/tierpsy/gui/Summarizer_ui.py +++ b/tierpsy/gui/Summarizer_ui.py @@ -11,7 +11,7 @@ class Ui_Summarizer(object): def setupUi(self, Summarizer): Summarizer.setObjectName("Summarizer") - Summarizer.resize(632, 514) + Summarizer.resize(632, 537) self.centralwidget = QtWidgets.QWidget(Summarizer) self.centralwidget.setObjectName("centralwidget") self.verticalLayout_3 = QtWidgets.QVBoxLayout(self.centralwidget) @@ -127,22 +127,28 @@ def setupUi(self, Summarizer): self.SelectByKeywords.setObjectName("SelectByKeywords") self.gridLayout_4 = QtWidgets.QGridLayout(self.SelectByKeywords) self.gridLayout_4.setObjectName("gridLayout_4") - self.label_10 = QtWidgets.QLabel(self.SelectByKeywords) - self.label_10.setObjectName("label_10") - self.gridLayout_4.addWidget(self.label_10, 0, 0, 1, 2) self.label_9 = QtWidgets.QLabel(self.SelectByKeywords) self.label_9.setObjectName("label_9") self.gridLayout_4.addWidget(self.label_9, 2, 0, 1, 1) + self.p_keywords_include = QtWidgets.QLineEdit(self.SelectByKeywords) + self.p_keywords_include.setObjectName("p_keywords_include") + self.gridLayout_4.addWidget(self.p_keywords_include, 1, 1, 1, 1) self.p_keywords_exclude = QtWidgets.QLineEdit(self.SelectByKeywords) self.p_keywords_exclude.setObjectName("p_keywords_exclude") self.gridLayout_4.addWidget(self.p_keywords_exclude, 2, 1, 1, 1) self.label_8 = QtWidgets.QLabel(self.SelectByKeywords) self.label_8.setObjectName("label_8") self.gridLayout_4.addWidget(self.label_8, 1, 0, 1, 1) - self.p_keywords_include = QtWidgets.QLineEdit(self.SelectByKeywords) - self.p_keywords_include.setObjectName("p_keywords_include") - self.gridLayout_4.addWidget(self.p_keywords_include, 1, 1, 1, 1) + self.label_10 = QtWidgets.QLabel(self.SelectByKeywords) + self.label_10.setObjectName("label_10") + self.gridLayout_4.addWidget(self.label_10, 0, 0, 1, 2) self.verticalLayout_3.addWidget(self.SelectByKeywords) + self.p_abbreviate_features = QtWidgets.QCheckBox(self.centralwidget) + self.p_abbreviate_features.setObjectName("p_abbreviate_features") + self.verticalLayout_3.addWidget(self.p_abbreviate_features) + self.p_dorsal_side_known = QtWidgets.QCheckBox(self.centralwidget) + self.p_dorsal_side_known.setObjectName("p_dorsal_side_known") + self.verticalLayout_3.addWidget(self.p_dorsal_side_known) self.pushButton_start = QtWidgets.QPushButton(self.centralwidget) self.pushButton_start.setObjectName("pushButton_start") self.verticalLayout_3.addWidget(self.pushButton_start) @@ -176,8 +182,10 @@ def retranslateUi(self, Summarizer): self.p_time_units.setItemText(1, _translate("Summarizer", "seconds")) self.p_time_windows.setText(_translate("Summarizer", "0:end")) self.label_6.setText(_translate("Summarizer", "Time windows")) - self.label_10.setText(_translate("Summarizer", "Select features by keywords :")) self.label_9.setText(_translate("Summarizer", "Exclude keywords")) self.label_8.setText(_translate("Summarizer", "Include keywords")) + self.label_10.setText(_translate("Summarizer", "Select features by keywords :")) + self.p_abbreviate_features.setText(_translate("Summarizer", "Abbreviate feature names")) + self.p_dorsal_side_known.setText(_translate("Summarizer", "Dorsal side annotated")) self.pushButton_start.setText(_translate("Summarizer", "START")) diff --git a/tierpsy/helper/misc/misc.py b/tierpsy/helper/misc/misc.py index 8b2eb4a8..86791af9 100755 --- a/tierpsy/helper/misc/misc.py +++ b/tierpsy/helper/misc/misc.py @@ -38,7 +38,7 @@ def get_local_or_sys_path(file_name): FFPROBE_CMD = get_local_or_sys_path('ffprobe.exe') else: FFPROBE_CMD = get_local_or_sys_path('ffprobe') -except FileNotFoundError: +except FileNotFoundError: FFPROBE_CMD = '' warnings.warn('ffprobe do not found. This might cause problems while extracting the raw videos timestamps.') @@ -79,5 +79,3 @@ def read(self): except Empty: line = None return line - - diff --git a/tierpsy/helper/params/docs_summarizer_param.py b/tierpsy/helper/params/docs_summarizer_param.py index db0e18a8..cda021a9 100644 --- a/tierpsy/helper/params/docs_summarizer_param.py +++ b/tierpsy/helper/params/docs_summarizer_param.py @@ -12,18 +12,18 @@ } dflt_args_list = [ - ('root_dir', - '', + ('root_dir', + '', 'Root directory where the features files are located and the results are going to be saved.' ), - ('feature_type', - 'tierpsy', + ('feature_type', + 'tierpsy', ''' Type of feature file to be used. Either the original OpenWorm features or the new Tierpsy Features. ''' ), - ('summary_type', - 'plate', + ('summary_type', + 'plate', ''' Indicates if the summary is going to be done over each individual plate, each individual trajectory or is going to be a data augmentation by randomingly sampling over a subset of the plate trajectories. @@ -37,7 +37,7 @@ '0:end', ''' Define time windows to extract features from the parts of the video included in each window. - Each window must be defined by the start_time and the end_time connected by \':\' (start_time:end_time). + Each window must be defined by the start_time and the end_time connected by \':\' (start_time:end_time). Different windows must be separated by \',\' (start_time_1:end_time_1, start_time_2:end_time_2). A sequence of equally sized windows can be defined using the format \'start_time:end_time:step'\. Attention: the start_time is included in the window, but the end_time is not included. @@ -53,6 +53,18 @@ Get a pre-selected subset of tierpsy features or select features by keywords. ''' ), + ('abbreviate_features', + False, + 'Shorten the feature names so that they are compatible with MATLAB' + ), + ('dorsal_side_known', + False, + ''' + Dorsal side of worm has been annotated on videos. If dorsal side is not + known then only absolute values for d/v signed features are returned. + (nb. d/v features are not included in the 2k, 512, 256, 16 and 8 sets anyway) + ''', + ), ('keywords_include', '', ''' @@ -71,7 +83,7 @@ Number of times each subsampling is going to be repeated (only for plate_augmentation). ''' ), - + ('frac_worms_to_keep', 0.8, 'Fraction of the total number trajectories that is going to be keep for each subsampling.' diff --git a/tierpsy/summary/collect.py b/tierpsy/summary/collect.py index 2b0490ce..36ad577a 100644 --- a/tierpsy/summary/collect.py +++ b/tierpsy/summary/collect.py @@ -9,46 +9,46 @@ import datetime import tables import pandas as pd -import multiprocessing as mp from tierpsy.helper.misc import TimeCounter, print_flush from tierpsy.summary.process_ow import ow_plate_summary, ow_trajectories_summary, ow_plate_summary_augmented from tierpsy.summary.process_tierpsy import tierpsy_plate_summary, tierpsy_trajectories_summary, tierpsy_plate_summary_augmented from tierpsy import AUX_FILES_DIR +from tierpsy.summary.helper import get_featsum_headers,get_fnamesum_headers -feature_files_ext = {'openworm' : ('_features.hdf5', '_feat_manual.hdf5'), +feature_files_ext = {'openworm' : ('_features.hdf5', '_feat_manual.hdf5'), 'tierpsy' : ('_featuresN.hdf5', '_featuresN.hdf5') } FEAT_SET_DIR = os.path.join(AUX_FILES_DIR,'feat_sets') -feature_sets_filenames = {'tierpsy_8' : 'tierpsy_8.csv', 'tierpsy_16' : 'tierpsy_16.csv', - 'tierpsy_256' : 'tierpsy_256.csv', +feature_sets_filenames = {'tierpsy_8' : 'tierpsy_8.csv', 'tierpsy_16' : 'tierpsy_16.csv', + 'tierpsy_256' : 'tierpsy_256.csv', 'tierpsy_2k' : 'top2k_tierpsy_no_blob_no_eigen_only_abs_no_norm.csv' } valid_feature_types = list(feature_files_ext.keys()) valid_summary_types = ['plate', 'trajectory', 'plate_augmented'] valid_time_windows_connector = ':' valid_time_windows_separator = ',' -time_windows_format_explain = 'Each time window must be defined by the start time and the end time connected by \'-\' (start_time-end_time). Different windows must be separated by {}. A sequence of equally sized windows can be defined with the format start_time:end_time:step.'.format(valid_time_windows_separator) +time_windows_format_explain = 'Each time window must be defined by the start time and the end time connected by \'-\' (start_time-end_time). Different windows must be separated by {}. A sequence of equally sized windows can be defined with the format start_time:end_time:step.'.format(valid_time_windows_separator) feat_df_id_cols = ['file_id','well_name'] def check_in_list(x, list_of_x, x_name): if not x in list_of_x: raise ValueError('{} invalid {}. Valid options {}.'.format(x, x_name, list_of_x)) - + def get_summary_func(feature_type, summary_type, time_windows_ints, time_units, is_manual_index, **fold_args): """ Chooses the function used for the extraction of feature summaries based on the input from the GUI """ - if feature_type == 'tierpsy': + if feature_type == 'tierpsy': if summary_type == 'plate': func = partial(tierpsy_plate_summary, time_windows=time_windows_ints, time_units=time_units, is_manual_index=is_manual_index) elif summary_type == 'trajectory': func = partial(tierpsy_trajectories_summary, time_windows=time_windows_ints, time_units=time_units, is_manual_index=is_manual_index) elif summary_type == 'plate_augmented': func = partial(tierpsy_plate_summary_augmented, time_windows=time_windows_ints, time_units=time_units, is_manual_index=is_manual_index, **fold_args) - + elif feature_type == 'openworm': if summary_type == 'plate': func = ow_plate_summary @@ -69,12 +69,12 @@ def time_windows_parser(time_windows): if valid_time_windows_connector not in time_windows: raise ValueError('Invalid format of time windows: '+time_windows_format_explain) return - + # Remove spaces and replace end with -1 windows = time_windows.replace(' ','').replace('end','-1') # Split at ',' to separate time windows, then split each non-empty time window at '-' or ':' windows = [x.split(valid_time_windows_connector) for x in windows.split(valid_time_windows_separator) if x] - + # Convert to integers try: windows = [[int(x) for x in wdw] for wdw in windows] @@ -100,7 +100,7 @@ def time_windows_parser(time_windows): fin_windows.append(window) else: ValueError('Invalid format of time windows: '+time_windows_format_explain) - + return fin_windows def keywords_parser(keywords): @@ -111,7 +111,7 @@ def keywords_parser(keywords): kwrds = keywords.replace(' ','') # Split at ',' to separate time windows, then keep non-empty words kwrds = [x for x in kwrds.split(',') if x] - + if kwrds: return kwrds else: @@ -120,7 +120,7 @@ def keywords_parser(keywords): def feat_set_parser(select_feat): """ EM : gets the full path of the file containing the selected feature set. - """ + """ if select_feat in feature_sets_filenames.keys(): feat_set_file = os.path.join(FEAT_SET_DIR,feature_sets_filenames[select_feat]) selected_feat = pd.read_csv(feat_set_file, header=None, index_col=None) @@ -129,24 +129,19 @@ def feat_set_parser(select_feat): selected_feat = None return selected_feat - -def make_df_filenames(fnames,time_windows_ints,time_units): +def make_df_filenames(fnames,time_windows_ints): """ EM : Create dataframe with filename summaries and time window info for every time window """ dd = tuple(zip(*enumerate(sorted(fnames)))) df_files = [pd.DataFrame({'file_id' : dd[0], 'file_name' : dd[1]}) for x in range(len(time_windows_ints))] - for iwin in range(len(time_windows_ints)): + for iwin in range(len(time_windows_ints)): df_files[iwin]['is_good'] = False - df_files[iwin]['window_id'] = iwin - df_files[iwin]['start_time'] = time_windows_ints[iwin][0] - df_files[iwin]['end_time'] = time_windows_ints[iwin][1] - df_files[iwin]['time_units'] = time_units return df_files def select_features(win_summaries,keywords_in,keywords_ex,selected_feat): id_cols = [col for col in feat_df_id_cols if col in win_summaries.columns] - + if not win_summaries.empty: if selected_feat is not None: win_summaries = win_summaries[id_cols+selected_feat] @@ -157,119 +152,133 @@ def select_features(win_summaries,keywords_in,keywords_ex,selected_feat): filter_col = [x for x in win_summaries.columns if any(key in x for key in keywords_ex)] win_summaries = win_summaries[win_summaries.columns.drop(filter_col)] return win_summaries - -def process_helper(dat_in, summary_func, time_windows_ints, time_units): - ifile, row = dat_in - fname = row['file_name'] - try: - df_list = summary_func(fname, time_windows_ints, time_units) - except (AttributeError, IOError, KeyError, tables.exceptions.HDF5ExtError, tables.exceptions.NoSuchNodeError): - - df_list = [] - return ifile, df_list - -def calculate_summaries(root_dir, - feature_type, - summary_type, - is_manual_index, - time_windows, - time_units, - n_processes = 1, - _is_debug = False, - **fold_args - ): +def shorten_feature_names(feat_summary): + """IB: shortens the feature names so that they are MATLAB compatible. + Does not change the feature names in the featuresN.hdf5. + Input + feat_summary = dataframe of features that are to be exported as a features + and values to be exported to summary file + Output + feat_summary = dataframe with feature names abbreviated + """ + + replace_vel = [(ft, ft.replace('velocity', 'vel')) for ft + in feat_summary.columns] + replace_ang = [(ft[0], ft[1].replace('angular', 'ang')) for ft + in replace_vel] + replace_rel = [(ft[0], ft[1].replace('relative', 'rel')) for ft + in replace_ang] + + renamed_feats = {k: v for (k, v) in replace_rel} + feat_summary.rename(columns=renamed_feats, inplace=True) + + return feat_summary + +def abs_features_only(feat_summary): + """ IB: drops the absolute features + Input: + feat_summary = dataframe of features + + Output: + feat_summary = dataframe with all d/v signed features removed + """ + + absft = [ft for ft in feat_summary.columns if '_abs' in ft] + ventr = [ft.replace('_abs', '') for ft in absft] + + feats_to_drop = list(set(ventr).intersection(feat_summary.columns)) + + if len(feats_to_drop)>0: + feat_summary.drop(columns=feats_to_drop, inplace=True) + return feat_summary + + else: + return feat_summary + + +def calculate_summaries(root_dir, feature_type, summary_type, is_manual_index, time_windows, time_units, + select_feat, keywords_include, keywords_exclude, abbreviate_features, dorsal_side_known, _is_debug = False, **fold_args): """ - Gets input from the GUI, calls the function that chooses the type of summary + Gets input from the GUI, calls the function that chooses the type of summary and runs the summary calculation for each file in the root_dir. """ save_base_name = 'summary_{}_{}'.format(feature_type, summary_type) if is_manual_index: save_base_name += '_manual' save_base_name += '_' + datetime.datetime.now().strftime('%Y%m%d_%H%M%S') - + #check the options are valid check_in_list(feature_type, valid_feature_types, 'feature_type') check_in_list(summary_type, valid_summary_types, 'summary_type') - + # EM : convert time windows to list of integers in frame number units time_windows_ints = time_windows_parser(time_windows) - # EM : get list of keywords to include and to exclude + # EM : get list of keywords to include and to exclude # TODO: catch conflicts keywords_in = keywords_parser(keywords_include) keywords_ex = keywords_parser(keywords_exclude) - + # EM : get full path to feature set file selected_feat = feat_set_parser(select_feat) - + #get summary function # INPUT time windows time units here summary_func = get_summary_func(feature_type, summary_type, time_windows_ints, time_units, is_manual_index, **fold_args) - + #get extension of results file possible_ext = feature_files_ext[feature_type] ext = possible_ext[1] if is_manual_index else possible_ext[0] - - fnames = glob.glob(os.path.join(root_dir, '**', '*' + ext), recursive=True) - - if len(fnames) == 0: + fnames = glob.glob(os.path.join(root_dir, '**', '*' + ext), recursive=True) + if not fnames: print_flush('No valid files found. Nothing to do here.') return None,None - + # EM :Make df_files list with one features_summaries dataframe per time window - df_files = make_df_filenames(fnames,time_windows_ints,time_units) - + df_files = make_df_filenames(fnames,time_windows_ints) + progress_timer = TimeCounter('') def _displayProgress(n): args = (n + 1, len(df_files[0]), progress_timer.get_time_str()) dd = "Extracting features summary. File {} of {} done. Total time: {}".format(*args) print_flush(dd) + _displayProgress(-1) - - # EM :Make all_summaries list with one element per time window. Each element contains + # EM :Make all_summaries list with one element per time window. Each element contains # the extracted feature summaries from all the files for the given time window. - all_summaries = [[] for x in range(len(time_windows_ints))] - - #i need to use partial and redifine this otherwise multiprocessing since it will not be pickable - _process_row = partial(process_helper, - summary_func=summary_func, - time_windows_ints=time_windows_ints, - time_units=time_units) - - - data2process = [x for x in df_files[0].iterrows()] - - - - n_processes = max(n_processes, 1) - if n_processes <= 1: - gen = map(_process_row, data2process) - else: - p = mp.Pool(n_processes) - gen = p.imap(_process_row, data2process) - + for ifile, row in df_files[0].iterrows(): + fname = row['file_name'] - for ii, (ifile, df_list) in enumerate(gen): - #reformat the outputs and remove any failed - for iwin, df in enumerate(df_list): - df.insert(0, 'file_id', ifile) - all_summaries[iwin].append(df) - if not df.empty: - df_files[iwin].loc[ifile, 'is_good'] = True - _displayProgress(ii + 1) - + df_list = summary_func(fname) + for iwin,df in enumerate(df_list): + try: + df.insert(0, 'file_id', ifile) + all_summaries[iwin].append(df) + except (AttributeError, IOError, KeyError, tables.exceptions.HDF5ExtError, tables.exceptions.NoSuchNodeError): + continue + else: + if not df.empty: + df_files[iwin].loc[ifile, 'is_good'] = True + _displayProgress(ifile) - # EM : Concatenate summaries for each window into one dataframe and select features for iwin in range(len(time_windows_ints)): all_summaries[iwin] = pd.concat(all_summaries[iwin], ignore_index=True, sort=False) all_summaries[iwin] = select_features(all_summaries[iwin],keywords_in,keywords_ex,selected_feat) - + + #IB : add in the option to abbreviate features + if abbreviate_features: + all_summaries[iwin] = shorten_feature_names(all_summaries[iwin]) + #IB : add in removal of signed features + + if not dorsal_side_known: + all_summaries[iwin] = abs_features_only(all_summaries[iwin]) + # EM : Save results if select_feat != 'all': win_save_base_name = save_base_name.replace('tierpsy',select_feat+'_tierpsy') @@ -281,27 +290,38 @@ def _displayProgress(n): f1 = os.path.join(root_dir, 'filenames_{}.csv'.format(win_save_base_name)) f2 = os.path.join(root_dir,'features_{}.csv'.format(win_save_base_name)) - - df_files[iwin].to_csv(f1, index=False) - all_summaries[iwin].to_csv(f2, index=False) - + + fnamesum_headers = get_fnamesum_headers( + f2,feature_type,summary_type,iwin,time_windows_ints[iwin], + time_units,len(time_windows_ints),select_feat) + featsum_headers = get_featsum_headers(f1) + + with open(f1,'w') as fid: + fid.write(fnamesum_headers) + df_files[iwin].to_csv(fid, index=False) + with open(f2,'w') as fid: + fid.write(featsum_headers) + all_summaries[iwin].to_csv(fid, index=False) + out = '****************************' out += '\nFINISHED. Created Files:\n-> {}\n-> {}'.format(f1,f2) - + print_flush(out) - - + + return df_files, all_summaries if __name__ == '__main__': - - root_dir = '/Users/em812/Documents/OneDrive - Imperial College London/Eleni/Tierpsy_GUI/test_results_2' + + root_dir = '/Users/ibarlow/Desktop/test_summarizer' is_manual_index = False -# feature_type = 'tierpsy' - feature_type = 'openworm' + feature_type = 'tierpsy' + # feature_type = 'openworm' summary_type = 'plate_augmented' # summary_type = 'plate' #summary_type = 'trajectory' + abbreviate_features = True + dorsal_sign_known = False # Luigi ## root_dir = '/Users/em812/Documents/OneDrive - Imperial College London/Eleni/Tierpsy_GUI/test_results_2' @@ -312,24 +332,23 @@ def _displayProgress(n): # #summary_type = 'plate_augmented' ## summary_type = 'plate' # summary_type = 'trajectory' - + fold_args = dict( - n_folds = 2, + n_folds = 2, frac_worms_to_keep = 0.8, time_sample_seconds = 10*60 ) - - time_windows = '0-end,100000-101000' '0:end:1000' #'0:end' # time_windows = '0:60,480:540' + + time_windows = '0:50'#'-end,100000-101000' '0:end:1000' #'0:end' # time_windows = '0:60,480:540' time_units = 'frame numbers' - select_feat = 'all' #'tierpsy_2k' + select_feat = 'tierpsy_2k' #'all' # #'tierpsy_256'# #'tierpsy_2k' keywords_include = '' keywords_exclude = '' #'curvature,velocity,norm,abs' - - df_files, all_summaries = calculate_summaries(root_dir, feature_type, summary_type, is_manual_index, time_windows, time_units, select_feat, keywords_include, keywords_exclude, **fold_args) - + + df_files, all_summaries = calculate_summaries(root_dir, feature_type, summary_type, is_manual_index, time_windows, time_units, select_feat, keywords_include, keywords_exclude, abbreviate_features, dorsal_sign_known, **fold_args) + # Luigi -# df_files, all_summaries = calculate_summaries(root_dir, feature_type, -# summary_type, is_manual_index, -# time_windows, time_units, +# df_files, all_summaries = calculate_summaries(root_dir, feature_type, +# summary_type, is_manual_index, +# time_windows, time_units, # **fold_args) - diff --git a/tierpsy/summary/helper.py b/tierpsy/summary/helper.py index 8d27319d..33f5cb6f 100644 --- a/tierpsy/summary/helper.py +++ b/tierpsy/summary/helper.py @@ -11,6 +11,7 @@ import random import math import pdb +import tables,json fold_args_dflt = {'n_folds' : 5, 'frac_worms_to_keep' : 0.8, @@ -91,3 +92,42 @@ def augment_data(df, fold_masks.append(good) return fold_masks + +def read_package_version(fname, + provenance_step, + pkg_name): + + fid = tables.open_file(fname, mode='r') + provenance_tracking = fid.get_node('/provenance_tracking/' + provenance_step).read() + provenance_tracking = json.loads(provenance_tracking.decode('utf-8')) + version = provenance_tracking['pkgs_versions'][pkg_name] + + return version + +def get_featsum_headers(fnamesum_fname): + header = ','.join(['# FILENAMES SUMMARY FILE', fnamesum_fname]) + '\n' + return header + +def get_fnamesum_headers(f2,feature_type,summary_type,iwin, + time_window,time_units,n_windows, + select_feat): + from tierpsy import __version__ as version + + if (n_windows==1 and time_window==[0,-1]): + header = '\n'.join([ + ','.join(['# FEATURE SUMMARIES FILE',f2]), + ','.join(['# TIERPSY_VERSION',version]), + ','.join(['# SUMMARY_TYPE','{}_{}'.format(feature_type, summary_type)]), + ','.join(['# SELECTED FEATURES',select_feat]) + ]) + '\n' + else: + header = '\n'.join([ + ','.join(['# FEATURE SUMMARIES FILE',f2]), + ','.join(['# TIERPSY_VERSION',version]), + ','.join(['# SUMMARY_TYPE','{}_{}'.format(feature_type, summary_type)]), + ','.join(['# SELECTED FEATURES',select_feat]), + ','.join(['# TIME WINDOW ID',str(iwin)]), + ','.join(['# TIME WINDOW START END',str(time_window[0]),str(time_window[1])]), + ','.join(['# TIME UNITS',time_units]) + ]) + '\n' + return header diff --git a/tierpsy/summary/process_tierpsy.py b/tierpsy/summary/process_tierpsy.py index b3b1d4ba..ad3caa09 100644 --- a/tierpsy/summary/process_tierpsy.py +++ b/tierpsy/summary/process_tierpsy.py @@ -22,7 +22,7 @@ def time_to_frame_nb(time_windows,time_units,fps,timestamp,fname): """ if timestamp.empty: return - + from copy import deepcopy time_windows_frames = deepcopy(time_windows) if time_units == 'seconds': @@ -31,24 +31,24 @@ def time_to_frame_nb(time_windows,time_units,fps,timestamp,fname): for ilim in range(2): if time_windows_frames[iwin][ilim]!=-1: time_windows_frames[iwin][ilim] = round(time_windows_frames[iwin][ilim]*fps) - + last_frame = timestamp.sort_values().iloc[-1] - for iwin in range(len(time_windows_frames)): + for iwin in range(len(time_windows_frames)): # If a window ends with -1, replace with the frame number of the last frame (or the start frame of the window+1 if window out of bounds) if time_windows_frames[iwin][1]==-1: time_windows_frames[iwin][1] = max(last_frame+1,time_windows_frames[iwin][0]) - + # If a window is out of bounds, print warning if time_windows_frames[iwin][0]>last_frame: print_flush('Warning: The start time of window {}/{} is out of bounds of file \'{}\'.'.format(iwin+1,len(time_windows_frames),fname)) - + return time_windows_frames def no_fps(time_units,fps): if time_units=='seconds' and fps==-1: print_flush( """ - Warning: The time windows were defined in seconds, but fps for file \'{}\' is unknown. + Warning: The time windows were defined in seconds, but fps for file \'{}\' is unknown. Define time windows in frame numbers instead. """.format(fname) ) @@ -64,17 +64,17 @@ def read_data(fname, time_windows, time_units, fps, is_manual_index): blob_features_list: list of blob_features for each time window (length of lists = number of windows) """ # EM: If time_units=seconds and fps is not defined, then return None with warning of no fps. - # Make this check here, to avoid wasting time reading the file + # Make this check here, to avoid wasting time reading the file if no_fps(time_units,fps): return - - with pd.HDFStore(fname, 'r') as fid: + + with pd.HDFStore(fname, 'r') as fid: timeseries_data = fid['/timeseries_data'] blob_features = fid['/blob_features'] if timeseries_data.empty: #no data, nothing to do here return - + if is_manual_index: #keep only data labeled as worm or worm clusters valid_labels = [WLAB[x] for x in ['WORM', 'WORMS']] @@ -82,21 +82,21 @@ def read_data(fname, time_windows, time_units, fps, is_manual_index): if not 'worm_index_manual' in trajectories_data: #no manual index, nothing to do here return - + good = trajectories_data['worm_label'].isin(valid_labels) good = good & (trajectories_data['skeleton_id'] >= 0) skel_id = trajectories_data['skeleton_id'][good] - + timeseries_data = timeseries_data.loc[skel_id] timeseries_data['worm_index'] = trajectories_data['worm_index_manual'][good].values timeseries_data = timeseries_data.reset_index(drop=True) - + blob_features = blob_features.loc[skel_id].reset_index(drop=True) - + # convert time windows to frame numbers for the given file time_windows_frames = time_to_frame_nb(time_windows,time_units,fps,timeseries_data['timestamp'],fname) - - #extract the timeseries_data and blob_features corresponding to each + + #extract the timeseries_data and blob_features corresponding to each #time window and store them in a list (length of lists = number of windows) timeseries_data_list = [] blob_features_list = [] @@ -106,28 +106,28 @@ def read_data(fname, time_windows, time_units, fps, is_manual_index): blob_features_list.append(blob_features.iloc[in_window.values].reset_index(drop=True)) return timeseries_data_list, blob_features_list -#%% +#%% def tierpsy_plate_summary(fname, time_windows, time_units, is_manual_index = False, delta_time = 1/3): """ - Calculate the plate summaries for a given file fname, within a given time window - (units of start time and end time are in frame numbers). + Calculate the plate summaries for a given file fname, within a given time window + (units of start time and end time are in frame numbers). """ fps = read_fps(fname) data_in = read_data(fname, time_windows, time_units, fps, is_manual_index) - - # if manual annotation was chosen and the trajectories_data does not contain + + # if manual annotation was chosen and the trajectories_data does not contain # worm_index_manual, then data_in is None # if time_windows in seconds and fps is not defined (fps=-1), then data_in is None if data_in is None: return [pd.DataFrame() for iwin in range(len(time_windows))] - + timeseries_data, blob_features = data_in - + # was the fov split in wells? only use the first window to detect that, # and to extract the list of well names is_fov_tosplit = was_fov_split(timeseries_data[0]) # is_fov_tosplit = False - + # initialize list of plate summaries for all time windows plate_feats_list = [] for iwin,window in enumerate(time_windows): @@ -135,8 +135,8 @@ def tierpsy_plate_summary(fname, time_windows, time_units, is_manual_index = Fal plate_feats = get_summary_stats(timeseries_data[iwin], fps, blob_features[iwin], delta_time) plate_feats_list.append(pd.DataFrame(plate_feats).T) else: - # get list of well names in this time window - # (maybe some wells looked empty during a whole window, + # get list of well names in this time window + # (maybe some wells looked empty during a whole window, # this prevents errors later on) well_names_list = list(set(timeseries_data[iwin]['well_name']) - set(['n/a'])) # create a list of well-specific, one-line long dataframes @@ -144,13 +144,13 @@ def tierpsy_plate_summary(fname, time_windows, time_units, is_manual_index = Fal for well_name in well_names_list: # find entries in timeseries_data[iwin] belonging to the right well idx_well = timeseries_data[iwin]['well_name'] == well_name - well_feats = get_summary_stats(timeseries_data[iwin][idx_well].reset_index(), - fps, - blob_features[iwin][idx_well].reset_index(), + well_feats = get_summary_stats(timeseries_data[iwin][idx_well].reset_index(), + fps, + blob_features[iwin][idx_well].reset_index(), delta_time) # first prepend the well_name_s to the well_feats series, # then transpose it so it is a single-row dataframe, - # and append it to the well_feats_list + # and append it to the well_feats_list well_name_s = pd.Series({'well_name':well_name}) well_feats_list.append(pd.DataFrame(pd.concat([well_name_s,well_feats])).T) # check: did we find any well? @@ -158,7 +158,7 @@ def tierpsy_plate_summary(fname, time_windows, time_units, is_manual_index = Fal plate_feats_list.append(pd.DataFrame()) else: # now concatenate all the single-row df in well_feats_list in a single df - # and append it to the growing list (1 entry = 1 window) + # and append it to the growing list (1 entry = 1 window) plate_feats = pd.concat(well_feats_list, ignore_index=True, sort=False) plate_feats_list.append(plate_feats) @@ -166,15 +166,15 @@ def tierpsy_plate_summary(fname, time_windows, time_units, is_manual_index = Fal def tierpsy_trajectories_summary(fname, time_windows, time_units, is_manual_index = False, delta_time = 1/3): """ - Calculate the trajectory summaries for a given file fname, within a given time window - (units of start time and end time are in frame numbers). + Calculate the trajectory summaries for a given file fname, within a given time window + (units of start time and end time are in frame numbers). """ fps = read_fps(fname) data_in = read_data(fname, time_windows, time_units, fps, is_manual_index) if data_in is None: return [pd.DataFrame() for iwin in range(len(time_windows))] timeseries_data, blob_features = data_in - + # initialize list of summaries for all time windows all_summaries_list = [] # loop over time windows @@ -187,25 +187,25 @@ def tierpsy_trajectories_summary(fname, time_windows, time_units, is_manual_inde # loop over worm indexes (individual trajectories) for w_ind, w_ts_data in timeseries_data[iwin].groupby('worm_index'): w_blobs = blob_features[iwin].loc[w_ts_data.index] - + w_ts_data = w_ts_data.reset_index(drop=True) w_blobs = w_blobs.reset_index(drop=True) - + worm_feats = get_summary_stats(w_ts_data, fps, w_blobs, delta_time) # returns empty dataframe when w_ts_data is empty worm_feats = pd.DataFrame(worm_feats).T worm_feats = add_trajectory_info(worm_feats, w_ind, w_ts_data, fps) - + all_summary.append(worm_feats) # concatenate all trajectories in given time window into one dataframe all_summary = pd.concat(all_summary, ignore_index=True, sort=False) - + # add dataframe to the list of summaries for all time windows all_summaries_list.append(all_summary) - + return all_summaries_list #%% - + def tierpsy_plate_summary_augmented(fname, time_windows, time_units, is_manual_index = False, delta_time = 1/3, **fold_args): fps = read_fps(fname) data_in = read_data(fname, time_windows, time_units, fps, is_manual_index) @@ -215,7 +215,7 @@ def tierpsy_plate_summary_augmented(fname, time_windows, time_units, is_manual_i # initialize list of summaries for all time windows all_summaries_list = [] - + # loop over time windows for iwin,window in enumerate(time_windows): if timeseries_data[iwin].empty: @@ -226,24 +226,24 @@ def tierpsy_plate_summary_augmented(fname, time_windows, time_units, is_manual_i all_summary = [] # loop over folds for i_fold, ind_fold in enumerate(fold_index): - - + + timeseries_data_r = timeseries_data[iwin][ind_fold].reset_index(drop=True) blob_features_r = blob_features[iwin][ind_fold].reset_index(drop=True) - - + + plate_feats = get_summary_stats(timeseries_data_r, fps, blob_features_r, delta_time) plate_feats = pd.DataFrame(plate_feats).T plate_feats.insert(0, 'i_fold', i_fold) - + all_summary.append(plate_feats) - + # concatenate all folds in given time window into one dataframe all_summary = pd.concat(all_summary, ignore_index=True, sort=False) - + # add dataframe to the list of summaries for all time windows all_summaries_list.append(all_summary) - + return all_summaries_list @@ -254,13 +254,13 @@ def tierpsy_plate_summary_augmented(fname, time_windows, time_units, is_manual_i fname='/Users/lferiani/Desktop/Data_FOVsplitter/evgeny/Results/20190808_subset/evgeny_plate01_r1_20190808_114758.22956805/metadata_featuresN.hdf5' is_manual_index = False - + fold_args = dict( - n_folds = 2, + n_folds = 2, frac_worms_to_keep = 0.8, time_sample_seconds = 10*60 ) - + # time_windows = [[0,10000],[10000,15000],[10000000,-1]] # time_units = 'frameNb' @@ -270,8 +270,3 @@ def tierpsy_plate_summary_augmented(fname, time_windows, time_units, is_manual_i # summary = tierpsy_plate_summary(fname,time_windows,time_units) summary = tierpsy_trajectories_summary(fname,time_windows,time_units) # summary = tierpsy_plate_summary_augmented(fname,time_windows,time_units,is_manual_index=False,delta_time=1/3,**fold_args) - - - - - \ No newline at end of file