Remove more "one file one function" stuff #8

ManifoldEM/CC/LoadPrDPsiMoviesMasked.py

@@ -3,7 +3,8 @@
 
 import numpy as np
 
-from ManifoldEM import annularMask, myio, p, projectMask
+from ManifoldEM import myio, p, projectMask
+from ManifoldEM.core import annularMask
 '''		
 Copyright (c) Columbia University Suvrajit Maji 2020		
 Modified:Sept 17,2021
@@ -23,7 +24,7 @@ def getMask2D(prD, maskType, radius):
             N2 = radius  # also includes radius = 0
         if prD == 0:
             print('Annular mask radius: {} pixels'.format(N2))
-        mask = annularMask.op(0, N2, N, N)
+        mask = annularMask(0, N2, N, N)
 
     elif maskType == 'volumetric':  #3d volume mask from user-input
         dist_file = p.get_dist_file(prD)

ManifoldEM/NLSA.py

@@ -5,8 +5,8 @@
 from scipy.fftpack import fft2
 from scipy.fftpack import ifft2
 
-from ManifoldEM import DMembeddingII, myio, get_wiener
-from ManifoldEM.core import L2_distance, svdRF
+from ManifoldEM import DMembeddingII, myio
+from ManifoldEM.core import L2_distance, svdRF, get_wiener
 from ManifoldEM.fit_1D_open_manifold_3D import fit_1D_open_manifold_3D
 '''
 Copyright (c) UWM, Ali Dashti 2016 (original matlab version)
@@ -48,7 +48,7 @@ def op(NLSAPar, DD, posPath, posPsi1, imgAll, msk2, CTF, ExtPar):  #pass the msk
     if 'prD' in ExtPar:
         IMG1 = imgAll[posPath[posPsi1], :, :]
         # Wiener filtering
-        wiener_dom, CTF1 = get_wiener.op(CTF, posPath, posPsi1, ConOrder, num)
+        wiener_dom, CTF1 = get_wiener(CTF, posPath, posPsi1, ConOrder, num)
     elif 'cuti' in ExtPar:
         IMG1 = imgAll[posPsi1, :, :]
 

ManifoldEM/PrepareOutputS2.py

@@ -5,18 +5,14 @@
 import numpy as np
 
 from ManifoldEM import writeRelionS2, p, myio, set_params
-''' %Version V 1.2
-    % Copyright (c) UWM, Ali Dashti 2016 (matlab version)
-    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    %This script prepares the image stacks and orientations for 3D reconstruction.
-    Copyright (c) Columbia Univ Hstau Liao 2018 (python version)
-    Copyright (c) Columbia University Suvrajit Maji 2020 (python version)
-'''
-
 
 def op(*argv):
+    """This script prepares the image stacks and orientations for 3D reconstruction."""
+    # Copyright (c) UWM, Ali Dashti 2016 (matlab version)
+    # Copyright (c) Columbia Univ Hstau Liao 2018 (python version)
+    # Copyright (c) Columbia University Suvrajit Maji 2020 (python version)
+
     set_params.op(1)
-    #ComputeEnergy1D.op() # this is s a repeat , commented out
     print("Writing output files...")
 
     data = myio.fin1(p.CC_file)

ManifoldEM/S2tessellation.py

@@ -5,7 +5,7 @@
 from mpl_toolkits.mplot3d import Axes3D
 from sklearn.neighbors import NearestNeighbors
 
-from ManifoldEM import distribute3Sphere
+from ManifoldEM.core import distribute3Sphere
 '''
 Copyright (c) UWM, Ali Dashti 2016 (original matlab version)
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -42,7 +42,7 @@ def _classS2(X, Q):
 def op(q, shAngWidth, PDsizeTh, visual, thres, *fig):
     nG = int(4 * np.pi / (shAngWidth**2))
     # reference angles
-    S20, it = distribute3Sphere.op(nG)
+    S20, it = distribute3Sphere(nG)
 
     S20 = S20.T
     # projection angles

ManifoldEM/core.py

@@ -138,7 +138,7 @@ def find_thres(logEps, D2):
         eps = np.exp(logEps[k])
         d = 1. / (2. * eps) * D2
         d = -d[d < thr]
-        Wij = np.exp(d)  #see Coifman 2008
+        Wij = np.exp(d)  # see Coifman 2008
         logSumWij[k] = np.log(np.sum(Wij))
 
     # curve fitting of a tanh():
@@ -151,8 +151,99 @@ def find_thres(logEps, D2):
         a0 = 1 * (np.random.rand(4, 1) - .5)
 
         residuals = logSumWij - fun(logEps, popt[0], popt[1], popt[2], popt[3])
-        ss_res = np.sum(residuals**2)  #residual sum of squares
-        ss_tot = np.sum((logSumWij - np.mean(logSumWij))**2)  #total sum of squares
-        R_squared = 1 - (ss_res / ss_tot)  #R**2-value
+        ss_res = np.sum(residuals**2)  # residual sum of squares
+        ss_tot = np.sum((logSumWij - np.mean(logSumWij))**2)  # total sum of squares
+        R_squared = 1 - (ss_res / ss_tot)  # R**2-value
 
     return (popt, logSumWij, resnorm, R_squared)
+
+
+def annularMask(a: float, b: float, N: int, M: int):
+    """
+    returns a N x M matrix with an annular (donut) mask of inner
+    radius a and outer radius b. Pixels outside the donut or inside the hole
+    are filled with 0, and the rest with 1.
+
+    The circles with radii a and b are centered on pixel (N/2,M/2).
+
+    Programmed December 2007, modified by Peter Schwander December 2008 (Python version by Hstau Liao 2018)
+    Copyright (c) Russell Fung 2007
+    """
+    aSq = a * a
+    bSq = b * b
+    mask = np.zeros((N, M))
+    for xx in range(N):
+        xDist = xx - N / 2 + 1
+        xDistSq = xDist * xDist
+        for yy in range(M):
+            yDist = yy - M / 2
+            yDistSq = yDist * yDist
+            rSq = xDistSq + yDistSq
+            mask[xx, yy] = (rSq >= aSq) & (rSq < bSq)
+
+    return mask
+
+
+def distribute3Sphere(numPts: int):
+    """
+    distributes numPts points roughly uniformly on a unit 3-sphere and
+    returns the coordinates in results. Number of iterations required is
+    returned in iter.
+
+    Algorithm adapted from L. Lovisolo and E.A.B. da Silva, Uniform
+    distribution of points on a hyper-sphere with applications to vector
+    bit-plane encoding, IEE Proc.-Vis. Image Signal Process., Vol. 148, No.
+    3, June 2001
+
+    Programmed February 2009
+    Copyright (c) Russell Fung 2009
+    Copyright (c) Columbia University Hstau Liao 2018 (python version)
+    """
+    maxIter = 100
+    K = numPts
+    A3 = 4 * np.pi  # surface area of a unit 3-sphere
+    delta = np.exp(np.log(A3 / K) / 2.)
+    results = np.zeros((2 * K, 3))
+    # algorithm sometimes returns more/ less points
+    it = 0
+    id = 0
+
+    while id != K and it < maxIter:
+        it = it + 1
+        id = 0
+        dw1 = delta
+        for w1 in np.arange(0.5 * dw1, np.pi, dw1):
+            cw1 = np.cos(w1)
+            sw1 = np.sin(w1)
+            x1 = cw1
+            dw2 = dw1 / sw1
+            for w2 in np.arange(0.5 * dw2, 2 * np.pi, dw2):
+                cw2 = np.cos(w2)
+                sw2 = np.sin(w2)
+                x2 = sw1 * cw2
+                x3 = sw1 * sw2
+
+                results[id, :] = np.hstack((x1, x2, x3))
+                id = id + 1
+
+        delta = delta * np.exp(np.log(float(id) / K) / 2.)
+
+    results = results[0:K, :]
+    return (results, it)
+
+
+def get_wiener(CTF, posPath, posPsi1, ConOrder, num):
+    # Copyright (c) UWM, Ali Dashti 2016 (original matlab version)
+    # Copyright (c) Columbia University Hstau Liao 2018 (python version)
+    dim = CTF.shape[1]
+    SNR = 5
+    CTF1 = CTF[posPath[posPsi1], :, :]
+    wiener_dom = np.zeros((num - ConOrder, dim, dim), dtype='float64')
+    for i in range(num - ConOrder):
+        for ii in range(ConOrder):
+            ind_CTF = ConOrder - ii + i
+            wiener_dom[i, :, :] = wiener_dom[i, :, :] + CTF1[ind_CTF, :, :]**2
+
+    wiener_dom = wiener_dom + 1. / SNR
+
+    return (wiener_dom, CTF1)

ManifoldEM/getDistanceCTF_local_Conj9combinedS2.py

@@ -57,7 +57,8 @@
 from scipy.fftpack import ifftshift, fft2, ifft2
 import matplotlib.pyplot as plt
 
-from ManifoldEM import ctemh_cryoFrank, myio, p, annularMask, projectMask
+from ManifoldEM import ctemh_cryoFrank, myio, p, projectMask
+from ManifoldEM.core import annularMask
 from ManifoldEM.quaternion import q2Spider
 from ManifoldEM.util import rotate_fill
 
@@ -241,7 +242,7 @@ def op(input_data, filterPar, imgFileName, sh, nStot, options):
     CTF = np.zeros((nS, N, N))  # each (i,:,:) is the CTF
     D = np.zeros((nS, nS))  # distances among the particles in the bin
 
-    msk = annularMask.op(0, N2, N, N)
+    msk = annularMask(0, N2, N, N)
 
     # read images with conjugates
     imgLabels = np.zeros(nS, dtype=int)

ManifoldEM/psiAnalysis.py

@@ -72,7 +72,7 @@ def op(*argv):
         offset = np.count_nonzero(fin_PDs == 1)
         progress3.emit(int((offset / float((p.numberofJobs) * p.num_psis)) * 100))
 
-    print("Processing {} projection directions.".format(len(input_data)))
+    print(f"Processing {len(input_data)} projection directions.")
 
     if p.ncpu == 1:  # avoids the multiprocessing package
         for i in range(len(input_data)):