diff --git a/src/supramolsim/analysis/sweep.py b/src/supramolsim/analysis/sweep.py index f3691cdf..81eeffce 100644 --- a/src/supramolsim/analysis/sweep.py +++ b/src/supramolsim/analysis/sweep.py @@ -64,6 +64,7 @@ def sweep_vasmples( "1XI5", ] if probes is None: + print("Probe list is None. Using default probe") probes = [ default_probe, ] @@ -89,54 +90,56 @@ def sweep_vasmples( for struct in structures: experiment.structure_id = struct experiment._build_structure() - for rep in range(repetitions): - probe_n = 0 - for probe in probes: - print(f"probe: {probe}") - # probe_param_n = 0 - for probe_param_n, p_param in probe_parameters.items(): - # copy p_param? - experiment.remove_probes() - if p_param is None: - experiment.add_probe( - probe_template=probe, - ) - #experiment.probe_parameters[probe] = default_params + #for rep in range(repetitions): + probe_n = 0 + for probe in probes: + print(f"probe: {probe}") + # probe_param_n = 0 + for probe_param_n, p_param in probe_parameters.items(): + # copy p_param? + experiment.remove_probes() + if p_param is None: + experiment.add_probe( + probe_template=probe, + ) + #experiment.probe_parameters[probe] = default_params + else: + p_param_copy = copy.deepcopy(p_param) + p_param_copy["fluorophore_id"] = default_fluorophore + experiment.add_probe( + probe_template=probe, + **p_param_copy + ) + for defect_n, defects_pars in particle_defects.items(): + particle_defects_copy = copy.deepcopy(defects_pars) + for d_key, d_val in particle_defects_copy.items(): + if d_key=="defect_small_cluster": + experiment.defect_eps["eps1"] = d_val + if d_key=="defect_large_cluster": + experiment.defect_eps["eps2"] = d_val + if d_key=="defect": + experiment.defect_eps["defect"] = d_val + # print(experiment.defect_eps) + print(experiment.probe_parameters) + experiment._build_particle(keep=True) + if experiment.generators_status("particle"): + if len(experiment.particle.emitters) == 0: + print(f"Skipping {probe}. No emitters were generated") + break else: - p_param_copy = copy.deepcopy(p_param) - p_param_copy["fluorophore_id"] = default_fluorophore - experiment.add_probe( - probe_template=probe, - **p_param_copy + break + #vsample_n = 0 + for vsample_n, vsample_pars in virtual_samples.items(): + _exported_field = None + # combination += str(vsample_n) + experiment.set_virtualsample_params( + **vsample_pars ) - for defect_n, defects_pars in particle_defects.items(): - particle_defects_copy = copy.deepcopy(defects_pars) - for d_key, d_val in particle_defects_copy.items(): - if d_key=="defect_small_cluster": - experiment.defect_eps["eps1"] = d_val - if d_key=="defect_large_cluster": - experiment.defect_eps["eps2"] = d_val - if d_key=="defect": - experiment.defect_eps["defect"] = d_val - # print(experiment.defect_eps) - print(experiment.probe_parameters) - experiment._build_particle(keep=True) - if experiment.generators_status("particle"): - if len(experiment.particle.emitters) == 0: - print(f"Skipping {probe}. No emitters were generated") - break - else: - break - #vsample_n = 0 - for vsample_n, vsample_pars in virtual_samples.items(): - _exported_field = None - # combination += str(vsample_n) - experiment.set_virtualsample_params( - **vsample_pars - ) - #_exported_field = experiment._build_coordinate_field( - # keep=False, use_self_particle=True, **vsample_pars - #) + #_exported_field = experiment._build_coordinate_field( + # keep=False, use_self_particle=True, **vsample_pars + #) + for rep in range(repetitions): + experiment.clear_virtual_sample() experiment.build(modules=["coordinate_field"], use_self_particle=True) _exported_field = experiment.coordinate_field.export_field() combination_n = ( @@ -159,11 +162,11 @@ def sweep_vasmples( vsample_params[combination_n] = _parameters vsample_outputs[combination_n] = [] vsample_outputs[combination_n].append(_exported_field) - # - # - #vsample_n += 1 - # defect_n += 1 - probe_n += 1 + # + # + #vsample_n += 1 + # defect_n += 1 + probe_n += 1 return experiment, vsample_outputs, vsample_params diff --git a/src/supramolsim/experiments.py b/src/supramolsim/experiments.py index ace99eb5..746fdc79 100644 --- a/src/supramolsim/experiments.py +++ b/src/supramolsim/experiments.py @@ -483,6 +483,9 @@ def _build_particle(self, lab_eff=1.0, defect_build=None, keep=False): deg_dissasembly=defect, ) else: + particle.add_defects( + deg_dissasembly=0, + ) print("Particle without defects") if keep: self.particle = particle @@ -942,6 +945,8 @@ def add_probe( probe_configuration["conjugation_target_info"] = ( probe_conjugation_target_info ) + probe_configuration["conjugation_sites"]["target"]["type"]= probe_conjugation_target_info["type"] + probe_configuration["conjugation_sites"]["target"]["value"]= probe_conjugation_target_info["value"] if probe_conjugation_efficiency is not None: probe_configuration["conjugation_efficiency"] = probe_conjugation_efficiency if probe_seconday_epitope is not None: @@ -956,6 +961,7 @@ def add_probe( probe_configuration["as_linker"] = False if probe_steric_hindrance is not None: probe_configuration["distance_between_epitope"] = probe_steric_hindrance + probe_configuration["binding"]["distance"]["between_targets"] = probe_steric_hindrance self.probe_parameters[probe_name] = probe_configuration self._update_probes() @@ -1138,7 +1144,8 @@ def current_settings(self, as_string=True, newline="
", modalities_acq_params def generate_virtual_sample( structure: str = "1XI5", - probe_name: str = None, + probe_template: str = "NHS_ester", + probe_name:str = None, probe_target_type: str = None, probe_target_value: str = None, probe_distance_to_epitope: float = None, @@ -1148,7 +1155,7 @@ def generate_virtual_sample( probe_conjugation_target_info=None, probe_conjugation_efficiency: float = None, probe_seconday_epitope=None, - probe_wobbling=False, + probe_wobble_theta=None, labelling_efficiency: float = 1.0, defect_small_cluster: float = None, defect_large_cluster: float = None, @@ -1160,6 +1167,7 @@ def generate_virtual_sample( random_orientations=False, random_placing=False, clear_probes=False, + clear_experiment = False, **kwargs, ): """ @@ -1223,19 +1231,24 @@ def generate_virtual_sample( - ExperimentParametrisation: The experiment containing all modules that were generated to build the virtual sample, and the virtual sample module itself. This experiment can be further used and tweaked for subsequent analysis or branching workflows. """ myexperiment = ExperimentParametrisation() + if clear_experiment: + myexperiment.clear_experiment() # load default configuration for probe if not clear_probes: - if probe_name is None: - probe_name = "NHS_ester" + print(probe_template) probe_configuration_file = os.path.join( - myexperiment.configuration_path, "probes", probe_name + ".yaml" + myexperiment.configuration_path, "probes", probe_template + ".yaml" ) probe_configuration = load_yaml(probe_configuration_file) - probe_configuration["probe_name"] = probe_name + probe_configuration["probe_template"] = probe_template + if probe_name is None: + probe_name = probe_template + else: + probe_configuration["label_name"] = probe_name if probe_target_type and probe_target_value: - probe_configuration["target_info"] = dict( - type=probe_target_type, value=probe_target_value - ) + print(probe_target_type, probe_target_value) + probe_configuration["probe_target_type"] = probe_target_type + probe_configuration["probe_target_value"] = probe_target_value if probe_distance_to_epitope is not None: probe_configuration["distance_to_epitope"] = probe_distance_to_epitope if probe_fluorophore is not None: @@ -1254,8 +1267,8 @@ def generate_virtual_sample( probe_configuration["conjugation_efficiency"] = probe_conjugation_efficiency if probe_seconday_epitope is not None: probe_configuration["epitope_target_info"] = probe_seconday_epitope - if probe_wobbling: - probe_configuration["enable_wobble"] = probe_wobbling + if probe_wobble_theta is not None: + probe_configuration["probe_wobble_theta"] = probe_wobble_theta myexperiment.add_probe(**probe_configuration) # load default configuration for virtual sample virtual_sample_template = os.path.join( diff --git a/src/supramolsim/generate/coordinates_field.py b/src/supramolsim/generate/coordinates_field.py index 532eab6f..1c824266 100644 --- a/src/supramolsim/generate/coordinates_field.py +++ b/src/supramolsim/generate/coordinates_field.py @@ -112,7 +112,9 @@ def init_from_file(self, field_yaml): # set molecules params self.set_molecules_params(**molecules) - def create_minimal_field(self, nmolecules=1, random_placing=False, random_orientations=False, **kwargs): + def create_minimal_field( + self, nmolecules=1, random_placing=False, random_orientations=False, **kwargs + ): """ Create a minimal field with a specified number of molecules and placement options. @@ -158,7 +160,7 @@ def create_minimal_field(self, nmolecules=1, random_placing=False, random_orient } else: # if none of the above, this is the case of a single particle in the center - #point = self.get_field_param("absolute_reference_point") + # point = self.get_field_param("absolute_reference_point") self.set_molecule_param("nMolecules", 1) if random_placing: self.random_placing = True @@ -169,8 +171,6 @@ def create_minimal_field(self, nmolecules=1, random_placing=False, random_orient self._set_fluo_plotting_params(fluo_name) self.random_orientations = random_orientations - - def calculate_absolute_reference(self): """ Calculate and set the absolute reference point based on the relative reference and dimension sizes. @@ -296,7 +296,9 @@ def generate_random_positions(self): # print(f"Total positions: {npositions}") # this can only work after the size of field has been established if self.molecules_params["minimal_distance"] is not None: - print(f"distributing with minimal distance: {self.molecules_params['minimal_distance']}") + print( + f"distributing with minimal distance: {self.molecules_params['minimal_distance']}" + ) self._random_pos_minimal_dist(npositions) else: print("Generating unconstrained random positions") @@ -318,7 +320,7 @@ def generate_random_orientations(self): orientations.append(np.array(sampl.sample_spherical_normalised(1, ndim=3))) self.set_molecule_param("orientations", orientations) - def generate_global_orientation(self, global_orientation = None): + def generate_global_orientation(self, global_orientation=None): """ Set a global orientation for all molecules. @@ -372,19 +374,19 @@ def _random_pos_minimal_dist(self, n): # print(np.linalg.norm(new - selected_positions[pos])) if np.linalg.norm(new - selected_positions[pos]) < minimal_distance: is_available = 0 - #break + # break if is_available: selected_positions.append(new) for abs_pos in selected_positions: rel_pos = abs_pos - rel_pos[0] = abs_pos[0]/maxabs_posx - rel_pos[1] = abs_pos[1]/maxabs_posy - rel_pos[2] = abs_pos[2]/maxabs_posz + rel_pos[0] = abs_pos[0] / maxabs_posx + rel_pos[1] = abs_pos[1] / maxabs_posy + rel_pos[2] = abs_pos[2] / maxabs_posz selected_relative.append(rel_pos) # print(f"end flag : {flag}") self.molecules_params["relative_positions"] = selected_relative - #self.set_molecule_param("absolute_positions", selected_positions) - #self.absolute_pos = selected_positions + # self.set_molecule_param("absolute_positions", selected_positions) + # self.absolute_pos = selected_positions def _calculate_absolute_position(self, relative_pos): fieldsizes = self.get_field_param("dimension_sizes") @@ -444,9 +446,7 @@ def create_molecules_from_InstanceObject(self, InstancePrototype: LabeledInstanc particle_copy.scale_coordinates_system(self.get_field_param("scale")) self.molecules_default_orientation = particle_copy.get_axis() if self.molecules_params["minimal_distance"] is None: - self._set_molecule_minimal_distance( - dist=particle_copy.radial_hindance - ) + self._set_molecule_minimal_distance(dist=particle_copy.radial_hindance) if self.random_placing: self.generate_random_positions() self._gen_abs_from_rel_positions() @@ -468,7 +468,6 @@ def create_molecules_from_InstanceObject(self, InstancePrototype: LabeledInstanc self.generate_random_orientations() # self.relabel_molecules() self.relabel_molecules() - def _create_instances_from_pdb(self, cif_f, structure_id, label_file): # this is exactly the same as the high level function @@ -513,7 +512,7 @@ def relabel_molecules(self): for mol in self.molecules: mol.generate_instance() # mol.show_instance() - #mol.scale_coordinates_system(self.get_field_param("scale")) + # mol.scale_coordinates_system(self.get_field_param("scale")) def relocate_molecules(self): """ @@ -578,9 +577,11 @@ def construct_static_field(self, relocate=True, reorient=True): self._construct_channels_by_fluorophores() else: default_fluorophore = "AF647" - #self.get_molecule_param("absolute_positions") + # self.get_molecule_param("absolute_positions") self.fluorophre_emitters = {} - self.fluorophre_emitters[default_fluorophore] = self.get_molecule_param("absolute_positions") + self.fluorophre_emitters[default_fluorophore] = self.get_molecule_param( + "absolute_positions" + ) def _construct_channels_by_fluorophores(self): """ @@ -639,7 +640,18 @@ def get_emitters_by_fluorophores(self): return dict(self.fluorophre_emitters) # methods for visualisation - def show_field(self, fluo_type="all", view_init=[30, 0, 0], initial_pos=True, return_fig=False): + def show_field( + self, + fluo_type="all", + view_init=[30, 0, 0], + initial_pos=True, + return_fig=False, + axesoff=False, + axis_object=None, + zoom_in=None, + emitters_plotsize=None, + **kwargs, + ): """ Visualize the field and emitters in 3D. @@ -659,14 +671,28 @@ def show_field(self, fluo_type="all", view_init=[30, 0, 0], initial_pos=True, re matplotlib.figure.Figure or None The figure if return_fig is True, otherwise None. """ - dimension_sizes = self.get_field_param("dimension_sizes") - xx, yy = np.meshgrid( - range(int(dimension_sizes[0])), range(int(dimension_sizes[1])) - ) - zz = (yy) * 0 - fig = plt.figure() - ax = fig.add_subplot(111, projection="3d") - ax.plot_surface(xx, yy, zz, alpha=0.2) + dimension_sizes = copy.copy(self.get_field_param("dimension_sizes")) + if zoom_in and zoom_in > 0.2 and zoom_in < 1: + x_to_remove = dimension_sizes[0] * zoom_in + x_size = dimension_sizes[0] - x_to_remove + y_to_remove = dimension_sizes[1] * zoom_in + y_size = dimension_sizes[1] - y_to_remove + xx, yy = np.meshgrid(range(int(x_size)), range(int(y_size))) + xx = xx + (x_to_remove / 2) + yy = yy + (y_to_remove / 2) + zz = (yy) * 0 + else: + xx, yy = np.meshgrid( + range(int(dimension_sizes[0])), range(int(dimension_sizes[1])) + ) + zz = (yy) * 0 + if axis_object is not None: + ax = axis_object + ax.plot_surface(xx, yy, zz, alpha=0.2) + else: + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + ax.plot_surface(xx, yy, zz, alpha=0.2) if initial_pos: if self.molecules_params["absolute_positions"] is not None: @@ -683,6 +709,8 @@ def show_field(self, fluo_type="all", view_init=[30, 0, 0], initial_pos=True, re print("Showing all fluorophores") for fname in self.fluorophre_emitters.keys(): print(fname) + if emitters_plotsize is not None: + self.plotting_params[fname]["plotsize"] = emitters_plotsize add_ax_scatter( ax, format_coordinates( @@ -690,6 +718,8 @@ def show_field(self, fluo_type="all", view_init=[30, 0, 0], initial_pos=True, re ), ) else: + if emitters_plotsize is not None: + self.plotting_params[fname]["plotsize"] = emitters_plotsize add_ax_scatter( ax, format_coordinates( @@ -701,10 +731,19 @@ def show_field(self, fluo_type="all", view_init=[30, 0, 0], initial_pos=True, re [ub - lb for lb, ub in (getattr(ax, f"get_{a}lim")() for a in "xyz")] ) ax.view_init(elev=view_init[0], azim=view_init[1], roll=view_init[2]) - if return_fig: - return fig + if axesoff: + ax.set_axis_off() + else: + ax.set_xlabel("X (Å)") + ax.set_ylabel("Y (Å)") + ax.set_zlabel("Z (Å)") + if axis_object is not None: + return ax else: - fig.show() + if return_fig: + return fig + else: + fig.show() def expand_isotropically(self, factor): """ @@ -772,7 +811,13 @@ def export_field(self): return export_field -def create_min_field(number_of_particles=1, random_placing=False, random_orientations=False, prints=False, **kwargs): +def create_min_field( + number_of_particles=1, + random_placing=False, + random_orientations=False, + prints=False, + **kwargs, +): """ Create a minimal field with a specified number of particles. @@ -797,19 +842,21 @@ def create_min_field(number_of_particles=1, random_placing=False, random_orienta if prints: print("Initialising default field") coordinates_field = Field() - #if molecule_pars: + # if molecule_pars: # for key, value in molecule_pars.items(): # coordinates_field.molecules_params[key] = value coordinates_field.create_minimal_field( nmolecules=number_of_particles, - random_placing=random_placing, + random_placing=random_placing, random_orientations=random_orientations, - **kwargs + **kwargs, ) return coordinates_field -def gen_positions_from_image(img, mode="mask", pixelsize = None, min_distance = None, **kwargs): +def gen_positions_from_image( + img, mode="mask", pixelsize=None, min_distance=None, **kwargs +): """ Generate relative positions from an image using either a mask or local maxima. @@ -834,7 +881,7 @@ def gen_positions_from_image(img, mode="mask", pixelsize = None, min_distance = if min_distance is None: min_distance = 1 else: - min_dist_pixels = min_distance/pixelsize + min_dist_pixels = min_distance / pixelsize min_distance = math.ceil(min_dist_pixels) npixels = list(img.shape) image_physical_size = np.zeros(shape=(2)) @@ -846,9 +893,7 @@ def gen_positions_from_image(img, mode="mask", pixelsize = None, min_distance = else: npositions = kwargs["npositions"] pixel_positions = sampl.get_random_pixels( - img, - num_pixels=npositions, - min_distance=min_distance + img, num_pixels=npositions, min_distance=min_distance ) elif mode == "localmaxima": if "background" not in kwargs.keys(): @@ -864,14 +909,13 @@ def gen_positions_from_image(img, mode="mask", pixelsize = None, min_distance = else: threshold = kwargs["threshold"] pixel_positions, img_processed = metrics.local_maxima_positions( - img, - min_distance=min_distance, - threshold=threshold, - sigma=sigma, - background=background) - xyz_relative = metrics.pixel_positions_to_relative( - pixel_positions, - image_sizes=image_physical_size, - pixelsize=pixelsize + img, + min_distance=min_distance, + threshold=threshold, + sigma=sigma, + background=background, ) - return xyz_relative, image_physical_size \ No newline at end of file + xyz_relative = metrics.pixel_positions_to_relative( + pixel_positions, image_sizes=image_physical_size, pixelsize=pixelsize + ) + return xyz_relative, image_physical_size diff --git a/src/supramolsim/generate/imaging.py b/src/supramolsim/generate/imaging.py index 2b0a2156..1bbacc99 100644 --- a/src/supramolsim/generate/imaging.py +++ b/src/supramolsim/generate/imaging.py @@ -657,6 +657,14 @@ def generate_imaging( asframes=False, ) beads = None + elif convolution_type == "raw_volume_activation": + images = conv.generate_frames_volume_convolution( + **field_data, + **psf_data, + asframes=False, + activation_only=True + ) + beads = None else: images, beads = self.images_by_convolutions( convolution_type, **simparams @@ -669,7 +677,7 @@ def generate_imaging( # wrap up the images from a single fluorophore in the channel output_per_fluoname[fluo] = dict(images=images, beads=beads) - if convolution_type != "raw_volume": + if convolution_type not in ["raw_volume", "raw_volume_activation"]: timeseries, beadstack = self._add_fluorophore_signals(output_per_fluoname) # # # Up to here only the photon information on arrival if noise: @@ -1068,6 +1076,7 @@ def show_field( initial_pos=True, reference_pt=False, axesoff=False, + emitters_plotsize=None, return_fig = False ): """ @@ -1110,7 +1119,7 @@ def show_field( zz = yy * 0 + (z_focus * factor) fig = plt.figure() ax = fig.add_subplot(111, projection="3d") - ax.plot_surface(xx, yy, zz, alpha=0.2, cmap="plasma") + ax.plot_surface(xx, yy, zz, alpha=0.2) # show ROI reference point if reference_pt: ref_pt = self.get_absoulte_reference_point() * factor @@ -1119,6 +1128,8 @@ def show_field( if self.emitters_by_fluorophore is not None: for fname, coords in self.emitters_by_fluorophore.items(): # print(fname) + if emitters_plotsize is not None: + self.plotting_params[fname]["plotsize"] = emitters_plotsize add_ax_scatter( ax, format_coordinates(coords * factor, **self.plotting_params[fname]), diff --git a/src/supramolsim/generate/labelled_instance.py b/src/supramolsim/generate/labelled_instance.py index df658d20..fa80f0f4 100644 --- a/src/supramolsim/generate/labelled_instance.py +++ b/src/supramolsim/generate/labelled_instance.py @@ -10,7 +10,7 @@ transform_displace_set, ) from ..utils.data_format.visualisation import format_coordinates, set_colorplot -from ..utils.visualisation.matplotlib_plots import add_ax_scatter, draw1nomral_segment +from ..utils.visualisation.matplotlib_plots import add_ax_scatter, draw1nomral_segment, draw_nomral_segments from ..utils.transform.cif_builder import create_instance_label from ..utils.transform.defects import xmersubset_byclustering from ..utils.data_format.structural_format import builder_format @@ -503,18 +503,23 @@ def add_defects( minsamples : int, optional Minimum samples for clustering. Default is 1. """ - d_cluster_params = dict( - eps1=eps1, - minsamples1=minsamples, - eps2=xmer_neigh_distance, - minsamples2=minsamples, - ) - self.defects_params["d_cluster_params"] = d_cluster_params - self.defects_params["deg_dissasembly"] = deg_dissasembly - self.defects_params["xmer_neigh_distance"] = xmer_neigh_distance - self.defects_params["fracture"] = fracture - self.defects = True - self.generate_instance() + if deg_dissasembly == 0: + self.defects = False + self.defects_target_normals = None + self.generate_instance() + else: + d_cluster_params = dict( + eps1=eps1, + minsamples1=minsamples, + eps2=xmer_neigh_distance, + minsamples2=minsamples, + ) + self.defects_params["d_cluster_params"] = d_cluster_params + self.defects_params["deg_dissasembly"] = deg_dissasembly + self.defects_params["xmer_neigh_distance"] = xmer_neigh_distance + self.defects_params["fracture"] = fracture + self.defects = True + self.generate_instance() def _model_defects(self, target_normals_dictionary): """ @@ -684,6 +689,11 @@ def transform_translate(self, newcenter): self.emitters[labeltype], _ = transform_displace_set( self.get_emitter_by_target(labeltype), self.get_ref_point(), nref ) + # translate sources as well + self.source["targets"][labeltype]["coordinates"], _ = transform_displace_set( + self._get_source_coords_normals(labeltype)["coordinates"], self.get_ref_point(), nref + ) + # then replace reference point self._set_ref_point(nref) self.axis["pivot"] = nref @@ -753,6 +763,8 @@ def scale_coordinates_system(self, new_scale: float): probe_scaling_factor = self.secondary[labeltype]["scale"] / new_scale if self.secondary[labeltype]["emitters"] is not None: self.secondary[labeltype]["emitters"] = self.secondary[labeltype]["emitters"].astype('float64') * probe_scaling_factor + if self.secondary[labeltype]["minimal_distance"] is not None: + self.secondary[labeltype]["minimal_distance"] *= probe_scaling_factor if "coordinates" in self.secondary[labeltype].keys(): self.secondary[labeltype]["coordinates"] = self.secondary[labeltype]["coordinates"].astype('float64') * probe_scaling_factor if self.secondary[labeltype]["binding"]["distance"]["to_target"] is not None: @@ -760,6 +772,10 @@ def scale_coordinates_system(self, new_scale: float): if self.secondary[labeltype]["binding"]["distance"]["between_targets"] is not None: self.secondary[labeltype]["binding"]["distance"]["between_targets"] *= probe_scaling_factor self.secondary[labeltype]["scale"] = new_scale + if self.defects: + self.defects_params["d_cluster_params"]["eps1"] *= scaling_factor + self.defects_params["d_cluster_params"]["eps2"] *= scaling_factor + self.defects_params["xmer_neigh_distance"] *= scaling_factor # methods to get emitters by target name @@ -818,6 +834,7 @@ def show_probe(self, xlims = [-100,100], ylims = [-100,100], zlims = None, + central_axis=True, **kwargs): """ Visualize the probe structure in 3D. @@ -851,7 +868,7 @@ def show_probe(self, probe_name = self.emitters[first_probe] fig = plt.figure() ax = fig.add_subplot(111, projection="3d") - probe_plotting_params = self._get_label_plotting_params(probe_name) + #probe_plotting_params = self._get_label_plotting_params(probe_name) total_coordinates = copy.copy(self.labels[probe_name]["emitters"]) total_number_coordinates = total_coordinates.shape[0] center = np.mean(total_coordinates, axis=0) @@ -864,19 +881,15 @@ def show_probe(self, format_coordinates( centered_emitters, plotmarker="o", - plotcolour="#984ea3", - plotsize=20 + plotsize=20, + **kwargs ), ) - add_ax_scatter( - ax, - format_coordinates( - centered_axis, - plotmarker="x", - plotcolour="k", - plotsize=20 - ), - ) + if central_axis: + axis_segment = dict() + axis_segment["pivot"] = centered_axis[0] + axis_segment["direction"] = centered_axis[1] - centered_axis[0] + draw1nomral_segment(axis_segment, ax, lenght=100, colors=["g", "y"]) ax.set_xlim(xlims) ax.set_ylim(ylims) if zlims is None: @@ -891,9 +904,9 @@ def show_probe(self, if axesoff: ax.set_axis_off() else: - ax.set_xlabel("X (Angstroms)") - ax.set_ylabel("Y (Angstroms)") - ax.set_zlabel("Z (Angstroms)") + ax.set_xlabel("X (Å)") + ax.set_ylabel("Y (Å)") + ax.set_zlabel("Z (Å)") if return_plot: return ax else: @@ -910,8 +923,8 @@ def show_instance( show_axis=False, with_sources=False, return_plot=False, - source_size=1, - emitter_plotsize=1, + source_size=None, + emitter_plotsize=None, ): """ Visualize the labelled instance in 3D. @@ -947,7 +960,8 @@ def show_instance( ax = fig.add_subplot(111, projection="3d") for labs in self.labelnames: lab_plotparams = self._get_label_plotting_params(labs) - lab_plotparams["plotsize"] = emitter_plotsize + if emitter_plotsize is not None: + lab_plotparams["plotsize"] = emitter_plotsize add_ax_scatter( ax, format_coordinates( @@ -955,6 +969,10 @@ def show_instance( ), ) if with_sources: + if source_size is not None: + source_size=source_size + else: + source_size=1 add_ax_scatter( ax, format_coordinates( @@ -998,9 +1016,12 @@ def gen_axis_plot( axesoff=True, show_axis=False, with_sources=False, - source_plotsize=1, + source_plotsize=None, axis_object=None, - emitter_plotsize=1, + emitter_plotsize=None, + source_plotcolour="#bbbbbb", + with_normals=False, + **kwargs ): """ Generate a 3D axis plot for the labelled instance. @@ -1030,7 +1051,8 @@ def gen_axis_plot( if labelnames == "All": for labs in self.labelnames: lab_plotparams = self._get_label_plotting_params(labs) - lab_plotparams["plotsize"] = emitter_plotsize + if emitter_plotsize is not None: + lab_plotparams["plotsize"] = emitter_plotsize add_ax_scatter( axis_object, format_coordinates(self.emitters[labs], **lab_plotparams), @@ -1038,6 +1060,10 @@ def gen_axis_plot( zlims[0] = np.min(self.emitters[labs]) zlims[1] = np.max(self.emitters[labs]) if with_sources: + if source_plotsize is not None: + source_plotsize = source_plotsize + else: + source_plotsize = 1 if self.defects_target_normals is not None: add_ax_scatter( axis_object, @@ -1049,7 +1075,7 @@ def gen_axis_plot( axis_object, format_coordinates( self._get_source_coords_normals(labs)["coordinates"], - plotcolour="#bbbbbb", + plotcolour=source_plotcolour, plotalpha=0.5, plotsize=source_plotsize, ), @@ -1062,6 +1088,13 @@ def gen_axis_plot( plotsize=source_plotsize, ), ) + if with_normals: + draw_nomral_segments( + [self.source["targets"][labs]["normals"], self.source["targets"][labs]["coordinates"]], + axis_object, + colors=["grey", "green"], + **kwargs + ) # add_ax_scatter(ax, format_coordinates(self.emitters[labs])) if reference_point: diff --git a/src/supramolsim/generate/molecular_structure.py b/src/supramolsim/generate/molecular_structure.py index e9fe3f14..73354305 100644 --- a/src/supramolsim/generate/molecular_structure.py +++ b/src/supramolsim/generate/molecular_structure.py @@ -309,6 +309,7 @@ def generate_assembly_reference_point(self): """ # this function could use a small portion of the atoms defined only allatoms_defined = self.get_atoms_infile() # Parse all atoms in file + self.atoms_in_file = allatoms_defined if self.assembly_operations is None: self.generate_assemmbly_operations() if self.assymetric_defined: @@ -767,7 +768,12 @@ def show_target_labels( axesoff=True, show_axis=True, with_normals=False, - return_plot=False + return_plot=False, + axis_object=None, + target_size = None, + target_plotcolour = None, + atoms_size = None, + atoms_alpha = None, ): """ Visualize target labels and optionally assembly atoms and reference point. @@ -798,10 +804,17 @@ def show_target_labels( matplotlib.figure.Figure or None The figure if return_plot is True, otherwise None. """ - if labelnames is None: + if axis_object is not None: + ax = axis_object + else: fig = plt.figure() ax = fig.add_subplot(111, projection="3d") + if labelnames is None: for trgt in list(self.label_targets.keys()): + if target_size is not None: + self.plotting_params[trgt]["plotsize"] = target_size + if target_plotcolour is not None: + self.plotting_params[trgt]["plotcolour"] = target_plotcolour if with_normals: draw_nomral_segments( [self.get_target_normals(trgt), self.get_target_coords(trgt)], @@ -820,6 +833,10 @@ def show_target_labels( atoms_subset = cif_builder.array_coords_subset( self.assembly_atoms, assembly_fraction ) + if atoms_size is not None: + self.plotting_params["assemblyatoms"]["plotsize"] = atoms_size + if atoms_alpha is not None: + self.plotting_params["assemblyatoms"]["plotalpha"] = atoms_alpha add_ax_scatter( ax, format_coordinates( @@ -846,11 +863,14 @@ def show_target_labels( ) if axesoff: ax.set_axis_off() - if return_plot: - plt.close() - return fig + if axis_object is not None: + return ax else: - fig.show + if return_plot: + plt.close() + return fig + else: + fig.show def show_assembly_atoms( self, assembly_fraction=0.01, view_init=[0, 0, 0], axesoff=True,return_plot=False diff --git a/src/supramolsim/utils/transform/image_convolution.py b/src/supramolsim/utils/transform/image_convolution.py index 9c60e0de..ad0d646d 100644 --- a/src/supramolsim/utils/transform/image_convolution.py +++ b/src/supramolsim/utils/transform/image_convolution.py @@ -77,6 +77,7 @@ def frame_by_volume_convolution( kernel3D, zfocus_slice, projection_depth=2, + activation_only=False, asframe=True ): """ @@ -93,6 +94,8 @@ def frame_by_volume_convolution( intensity_voxel = voxelize_active_fluorophores_withrange( coordinates, photon_vector, ranges, bin_pixelsize=bin_size ) + if activation_only: + return intensity_voxel if np.sum(photon_vector) > 0: convolved_intensity = convolve3D(intensity_voxel, kernel3D) else: @@ -132,6 +135,7 @@ def generate_frames_volume_convolution( psf_focus_slice: int, psf_pixelsizeXY: float, asframes=True, + activation_only=False, psf_projection_depth=1, **kwargs, ): @@ -202,7 +206,8 @@ def generate_frames_volume_convolution( psf_array, psf_focus_slice, projection_depth=psf_projection_depth, - asframe=False + asframe=asframes, + activation_only=activation_only ) volume_frames.append(v_frame) return volume_frames diff --git a/src/supramolsim/utils/visualisation/matplotlib_plots.py b/src/supramolsim/utils/visualisation/matplotlib_plots.py index ccf99b79..9836ab15 100644 --- a/src/supramolsim/utils/visualisation/matplotlib_plots.py +++ b/src/supramolsim/utils/visualisation/matplotlib_plots.py @@ -114,7 +114,7 @@ def draw1nomral_segment(points_normal, figure, lenght=100, colors=["g", "y"]): figure.scatter(ends[0], ends[1], ends[2], color=colors[1], marker="o") -def draw_nomral_segments(points_normal, figure, lenght=100, colors=["g", "y"]): +def draw_nomral_segments(points_normal, figure, lenght=100, colors=["g", "y"], **kwargs): # points_normals is a list of 2 elements # first element are the normals # second element are ponts in space @@ -138,15 +138,27 @@ def draw_nomral_segments(points_normal, figure, lenght=100, colors=["g", "y"]): return figure -def stack_projection(stack, angle=45, axes=(1,2), method = "max"): +def stack_projection(stack, angle=45, axes=(1,2), method = "sd", zdepth = False): projection = None zstack = copy.copy(stack) rotated = rotate(zstack, angle=angle, axes=axes, reshape=False, order=3) # order=3: cubic interpolation + if zdepth: + nslices = rotated.shape[2] + for i in range(1, nslices): + #intensity_factor = np.log(1.8 + ((nslices - i)/(nslices))) + intensity_factor = 0.6 * np.exp(((nslices - i)/(nslices)) * (1/(i+1))) + rotated[:,:,i] = rotated[:,:,i]*intensity_factor if method == "max": - projection = np.max(rotated, axis=2) + projection = np.max(rotated, axis=2) + elif method == "sd": + projection = np.std(rotated, axis=2) + elif method == "mean": + projection = np.mean(rotated, axis=2) + elif method == "sum": + projection = np.sum(rotated, axis=2) return projection -def plot_projection(stack, angle=45, plane="XY", method = "max"): +def plot_projection(stack, angle=45, plane="XY", method = "max", zdepth=False): if plane == "XY": axes = (1, 2) elif plane == "XZ": @@ -155,7 +167,7 @@ def plot_projection(stack, angle=45, plane="XY", method = "max"): axes = (0, 1) else: raise ValueError("Invalid plane. Choose from 'XY', 'XZ', or 'YZ'.") - projection = stack_projection(stack, angle=angle, axes=axes, method=method) + projection = stack_projection(stack, angle=angle, axes=axes, method=method, zdepth=zdepth) plt.imshow(projection, cmap='gray') plt.axis('off') plt.show() \ No newline at end of file diff --git a/src/supramolsim/workflows.py b/src/supramolsim/workflows.py index 4640c976..4a7e20ad 100644 --- a/src/supramolsim/workflows.py +++ b/src/supramolsim/workflows.py @@ -173,7 +173,7 @@ def probe_model( # print(structural_model.label_targets) - return structural_model, target_sites, anchor_point, direction_point, probe_epitope + return structural_model, target_sites, np.squeeze(anchor_point, axis=None), np.squeeze(direction_point, axis=None), probe_epitope def particle_from_structure( @@ -225,7 +225,10 @@ def particle_from_structure( ) if anchor_point.shape == (3,): print("setting new axis") - label_object.set_axis(pivot=anchor_point, direction=direction_point) + print(anchor_point, direction_point) + direction_vector = direction_point - anchor_point + label_object.set_axis(pivot=anchor_point, direction=direction_vector) + print(f'label object axis: {label_object.params["axis"]}') if ( probe_epitope["coordinates"] is not None and label_params["as_linker"] @@ -240,7 +243,7 @@ def particle_from_structure( label_object.set_emitters(probe_emitter_sites) label_params["coordinates"] = label_object.gen_labeling_entity() print(label_params["coordinates"]) - label_params_list.append(label_params) + structure.add_label(label_object) if label_params["binding"]["distance"]["to_target"] is not None: print("Label is indirect label") @@ -250,6 +253,7 @@ def particle_from_structure( if label_params["binding"]["distance"]["between_targets"] == "estimate": distances = pdist(label_params["coordinates"]) label_params["binding"]["distance"]["between_targets"] = np.max(distances) + label_params_list.append(label_params) inst_builder = structure.create_instance_builder() particle = labinstance.create_particle( source_builder=inst_builder, label_params_list=label_params_list diff --git a/tests/test_feature_defects.py b/tests/test_feature_defects.py index 6fa8a5fa..5c993d67 100644 --- a/tests/test_feature_defects.py +++ b/tests/test_feature_defects.py @@ -13,4 +13,7 @@ def test_add_defects(defect, experiment_7r5k_base): xmer_neigh_distance=600, deg_dissasembly=defect, ) - assert experiment_7r5k_base.particle.defects_target_normals is not None + if defect == 0: + assert experiment_7r5k_base.particle.defects_target_normals is None + else: + assert experiment_7r5k_base.particle.defects_target_normals is not None