[fix/doc]add visualizing radonfit membrane analysis jupyter notebook …

…and fix some small bugs

.github/workflows/ci.yml

@@ -25,4 +25,5 @@ jobs:
           restore-keys: |
             mkdocs-material-
       - run: pip install mkdocs-material
+      - run: pip install mkdocs-jupyter
       - run: mkdocs gh-deploy --force

build/lib/memxterminator/__init__.py

@@ -1,3 +1,3 @@
 __author__ = 'Zhen Huang'
 __email__ = '[email protected]'
-__version__ = '1.0.0'
+__version__ = '1.2.1'

build/lib/memxterminator/radonfit/lib/calculate_curve.py

@@ -10,16 +10,7 @@
 
 class Curve:
     def __init__(self, output_filename, i, image, y0, x0, theta, kappa):
-        # self.radonanalyze = RadonAnalyzer(image, thr=thr)
-        # self.centerfit = Template_centerfitting(3, 3, image, thr=thr)
         df_star = readstar(output_filename)
-        # self.y0, self.x0 = df_star.loc[0, 'rlnCenterX'], df_star.loc[0, 'rlnCenterY'] # center of membrane
-        # if mode == 0:
-        #     self.theta = df_star.loc[0, 'rlnAngleTheta'] * np.pi / 180
-        # elif mode == 1:
-        #     self.theta = theta
-        #     # self.y0, self.x0 = self.centerfit.centerfinder() # center of membrane
-        #     self.y0, self.x0 = y0, x0
         self.y0 = y0
         self.x0 = x0
         self.theta = theta
@@ -80,17 +71,16 @@ def compute_line_dist(self):
             return self.distance_matrix
     def generate_membrane(self):
         self.distance_matrix = self.compute()
+        if self.kappa >= 0:
+            self.distance_matrix = -self.distance_matrix
         self.simulate_membrane = gaussian2_gpu(self.distance_matrix, self.membrane_distance, self.sigma1, self.sigma2)
         return self.simulate_membrane
 
 class Curvefitting:
     def __init__(self, output_filename, i, image, kappa_start, kappa_end, kappa_step):
-        # self.thr = thr
         df_star = readstar(output_filename)
         self.output_filename = output_filename
         self.i = i
-        # radonanalyze = RadonAnalyzer(image, thr=thr)
-        # centerfit = Template_centerfitting(3, 3, image, thr=thr)
         num = int((kappa_end - kappa_start) / kappa_step) + 1
         self.kappa_lst = np.linspace(kappa_start, kappa_end, num)
         self.corr_lst = []
@@ -110,7 +100,6 @@ def fit_curve(self):
         i = 0
         print('>>>Start curve fitting...')
         for kappa in self.kappa_lst:
-            # print(f'=====iter: {i}=====')
             self.simulated_membrane  = self.generate_membrane(kappa)
             self.simulated_membrane = self.simulated_membrane.astype(cp.float64)
             self.gray_image = self.gray_image.astype(cp.float64)
@@ -129,24 +118,4 @@ def fit_curve(self):
         self.best_kappa = self.kappa_lst[self.corr_max_index]
         print('Best kappa:', self.best_kappa)
         self.mem_best = self.generate_membrane(self.best_kappa)
-        return self.best_kappa
-    def fit_curve_visualize(self):
-        fig, axes = plt.subplots(1, 3, figsize=(10, 5))
-        axes[0].plot(self.kappa_lst, self.corr_lst, 'x', color = 'red')
-        axes[2].imshow(self.mem_best, cmap='gray', origin='lower')
-        axes[1].imshow(self.gray_image, cmap='gray', origin='lower')
-        plt.show()
-
-if __name__ == '__main__':
-    # gene_curve = Curve('image2.mrc', section=1, thr=0.5, kappa=0)
-    # membrane = gene_curve.generate_membrane()
-    # plt.imshow(membrane, cmap='gray', origin='lower')
-    # plt.show()
-    image = readmrc('neuron_templates.mrc', section=3, mode='gpu')
-    image = zoom(image, 256/64)
-    curve_fitting = Curvefitting(image, -0.01, 0.01, 21)
-    # curve_fitting.fit_curve()
-    # curve_fitting.fit_curve_visualize()
-    # mem_sub = curve_fitting.mem_subtract()
-    # plt.imshow(mem_sub, cmap='gray', origin='lower')
-    # plt.show()
+        return self.best_kappa

build/lib/memxterminator/radonfit/lib/generate_gaussian_template.py

@@ -32,7 +32,7 @@ def gaussian2_gpu(x, membrane_dist, std1, std2):
 def generate_template(size, theta, membrane_dist, std1, std2, mode='cpu'):
     if mode == 'cpu':
         distance_matrix = np.zeros((size, size))
-        theta = (theta + 90) * np.pi / 180
+        theta = (theta - 90) * np.pi / 180
         x = range(size)
         y = range(size)
         x0 = size / 2
@@ -47,7 +47,8 @@ def generate_template(size, theta, membrane_dist, std1, std2, mode='cpu'):
             for i in x:
                 for j in y:
                     distance_matrix[i, j] = (k * i - j + b) / np.sqrt(k**2 + 1)
-        template = gaussian2(distance_matrix, membrane_dist, std1, std2)
+        # add a minus sign due to some annoying reasons
+        template = gaussian2(-distance_matrix, membrane_dist, std1, std2)
         return template
     elif mode == 'gpu':
         distance_matrix = cp.zeros((size, size))
@@ -56,7 +57,7 @@ def generate_template(size, theta, membrane_dist, std1, std2, mode='cpu'):
         y = range(size)
         x0 = size / 2
         y0 = size / 2
-        if cp.isclose(theta, cp.pi/2) or np.isclose(theta, 3*cp.pi/2):
+        if cp.isclose(theta, cp.pi/2) or cp.isclose(theta, 3*cp.pi/2):
             for i in x:
                 for j in y:
                     distance_matrix[i, j] = i - x0
@@ -66,15 +67,6 @@ def generate_template(size, theta, membrane_dist, std1, std2, mode='cpu'):
             for i in x:
                 for j in y:
                     distance_matrix[i, j] = (k * i - j + b) / cp.sqrt(k**2 + 1)
-        template = gaussian2_gpu(distance_matrix, membrane_dist, std1, std2)
-        return template
-
-if __name__ == '__main__':
-    x = np.linspace(-100, 100, 200)
-    g = gaussian2(x,12,4,4)
-    alpha = 30
-    size = 32
-    template_rotated = generate_template(size, alpha, 12, 4, 2)
-    plt.imshow(template_rotated, cmap='gray', origin='lower')
-    plt.show()
-
+        # add a minus sign due to some annoying reasons
+        template = gaussian2_gpu(-distance_matrix, membrane_dist, std1, std2)
+        return template

build/lib/memxterminator/radonfit/lib/generate_membrane_mask.py

@@ -9,21 +9,14 @@
 
 class mem_mask:
     def __init__(self,output_filename, i, image, edge_sigma):
-        # fit_curve = Curvefitting(image, thr)
-        # self.gray_image = fit_curve.gray_image
         self.i = i
         self.output_filename = output_filename
         df_star = readstar(output_filename)
         self.gray_image = image
-        # self.thr = thr
         self.edge_sigma = edge_sigma
-        # self.y0, self.x0 = fit_curve.yc, fit_curve.xc
         self.x0, self.y0 = df_star.loc[i, 'rlnCenterY'], df_star.loc[i, 'rlnCenterX']
-        # self.theta = fit_curve.theta
         self.theta = df_star.loc[i, 'rlnAngleTheta'] * np.pi / 180
-        # self.membrane_distance = fit_curve.membrane_distance
         self.membrane_distance = df_star.loc[i, 'rlnMembraneDistance']
-        # self.kappa = fit_curve.best_kappa
         self.kappa = df_star.loc[i, 'rlnCurveKappa']
         self.edge_sigma = edge_sigma
         self.image_size = image.shape[0]
@@ -34,19 +27,7 @@ def compute_dist_mat(self):
         curve_generator = Curve(self.output_filename, self.i, self.gray_image, self.y0, self.x0, self.theta, self.kappa)
         self.distance_mat = curve_generator.compute()
         return self.distance_mat
-    # def generate_mem_mask(self):
-    #     self.distance_mat = self.compute_dist_mat()
-    #     for i in range(self.image_size):
-    #         for j in range(self.image_size):
-    #             if abs(self.distance_mat[i, j]) <= self.membrane_distance:
-    #                 self.membrane_mask[i, j] = 1
-    #             else:
-    #                 gray_value = np.exp(-(abs(self.distance_mat[i, j])-self.membrane_distance)**2 / (2 * self.edge_sigma**2))
-    #                 if gray_value < 0.001:
-    #                     self.membrane_mask[i, j] = 0
-    #                 else:
-    #                     self.membrane_mask[i, j] = gray_value
-    #     return self.membrane_mask
+
     def generate_mem_mask(self):
         print('>>>Start generating membrane mask...')
         distance_mat_gpu = self.compute_dist_mat()
@@ -59,10 +40,7 @@ def generate_mem_mask(self):
         membrane_mask_gpu[mask_outside_distance & ~mask_small_gray_value] = gray_value[mask_outside_distance & ~mask_small_gray_value]
         self.membrane_mask = membrane_mask_gpu.get()
         return self.membrane_mask
-    # def apply_mem_mask(self):
-    #     self.membrane_mask = self.generate_mem_mask()
-    #     self.masked_image = self.gray_image * self.membrane_mask
-    #     return self.masked_image
+
     def visualize_mem_mask(self):
         self.masked_image = self.apply_mem_mask()
         fig, ax = plt.subplots(1,3,figsize=(15,5))
@@ -72,14 +50,4 @@ def visualize_mem_mask(self):
         plt.show()
     def save_mem_mask(self, filename):
         self.membrane_mask = self.membrane_mask.astype('float32')
-        savemrc(self.membrane_mask, filename)
-
-
-if __name__ == '__main__':
-    image = readmrc('neuron_templates.mrc', section=3, mode='gpu')
-    image = zoom(image, 256/64)
-    mem_masker = mem_mask(image, edge_sigma=3)
-    mem_masker.generate_mem_mask()
-    # mem_masker.visualize_mem_mask()
-    # masked_mem = mem_masker.apply_mem_mask()
-    mem_masker.save_mem_mask('neuron_templates_mask_3.mrc')
+        savemrc(self.membrane_mask, filename)

build/lib/memxterminator/radonfit/lib/mem_average.py

@@ -85,9 +85,7 @@ def kappa_templates_generator(self, kappa_start=-0.01, kappa_end=0.01, kappa_num
             interpolated_values = f(distance_mat_temp_cpu)
             membrane_average_2d_gpu = cp.zeros((self.image_size, self.image_size))
             membrane_average_2d_gpu[mask_within_range] = cp.array(interpolated_values)
-            # membrane_average_2d_gpu[mask_within_range] = f(distance_mat_temp_gpu[mask_within_range])
             gray_value = cp.exp(-(cp.abs(distance_mat_temp_gpu) - self.membrane_distance)**2 / (2 * self.sigma**2))
-            # gray_value = cp.exp(-(distance_mat_temp_gpu - self.membrane_distance)**2 / (2 * self.sigma**2))
             update_mask = (cp.abs(distance_mat_temp_gpu) > self.membrane_distance) & (gray_value >= 0.001)
 
             membrane_average_2d_gpu[update_mask] *= gray_value[update_mask]

build/lib/memxterminator/radonfit/lib/radonanalyser.py

@@ -12,10 +12,6 @@ def __init__(self, output_filename, i, image, crop_rate, thr, theta_start=0, the
         self.gray_image = image
         self.theta_start = theta_start
         self.theta_end = theta_end
-        # kernel = gaussian_kernel(size=5, sigma=20)
-        # print('filter')
-        # self.gray_image = create_gaussian_low_pass_filter(self.gray_image, cutoff_frequency=20)
-        # self.gray_image = convolve(self.gray_image, kernel)
         self.image_size = self.gray_image.shape[0]
         self.crop_rate = crop_rate
         radius = (self.image_size/2) * self.crop_rate
@@ -62,19 +58,6 @@ def find_peaks_2d(self, image, neighborhood_size, threshold):
     def find_2_peaks(self):
         theta_start = self.theta_start
         theta_end = self.theta_end
-        # for dl in range(0,180):
-        #     theta, projection = self.radon_transform(theta_start+dl, theta_end+dl)
-        #     peaks = self.find_peaks_2d(cp.flipud(projection), 7, self.threshold * np.max(projection))
-        #     if len(peaks[0]) == 2:
-        #         if abs(peaks[1][0] - peaks[1][1]) < 1:
-        #             peaks[1][0] = peaks[1][0] + dl
-        #             peaks[1][1] = peaks[1][1] + dl
-        #             break
-        #         else:
-        #             pass
-        #     else:
-        #         continue
-        # return theta, projection, peaks
         theta, projection = self.radon_transform(theta_start, theta_end)
         peaks = self.find_peaks_2d(np.flipud(projection), 7, self.threshold * np.max(projection))
         print('peaks:', peaks)
@@ -141,9 +124,3 @@ def visualize_analyze(self):
         plt.tight_layout()
         plt.show()
 
-
-# if __name__ == '__main__':
-#     image = readmrc('J320/templates_selected.mrc', section=3, mode='gpu')
-#     analyzer = RadonAnalyzer(0, image, crop_rate=0.9, thr=0.8)
-#     analyzer.visualize_analyze()
-

build/lib/memxterminator/radonfit/lib/template_centerfitting.py

@@ -12,16 +12,10 @@ def __init__(self, output_filename, i, sigma1, sigma2, image, crop_rate, thr, th
         df_star = readstar(output_filename)
         self.output_filename = output_filename
         radonanalyze = RadonAnalyzer(output_filename, i, image, crop_rate=crop_rate, thr=thr, theta_start=theta_start, theta_end=theta_end)
-        # self.gray_image = self.radonanalyze.gray_image
         self.i = i
         self.gray_image = image
-        # self.image_size = self.radonanalyze.image_size
         self.image_size = self.gray_image.shape[0]
-        # self.peaks = self.radonanalyze.peaks
-        # self.theta = self.radonanalyze.average_theta
         self.theta = df_star.loc[self.i, 'rlnAngleTheta']
-        # print('theta:', self.theta)
-        # self.membrane_distance = self.radonanalyze.membrane_distance
         self.membrane_distance = df_star.loc[self.i, 'rlnMembraneDistance']
         self.template_size = template_size
         self.sigma1 = sigma1
@@ -53,29 +47,18 @@ def centerfinder(self):
             correlation_result = correlate2d(region, template, mode='full')
             corr_score = np.max(correlation_result)
             self.corr_lst.append(corr_score)
-        # self.corr_lst = cp.asarray(self.corr_lst)
         self.corr_max = max(self.corr_lst)
         self.corr_max_index = self.corr_lst.index(self.corr_max)
         self.yc = self.x0[self.corr_max_index]
         self.xc = self.y0[self.corr_max_index]
         print('membrane center y:', self.yc)
         print('membrane center x:', self.xc)
-        # print('corr_max_index:', self.corr_max_index)
 
         df_star = readstar(self.output_filename)
         df_star.loc[self.i, 'rlnCenterX'] = self.xc
         df_star.loc[self.i, 'rlnCenterY'] = self.yc
         write_star(df_star, self.output_filename)
         return self.xc, self.yc
-
-    def visualize_center(self):
-        fig, axes = plt.subplots(1, 4, figsize=(10, 5))
-        axes[0].plot(self.corr_lst, 'x', color = 'red')
-        axes[1].imshow(self.template, cmap='gray', origin='lower')
-        axes[2].imshow(self.get_region(self.corr_max_index), cmap='gray', origin='lower')
-        axes[3].imshow(self.gray_image, cmap='gray', origin='lower')
-        axes[3].plot(self.xc, self.yc, 'x', color = 'red')
-        plt.show()
 
     def fit_sigma(self):
         self.cc_region = self.get_region(self.corr_max_index)
@@ -87,11 +70,9 @@ def fit_sigma(self):
         print('>>>Start fitting sigma...')
         for i in corr_sigma:
             sigma_temp = i
-            # print('sigma_temp:', sigma_temp)
             template_n = generate_template(self.template_size, self.theta+90, self.membrane_distance, sigma_temp, self.sigma2, mode='gpu')
             template_n = (template_n - cp.mean(template_n)) / cp.std(template_n)
-            corr_result = correlate2d(self.cc_region, template_n, mode='full')
-            # print('corr_result:', np.nanmax(corr_result))
+            corr_result = correlate2d(self.cc_region, template_n, mode='same')
             self.corr_sigma_score1.append(cp.max(corr_result))
         corr_best_sigma_score = max(self.corr_sigma_score1)
         corr_best_sigma_index = self.corr_sigma_score1.index(corr_best_sigma_score)
@@ -101,53 +82,17 @@ def fit_sigma(self):
         corr_sigma = cp.linspace(0.1, self.sigma_range*self.sigma2, number2)
         for i in corr_sigma:
             sigma_temp = i
-            # print('sigma_temp:', sigma_temp)
             template_n_n = generate_template(self.template_size, self.theta+90, self.membrane_distance, best_sigma1, sigma_temp, mode='gpu')
             template_n_n = (template_n - cp.mean(template_n_n)) / cp.std(template_n_n)
-            corr_result = correlate2d(self.cc_region, template_n_n, mode='full')
-            # print('corr_result:', np.nanmax(corr_result))
+            corr_result = correlate2d(self.cc_region, template_n_n, mode='same')
             self.corr_sigma_score2.append(cp.max(corr_result))
         corr_best_sigma_score = max(self.corr_sigma_score2)
         corr_best_sigma_index = self.corr_sigma_score2.index(corr_best_sigma_score)
         best_sigma2 = corr_sigma[corr_best_sigma_index]
         print('best sigma2:', best_sigma2)        
-        # self.template = generate_template(self.template_size, self.theta+90, 12, best_sigma1, best_sigma2)
         df_star = readstar(self.output_filename)
         df_star.loc[self.i, 'rlnSigma1'] = best_sigma1
         df_star.loc[self.i, 'rlnSigma2'] = best_sigma2
         write_star(df_star, self.output_filename)
         return best_sigma1, best_sigma2
-
-    def fit_sigma_visualize(self):
-        fig, axes = plt.subplots(1, 4, figsize=(10, 5))
-        axes[0].imshow(self.cc_region, cmap='gray', origin='lower')
-        axes[1].imshow(self.template, cmap='gray', origin='lower')
-        axes[2].plot(self.corr_sigma_score1, 'x')
-        axes[3].plot(self.corr_sigma_score2, 'x')
-        plt.show()
-    # def refine_center(self):
-    #     self.centerfinder()
-    #     bestsigma1, bestsigma2 = self.fit_sigma()
-    #     self.template = generate_template(self.template_size, self.theta+90, self.membrane_distance, bestsigma1, bestsigma2)
-    #     self.centerfinder()
-    #     self.visualize_center()
-    #     return self.xc, self.yc
-    # def fit_2_sigma(self):
-    #     bestsigma1, bestsigma2 = self.fit_sigma()
-    #     return bestsigma1, bestsigma2
 
-if __name__ == '__main__':
-    image = readmrc('neuron_templates.mrc', section=3, mode='gpu')
-    image = zoom(image, 256/64)
-    analyzer = Template_centerfitting(6, 6, image, crop_rate=1, thr=0.4, template_size=80)
-    analyzer.centerfinder()
-    # analyzer.visualize_center()
-    analyzer.fit_sigma()
-    # analyzer.fit_sigma_visualize()
-    # # analyzer.refine_center()
-    # best_sigma1, best_sigma2 = analyzer.fit_sigma()
-    # print('best_sigma1:', best_sigma1)
-    # print('best_sigma2:', best_sigma2)
-    # fitted_analyzer = Template_centerfitting(best_sigma1, best_sigma2, image, thr=0.5)
-    # fitted_analyzer.centerfinder()
-    # fitted_analyzer.visualize_center()

mkdocs.yml

@@ -42,6 +42,7 @@ theme:
     - announce.dismiss
     - toc.follow
     - content.code.copy
+    - content.code.annotate
     - navigation.top
     - search.suggest
     - search.highlight
@@ -52,6 +53,9 @@ site_dir: ./wiki/site
 
 plugins:
   - search
+  - mkdocs-jupyter:
+        include_source: True
+        kernel_name: python3
 
 markdown_extensions:
   # Python Markdown
@@ -119,5 +123,6 @@ nav:
     - Micrograph Membrane Subtraction: ./tutorials/micrograph-membrane-subtraction.md
   - Frequently Asked Questions: ./tutorials/faq.md
   - Reference Pages:
+    - "Visualize Radonfit (.ipynb)": ./tutorials/reference/radonfit-mem-analysis-visualizer.ipynb
     - "Conventions": ./tutorials/reference/conventions.md
     - "API Reference": ./tutorials/reference/api.md