now avoid crashing when invoked on empty arrays (length 0)

setup.py

@@ -38,9 +38,9 @@ def Superpose3D(X,    # <-- Nx3 array of coords for the "frozen" point cloud
 ```
 ...where:
 ```
-    R  = a rotation matrix    (a 3x3 numpy array whose determinant = 1),
-    T  = a translation vector (a 1-D numpy array containing x,y,z displacements),
-    c  = a scalar             (a number, 1 by default)
+    R = a rotation matrix    (a 3x3 numpy array representing the rotation. |R|=1)
+    T = a translation vector (a 1-D numpy array containing x,y,z displacements)
+    c = a scalar             (a number, 1 by default)
 ```
 This function returns a 4-tuple containing the optimal values of:
 ```
@@ -76,9 +76,9 @@ def Superpose3D(X,    # <-- Nx3 array of coords for the "frozen" point cloud
 
   url='https://github.com/jewettaij/superpose3d',
 
-  download_url='https://github.com/jewettaij/superpose3d/archive/v1.0.0.zip',
+  download_url='https://github.com/jewettaij/superpose3d/archive/v1.0.1.zip',
 
-  version='1.0.0',
+  version='1.0.1',
 
   install_requires=[
       'numpy',

superpose3d/__init__.py

@@ -45,14 +45,15 @@ def Superpose3D(aaXf_orig,   # <-- coordinates for the "frozen" object
             aCenter_f[d] += aaXf_orig[n][d]*aWeights[n]
             aCenter_m[d] += aaXm_orig[n][d]*aWeights[n]
         sum_weights += aWeights[n]
-    for d in range(0, 3):
-        aCenter_f[d] /= sum_weights
-        aCenter_m[d] /= sum_weights
+    if sum_weights != 0:
+        for d in range(0, 3):
+            aCenter_f[d] /= sum_weights
+            aCenter_m[d] /= sum_weights
 
     # Subtract the centers-of-mass from the original coordinates for each object
     aaXf = np.empty((N,3))
     aaXm = np.empty((N,3))
-    aaXf[0][0] = 0.0
+
     for n in range(0, N):
         for d in range(0, 3):
             aaXf[n][d] = aaXf_orig[n][d] - aCenter_f[d]
@@ -97,10 +98,10 @@ def Superpose3D(aaXf_orig,   # <-- coordinates for the "frozen" object
     P[3][3] = 0.0
 
     p = np.zeros(4)
-    # default values (in case the caller supplies nonsensical because N<2)
+    # default values for p and pPp
     p[3] = 1.0
     pPp = 0.0
-    singular = (N < 2)
+    singular = (N < 2);  # (it doesn't make sense to rotate a single point)
 
     try:
         #http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.eigh.html
@@ -109,8 +110,8 @@ def Superpose3D(aaXf_orig,   # <-- coordinates for the "frozen" object
     except LinAlgError:
         singular = True  # (I have never seen this happen.)
 
-    if (not singular):
-
+
+    if (not singular):  # (don't crash if the caller supplies nonsensical input)
         eval_max = aEigenvals[0]
         i_eval_max = 0
         for i in range(1,4):
@@ -166,7 +167,10 @@ def Superpose3D(aaXf_orig,   # <-- coordinates for the "frozen" object
     sum_sqr_dist = E0 - c*2.0*pPp
     if sum_sqr_dist < 0.0: #(edge case due to rounding error)
         sum_sqr_dist = 0.0
-    rmsd = sqrt(sum_sqr_dist/sum_weights)
+
+    rmsd = 0.0
+    if sum_weights != 0.0:
+        rmsd = sqrt(sum_sqr_dist/sum_weights)
 
     # Lastly, calculate the translational offset:
     # Recall that:

test_superpose3d.py

@@ -21,14 +21,16 @@ def test_superpose3d():
     # now test weights, rescale, and quaternions
     w = [1.0, 1.0, 1.0, 1.0]
     result = Superpose3D(X, xscshift, w, True)
-
     # Does the RMSD returned in result[0] match the RMSD calculated manually?
     R = np.matrix(result[1])              # rotation matrix
     T = np.matrix(result[2]).transpose()  # translation vector (3x1 matrix)
     c = result[3]                         # scalar
-    _x = np.matrix(xscshift).transpose()
-    _xprime = c*R*_x + T
-    xprime = np.array(_xprime.transpose()) # convert to length 3 numpy array
+    if len(X) > 0:
+        _x = np.matrix(xscshift).transpose()
+        _xprime = c*R*_x + T
+        xprime = np.array(_xprime.transpose()) # convert to length 3 numpy array
+    else:
+        xprime = np.array([])
 
     print('1st (frozen) point cloud:\n'+str(X))
     print('2nd (mobile) point cloud:\n'+str(xscshift))
@@ -46,7 +48,9 @@ def test_superpose3d():
                  (X[i][1] - xprime[i][1])**2 +
                  (X[i][2] - xprime[i][2])**2)
 
-    RMSD = sqrt(RMSD / len(X))
+    if len(X) > 0:
+        RMSD = sqrt(RMSD / len(X))
+
     assert(abs(RMSD - result[0]) < 1.0e-6)