--- a/python/cvutils.py	Wed Jul 12 18:00:53 2017 -0400
+++ b/python/cvutils.py	Thu Jul 13 00:31:54 2017 -0400
@@ -22,6 +22,7 @@
 from math import floor, log10, ceil
 
 from numpy import dot, array, append, float32, loadtxt, savetxt, append, zeros, ones, identity, abs as npabs, logical_and, unravel_index, sum as npsum, isnan, mgrid, median, floor as npfloor, ceil as npceil
+from numpy.linalg import inv
 from matplotlib.mlab import find
 from matplotlib.pyplot import imread, imsave
 
@@ -518,7 +519,7 @@
             out.write('{0} '.format(cvmat[i,j]))
         out.write('\n')
 
-def projectArray(homography, points, intrinsicCameraMatrix = None, distortionCoefficients = None):
+def projectArray(homography, points, intrinsicCameraMatrix = None, distortionCoefficients = None, newCameraMatrix = None):
     '''Returns the coordinates of the projected points through homography
     (format: array 2xN points)'''
     if points.shape[0] != 2:
@@ -532,10 +533,14 @@
         projected = augmentedPoints
 
     if intrinsicCameraMatrix is not None and distortionCoefficients is not None:
-        #projected[2,:] = 0
-        projected, jacobian = cv2.projectPoints(projected.T, (0,0,0), (0,0,0), intrinsicCameraMatrix, distortionCoefficients) # in: 3xN, out: 2x1xN
+        if newCameraMatrix is not None:
+            invNewCameraMatrix = inv(newCameraMatrix)
+            reducedPoints = dot(invNewCameraMatrix, projected)
+        else:
+            reducedPoints = projected
+        projected, jacobian = cv2.projectPoints(reducedPoints.T, (0.,0.,0.), (0.,0.,0.), intrinsicCameraMatrix, distortionCoefficients) # in: 3xN, out: 2x1xN
         projected = projected.reshape(-1,2).T
-    return projected[:2,:]
+        return projected[:2,:]
 
 def project(homography, p, intrinsicCameraMatrix = None, distortionCoefficients = None):
     '''Returns the coordinates of the projection of the point p with coordinates p[0], p[1]
@@ -551,7 +556,6 @@
 def invertHomography(homography):
     '''Returns an inverted homography
     Unnecessary for reprojection over camera image'''
-    from numpy.linalg import inv
     invH = inv(homography)
     invH /= invH[2,2]
     return invH

--- a/python/tests/cvutils.txt	Wed Jul 12 18:00:53 2017 -0400
+++ b/python/tests/cvutils.txt	Thu Jul 13 00:31:54 2017 -0400
@@ -1,11 +1,10 @@
 >>> import cv2, cvutils
->>> from numpy import array, round, ones
+>>> from numpy import array, round, ones, dot, linalg, absolute
 >>> img = cv2.imread("../samples/val-dor-117-111.png")
 >>> width = img.shape[1]
 >>> height = img.shape[0]
 >>> intrinsicCameraMatrix = array([[ 377.42,    0.  ,  639.12], [   0.  ,  378.43,  490.2 ], [   0.  ,    0.  ,    1.  ]])
->>> distortionCoefficients = array([-0.11759321, 0.0148536, 0.00030756, -0.00020578, -0.00091816])
->>> distortionCoefficients = array([-0.11759321, 0., 0., 0., 0.])
+>>> distortionCoefficients = array([-0.11759321, 0.0148536, 0.00030756, -0.00020578, -0.00091816])# distortionCoefficients = array([-0.11759321, 0., 0., 0., 0.])
 >>> multiplicationFactor = 1.31
 >>> [map1, map2] = cvutils.computeUndistortMaps(width, height, multiplicationFactor, intrinsicCameraMatrix, distortionCoefficients)
 >>> undistorted = cv2.remap(img, map1, map2, interpolation=cv2.INTER_LINEAR)
@@ -14,6 +13,17 @@
 >>> imgPoints = array([[[150.,170.],[220.,340.],[340.,440.],[401.,521.]]])
 >>> newCameraMatrix = cv2.getDefaultNewCameraMatrix(intrinsicCameraMatrix, (int(round(width*multiplicationFactor)), int(round(height*multiplicationFactor))), True)
 >>> undistortedPoints = cv2.undistortPoints(imgPoints, intrinsicCameraMatrix, distortionCoefficients, P = newCameraMatrix).reshape(-1, 2)
+>>> invNewCameraMatrix = linalg.inv(newCameraMatrix)
 >>> tmp = ones((imgPoints[0].shape[0], 3))
 >>> tmp[:,:2] = undistortedPoints
->>> origPoints = cv2.projectPoints(tmp, (0.,0.,0.), (0.,0.,0.), intrinsicCameraMatrix, distortionCoefficients)[0].reshape(-1,2)
+>>> reducedPoints = dot(invNewCameraMatrix, tmp.T).T
+>>> origPoints = cv2.projectPoints(reducedPoints, (0.,0.,0.), (0.,0.,0.), intrinsicCameraMatrix, distortionCoefficients)[0].reshape(-1,2)
+>>> (round(origPoints[1:,:]) == imgPoints[0][1:,:]).all()
+True
+>>> (absolute(origPoints[0,:]-imgPoints[0][0,:])).max() < 6.
+True
+>>> origPoints = cvutils.projectArray(None, undistortedPoints.T, intrinsicCameraMatrix, distortionCoefficients, newCameraMatrix).T
+>>> (round(origPoints[1:,:]) == imgPoints[0][1:,:]).all()
+True
+>>> (absolute(origPoints[0,:]-imgPoints[0][0,:])).max() < 6.
+True