From 4bd0edd472aae5c22f28fb7e8f35c6729500a19b Mon Sep 17 00:00:00 2001 From: Constantin Pape Date: Tue, 17 Oct 2023 15:18:12 +0200 Subject: [PATCH] Add changes from Sarah to tutorial code --- .../synth_sumb/tracing/mb_fils_network.py | 10 +- .../synth_sumb/tracing/mb_graph_mp.py | 437 ++++++------------ .../synth_sumb/tracing/mbu_picking.py | 19 +- 3 files changed, 152 insertions(+), 314 deletions(-) diff --git a/code/tutorials/synth_sumb/tracing/mb_fils_network.py b/code/tutorials/synth_sumb/tracing/mb_fils_network.py index 6bc8f86..c879cec 100644 --- a/code/tutorials/synth_sumb/tracing/mb_fils_network.py +++ b/code/tutorials/synth_sumb/tracing/mb_fils_network.py @@ -36,20 +36,20 @@ ####### Input data -ROOT_PATH = '../../../..' # Data path +ROOT_PATH = '/scratch/users/muth9/simsiam/particle_picking' # Data path # Input STAR file with the GraphMCF pickles -in_star = ROOT_PATH + '/data/tutorials/synth_sumb/graphs/test_1_seg_mb_graph.star' # The outuput of mb_graph.py +in_star = ROOT_PATH + '/data/graphs/PN_seg_mb_graph.star' # The outuput of mb_graph.py # Sources slice XML file -in_sources = ROOT_PATH + '/data/tutorials/synth_sumb/fils/in/mb_sources.xml' +in_sources = ROOT_PATH + '/pyseg-container/pyseg_system/data/tutorials/synth_sumb/tracing/fils/in/mb_sources.xml' # Targets slice XML file -in_targets = ROOT_PATH + '/data/tutorials/synth_sumb/fils/in/no_mb_targets.xml' +in_targets = ROOT_PATH + '/pyseg-container/pyseg_system/data/tutorials/synth_sumb/tracing/fils/in/no_mb_targets.xml' ####### Output data -out_dir = ROOT_PATH + '/data/tutorials/synth_sumb/fils/out' # '/fils/cyto' +out_dir = ROOT_PATH + '/data/fils' # '/fils/cyto' out_int_g = True ####### Filament geometrical parameters diff --git a/code/tutorials/synth_sumb/tracing/mb_graph_mp.py b/code/tutorials/synth_sumb/tracing/mb_graph_mp.py index 98838a8..efaf554 100644 --- a/code/tutorials/synth_sumb/tracing/mb_graph_mp.py +++ b/code/tutorials/synth_sumb/tracing/mb_graph_mp.py @@ -24,7 +24,6 @@ import os import numpy as np import multiprocessing as mp -import shutil try: import pickle as pickle @@ -45,23 +44,23 @@ ####### Input data -###ROOT_PATH = '/scratch1/projects/rubsak/tapu/Pyseg/pyseg-test/test-20823/' # Data path +ROOT_PATH = '/scratch/users/muth9/simsiam/particle_picking/data' # Data path -### Input STAR file with segmentations -##in_star = ROOT_PATH + 'segs/test_1_seg.star' -npr = 1 # number of parallel processes +# Input STAR file with segmentations +in_star = ROOT_PATH + '/test_time_limit_arsens_vol.star' +npr = 5 # number of parallel processes ####### Output data -output_dir = 'graphs' +output_dir = ROOT_PATH + '/graphs' ####### GraphMCF perameter -###res = 0.756 # nm/pix -s_sig = 1 # 1.5 +res = 1.048 # nm/pix +s_sig = 1.0 # 1.5 v_den = 0.0035 # 0.007 # 0.0025 # nm^3 ve_ratio = 4 # 2 -max_len = 20 # 15 # 30 # nm +max_len = 10 # 15 # 30 # nm ####### Advanced parameters @@ -83,93 +82,89 @@ e_mode = 'low' prop_topo = ps.globals.STR_FIELD_VALUE # ps.globals.STR_FIELD_VALUE_EQ # None is ps.globals.STR_FIELD_VALUE -def parse_command_line(): - """ - Parse the command line. Adapted from sxmask.py - - Arguments: - None: - - Returns: - Parsed arguments object - """ - - # Formatter displays default value - parser= argparse.ArgumentParser( - formatter_class=argparse.ArgumentDefaultsHelpFormatter, - ) - - required= parser.add_argument_group( - title="Required parameters", - description="These parameters are required.") - - required.add_argument('--inStar', help='Input star file.', required=True) - required.add_argument('--pixelSize', type=float, help='Resolution in nm/pix.', required=True) - - options= parser.add_argument_group( - title="Optional parameters", - description="These parameters can be overridden.") - - options.add_argument('--outDir', default=output_dir, help='Output files directory.') - options.add_argument('--nprocs', '-j', default=npr, type=int, help='Number of processors.') - options.add_argument('--sSig', type=float, default=s_sig, help='Sigma for Gaussian filtering.') - options.add_argument('--vDen', type=float, default=v_den, help='Vertex density within membranes in cubic nm.') - options.add_argument('--veRatio', type=float, default=ve_ratio, - help='Averaged ratio vertex/edge in the graph within membrane.') - options.add_argument('--maxLen', type=float, default=max_len, - help='Maximum euclidean shortest distance to membrane in nm.') - - advanced= parser.add_argument_group( - title="Advanced parameters", - description="See the PySeg documentation for a more detailed description of these parameters.") - - advanced.add_argument('--s_sig', default=s_sig, type=float, help='Sigma for Gaussian pre-processing.') - advanced.add_argument('--v_den', default=v_den, type=float, help='Target vertex density (membrane).') - advanced.add_argument('--ve_ratio', default=ve_ratio, type=float, help='Maximum ratio between graph vertices and edges.') - advanced.add_argument('--max_len', default=max_len, type=float, help='Maximum distance to the segmented membrane to process.') - advanced.add_argument('--csig', default=csig, type=float, help='DisPerSe persistence threshold.') - advanced.add_argument('--ang_rot', default=ang_rot, type=float, help='Missing wedge rotation angle.') - advanced.add_argument('--ang_tilt', default=ang_tilt, type=float, help='Missing wedge tilt angle.') - advanced.add_argument('--nstd', default=nstd, type=float, help='Sigma for contrast enhancement.') - advanced.add_argument('--smooth', default=smooth, type=float, help='Skeleton smoothing factor.') - advanced.add_argument('--mb_dst_off', default=mb_dst_off, type=float, help='Mask offset.') - advanced.add_argument('--DILATE_NITER', default=DILATE_NITER, type=float, help='Number of mask dilations.') - advanced.add_argument('--do_clahe', action="store_true", help='Flag to compute CLAHE.') - advanced.add_argument('--v_prop', default=v_prop, type=str, help='Vertex simplification property.') - advanced.add_argument('--e_prop', default=e_prop, type=str, help='Edge simplification property.') - advanced.add_argument('--v_mode', default=v_mode, type=str, help="If 'high' then vertices with highest property values are preserved.") - advanced.add_argument('--e_mode', default=e_mode, type=str, help="If 'high' then edges with highest property values are preserved.") - advanced.add_argument('--prop_topo', default=prop_topo, type=str, help='GraphMCF property.') - - return parser.parse_args() +######################################################################################## +# MAIN ROUTINE +######################################################################################## + +# Get them from the command line if they were passed through itparser = argparse.ArgumentParser() +parser = argparse.ArgumentParser() +parser.add_argument('--inStar', default=in_star, help='Input star file.') +parser.add_argument('--outDir', default=output_dir, help='Output files directory.') +parser.add_argument('--pixelSize', type=float, default=res, help='Resolution in nm/pix.') +parser.add_argument('--sSig', type=float, default=s_sig, help='Sigma for gaussian filtering.') +parser.add_argument('--vDen', type=float, default=v_den, help='Vertex density within membranes in cubic nm.') +parser.add_argument('--veRatio', type=float, default=ve_ratio, + help='Averaged ratio vertex/edge in the graph within membrane.') +parser.add_argument('--maxLen', type=float, default=max_len, + help='Maximum euclidean shortest distance to membrane in nm.') +parser.add_argument('-j', default=npr, type=int, help='Number of processors.') + +args = parser.parse_args() +in_star = args.inStar +output_dir = args.outDir +res = args.pixelSize +s_sig = args.sSig +v_den = args.vDen +ve_ratio = args.veRatio +max_len = args.maxLen +npr = args.j + +# Print initial message +print('Extracting GraphMCF and NetFilament objects from tomograms') +print('\tAuthor: ' + __author__) +print('\tDate: ' + time.strftime("%c") + '\n') +print('Options:') +# print '\tDisPerSe persistence threshold (nsig): ' + str(nsig) +print('\tSTAR file with the segmentations: ' + str(in_star)) +print('\tNumber of parallel processes: ' + str(npr)) +print('\tDisPerSe persistence threshold (csig): ' + str(csig)) +if ang_rot is not None: + print('Missing wedge edge compensation (rot, tilt): (' + str(ang_rot) + ', ' + str(ang_tilt) + ')') +print('\tSigma for gaussian pre-processing: ' + str(s_sig)) +print('\tSigma for contrast enhancement: ' + str(nstd)) +print('\tSkeleton smoothing factor: ' + str(smooth)) +print('\tData resolution: ' + str(res) + ' nm/pixel') +print('\tMask offset: ' + str(mb_dst_off) + ' nm') +print('\tOutput directory: ' + output_dir) +print('Graph density thresholds:') +if v_prop is None: + print('\tTarget vertex density (membrane) ' + str(v_den) + ' vertex/nm^3 for topological simplification') +else: + print('\tTarget vertex density (membrane) ' + + str(v_den) + ' vertex/nm^3 for property ' + v_prop + ' with mode ' + v_mode) +print('\tTarget edge/vertex ratio (non membrane) ' + str(ve_ratio) + ' for property ' + e_prop + ' with mode ' + e_mode) +if do_clahe: + print('\t-Computing CLAHE.') +print('') + +print('Paring input star file...') +star = ps.sub.Star() +star.load(in_star) +in_seg_l = star.get_column_data('_psSegImage') +in_tomo_l = star.get_column_data('_rlnImageName') +star.add_column('_psGhMCFPickle') ### Parallel worker -def pr_worker(pr_id, ids, q_pkls=None): +def pr_worker(pr_id, ids, q_pkls): + pkls_dic = dict() - - # Loop through IDs for row in ids: input_seg, input_tomo = in_seg_l[row], in_tomo_l[row] - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Sub-volume to process found: {input_tomo}', file=sys.stdout, flush=True) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Computing paths for: {input_tomo}', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Sub-volume to process found: ' + input_tomo) + print('\tP[' + str(pr_id) + '] Computing paths for ' + input_tomo + ' ...') path, stem_tomo = os.path.split(input_tomo) stem_pkl, _ = os.path.splitext(stem_tomo) - input_file = os.path.join( - output_dir, - f"{stem_pkl}_g{str(s_sig)}.fits" - ) + input_file = output_dir + '/' + stem_pkl + '_g' + str(s_sig) + '.fits' _, stem = os.path.split(input_file) stem, _ = os.path.splitext(stem) - ###print('\tP[' + str(pr_id) + '] Loading input data: ' + stem_tomo, file=sys.stdout, flush=True) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Loading input data: {stem_tomo}', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Loading input data: ' + stem_tomo) tomo = ps.disperse_io.load_tomo(input_tomo).astype(np.float32) segh = ps.disperse_io.load_tomo(input_seg) - ###print('{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[' + str(pr_id) + '] Computing distance, mask and segmentation tomograms...', file=sys.stdout, flush=True) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Computing distance, mask and segmentation tomograms...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Computing distance, mask and segmentation tomograms...') tomod = ps.disperse_io.seg_dist_trans(segh == SEG_MB) * res maskh = np.ones(shape=segh.shape, dtype=int) maskh[DILATE_NITER:-DILATE_NITER, DILATE_NITER:-DILATE_NITER, DILATE_NITER:-DILATE_NITER] = 0 @@ -177,11 +172,7 @@ def pr_worker(pr_id, ids, q_pkls=None): maskh += mask maskh += (segh == 0) mask = np.asarray(maskh > 0, dtype=float) - - input_msk = os.path.join( - output_dir, - stem + '_mask.fits' - ) + input_msk = output_dir + '/' + stem + '_mask.fits' ps.disperse_io.save_numpy(mask.transpose(), input_msk) if mb_dst_off > 0: seg = np.zeros(shape=segh.shape, dtype=segh.dtype) @@ -193,37 +184,22 @@ def pr_worker(pr_id, ids, q_pkls=None): seg = segh mask_den = np.asarray(tomod <= mb_dst_off, dtype=bool) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Smoothing input tomogram (s={str(s_sig)})...', file=sys.stdout, flush=True) - density = sp.ndimage.gaussian_filter(tomo, s_sig) # WAS sp.ndimage.filters.gaussian_filter(tomo, s_sig) + print('\tP[' + str(pr_id) + '] Smoothing input tomogram (s=' + str(s_sig) + ')...') + density = sp.ndimage.filters.gaussian_filter(tomo, s_sig) density = ps.globals.cont_en_std(density, nstd=nstd, lb=0, ub=1) - ps.disperse_io.save_numpy( - tomo, - os.path.join( - output_dir, - stem_pkl + '.vti' - ) - ) + ps.disperse_io.save_numpy(tomo, output_dir + '/' + stem_pkl + '.vti') ps.disperse_io.save_numpy(density.transpose(), input_file) - ps.disperse_io.save_numpy( - density, - os.path.join( - output_dir, - stem + '.vti' - ) - ) - - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Initializing DisPerSe...', file=sys.stdout, flush=True) - work_dir = os.path.join( - output_dir, - 'disperse_pr_' + str(pr_id) - ) - disperse = ps.disperse_io.DisPerSe(input_file, work_dir) # class defined in handler.py + ps.disperse_io.save_numpy(density, output_dir + '/' + stem + '.vti') + + print('\tP[' + str(pr_id) + '] Initializing DisPerSe...') + work_dir = output_dir + '/disperse_pr_' + str(pr_id) + disperse = ps.disperse_io.DisPerSe(input_file, work_dir) try: disperse.clean_work_dir() # except ps.pexceptions.PySegInputWarning as e: # print e.get_message() except Warning: - print('Jol!!!', file=sys.stdout, flush=True) + print('Jol!!!') # Manifolds for descending fields with the inverted image disperse.set_manifolds('J0a') @@ -231,85 +207,83 @@ def pr_worker(pr_id, ids, q_pkls=None): disperse.set_dump_arcs(-1) # disperse.set_nsig_cut(nsig) rcut = round(density[mask_den].std() * csig, 4) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Persistence cut thereshold set to: {str(rcut)} grey level', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Persistence cut thereshold set to: ' + str(rcut) + ' grey level') disperse.set_cut(rcut) disperse.set_mask(input_msk) disperse.set_smooth(smooth) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Running DisPerSe...', file=sys.stdout, flush=True) - return_code= disperse.mse(no_cut=False, inv=False) - + + print('\tP[' + str(pr_id) + '] Running DisPerSe...') + disperse.mse(no_cut=False, inv=False) skel = disperse.get_skel() - manifolds, return_code = disperse.get_manifolds(no_cut=False, inv=False) - ##print(f"243 manifolds : {type(manifolds)}", file=sys.stdout, flush=True) - ##print(f"244 return_code : {return_code}", file=sys.stdout, flush=True) + manifolds = disperse.get_manifolds(no_cut=False, inv=False) # Build the GraphMCF for the membrane - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Building MCF graph...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Building MCF graph...') graph = ps.mb.MbGraphMCF(skel, manifolds, density, seg) graph.set_resolution(res) graph.build_from_skel(basic_props=False) graph.filter_self_edges() graph.filter_repeated_edges() - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Filtering nodes close to mask border...', file=sys.stdout, flush=True) - mask = sp.ndimage.binary_dilation(mask, iterations=DILATE_NITER) # WAS sp.ndimage.morphology.binary_dilation(mask, iterations=DILATE_NITER) + print('\tP[' + str(pr_id) + '] Filtering nodes close to mask border...') + mask = sp.ndimage.morphology.binary_dilation(mask, iterations=DILATE_NITER) for v in graph.get_vertices_list(): x, y, z = graph.get_vertex_coords(v) if mask[int(round(x)), int(round(y)), int(round(z))]: graph.remove_vertex(v) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Building geometry...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Building geometry...') graph.build_vertex_geometry() if do_clahe: - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] CLAHE on filed_value_inv property...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] CLAHE on filed_value_inv property...') graph.compute_edges_length(ps.globals.SGT_EDGE_LENGTH, 1, 1, 1, False) graph.clahe_field_value(max_geo_dist=50, N=256, clip_f=100., s_max=4.) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Computing vertices and edges properties...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Computing vertices and edges properties...') graph.compute_vertices_dst() graph.compute_edge_filamentness() graph.add_prop_inv(prop_topo, edg=True) graph.compute_edge_affinity() - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Applying general thresholds...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Applying general thresholds...') if ang_rot is not None: - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Deleting edges in MW area...', file=sys.stdout, flush=True) + print('\tDeleting edges in MW area...') graph.filter_mw_edges(ang_rot, ang_tilt) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Computing graph global statistics (before simplification)...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Computing graph global statistics (before simplification)...') nvv, nev, nepv = graph.compute_global_stat(mask=mask_den) - print(f'\t{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[' + str(pr_id) + '] Vertex density: ' + str(round(nvv, 5)) + ' nm^3', file=sys.stdout, flush=True) - print(f'\t{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[' + str(pr_id) + '] Edge density: ' + str(round(nev, 5)) + ' nm^3', file=sys.stdout, flush=True) - print(f'\t{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[' + str(pr_id) + '] Edge/Vertex ratio: ' + str(round(nepv, 5)), file=sys.stdout, flush=True) + print('\t\t-P[' + str(pr_id) + '] Vertex density: ' + str(round(nvv, 5)) + ' nm^3') + print('\t\t-P[' + str(pr_id) + '] Edge density: ' + str(round(nev, 5)) + ' nm^3') + print('\t\t-P[' + str(pr_id) + '] Edge/Vertex ratio: ' + str(round(nepv, 5))) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Graph density simplification for vertices...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Graph density simplification for vertices...') if prop_topo != ps.globals.STR_FIELD_VALUE: print('\t\tProperty used: ' + prop_topo) graph.set_pair_prop(prop_topo) try: graph.graph_density_simp_ref(mask=np.asarray(mask_den, dtype=int), v_den=v_den, - v_prop=v_prop, v_mode=v_mode) + v_prop=v_prop, v_mode=v_mode) except ps.pexceptions.PySegInputWarning as e: - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] WARNING: graph density simplification failed:', file=sys.stdout, flush=True) + print('P[' + str(pr_id) + '] WARNING: graph density simplification failed:') print('\t-' + e.get_message()) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Graph density simplification for edges in the membrane...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Graph density simplification for edges in the membrane...') mask_mb = (seg == 1) * (~mask) nvv, nev, nepv = graph.compute_global_stat(mask=mask_mb) if nepv > ve_ratio: e_den = nvv * ve_ratio hold_e_prop = e_prop graph.graph_density_simp_ref(mask=np.asarray(mask_mb, dtype=int), e_den=e_den, - e_prop=hold_e_prop, e_mode=e_mode, fit=True) + e_prop=hold_e_prop, e_mode=e_mode, fit=True) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Computing graph global statistics (after simplification)...', file=sys.stdout, flush=True) + print('\tP[' + str(pr_id) + '] Computing graph global statistics (after simplification)...') nvv, _, _ = graph.compute_global_stat(mask=mask_den) _, nev_mb, nepv_mb = graph.compute_global_stat(mask=mask_mb) - print(f'\t\tProcess[' + str(pr_id) + '] Vertex density (membrane): ' + str(round(nvv, 5)) + ' nm^3', file=sys.stdout, flush=True) - print(f'\t\tProcess[' + str(pr_id) + '] Edge density (membrane):' + str(round(nev_mb, 5)) + ' nm^3', file=sys.stdout, flush=True) - print(f'\t\tProcess[' + str(pr_id) + '] Edge/Vertex ratio (membrane): ' + str(round(nepv_mb, 5)), file=sys.stdout, flush=True) + print('\t\t-P[' + str(pr_id) + '] Vertex density (membrane): ' + str(round(nvv, 5)) + ' nm^3') + print('\t\t-P[' + str(pr_id) + '] Edge density (membrane):' + str(round(nev_mb, 5)) + ' nm^3') + print('\t\t-P[' + str(pr_id) + '] Edge/Vertex ratio (membrane): ' + str(round(nepv_mb, 5))) - print('\nComputing graph properties (2)...', file=sys.stdout, flush=True) + print('\tComputing graph properties (2)...') graph.compute_mb_geo(update=True) graph.compute_mb_eu_dst() graph.compute_edge_curvatures() @@ -318,195 +292,48 @@ def pr_worker(pr_id, ids, q_pkls=None): graph.compute_edge_filamentness() graph.compute_edge_affinity() - print('\n{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Saving intermediate graphs...', file=sys.stdout, flush=True) - ps.disperse_io.save_vtp( - graph.get_vtp(av_mode=True, edges=True), - os.path.join( - output_dir, - stem + '_edges.vtp' - ) - ) - ps.disperse_io.save_vtp( - graph.get_vtp(av_mode=False, edges=True), - os.path.join( - output_dir, - stem + '_edges_2.vtp' - ) - ) + print('\tSaving intermediate graphs...') + ps.disperse_io.save_vtp(graph.get_vtp(av_mode=True, edges=True), + output_dir + '/' + stem + '_edges.vtp') + ps.disperse_io.save_vtp(graph.get_vtp(av_mode=False, edges=True), + output_dir + '/' + stem + '_edges_2.vtp') # ps.disperse_io.save_vtp(graph.get_scheme_vtp(nodes=True, edges=True), # output_dir + '/' + stem + '_sch.vtp') - out_pkl = os.path.join( - output_dir, - stem_pkl + '.pkl' - ) - print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Process[{str(pr_id)}] Pickling the graph as: ' + out_pkl, file=sys.stdout, flush=True) + out_pkl = output_dir + '/' + stem_pkl + '.pkl' + print('\tP[' + str(pr_id) + '] Pickling the graph as: ' + out_pkl) graph.pickle(out_pkl) # star.set_element('_psGhMCFPickle', row, out_pkl) q_pkls.put((row, out_pkl)) pkls_dic[row] = out_pkl - # End ID loop - - ###sys.exit(pr_id) - ##sys.exit(5) - if q_pkls : sys.exit(return_code) - -###def run(): -######################################################################################## -# MAIN ROUTINE -######################################################################################## - -args = parse_command_line() - -##print(args) -##exit() - -in_star = args.inStar -output_dir = args.outDir -res = args.pixelSize -s_sig = args.sSig -v_den = args.vDen -ve_ratio = args.veRatio -max_len = args.maxLen -npr = args.nprocs - -# Print initial message -print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )}') -print('Extracting GraphMCF and NetFilament objects from tomograms', file=sys.stdout, flush=True) -print('\tAuthor: ' + __author__, file=sys.stdout, flush=True) -print('\tDate: ' + time.strftime("%c") + '\n', file=sys.stdout, flush=True) -print('Options:', file=sys.stdout, flush=True) -# print '\tDisPerSe persistence threshold (nsig): ' + str(nsig) -print('\tSTAR file with the segmentations: ' + str(in_star), file=sys.stdout, flush=True) -print('\tNumber of parallel processes: ' + str(npr), file=sys.stdout, flush=True) -print('\tDisPerSe persistence threshold (csig): ' + str(csig), file=sys.stdout, flush=True) -if ang_rot is not None: - print('Missing wedge edge compensation (rot, tilt): (' + str(ang_rot) + ', ' + str(ang_tilt) + ')', file=sys.stdout, flush=True) -print('\tSigma for Gaussian pre-processing: ' + str(s_sig), file=sys.stdout, flush=True) -print('\tSigma for contrast enhancement: ' + str(nstd), file=sys.stdout, flush=True) -print('\tSkeleton smoothing factor: ' + str(smooth), file=sys.stdout, flush=True) -print('\tData resolution: ' + str(res) + ' nm/pixel', file=sys.stdout, flush=True) -print('\tMask offset: ' + str(mb_dst_off) + ' nm', file=sys.stdout, flush=True) -print('\tOutput directory: ' + output_dir, file=sys.stdout, flush=True) - -print('\nGraph density thresholds:', file=sys.stdout, flush=True) -if v_prop is None: - print('\tTarget vertex density (membrane) ' + str(v_den) + ' vertex/nm^3 for topological simplification', file=sys.stdout, flush=True) -else: - print('\tTarget vertex density (membrane) ' + - str(v_den) + ' vertex/nm^3 for property ' + v_prop + ' with mode ' + v_mode, file=sys.stdout, flush=True) -print('\tTarget edge/vertex ratio (non membrane) ' + str(ve_ratio) + ' for property ' + e_prop + ' with mode ' + e_mode, file=sys.stdout, flush=True) -if do_clahe: - print('\t-Computing CLAHE.', file=sys.stdout, flush=True) -print('', file=sys.stdout, flush=True) -print('Paring input star file...', file=sys.stdout, flush=True) -star = ps.sub.Star() -star.load(in_star) -in_seg_l = star.get_column_data('_psSegImage') -in_tomo_l = star.get_column_data('_rlnImageName') -star.add_column('_psGhMCFPickle') + sys.exit(pr_id) # Loop for processing the input data -print('Running main loop in parallel: \n', file=sys.stdout, flush=True) +print('Running main loop in parallel: ') q_pkls = mp.Queue() processes, pr_results = dict(), dict() spl_ids = np.array_split(list(range(star.get_nrows())), npr) - - -# Sanity checks -if not shutil.which(ps.globals.variables.DPS_MSE_CMD): - print("ERROR!! DisPerSE executable 'mse' not found! Exiting...", file=sys.stdout, flush=True) - exit() -if not os.path.isdir(output_dir) : os.makedirs(output_dir) # os.mkdir() can only operate one directory deep - - for pr_id in range(npr): - print(f'\tProcess[{str(pr_id)}] Starting...') pr = mp.Process(target=pr_worker, args=(pr_id, spl_ids[pr_id], q_pkls)) pr.start() processes[pr_id] = pr - ###print('process ' + str(pr_id) + ', Exit code: [' + str(pr.exitcode) + ']') +for pr_id, pr in zip(iter(processes.keys()), iter(processes.values())): pr.join() - #if pr_id != pr.exitcode: - #print('ERROR: the process ' + str(pr_id) + ' ended unsuccessfully [' + str(pr.exitcode) + ']') - #print('Unsuccessfully terminated. (' + time.strftime("%c") + ')') - #exit() - #else: - #pr.join() - #num_success += 1 - -#if num_success == 0 : - #print(f"333 num_success '{num_success}'") - #exit(3) - -num_success = 0 - -for pr_id in range(npr): - pr = processes[pr_id] - - ###if pr_id != pr.exitcode: - if pr.exitcode != 0: - print(f'\t\tERROR!! Process [{str(pr_id)}] ended unsuccessfully! Exit code: {str(pr.exitcode)}') - print(f'\tProcess[{str(pr_id)}] FAILED. ({time.strftime("%c")})') - ###exit() - else: - print(f'\tProcess[{str(pr_id)}] successful, exit code: {str(pr.exitcode)}') - num_success += 1 - -###exit(4) - -###if num_success == 0 : -print(f"333 num_success '{num_success}'") -###exit(3) - + if pr_id != pr.exitcode: + print('ERROR: the process ' + str(pr_id) + ' ended unsuccessfully [' + str(pr.exitcode) + ']') + print('Unsuccessfully terminated. (' + time.strftime("%c") + ')') + print('stuck in for loop') +print('not stuck in for loop') count, n_rows = 0, star.get_nrows() while count < n_rows: + print(f'in while loop: {count}') hold_out_pkl = q_pkls.get() star.set_element(key='_psGhMCFPickle', row=hold_out_pkl[0], val=hold_out_pkl[1]) count += 1 -out_star = os.path.join( - output_dir, - os.path.splitext(os.path.split(in_star)[1])[0] + '_mb_graph.star' - ) -print(f'{time.strftime( "%Y-%m-%d %H:%M:%S", time.localtime() )} Storing output STAR file in: ' + out_star, file=sys.stdout, flush=True) +out_star = output_dir + '/' + os.path.splitext(os.path.split(in_star)[1])[0] + '_mb_graph.star' +print('\tStoring output STAR file in: ' + out_star) star.store(out_star) -print('\nTerminated. (' + time.strftime("%c") + ')\n') - -#if __name__ == "__main__": - #run() - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +print('Terminated. (' + time.strftime("%c") + ')') diff --git a/code/tutorials/synth_sumb/tracing/mbu_picking.py b/code/tutorials/synth_sumb/tracing/mbu_picking.py index e80337e..517bcfd 100644 --- a/code/tutorials/synth_sumb/tracing/mbu_picking.py +++ b/code/tutorials/synth_sumb/tracing/mbu_picking.py @@ -17,6 +17,7 @@ __author__ = 'Antonio Martinez-Sanchez' +import argparse import operator import pyseg as ps from pyseg.globals.utils import coords_scale_supression @@ -26,18 +27,28 @@ ######################################################################################## -ROOT_PATH = '../../../..' +ROOT_PATH = '/scratch/users/muth9/simsiam/particle_picking' # Input STAR file -in_star = ROOT_PATH + '/data/tutorials/synth_sumb/fils/out/fil_mb_sources_to_no_mb_targets_net.star' +in_star = ROOT_PATH + '/data/fils/fil_mb_sources_to_no_mb_targets_net.star' ####### Output data -output_dir = ROOT_PATH + '/data/tutorials/synth_sumb/pick/out' +output_dir = ROOT_PATH + '/data/pick' ###### Slices settings file -slices_file = ROOT_PATH + '/data/tutorials/synth_sumb/pick/in/mb_cont_1.xml' +slices_file = ROOT_PATH + '/data/pick/mb_cont_1.xml' + +#Sarahs addition +parser = argparse.ArgumentParser() +parser.add_argument('--inStar', default=in_star, help='Input star file.') +parser.add_argument('--outDir', default=output_dir, help='Output subtomograms directory.') + +args = parser.parse_args() + +in_star = args.inStar +output_dir = args.outDir ###### Peaks configuration