Index: modules/imgproc/include/opencv2/imgproc/imgproc.hpp
===================================================================
--- modules/imgproc/include/opencv2/imgproc/imgproc.hpp	(revision 8691)
+++ modules/imgproc/include/opencv2/imgproc/imgproc.hpp	(working copy)
@@ -600,7 +600,7 @@
 //! computes PSNR image/video quality metric
 CV_EXPORTS_W double PSNR(InputArray src1, InputArray src2);
 
-CV_EXPORTS_W Point2d phaseCorrelate(InputArray src1, InputArray src2, InputArray window = noArray());
+CV_EXPORTS_W Point2d phaseCorrelate(InputArray src1, InputArray src2, InputArray window = noArray(), double* response=0);
 CV_EXPORTS_W void createHanningWindow(OutputArray dst, Size winSize, int type);
 
 //! type of the threshold operation
Index: modules/imgproc/doc/motion_analysis_and_object_tracking.rst
===================================================================
--- modules/imgproc/doc/motion_analysis_and_object_tracking.rst	(revision 8691)
+++ modules/imgproc/doc/motion_analysis_and_object_tracking.rst	(working copy)
@@ -151,11 +151,12 @@
 
 Calculates the cross-power spectrum of two supplied source arrays. The arrays are padded if needed with :ocv:func:`getOptimalDFTSize`.
 
-.. ocv:function:: Point2d phaseCorrelate(InputArray src1, InputArray src2, InputArray window = noArray())
+.. ocv:function:: Point2d phaseCorrelate(InputArray src1, InputArray src2, InputArray window = noArray(), double* response = 0)
 
     :param src1: Source floating point array (CV_32FC1 or CV_64FC1)
     :param src2: Source floating point array (CV_32FC1 or CV_64FC1)
     :param window: Floating point array with windowing coefficients to reduce edge effects (optional).
+    :param response: Signal power within the 5x5 centroid around the peak, between 0 and 1 (optional).
 
 Return value: detected phase shift (sub-pixel) between the two arrays.
 
@@ -165,31 +166,33 @@
 
 * Next it computes the forward DFTs of each source array:
 
-     .. math::
+.. math::
 
         \mathbf{G}_a = \mathcal{F}\{src_1\}, \; \mathbf{G}_b = \mathcal{F}\{src_2\}
 
   where
-  :math:`\mathcal{F}` is the forward DFT.
+..  math::`\mathcal{F}` is the forward DFT.
 
 * It then computes the cross-power spectrum of each frequency domain array:
 
-    .. math::
+.. math::
 
           R = \frac{ \mathbf{G}_a \mathbf{G}_b^*}{|\mathbf{G}_a \mathbf{G}_b^*|}
 
 * Next the cross-correlation is converted back into the time domain via the inverse DFT:
 
-    .. math::
+.. math::
 
           r = \mathcal{F}^{-1}\{R\}
 
 * Finally, it computes the peak location and computes a 5x5 weighted centroid around the peak to achieve sub-pixel accuracy.
 
-    .. math::
+.. math::
 
          (\Delta x, \Delta y) = \texttt{weightedCentroid} \{\arg \max_{(x, y)}\{r\}\}
 
+* If non-zero, the response parameter is computed as the sum of the elements of r within the 5x5 centroid around the peak location. It is normalized to a maximum of 1 (meaning there is a single peak) and will be smaller when there are multiple peaks.
+
 .. seealso::
     :ocv:func:`dft`,
     :ocv:func:`getOptimalDFTSize`,
Index: modules/imgproc/src/phasecorr.cpp
===================================================================
--- modules/imgproc/src/phasecorr.cpp	(revision 8691)
+++ modules/imgproc/src/phasecorr.cpp	(working copy)
@@ -406,7 +406,7 @@
     merge(planes, out);
 }
 
-static Point2d weightedCentroid(InputArray _src, cv::Point peakLocation, cv::Size weightBoxSize)
+static Point2d weightedCentroid(InputArray _src, cv::Point peakLocation, cv::Size weightBoxSize, double* response)
 {
     Mat src = _src.getMat();
 
@@ -475,6 +475,9 @@
         }
     }
 
+    if(response)
+        *response = sumIntensity;
+
     sumIntensity += DBL_EPSILON; // prevent div0 problems...
 
     centroid.x /= sumIntensity;
@@ -485,7 +488,7 @@
 
 }
 
-cv::Point2d cv::phaseCorrelate(InputArray _src1, InputArray _src2, InputArray _window)
+cv::Point2d cv::phaseCorrelate(InputArray _src1, InputArray _src2, InputArray _window, double* response)
 {
     Mat src1 = _src1.getMat();
     Mat src2 = _src2.getMat();
@@ -553,8 +556,12 @@
 
     // get the phase shift with sub-pixel accuracy, 5x5 window seems about right here...
     Point2d t;
-    t = weightedCentroid(C, peakLoc, Size(5, 5));
+    t = weightedCentroid(C, peakLoc, Size(5, 5), response);
 
+    // max response is M*N (not exactly, might be slightly larger due to rounding errors)
+    if(response)
+        *response /= M*N;
+
     // adjust shift relative to image center...
     Point2d center((double)padded1.cols / 2.0, (double)padded1.rows / 2.0);
 
