Fast Point Cloud Fusion of Incomplete SMLM Particles
Sabri Bolkar∗
Abstract
Single molecule localization microscopy (SMLM) allows fluorescent imaging of nano structures well
below the diffraction limit, but it suffers from low spatial resolution. Registration of identical particles offers a straightforward solution to increase the apparent resolution and SNR of localization
microscopy sets. In this technical report, we propose a fast and global joint registration framework for
SMLM particle averaging. We compare our results with the all-to-all registration framework proposedin [1].

Draft

The main motivation of this work is to develop a fast algorithm that will allow global registration of underlabeled 3d SMLM particles. Current state-of-the-art [1] is based on all-successive pairwise registration ofparticles that has quadratic computational complexity, hence it is not feasible to densely sample azimuthand elevation initialization angles in order to escape local optima points. The proposed joint registration
method is inspired by [2], but we put additional constraints on transformation parameters by prohibitingreflection during registration. We further enhance the idea by enabling multiple initializations that allow
us to reach the global maximizer. It is a work-in-progress with ongoing improvements that are summarizedat the end of the report.
In SMLM particle averaging, each particle can be treated as noisy and incomplete realizations of
the same underlying structure (i.e., binding sites). Furthermore, emission from a binding site generates
multiple localizations captured during image acquisition. It is possible to use the exactly same scenario
when modeling the registration algorithm. In the end, this will provide us a physically plausible particle
averaging framework with inherent localization clustering abilities.
Let assume a central Gaussian mixture model with K components G(x, µk , Σk ), k ∈ {1, K} having themeans µk and the covariance matrices Σk actually generated the observed set of M particles V : {vj }M
j=1 .
But these particles are incomplete realizations of the central GMM by rigid transformations (exceptreflection) T : {Rj , tj }M
j=1 . Now we can model probability distributions of each SMLM particle vj withtotal Nj localizations using the central GMM as following:
P (vj ) =Nj K
XXpk G((vji − tj )Rj−1 ; µk , Σk ).

(1)

i=1 k=1

where pk is the weight of each GMM component. This formulation exposes the unknown parametersKdenoted by θ that need to be estimated by our registration algorithm: θ = {{Rj , tj }M
j=1 , {µk , Σk , pk }k=1 }.
P This registration framework can be cast as a maximum likelihood problem where log likehood l(θ) =j ln(P (vj ; θ)) is maximized. But under-labeled characteristics of particles prevent direct estimationof parameters. In this work, expectation conditional maximization (ECM) [3] is used to iteratively
estimate each parameter by holding the others constant. This algorithm is a variant of famous expectation
maximization (EM) [4] that is known to handle missing data well. But, this would still not be sufficient
to obtain high quality particle averaging results. Due to lack of control at nano-level, SMLM particles
have rather arbitrary orientation and directly performing registration will be plagued by local optima. To
reach global maximizer, we also initialize each particle by rotating them with predefined angles uniformly
sampled from [0, 360], then choosing the parameters giving the best fit to the central GMM. In the below
sections, we first describe our algorithm then present the details about the implementation.

Theoretical Analysis
Let us explicitly write the log likelihood function l(θ) previously stated in the above paragraph:l(θ) =Nj
M XX

 X
X
Nj
M XK
ln P (vji ; θ) =ln
pk Gk (vji ; θ) .

j=1 i=1j=1 i=1∗ This work is carried out at TU Delft Quantitative Imaging Group

1k=1(2)

The main problem in (2) is that summation of Gaussian terms appears inside the logarithmic function
which in turn complicates analysis and prevents obtaining closed-form expressions for the parameters θ[5]. The trick applied in EM algorithm is to assume the existence of a hidden random variable Z : {zk }Kk=1 ,which is included in our formulation by defining P (V ) as the marginal distribution of the complete-data
distribution P (V, Z) (double summation is avoided for brevity):l(θ) =Xj,iXKln
P (vji , zh ; θ) .

(3)

h=1

Here, Z is an assignment variable representing which Gaussian component is responsible for which localizations. It takes values {0, 1} such that P (vji |z = k) = Gk (vji ; θ) and P (z = k) = pk (i.e., also denotedby P (zk ) for simplification) for corresponding kth assignment. Therefore, being aware of the joint distribution P (V, Z) simplifies the formulation by removing the summation appearing inside logarithmicfunction in (2). However, Z is actually not known to us a priori and we need additional considerations
to solve the problem.
At this point, we will exploit the Jensen’s inequality, f (E[x]) ≤ E[f (x)], which is defined for convex
functions f (x) : x ∈ R. Let us introduce a distribution, denoted by Q(Z), into our formulation in (3):XK
P (vji , zh ; θ),Q(zh )
lnQ(zh )j,iX

Draftl(θ) =(4)h=1

the expression in (4) can be expressed as expectation:
X  P (vji , zh ; θ) l(θ) =ln E() ,Q(zh )j,i(5)

now we can apply Jensen’s inequality to (5) and find a lower bound that can be used during maximization:X  P (vji , zh ; θ)  X
P (vji , zh ; θ)) ≤E(ln).
(6)
ln E(Q(zh )Q(zh )
j,ij,iAs can be seen in (6), now logarithmic function operates inside the summation. In a way, Jensen’sinequality allow us to carry out the optimization in a more straightforward manner. Further tightening
the inequality to an equality condition will result in quite convenient expression that we will use during
EM optimization. Equality holds if the expectation is performed over a constant valued function, thus
we can choose Q(Z) = PcP (V, Z; θ) where c is a constant real value. As both Q(z) and P (V, Z; θ) aredistributions, to satisfy h Q(zh ) = 1, Q(Z) can be chosen to be the normalized expression given by [6]:P (v, z; θ)Q(z) = P= P (z|v; θ).
z P (v, z; θ)(7)This expression gives the E-step of the EM algorithm where posterior distributions of Z given V
are computed. Then optimization is performed to find maximizers for each parameter in θ (M-step).
Considering Jensen’s inequality, here we actually optimize the lower bound instead of the originalformulation. It can be shown that both converges to the same local optima [5].
Now let us define the EM optimization as following:!
XP (vji , zh ; θ).
(8)
θ = argmaxθ
P (zh |vji ; θ)ln
P (zh |vji ; θ)i,j,hDue to the extended set of parameters, derivation of closed form expressions for all parameters θ =K{{Rj , tj }M
j=1 , {µk , Σk , pk }k=1 } is rather difficult using EM directly. Hence, as done in [2], we apply
Expectation Conditional Maximization algorithm where maximization steps of EM (i.e., M-steps) are
divided into smaller steps where closed form update expressions for each of the parameters obtained bytaking the rest of the parameters as they are. In the below section, we present details for each step in
our algorithmic implementation.
2

Algorithmic Analysis and Implementation
Theoretical analysis explained in the above paragraphs gives little intuition about how this framework
can actually be implemented. In this section, specific details related to the algorithm is provided. First,a pseudo-code of our iterative algorithm is shown below:

DraftK
Initialize θ = {{Rj , tj }M
j=1 , {µk , Σk , pk }k=1 };while iter< MaxIter (e.g., 100) do
if iter<FineIter (e.g., 20) then1) Coarse initialization: Rotate particles {Rja }M
j=1 using coarse sampled angles, a ∈ {1, 4};2) E step: Update P a (Z|V ; θ), ∀a;3) M-T step: Update {Rja , taj }, ∀a;4) Clamp step: When detected, clamp reflection in Rja by Rja ← I;
5) Transform particles using updated {Rja , taj }M
j=1 , ∀a;6) Max step: Choose the best Rj , tj and P (Z|V ) from {Rja , taj , P a (Z|V )} givingmax similarity(V, GM M ) ;7) M-GMM step: Update {µk , Σk };8) Continue with the particles transformed using the best parameters
else1)Fine initialization: Rotate particles {Rja }M
j=1 using fine sampled angles, a ∈ {1, 64};... (Repeat steps 2 to 8) ...;
end
end
Algorithm 1: Global joint SMLM particle averaging algorithm

The algorithm consist of alternating E and M steps of ECM algorithm until a predefined iteration isreached (i.e., 100). At the outset, we initialize all parameters. Rotation matrices are initialized by using
uniformly sampled angles from 0 to 360 degrees and Rja of jth particle is estimated at each initialization(a). Then, the one initialization with the highest Battacharyya similarity value with respect to the centralGMM is kept to be used for the next steps. Battacharyya similarity is defined exactly as in [1]. Here, theterm coarse to fine used to indicate that we actually increase number of initializations from 4 to 64 after
a predefined iteration passes (e.g., 20th iteration) to exploit roughly aligned particles for subsequent fineregistration. Central Gaussian components are initialized from one of the particles. The chosen particle isdown-sampled to ensure that the number of GMM components are lower than the number of localizations(i..e, 50-100 GMM means preferred). This scheme encourages clustering of localizations around central
GMM components.
E StepIn the E step, posterior distributions of Z given particles V and GMM centers, P (Z|V ; θ) are determined.
Briefly, it is a soft assignment step where kth GMM component’s probability of being assigned to thelocalization vji is estimated:Gjikajik = P (zji = k|vji ; θ) = P,k Gjik
where isotropic Guassian term Gjik with Σk = σk I is defined as:pk−||Rj vji + tj − µk ||2Gjik =exp.
σk2σk2(9)(10)

Note that we assumed spherical (i.e. isotropic) Gaussian mixture components in the implementation.
Previous works indicate no sign of significant improvement when using non-isotropic GMM mixtures [2].
M-T Step
In this step, transformation parameters for each particle vji are estimated fixing current ajik values and3

GMM components. This section simply formulates the problem with least-squares minimization andupdates parameters using the work of Umeyana [7]. The {Rj , tj }M
j=1 updates are computed given the soft
assignments of each localization to the Gaussian components (ajik ). The minimization can be stated asfollowing:minRj ,tj ||(Wj + tj eT − M )Λ||2F
RT R = I, det(Rj ) = 1
P
1where Λj ∈ IRKxK is a diagonal matrix with entries as λjkk = ( i ajik /σk ) 2 , M indicates vectorized
GMM centers M = [µ1 , .., µK ] ∈ IRdxK (d: dimension) and Wj = [wj1 , .., wjK ] ∈ IRdxK such that eachwjk defined as:
Pi ajik vji
wjk = P.
(11)
i ajikThis minimization problem has closed form solutions for both rotation and translation using singularvalue decomposition (see [7] for details).

Draft

Clamp Step
It should be noted that M-T step can be improved by simply using a more modern technique for estimation of point-to-point transformation parameters. The main reason of this claim is that, for under-labeled
particles, Umeyana’s method fails to correct for reflection. Therefore, we include a reflection clamping
step when estimating rotation parameters. In this step, reflection is checked and equalized to identity
matrix when detected.
Max Step
In this step, the parameters belonging to the highest Bhattacharyya similarity computed between thetransformed particle and the central GMM are preserved for the next iteration. This operation is performed for all particles in the dataset.
M-GMM Step
As the last step in the ECM algorithm, GMM parameters {µk , σk }K
k=1 are estimated by from the besttransformation parameters and posterior distributions. The expression for the local optima points are
simply obtained by equalizing the respective derivatives to zero. It can be shown that actually the GMM
centers maximizing log likelihood are actually weighted average of the transformed particles:
Pji ajik (Rj vji + tj )P
µk =(12)ji ajikP2
ji ajik ||Rj vji + tj − µk ||2P
σk =(13)ji ajikIn the implementation, we ignored updates for weights (a.k.a priors) pk and fixed them to the uniform
distribution pk = 1/k. Previous works show no significant improvement for updates applied to the GMMweights.
Computational Complexity
Our joint registration framework has computational complexity of O(M N K) where N is the average
number of localizations per particle, M is the number of particles, and K is the number of GMM
components. On the other hand, all-to-all framework has complexity of O((M N )2 ), which is order ofmagnitude slower comparing to the proposed algorithm.

Results
Particle averaging results for simulation and experimental datasets are shown in Figure 2 and Figure 3
respectively. We compare our reconstruction outcomes against the all-to-all framework [1].
4

50 Dol Simulation

30 Dol Simulation

After bootstrapping

Main algorithm

After bootstrapping

Draft

All-to-all

Global Joint

Main algorithm

Figure 1: Qualitative comparison of Global Joint and All-to-all registration for simulation set.
50 Dol Experimental

30 Dol Experimental

After bootstrapping

Main algorithm

After bootstrappingAll-to-all

Global Joint

Main algorithm

Figure 2: Qualitative comparison of Global Joint and All-to-all registration for experimental set.

Future Work
Global joint SMLM particle averaging is a work-in-progress with many aspects to include into the framework:• Include individual localization uncertainties. At the moment, algorithm assumes fixed uncertainties
for all localizations.
• Implement a better formulation for the M-T step instead of Umeya’s solution developed in 1991.
• Improve global optimization framework by using Bayesian optimization rather then uniform sampling of rotation space.

5• Implement a consistency check block such that transformed particles with less similar to the central
GMM are re-processed at later iterations.
• Enable paralleled GPU computation when estimating transformation parameters as the algorithm
steps are independent for each particle (except M-GMM step).

ReferencesDraft[1] Hamidreza Heydarian, Florian Schueder, Maximilian T Strauss, Ben Van Werkhoven, MohamadrezaFazel, Keith A Lidke, Ralf Jungmann, Sjoerd Stallinga, and Bernd Rieger, “Template-free 2d particle
fusion in localization microscopy,” Nature methods, vol. 15, no. 10, pp. 781, 2018.
[2] Georgios Dimitrios Evangelidis and Radu Horaud, “Joint alignment of multiple point sets with batchand incremental expectation-maximization,” IEEE transactions on pattern analysis and machine
intelligence, vol. 40, no. 6, pp. 1397–1410, 2018.
[3] Xiao-Li Meng and Donald B Rubin, “Maximum likelihood estimation via the ecm algorithm: Ageneral framework,” Biometrika, vol. 80, no. 2, pp. 267–278, 1993.
[4] Arthur P Dempster, Nan M Laird, and Donald B Rubin, “Maximum likelihood from incomplete data
via the em algorithm,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 39,no. 1, pp. 1–22, 1977.
[5] Christopher M Bishop, Pattern recognition and machine learning, springer, 2006.
[6] Andrew Ng, “Cs229 lecture notes,” CS229 Lecture notes, vol. 1, no. 1, pp. 1–3, 2000.
[7] Shinji Umeyama, “Least-squares estimation of transformation parameters between two point patterns,” IEEE Transactions on Pattern Analysis & Machine Intelligence, , no. 4, pp. 376–380, 1991.

6