Merge pull request #23664 from asmorkalov:as/solveLP_infeasible

Added final constrants check to solveLP to filter out flating-point numeric issues

modules/core/include/opencv2/core/optim.hpp

@@ -256,6 +256,7 @@ public:
 //! return codes for cv::solveLP() function
 enum SolveLPResult
 {
+    SOLVELP_LOST   = -3, //!< problem is feasible, but solver lost solution due to floating-point arithmetic errors
     SOLVELP_UNBOUNDED    = -2, //!< problem is unbounded (target function can achieve arbitrary high values)
     SOLVELP_UNFEASIBLE    = -1, //!< problem is unfeasible (there are no points that satisfy all the constraints imposed)
     SOLVELP_SINGLE    = 0, //!< there is only one maximum for target function
@@ -291,9 +292,13 @@ in the latter case it is understood to correspond to \f$c^T\f$.
 and the remaining to \f$A\f$. It should contain 32- or 64-bit floating point numbers.
 @param z The solution will be returned here as a column-vector - it corresponds to \f$c\f$ in the
 formulation above. It will contain 64-bit floating point numbers.
+@param constr_eps allowed numeric disparity for constraints
 @return One of cv::SolveLPResult
  */
-CV_EXPORTS_W int solveLP(const Mat& Func, const Mat& Constr, Mat& z);
+CV_EXPORTS_W int solveLP(const Mat& Func, const Mat& Constr, Mat& z, double constr_eps);
+
+/** @overload */
+CV_EXPORTS int solveLP(const Mat& Func, const Mat& Constr, Mat& z);
 
 //! @}
 

modules/core/src/lpsolver.cpp

@@ -90,7 +90,7 @@ static void swap_columns(Mat_<double>& A,int col1,int col2);
 #define SWAP(type,a,b) {type tmp=(a);(a)=(b);(b)=tmp;}
 
 //return codes:-2 (no_sol - unbdd),-1(no_sol - unfsbl), 0(single_sol), 1(multiple_sol=>least_l2_norm)
-int solveLP(const Mat& Func, const Mat& Constr, Mat& z){
+int solveLP(const Mat& Func, const Mat& Constr, Mat& z, double constr_eps){
     dprintf(("call to solveLP\n"));
 
     //sanity check (size, type, no. of channels)
@@ -140,9 +140,24 @@ int solveLP(const Mat& Func, const Mat& Constr, Mat& z){
         }
     }
 
+    //check constraints feasibility
+    Mat prod = Constr(Rect(0, 0, Constr.cols - 1, Constr.rows)) * z;
+    Mat constr_check = Constr.col(Constr.cols - 1) - prod;
+    double min_value = 0.0;
+    minMaxIdx(constr_check, &min_value);
+    if (min_value < -constr_eps)
+    {
+        return SOLVELP_LOST;
+    }
+
     return res;
 }
 
+CV_EXPORTS_W int solveLP(const Mat& Func, const Mat& Constr, Mat& z)
+{
+    return solveLP(Func, Constr, z, 1e-12);
+}
+
 static int initialize_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>& N,vector<int>& B,vector<unsigned int>& indexToRow){
     N.resize(c.cols);
     N[0]=0;
@@ -255,7 +270,7 @@ static int inner_simplex(Mat_<double>& c, Mat_<double>& b,double& v,vector<int>&
         dprintf(("iteration #%d\n",count));
         count++;
 
-        static MatIterator_<double> pos_ptr;
+        MatIterator_<double> pos_ptr;
         int e=-1,pos_ctr=0,min_var=INT_MAX;
         bool all_nonzero=true;
         for(pos_ptr=c.begin();pos_ptr!=c.end();pos_ptr++,pos_ctr++){

modules/core/test/test_lpsolver.cpp

@@ -151,4 +151,18 @@ TEST(Core_LPSolver, issue_12337)
     //need to update interface: EXPECT_ANY_THROW(Mat1b z_8u; cv::solveLP(A, B, z_8u));
 }
 
+// NOTE: Test parameters found experimentally to get numerically inaccurate result.
+// The test behaviour may change after algorithm tuning and may removed.
+TEST(Core_LPSolver, issue_12343)
+{
+    Mat A = (cv::Mat_<double>(4, 1) << 3., 3., 3., 4.);
+    Mat B = (cv::Mat_<double>(4, 5) << 0., 1., 4., 4., 3.,
+                                       3., 1., 2., 2., 3.,
+                                       4., 4., 0., 1., 4.,
+                                       4., 0., 4., 1., 4.);
+    Mat z;
+    int result = cv::solveLP(A, B, z);
+    EXPECT_EQ(SOLVELP_LOST, result);
+}
+
 }} // namespace