debug the script for plot and clean the particles for standard template matching

test/s210/tm/plot_and_clean_results_reomve_clash.py

@@ -6,6 +6,7 @@ from file_handler import read_mrc, write_mrc
 from scipy.spatial import cKDTree
 from tqdm import tqdm
 import os
+from config import shrinkage_factor,tomogram
 
 current_dir = os.getcwd()
 split_path = current_dir.split("/")
@@ -18,7 +19,7 @@ coords = df[["x", "y", "z"]].to_numpy()
 rots = df[["phi", "theta", "psi"]].to_numpy()
 ccc_values = df["ccc"].to_numpy()
 # Distance cutoff
-distance_cutoff = 5
+distance_cutoff = 30
 
 # Initialize KDTree
 tree = cKDTree(coords)
@@ -54,18 +55,19 @@ filtered_df = pd.DataFrame(
     }
 )
 print("Filtered particle number: ", len(filtered_df))
-vol_array = np.zeros([40, 40, 40])
-filtered_df.to_csv(tomoID + ".tm_cut" + str(distance_cutoff) + ".csv", index=False)
+tomo = read_mrc(tomogram)
+vol_array = np.zeros(tomo.shape)
+filtered_df.to_csv('results/'+tomoID + ".tm_cut" + str(distance_cutoff) + ".csv", index=False)
 
-template = read_mrc("inputs/mono-8Apx-lps30-box30-rot12-core.mrc")
+template = read_mrc("inputs/emd_8799_Yeast80S_ribosome_pix8.24_vol40_low40.mrc")
 df2 = pd.DataFrame()
-ccc_list = [0.3, 0.2]
+ccc_list = [0.4, 0.33]
 
 print("start mapping")
 for i in tqdm(range(len(filtered_df))):
     template_rot = rotate(template, filtered_rots[i])
     vol_array, ccc = try_add_obj(
-        vol_array, template_rot, filtered_coords[i], 1, filtered_ccc[i]
+        vol_array, template_rot, filtered_coords[i], shrinkage_factor, filtered_ccc[i]
     )
     if ccc > 0:
         # save the dataframe row i to df2
@@ -73,11 +75,11 @@ for i in tqdm(range(len(filtered_df))):
 
         if ccc < ccc_list[0]:
             print("with " + str(ccc) + "got paritcle " + str(len(df2)))
-            df2.to_csv(tomoID + ".test1.tm_" + str(ccc_list[0]) + ".csv", index=False)
+            df2.to_csv('results/'+tomoID + ".test1.tm_" + str(ccc_list[0]) + ".csv", index=False)
 
             write_mrc(
                 np.float32(vol_array),
-                tomoID + ".test1.tm_" + str(ccc_list[0]) + "_rc.mrc",
+                'results/'+tomoID + ".test1.tm_" + str(ccc_list[0]) + "_rc.mrc",
             )
             ccc_list.pop(0)
         if not ccc_list:

test/s68/tm/plot_and_clean_results_reomve_clash.py

@@ -6,6 +6,7 @@ from file_handler import read_mrc, write_mrc
 from scipy.spatial import cKDTree
 from tqdm import tqdm
 import os
+from config import shrinkage_factor,tomogram
 
 current_dir = os.getcwd()
 split_path = current_dir.split("/")
@@ -18,7 +19,7 @@ coords = df[["x", "y", "z"]].to_numpy()
 rots = df[["phi", "theta", "psi"]].to_numpy()
 ccc_values = df["ccc"].to_numpy()
 # Distance cutoff
-distance_cutoff = 5
+distance_cutoff = 30
 
 # Initialize KDTree
 tree = cKDTree(coords)
@@ -54,18 +55,19 @@ filtered_df = pd.DataFrame(
     }
 )
 print("Filtered particle number: ", len(filtered_df))
-vol_array = np.zeros([40, 40, 40])
-filtered_df.to_csv(tomoID + ".tm_cut" + str(distance_cutoff) + ".csv", index=False)
+tomo = read_mrc(tomogram)
+vol_array = np.zeros(tomo.shape)
+filtered_df.to_csv('results/'+tomoID + ".tm_cut" + str(distance_cutoff) + ".csv", index=False)
 
-template = read_mrc("inputs/mono-8Apx-lps30-box30-rot12-core.mrc")
+template = read_mrc("inputs/emd_8799_Yeast80S_ribosome_pix8.24_vol40_low40.mrc")
 df2 = pd.DataFrame()
-ccc_list = [0.3, 0.2]
+ccc_list = [0.4, 0.33]
 
 print("start mapping")
 for i in tqdm(range(len(filtered_df))):
     template_rot = rotate(template, filtered_rots[i])
     vol_array, ccc = try_add_obj(
-        vol_array, template_rot, filtered_coords[i], 1, filtered_ccc[i]
+        vol_array, template_rot, filtered_coords[i], shrinkage_factor, filtered_ccc[i]
     )
     if ccc > 0:
         # save the dataframe row i to df2
@@ -73,11 +75,11 @@ for i in tqdm(range(len(filtered_df))):
 
         if ccc < ccc_list[0]:
             print("with " + str(ccc) + "got paritcle " + str(len(df2)))
-            df2.to_csv(tomoID + ".test1.tm_" + str(ccc_list[0]) + ".csv", index=False)
+            df2.to_csv('results/'+tomoID + ".test1.tm_" + str(ccc_list[0]) + ".csv", index=False)
 
             write_mrc(
                 np.float32(vol_array),
-                tomoID + ".test1.tm_" + str(ccc_list[0]) + "_rc.mrc",
+                'results/'+tomoID + ".test1.tm_" + str(ccc_list[0]) + "_rc.mrc",
             )
             ccc_list.pop(0)
         if not ccc_list:

tools/counter_level.py

@@ -7,8 +7,8 @@ current_dir = os.getcwd()
 split_path = current_dir.split("/")
 tomoID = split_path[-2]
 
-vol = sys.argv[1]
-ctf = sys.argv[2]
+vol = read_mrc(sys.argv[1])
+ctf = read_mrc(sys.argv[2])
 
 rescaled_ctf = rescale_ctf_volume(ctf,ctf.shape,vol.shape)
 ctf_vol = apply_ctf(vol,rescaled_ctf)

tools/plot_and_clean_results_reomve_clash.py

+from utils import rotate
+from utils import try_add_obj
+import pandas as pd
+import numpy as np
+from file_handler import read_mrc, write_mrc
+from scipy.spatial import cKDTree
+from tqdm import tqdm
+import os
+from config import shrinkage_factor,tomogram
+
+current_dir = os.getcwd()
+split_path = current_dir.split("/")
+tomoID = split_path[-2]
+
+df = pd.read_csv("results/" + tomoID + ".test1.tm.csv")
+# df = df[(df["x"] > 20) & (df["y"] > 20) & (df["z"] > 10)]
+# df = df[(df["x"] < 276) & (df["y"] < 276) & (df["z"] < 110)]
+coords = df[["x", "y", "z"]].to_numpy()
+rots = df[["phi", "theta", "psi"]].to_numpy()
+ccc_values = df["ccc"].to_numpy()
+# Distance cutoff
+distance_cutoff = 30
+
+# Initialize KDTree
+tree = cKDTree(coords)
+
+# Array to keep track of points to keep
+keep = np.ones(len(coords), dtype=bool)
+print("Total partilce number: ", len(df))
+# Filter coordinates and rotations
+for i in tqdm(range(len(coords))):
+    if keep[i]:
+        # Find all points within the distance cutoff
+        indices = tree.query_ball_point(coords[i], distance_cutoff)
+        # Set keep to False for all points within the cutoff distance, except the current point
+        for index in indices:
+            if index > i:
+                keep[index] = False
+
+
+# Filter the arrays
+filtered_coords = coords[keep]
+filtered_rots = rots[keep]
+filtered_ccc = ccc_values[keep]
+# Combine filtered results into a new DataFrame
+filtered_df = pd.DataFrame(
+    {
+        "x": filtered_coords[:, 0],
+        "y": filtered_coords[:, 1],
+        "z": filtered_coords[:, 2],
+        "phi": filtered_rots[:, 0],
+        "theta": filtered_rots[:, 1],
+        "psi": filtered_rots[:, 2],
+        "ccc": filtered_ccc,
+    }
+)
+print("Filtered particle number: ", len(filtered_df))
+tomo = read_mrc(tomogram)
+vol_array = np.zeros(tomo.shape)
+filtered_df.to_csv('results/'+tomoID + ".tm_cut" + str(distance_cutoff) + ".csv", index=False)
+
+template = read_mrc("inputs/emd_8799_Yeast80S_ribosome_pix8.24_vol40_low40.mrc")
+df2 = pd.DataFrame()
+ccc_list = [0.4, 0.33]
+
+print("start mapping")
+for i in tqdm(range(len(filtered_df))):
+    template_rot = rotate(template, filtered_rots[i])
+    vol_array, ccc = try_add_obj(
+        vol_array, template_rot, filtered_coords[i], shrinkage_factor, filtered_ccc[i]
+    )
+    if ccc > 0:
+        # save the dataframe row i to df2
+        df2 = pd.concat([df2, pd.DataFrame(filtered_df.iloc[i]).T], ignore_index=1)
+
+        if ccc < ccc_list[0]:
+            print("with " + str(ccc) + "got paritcle " + str(len(df2)))
+            df2.to_csv('results/'+tomoID + ".test1.tm_" + str(ccc_list[0]) + ".csv", index=False)
+
+            write_mrc(
+                np.float32(vol_array),
+                'results/'+tomoID + ".test1.tm_" + str(ccc_list[0]) + "_rc.mrc",
+            )
+            ccc_list.pop(0)
+        if not ccc_list:
+            print("Done!")
+            break