diff --git a/ManifoldEM/CC/LoadPrDPsiMoviesMasked.py b/ManifoldEM/CC/LoadPrDPsiMoviesMasked.py index 8368093..73a3f45 100644 --- a/ManifoldEM/CC/LoadPrDPsiMoviesMasked.py +++ b/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) diff --git a/ManifoldEM/NLSA.py b/ManifoldEM/NLSA.py index f42b41a..5a102aa 100644 --- a/ManifoldEM/NLSA.py +++ b/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, :, :] diff --git a/ManifoldEM/PrepareOutputS2.py b/ManifoldEM/PrepareOutputS2.py index e12eb3f..4010518 100644 --- a/ManifoldEM/PrepareOutputS2.py +++ b/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) diff --git a/ManifoldEM/S2tessellation.py b/ManifoldEM/S2tessellation.py index 94756cd..774378b 100644 --- a/ManifoldEM/S2tessellation.py +++ b/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 diff --git a/ManifoldEM/annularMask.py b/ManifoldEM/annularMask.py deleted file mode 100644 index ba6b8eb..0000000 --- a/ManifoldEM/annularMask.py +++ /dev/null @@ -1,32 +0,0 @@ -''' -Copyright (c) Columbia University Hstau Liao 2019 -''' -"""function mask = annularMask(a,b,N,M) -% mask = annularMask(a,b,N,M) -% -% 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 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -""" -import numpy as np - -def op(a, b, N, M): - 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 diff --git a/ManifoldEM/core.py b/ManifoldEM/core.py index 448fb88..363628e 100644 --- a/ManifoldEM/core.py +++ b/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) diff --git a/ManifoldEM/distribute3Sphere.py b/ManifoldEM/distribute3Sphere.py deleted file mode 100644 index c983b56..0000000 --- a/ManifoldEM/distribute3Sphere.py +++ /dev/null @@ -1,54 +0,0 @@ -""" function [results, iter] = distribute3Sphere(numPts) -% [results, iter] = distribute3Sphere(numPts) -% 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) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -""" - -import numpy as np -import math - - -def op(numPts): - maxIter = 100 - K = numPts - A3 = 4 * np.pi # surface area of a unit 3-sphere - delta = math.exp(math.log(A3 / K) / 2.) - results = np.zeros((2 * K, 3)) - # algorithm sometimes returns more/ less points - iter = 0 - id = 0 - - while id != K and iter < maxIter: - iter = iter + 1 - id = 0 - dw1 = delta - for w1 in np.arange(0.5 * dw1, np.pi, dw1): - cw1 = math.cos(w1) - sw1 = math.sin(w1) - x1 = cw1 - dw2 = dw1 / sw1 - for w2 in np.arange(0.5 * dw2, 2 * np.pi, dw2): - cw2 = math.cos(w2) - sw2 = math.sin(w2) - x2 = sw1 * cw2 - x3 = sw1 * sw2 - - results[id, :] = np.hstack((x1, x2, x3)) - id = id + 1 - - delta = delta * math.exp(math.log(float(id) / K) / 2.) - - results = results[0:K, :] - return (results, iter) diff --git a/ManifoldEM/getDistanceCTF_local_Conj9combinedS2.py b/ManifoldEM/getDistanceCTF_local_Conj9combinedS2.py index ad17b7a..7b739d4 100644 --- a/ManifoldEM/getDistanceCTF_local_Conj9combinedS2.py +++ b/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) diff --git a/ManifoldEM/get_wiener.py b/ManifoldEM/get_wiener.py deleted file mode 100644 index 09c4c2b..0000000 --- a/ManifoldEM/get_wiener.py +++ /dev/null @@ -1,22 +0,0 @@ -''' -Copyright (c) UWM, Ali Dashti 2016 (original matlab version) -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Copyright (c) Columbia University Hstau Liao 2018 (python version) -''' - -import numpy as np - - -def op(CTF, posPath, posPsi1, ConOrder, num): - 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) diff --git a/ManifoldEM/psiAnalysis.py b/ManifoldEM/psiAnalysis.py index 75437de..0f5fe04 100644 --- a/ManifoldEM/psiAnalysis.py +++ b/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)):