Merge branch 'master' into gpu-cuda-rename

Conflicts:
	modules/cudaoptflow/perf/perf_optflow.cpp
	modules/cudaoptflow/src/tvl1flow.cpp
	samples/gpu/stereo_multi.cpp

3rdparty/libtiff/CMakeLists.txt

@@ -17,7 +17,7 @@ check_include_file(string.h HAVE_STRING_H)
 check_include_file(sys/types.h HAVE_SYS_TYPES_H)
 check_include_file(unistd.h HAVE_UNISTD_H)
 
-if(WIN32)
+if(WIN32 AND NOT HAVE_WINRT)
   set(USE_WIN32_FILEIO 1)
 endif()
 
@@ -79,14 +79,12 @@ set(lib_srcs
     "${CMAKE_CURRENT_BINARY_DIR}/tif_config.h"
     )
 
-if(UNIX)
+if(WIN32 AND NOT HAVE_WINRT)
+  list(APPEND lib_srcs tif_win32.c)
+else()
   list(APPEND lib_srcs tif_unix.c)
 endif()
-
-if(WIN32)
-  list(APPEND lib_srcs tif_win32.c)
-endif(WIN32)
-
+  
 ocv_warnings_disable(CMAKE_C_FLAGS -Wno-unused-but-set-variable -Wmissing-prototypes -Wmissing-declarations -Wundef -Wunused -Wsign-compare
                                    -Wcast-align -Wshadow -Wno-maybe-uninitialized -Wno-pointer-to-int-cast -Wno-int-to-pointer-cast)
 ocv_warnings_disable(CMAKE_C_FLAGS -Wunused-parameter) # clang

cmake/OpenCVCRTLinkage.cmake

@@ -36,6 +36,9 @@ endif()
 
 if (HAVE_WINRT)
   add_definitions(/DWINVER=0x0602 /DNTDDI_VERSION=NTDDI_WIN8 /D_WIN32_WINNT=0x0602)
+  set(CMAKE_SHARED_LINKER_FLAGS "${CMAKE_SHARED_LINKER_FLAGS} /appcontainer")
+  set(CMAKE_MODULE_LINKER_FLAGS "${CMAKE_MODULE_LINKER_FLAGS} /appcontainer")
+  set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /appcontainer")
 endif()
 
 if(NOT BUILD_SHARED_LIBS AND BUILD_WITH_STATIC_CRT)

doc/tutorials/core/mat_the_basic_image_container/mat_the_basic_image_container.rst

@@ -143,7 +143,7 @@ Although *Mat* works really well as an image container, it is also a general mat
 
     You cannot initialize the matrix values with this construction. It will only reallocate its matrix data memory if the new size will not fit into the old one.
 
-   + MATLAB style initializer: :basicstructures:`zeros() <mat-zeros>`, :basicstructures:`ones() <mat-ones>`, ::basicstructures:`eyes() <mat-eye>`. Specify size and data type to use:
+   + MATLAB style initializer: :basicstructures:`zeros() <mat-zeros>`, :basicstructures:`ones() <mat-ones>`, :basicstructures:`eye() <mat-eye>`. Specify size and data type to use:
 
      .. literalinclude:: ../../../../samples/cpp/tutorial_code/core/mat_the_basic_image_container/mat_the_basic_image_container.cpp
         :language: cpp

modules/bioinspired/src/precomp.cpp

-/*M///////////////////////////////////////////////////////////////////////////////////////
-//
-//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
-//
-//  By downloading, copying, installing or using the software you agree to this license.
-//  If you do not agree to this license, do not download, install,
-//  copy or use the software.
-//
-//
-//                        Intel License Agreement
-//                For Open Source Computer Vision Library
-//
-// Copyright (C) 2000, Intel Corporation, all rights reserved.
-// Third party copyrights are property of their respective owners.
-//
-// Redistribution and use in source and binary forms, with or without modification,
-// are permitted provided that the following conditions are met:
-//
-//   * Redistribution's of source code must retain the above copyright notice,
-//     this list of conditions and the following disclaimer.
-//
-//   * Redistribution's in binary form must reproduce the above copyright notice,
-//     this list of conditions and the following disclaimer in the documentation
-//     and/or other materials provided with the distribution.
-//
-//   * The name of Intel Corporation may not be used to endorse or promote products
-//     derived from this software without specific prior written permission.
-//
-// This software is provided by the copyright holders and contributors "as is" and
-// any express or implied warranties, including, but not limited to, the implied
-// warranties of merchantability and fitness for a particular purpose are disclaimed.
-// In no event shall the Intel Corporation or contributors be liable for any direct,
-// indirect, incidental, special, exemplary, or consequential damages
-// (including, but not limited to, procurement of substitute goods or services;
-// loss of use, data, or profits; or business interruption) however caused
-// and on any theory of liability, whether in contract, strict liability,
-// or tort (including negligence or otherwise) arising in any way out of
-// the use of this software, even if advised of the possibility of such damage.
-//
-//M*/
-
-#include "precomp.hpp"
-
-/* End of file. */

modules/bioinspired/src/retina.cpp

@@ -628,6 +628,7 @@ void RetinaImpl::_init(const cv::Size inputSz, const bool colorMode, int colorSa
            delete _retinaFilter;
     _retinaFilter = new RetinaFilter(inputSz.height, inputSz.width, colorMode, colorSamplingMethod, useRetinaLogSampling, reductionFactor, samplingStrenght);
 
+    _retinaParameters.OPLandIplParvo.colorMode = colorMode;
     // prepare the default parameter XML file with default setup
     setup(_retinaParameters);
 

modules/bioinspired/test/test_precomp.cpp

-#include "test_precomp.hpp"

modules/core/include/opencv2/core.hpp

@@ -670,6 +670,10 @@ public:
     //! reconstructs the original vector from the projection
     void backProject(InputArray vec, OutputArray result) const;
 
+    //! write and load PCA matrix
+    void write(FileStorage& fs ) const;
+    void read(const FileNode& fs);
+
     Mat eigenvectors; //!< eigenvectors of the covariation matrix
     Mat eigenvalues; //!< eigenvalues of the covariation matrix
     Mat mean; //!< mean value subtracted before the projection and added after the back projection

modules/core/perf/perf_stat.cpp

@@ -83,8 +83,8 @@ PERF_TEST_P(Size_MatType, meanStdDev_mask, TYPICAL_MATS)
 
     TEST_CYCLE() meanStdDev(src, mean, dev, mask);
 
-    SANITY_CHECK(mean, 1e-6);
-    SANITY_CHECK(dev, 1e-6);
+    SANITY_CHECK(mean, 1e-5);
+    SANITY_CHECK(dev, 1e-5);
 }
 
 PERF_TEST_P(Size_MatType, countNonZero, testing::Combine( testing::Values( TYPICAL_MAT_SIZES ), testing::Values( CV_8UC1, CV_8SC1, CV_16UC1, CV_16SC1, CV_32SC1, CV_32FC1, CV_64FC1 ) ))

modules/core/src/arithm.cpp

@@ -2153,10 +2153,30 @@ cmp_(const T* src1, size_t step1, const T* src2, size_t step2,
     }
 }
 
+#if ARITHM_USE_IPP
+inline static IppCmpOp convert_cmp(int _cmpop)
+{
+    return _cmpop == CMP_EQ ? ippCmpEq :
+        _cmpop == CMP_GT ? ippCmpGreater :
+        _cmpop == CMP_GE ? ippCmpGreaterEq :
+        _cmpop == CMP_LT ? ippCmpLess :
+        _cmpop == CMP_LE ? ippCmpLessEq :
+        (IppCmpOp)-1;
+}
+#endif
 
 static void cmp8u(const uchar* src1, size_t step1, const uchar* src2, size_t step2,
                   uchar* dst, size_t step, Size size, void* _cmpop)
 {
+#if ARITHM_USE_IPP
+    IppCmpOp op = convert_cmp(*(int *)_cmpop);
+    if( op  >= 0 )
+    {
+        fixSteps(size, sizeof(dst[0]), step1, step2, step);
+        if( ippiCompare_8u_C1R(src1, (int)step1, src2, (int)step2, dst, (int)step, (IppiSize&)size, op) >= 0 )
+            return;
+    }
+#endif
   //vz optimized  cmp_(src1, step1, src2, step2, dst, step, size, *(int*)_cmpop);
     int code = *(int*)_cmpop;
     step1 /= sizeof(src1[0]);
@@ -2231,12 +2251,30 @@ static void cmp8s(const schar* src1, size_t step1, const schar* src2, size_t ste
 static void cmp16u(const ushort* src1, size_t step1, const ushort* src2, size_t step2,
                   uchar* dst, size_t step, Size size, void* _cmpop)
 {
+#if ARITHM_USE_IPP
+    IppCmpOp op = convert_cmp(*(int *)_cmpop);
+    if( op  >= 0 )
+    {
+        fixSteps(size, sizeof(dst[0]), step1, step2, step);
+        if( ippiCompare_16u_C1R(src1, (int)step1, src2, (int)step2, dst, (int)step, (IppiSize&)size, op) >= 0 )
+            return;
+    }
+#endif
     cmp_(src1, step1, src2, step2, dst, step, size, *(int*)_cmpop);
 }
 
 static void cmp16s(const short* src1, size_t step1, const short* src2, size_t step2,
                   uchar* dst, size_t step, Size size, void* _cmpop)
 {
+#if ARITHM_USE_IPP
+    IppCmpOp op = convert_cmp(*(int *)_cmpop);
+    if( op  > 0 )
+    {
+        fixSteps(size, sizeof(dst[0]), step1, step2, step);
+        if( ippiCompare_16s_C1R(src1, (int)step1, src2, (int)step2, dst, (int)step, (IppiSize&)size, op) >= 0 )
+            return;
+    }
+#endif
    //vz optimized cmp_(src1, step1, src2, step2, dst, step, size, *(int*)_cmpop);
 
     int code = *(int*)_cmpop;
@@ -2334,6 +2372,15 @@ static void cmp32s(const int* src1, size_t step1, const int* src2, size_t step2,
 static void cmp32f(const float* src1, size_t step1, const float* src2, size_t step2,
                   uchar* dst, size_t step, Size size, void* _cmpop)
 {
+#if ARITHM_USE_IPP
+    IppCmpOp op = convert_cmp(*(int *)_cmpop);
+    if( op  >= 0 )
+    {
+        fixSteps(size, sizeof(dst[0]), step1, step2, step);
+        if( ippiCompare_32f_C1R(src1, (int)step1, src2, (int)step2, dst, (int)step, (IppiSize&)size, op) >= 0 )
+            return;
+    }
+#endif
     cmp_(src1, step1, src2, step2, dst, step, size, *(int*)_cmpop);
 }
 

modules/core/src/dxt.cpp

@@ -50,6 +50,13 @@ namespace cv
 # pragma warning(disable: 4748)
 #endif
 
+#if defined HAVE_IPP && IPP_VERSION_MAJOR >= 7
+#define USE_IPP_DFT 1
+#else
+#undef USE_IPP_DFT
+#endif
+
+
 /****************************************************************************************\
                                Discrete Fourier Transform
 \****************************************************************************************/
@@ -455,7 +462,7 @@ template<> struct DFT_VecR4<float>
 
 #endif
 
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
 static void ippsDFTFwd_CToC( const Complex<float>* src, Complex<float>* dst,
                              const void* spec, uchar* buf)
 {
@@ -517,7 +524,7 @@ DFT( const Complex<T>* src, Complex<T>* dst, int n,
      int nf, const int* factors, const int* itab,
      const Complex<T>* wave, int tab_size,
      const void*
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
      spec
 #endif
      , Complex<T>* buf,
@@ -537,7 +544,7 @@ DFT( const Complex<T>* src, Complex<T>* dst, int n,
     T scale = (T)_scale;
     int tab_step;
 
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
     if( spec )
     {
         if( !inv )
@@ -957,7 +964,7 @@ DFT( const Complex<T>* src, Complex<T>* dst, int n,
 template<typename T> static void
 RealDFT( const T* src, T* dst, int n, int nf, int* factors, const int* itab,
          const Complex<T>* wave, int tab_size, const void*
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
          spec
 #endif
          ,
@@ -968,11 +975,18 @@ RealDFT( const T* src, T* dst, int n, int nf, int* factors, const int* itab,
     int j, n2 = n >> 1;
     dst += complex_output;
 
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
     if( spec )
     {
         ippsDFTFwd_RToPack( src, dst, spec, (uchar*)buf );
-        goto finalize;
+        if( complex_output )
+        {
+            dst[-1] = dst[0];
+            dst[0] = 0;
+            if( (n & 1) == 0 )
+                dst[n] = 0;
+        }
+        return;
     }
 #endif
     assert( tab_size == n );
@@ -1056,15 +1070,11 @@ RealDFT( const T* src, T* dst, int n, int nf, int* factors, const int* itab,
         }
     }
 
-#ifdef HAVE_IPP
-finalize:
-#endif
     if( complex_output && (n & 1) == 0 )
     {
         dst[-1] = dst[0];
         dst[0] = 0;
-        if( (n & 1) == 0 )
-            dst[n] = 0;
+        dst[n] = 0;
     }
 }
 
@@ -1076,7 +1086,7 @@ template<typename T> static void
 CCSIDFT( const T* src, T* dst, int n, int nf, int* factors, const int* itab,
          const Complex<T>* wave, int tab_size,
          const void*
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
          spec
 #endif
          , Complex<T>* buf,
@@ -1097,7 +1107,7 @@ CCSIDFT( const T* src, T* dst, int n, int nf, int* factors, const int* itab,
         ((T*)src)[1] = src[0];
         src++;
     }
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
     if( spec )
     {
         ippsDFTInv_PackToR( src, dst, spec, (uchar*)buf );
@@ -1225,7 +1235,7 @@ CCSIDFT( const T* src, T* dst, int n, int nf, int* factors, const int* itab,
         }
     }
 
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
 finalize:
 #endif
     if( complex_input )
@@ -1458,7 +1468,7 @@ static void CCSIDFT_64f( const double* src, double* dst, int n, int nf, int* fac
 
 }
 
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
 typedef IppStatus (CV_STDCALL* IppDFTGetSizeFunc)(int, int, IppHintAlgorithm, int*, int*, int*);
 typedef IppStatus (CV_STDCALL* IppDFTInitFunc)(int, int, IppHintAlgorithm, void*, uchar*);
 #endif
@@ -1486,7 +1496,7 @@ void cv::dft( InputArray _src0, OutputArray _dst, int flags, int nonzero_rows )
     int elem_size = (int)src.elemSize1(), complex_elem_size = elem_size*2;
     int factors[34];
     bool inplace_transform = false;
-#ifdef HAVE_IPP
+#ifdef USE_IPP_DFT
     AutoBuffer<uchar> ippbuf;
     int ipp_norm_flag = !(flags & DFT_SCALE) ? 8 : inv ? 2 : 1;
 #endif
@@ -1546,12 +1556,8 @@ void cv::dft( InputArray _src0, OutputArray _dst, int flags, int nonzero_rows )
         }
 
         spec = 0;
-#ifdef HAVE_IPP
-        if(
-#if IPP_VERSION_MAJOR >= 7
-           depth == CV_32F && // IPP 7.x and 8.0 have bug somewhere in double-precision DFT
-#endif
-           len*count >= 64 ) // use IPP DFT if available
+#ifdef USE_IPP_DFT
+        if( len*count >= 64 ) // use IPP DFT if available
         {
             int specsize=0, initsize=0, worksize=0;
             IppDFTGetSizeFunc getSizeFunc = 0;

modules/core/src/glob.cpp

@@ -91,6 +91,7 @@ namespace
         if(dir->handle == INVALID_HANDLE_VALUE)
         {
             /*closedir will do all cleanup*/
+            delete dir;
             return 0;
         }
         return dir;
@@ -140,6 +141,7 @@ static bool isDir(const cv::String& path, DIR* dir)
 {
 #if defined WIN32 || defined _WIN32 || defined WINCE
     DWORD attributes;
+    BOOL status = TRUE;
     if (dir)
         attributes = dir->data.dwFileAttributes;
     else
@@ -149,14 +151,14 @@ static bool isDir(const cv::String& path, DIR* dir)
         wchar_t wpath[MAX_PATH];
         size_t copied = mbstowcs(wpath, path.c_str(), MAX_PATH);
         CV_Assert((copied != MAX_PATH) && (copied != (size_t)-1));
-        ::GetFileAttributesExW(wpath, GetFileExInfoStandard, &all_attrs);
+        status = ::GetFileAttributesExW(wpath, GetFileExInfoStandard, &all_attrs);
 #else
-        ::GetFileAttributesExA(path.c_str(), GetFileExInfoStandard, &all_attrs);
+        status = ::GetFileAttributesExA(path.c_str(), GetFileExInfoStandard, &all_attrs);
 #endif
         attributes = all_attrs.dwFileAttributes;
     }
 
-    return (attributes != INVALID_FILE_ATTRIBUTES) && ((attributes & FILE_ATTRIBUTE_DIRECTORY) != 0);
+    return status && ((attributes & FILE_ATTRIBUTE_DIRECTORY) != 0);
 #else
     (void)dir;
     struct stat stat_buf;

modules/core/src/matmul.cpp

@@ -2911,6 +2911,27 @@ PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, int maxComp
     return *this;
 }
 
+void PCA::write(FileStorage& fs ) const
+{
+    CV_Assert( fs.isOpened() );
+
+    fs << "name" << "PCA";
+    fs << "vectors" << eigenvectors;
+    fs << "values" << eigenvalues;
+    fs << "mean" << mean;
+}
+
+void PCA::read(const FileNode& fs)
+{
+    CV_Assert( !fs.empty() );
+    String name = (String)fs["name"];
+    CV_Assert( name == "PCA" );
+
+    cv::read(fs["vectors"], eigenvectors);
+    cv::read(fs["values"], eigenvalues);
+    cv::read(fs["mean"], mean);
+}
+
 template <typename T>
 int computeCumulativeEnergy(const Mat& eigenvalues, double retainedVariance)
 {

modules/core/test/test_mat.cpp

@@ -510,6 +510,32 @@ protected:
             return;
         }
     #endif
+        // Test read and write
+        FileStorage fs( "PCA_store.yml", FileStorage::WRITE );
+        rPCA.write( fs );
+        fs.release();
+
+        PCA lPCA;
+        fs.open( "PCA_store.yml", FileStorage::READ );
+        lPCA.read( fs.root() );
+        err = norm( rPCA.eigenvectors, lPCA.eigenvectors, CV_RELATIVE_L2 );
+        if( err > 0 )
+        {
+            ts->printf( cvtest::TS::LOG, "bad accuracy of write/load functions (YML); err = %f\n", err );
+            ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY );
+        }
+        err = norm( rPCA.eigenvalues, lPCA.eigenvalues, CV_RELATIVE_L2 );
+        if( err > 0 )
+        {
+            ts->printf( cvtest::TS::LOG, "bad accuracy of write/load functions (YML); err = %f\n", err );
+            ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY );
+        }
+        err = norm( rPCA.mean, lPCA.mean, CV_RELATIVE_L2 );
+        if( err > 0 )
+        {
+            ts->printf( cvtest::TS::LOG, "bad accuracy of write/load functions (YML); err = %f\n", err );
+            ts->set_failed_test_info( cvtest::TS::FAIL_BAD_ACCURACY );
+        }
     }
 };
 

modules/cudaoptflow/perf/perf_optflow.cpp

@@ -368,8 +368,8 @@ PERF_TEST_P(ImagePair, OpticalFlowDual_TVL1,
 
         TEST_CYCLE() d_alg(d_frame0, d_frame1, u, v);
 
-        CUDA_SANITY_CHECK(u, 1e-2);
-        CUDA_SANITY_CHECK(v, 1e-2);
+        CUDA_SANITY_CHECK(u, 1e-1);
+        CUDA_SANITY_CHECK(v, 1e-1);
     }
     else
     {

modules/cudaoptflow/src/cuda/tvl1flow.cu

@@ -211,7 +211,7 @@ namespace tvl1flow
                               const PtrStepf grad, const PtrStepf rho_c,
                               const PtrStepf p11, const PtrStepf p12, const PtrStepf p21, const PtrStepf p22,
                               PtrStepf u1, PtrStepf u2, PtrStepf error,
-                              const float l_t, const float theta)
+                              const float l_t, const float theta, const bool calcError)
     {
         const int x = blockIdx.x * blockDim.x + threadIdx.x;
         const int y = blockIdx.y * blockDim.y + threadIdx.y;
@@ -265,21 +265,24 @@ namespace tvl1flow
         u1(y, x) = u1NewVal;
         u2(y, x) = u2NewVal;
 
-        const float n1 = (u1OldVal - u1NewVal) * (u1OldVal - u1NewVal);
-        const float n2 = (u2OldVal - u2NewVal) * (u2OldVal - u2NewVal);
-        error(y, x) = n1 + n2;
+        if (calcError)
+        {
+            const float n1 = (u1OldVal - u1NewVal) * (u1OldVal - u1NewVal);
+            const float n2 = (u2OldVal - u2NewVal) * (u2OldVal - u2NewVal);
+            error(y, x) = n1 + n2;
+        }
     }
 
     void estimateU(PtrStepSzf I1wx, PtrStepSzf I1wy,
                    PtrStepSzf grad, PtrStepSzf rho_c,
                    PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22,
                    PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf error,
-                   float l_t, float theta)
+                   float l_t, float theta, bool calcError)
     {
         const dim3 block(32, 8);
         const dim3 grid(divUp(I1wx.cols, block.x), divUp(I1wx.rows, block.y));
 
-        estimateUKernel<<<grid, block>>>(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, u1, u2, error, l_t, theta);
+        estimateUKernel<<<grid, block>>>(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, u1, u2, error, l_t, theta, calcError);
         cudaSafeCall( cudaGetLastError() );
 
         cudaSafeCall( cudaDeviceSynchronize() );

modules/cudaoptflow/src/tvl1flow.cpp

@@ -173,7 +173,7 @@ namespace tvl1flow
                    PtrStepSzf grad, PtrStepSzf rho_c,
                    PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22,
                    PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf error,
-                   float l_t, float theta);
+                   float l_t, float theta, bool calcError);
     void estimateDualVariables(PtrStepSzf u1, PtrStepSzf u2, PtrStepSzf p11, PtrStepSzf p12, PtrStepSzf p21, PtrStepSzf p22, float taut);
 }
 
@@ -218,12 +218,24 @@ void cv::cuda::OpticalFlowDual_TVL1_CUDA::procOneScale(const GpuMat& I0, const G
         warpBackward(I0, I1, I1x, I1y, u1, u2, I1w, I1wx, I1wy, grad, rho_c);
 
         double error = std::numeric_limits<double>::max();
+        double prevError = 0.0;
         for (int n = 0; error > scaledEpsilon && n < iterations; ++n)
         {
-            estimateU(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, u1, u2, diff, l_t, static_cast<float>(theta));
+            // some tweaks to make sum operation less frequently
+            bool calcError = (epsilon > 0) && (n & 0x1) && (prevError < scaledEpsilon);
 
-            if (epsilon > 0)
+            estimateU(I1wx, I1wy, grad, rho_c, p11, p12, p21, p22, u1, u2, diff, l_t, static_cast<float>(theta), calcError);
+
+            if (calcError)
+            {
                 error = cuda::sum(diff, norm_buf)[0];
+                prevError = error;
+            }
+            else
+            {
+                error = std::numeric_limits<double>::max();
+                prevError -= scaledEpsilon;
+            }
 
             estimateDualVariables(u1, u2, p11, p12, p21, p22, taut);
         }