diff --git a/3rdparty/libwebp/CMakeLists.txt b/3rdparty/libwebp/CMakeLists.txt index 13593526cd..ff9a383229 100644 --- a/3rdparty/libwebp/CMakeLists.txt +++ b/3rdparty/libwebp/CMakeLists.txt @@ -21,7 +21,7 @@ if(ANDROID AND ARMEABI_V7A AND NOT NEON) endforeach() endif() -if(WIN32) +if(WIN32 AND (X86_64 OR X86)) foreach(file ${lib_srcs}) if("${file}" MATCHES "_avx2.c") if(MSVC) diff --git a/cmake/OpenCVDownload.cmake b/cmake/OpenCVDownload.cmake index 40e48d7465..0f8de153d8 100644 --- a/cmake/OpenCVDownload.cmake +++ b/cmake/OpenCVDownload.cmake @@ -282,11 +282,17 @@ ${OPENCV_DOWNLOAD_LOG} ocv_download_log("#mkdir \"${DL_DESTINATION_DIR}\"") file(MAKE_DIRECTORY "${DL_DESTINATION_DIR}") ocv_download_log("#unpack \"${DL_DESTINATION_DIR}\" \"${CACHE_CANDIDATE}\"") - execute_process(COMMAND "${CMAKE_COMMAND}" -E tar xzf "${CACHE_CANDIDATE}" - WORKING_DIRECTORY "${DL_DESTINATION_DIR}" - RESULT_VARIABLE res) - if(NOT res EQUAL 0) - message(FATAL_ERROR "${__msg_prefix}Unpack failed: ${res}") + if(CMAKE_VERSION VERSION_GREATER_EQUAL "3.18") + file(ARCHIVE_EXTRACT + INPUT "${CACHE_CANDIDATE}" + DESTINATION "${DL_DESTINATION_DIR}") + else() + execute_process(COMMAND "${CMAKE_COMMAND}" -E tar xzf "${CACHE_CANDIDATE}" + WORKING_DIRECTORY "${DL_DESTINATION_DIR}" + RESULT_VARIABLE res) + if(NOT res EQUAL 0) + message(FATAL_ERROR "${__msg_prefix}Unpack failed: ${res}") + endif() endif() else() ocv_download_log("#copy \"${COPY_DESTINATION}\" \"${CACHE_CANDIDATE}\"") diff --git a/cmake/OpenCVFindIPPIW.cmake b/cmake/OpenCVFindIPPIW.cmake index 1e065e10f6..c17d960fb3 100644 --- a/cmake/OpenCVFindIPPIW.cmake +++ b/cmake/OpenCVFindIPPIW.cmake @@ -43,7 +43,12 @@ file(TO_CMAKE_PATH "${IPPROOT}" IPPROOT) # This function detects Intel IPP Integration Wrappers version by analyzing .h file macro(ippiw_setup PATH BUILD) + set(IPP_NEW_LAYOUT 0) set(FILE "${PATH}/include/iw/iw_version.h") + if(NOT EXISTS "${FILE}") + set(IPP_NEW_LAYOUT 1) + set(FILE "${PATH}/include/ipp/iw/iw_version.h") + endif() if(${BUILD}) ippiw_debugmsg("Checking sources: ${PATH}") else() @@ -83,7 +88,9 @@ macro(ippiw_setup PATH BUILD) set(IPP_IW_LIBRARY ippiw) set(IPP_IW_INCLUDES "${IPP_IW_PATH}/include") set(IPP_IW_LIBRARIES ${IPP_IW_LIBRARY}) - execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different "${OpenCV_SOURCE_DIR}/3rdparty/ippicv/CMakeLists.txt" "${IPP_IW_PATH}/") + if(NOT EXISTS "${IPP_IW_PATH}/CMakeLists.txt") + file(COPY "${OpenCV_SOURCE_DIR}/3rdparty/ippicv/CMakeLists.txt" DESTINATION "${IPP_IW_PATH}") + endif() add_subdirectory("${IPP_IW_PATH}/" ${OpenCV_BINARY_DIR}/3rdparty/ippiw) set(HAVE_IPP_IW 1) @@ -97,11 +104,20 @@ macro(ippiw_setup PATH BUILD) endif() else() # check binaries - if(IPP_X64) - set(FILE "${PATH}/lib/intel64/${CMAKE_STATIC_LIBRARY_PREFIX}ipp_iw${CMAKE_STATIC_LIBRARY_SUFFIX}") + if (IPP_NEW_LAYOUT) + if(IPP_X64) + set(IW_LIB_DIR "${PATH}/lib") + else() + set(IW_LIB_DIR "${PATH}/lib32") + endif() else() - set(FILE "${PATH}/lib/ia32/${CMAKE_STATIC_LIBRARY_PREFIX}ipp_iw${CMAKE_STATIC_LIBRARY_SUFFIX}") + if(IPP_X64) + set(IW_LIB_DIR "${PATH}/lib/intel64") + else() + set(IW_LIB_DIR "${PATH}/lib/ia32") + endif() endif() + set(FILE "${IW_LIB_DIR}/${CMAKE_STATIC_LIBRARY_PREFIX}ipp_iw${CMAKE_STATIC_LIBRARY_SUFFIX}") if(EXISTS ${FILE}) ippiw_debugmsg("binaries\tyes (64=${IPP_X64})") set(IPP_IW_PATH "${PATH}") diff --git a/cmake/OpenCVModule.cmake b/cmake/OpenCVModule.cmake index dcbf46f05e..a171ebfb5c 100644 --- a/cmake/OpenCVModule.cmake +++ b/cmake/OpenCVModule.cmake @@ -991,7 +991,7 @@ macro(_ocv_create_module) ocv_target_link_libraries(${the_module} PUBLIC ${OPENCV_MODULE_${the_module}_DEPS_EXT} INTERFACE ${OPENCV_MODULE_${the_module}_DEPS_EXT} ) - ocv_target_link_libraries(${the_module} PRIVATE ${OPENCV_LINKER_LIBS} ${OPENCV_HAL_LINKER_LIBS} ${IPP_LIBS} ${ARGN}) + ocv_target_link_libraries(${the_module} PRIVATE ${OPENCV_LINKER_LIBS} ${OPENCV_HAL_LINKER_LIBS} ${ARGN}) if (NOT ENABLE_CUDA_FIRST_CLASS_LANGUAGE AND HAVE_CUDA) ocv_target_link_libraries(${the_module} PRIVATE ${CUDA_LIBRARIES} ${CUDA_npp_LIBRARY}) endif() diff --git a/doc/opencv.bib b/doc/opencv.bib index 316baface9..0bd5731ee0 100644 --- a/doc/opencv.bib +++ b/doc/opencv.bib @@ -769,6 +769,22 @@ number = {8}, publisher = {IOP Publishing Ltd} } +@article{Lourakis2009_sba, + author = {Lourakis, Manolis I. A. and Argyros, Antonis A.}, + title = {SBA: A Software Package for Generic Sparse Bundle Adjustment}, + year = {2009}, + month = mar, + journal = {ACM Transactions on Mathematical Software}, + volume = {36}, + number = {1}, + articleno = {2}, + pages = {2:1--2:30}, + numpages = {30}, + publisher = {Association for Computing Machinery}, + doi = {10.1145/1486525.1486527}, + url = {https://scispace.com/pdf/sba-a-software-package-for-generic-sparse-bundle-adjustment-1d4hp0z31z.pdf}, + month_numeric = {3} +} @article{LowIlie2003, author = {Kok-Lim Low, Adrian Ilie}, year = {2003}, @@ -1253,6 +1269,20 @@ publisher = {Taylor \& Francis}, url = {https://www.olivier-augereau.com/docs/2004JGraphToolsTelea.pdf} } +@incollection{Triggs2000_bundle_adjustment, + author = {Triggs, Bill and McLauchlan, Philip F. and Hartley, Richard I. and Fitzgibbon, Andrew W.}, + title = {Bundle Adjustment---A Modern Synthesis}, + booktitle = {Vision Algorithms: Theory and Practice}, + year = {2000}, + pages = {298--372}, + publisher = {Springer}, + series = {Lecture Notes in Computer Science}, + volume = {1883}, + editor = {Triggs, Bill and Zisserman, Andrew and Szeliski, Richard}, + doi = {10.1007/3-540-44480-7_21}, + isbn = {978-3-540-67973-8}, + url = {https://www.cs.jhu.edu/~misha/ReadingSeminar/Papers/Triggs00.pdf} +} @article{Tsai89, author = {R. Y. Tsai and R. K. Lenz}, journal = {IEEE Transactions on Robotics and Automation}, diff --git a/doc/py_tutorials/py_calib3d/py_calibration/py_calibration.markdown b/doc/py_tutorials/py_calib3d/py_calibration/py_calibration.markdown index a08a204d6c..401fc45dfb 100644 --- a/doc/py_tutorials/py_calib3d/py_calibration/py_calibration.markdown +++ b/doc/py_tutorials/py_calib3d/py_calibration/py_calibration.markdown @@ -206,17 +206,17 @@ Re-projection Error Re-projection error gives a good estimation of just how exact the found parameters are. The closer the re-projection error is to zero, the more accurate the parameters we found are. Given the intrinsic, distortion, rotation and translation matrices, we must first transform the object point to image point using **cv.projectPoints()**. Then, we can calculate -the absolute norm between what we got with our transformation and the corner finding algorithm. To -find the average error, we calculate the arithmetical mean of the errors calculated for all the -calibration images. +the norm between what we got with our transformation and the corner finding algorithm. To find the +RMSE (root mean squared error), we average the squared errors over all points and images, then take +the square root. @code{.py} mean_error = 0 for i in range(len(objpoints)): imgpoints2, _ = cv.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist) - error = cv.norm(imgpoints[i], imgpoints2, cv.NORM_L2)/len(imgpoints2) + error = cv.norm(imgpoints[i], imgpoints2, cv.NORM_L2SQR) / len(imgpoints2) mean_error += error -print( "total error: {}".format(mean_error/len(objpoints)) ) +print( "total error: {}".format(np.sqrt(mean_error/len(objpoints))) ) @endcode Exercises diff --git a/hal/ipp/CMakeLists.txt b/hal/ipp/CMakeLists.txt index 5ff630c8c0..345a7534b1 100644 --- a/hal/ipp/CMakeLists.txt +++ b/hal/ipp/CMakeLists.txt @@ -45,7 +45,7 @@ target_include_directories(ipphal PRIVATE ${IPP_INCLUDE_DIRS} ) -target_link_libraries(ipphal PUBLIC ${IPP_IW_LIBRARY} ${IPP_LIBRARIES}) +target_link_libraries(ipphal PUBLIC ${IPP_IW_LIBRARIES} ${IPP_LIBRARIES}) set_target_properties(ipphal PROPERTIES ARCHIVE_OUTPUT_DIRECTORY ${3P_LIBRARY_OUTPUT_PATH} DEBUG_POSTFIX "${OPENCV_DEBUG_POSTFIX}") diff --git a/modules/calib/include/opencv2/calib.hpp b/modules/calib/include/opencv2/calib.hpp index 5a06ad2af6..97f6249ca7 100644 --- a/modules/calib/include/opencv2/calib.hpp +++ b/modules/calib/include/opencv2/calib.hpp @@ -534,6 +534,7 @@ enum { CALIB_USE_INTRINSIC_GUESS = 0x00001, //!< Use user provided intrinsics as // for stereo rectification CALIB_ZERO_DISPARITY = 0x00400, //!< Deprecated synonim of @ref STEREO_ZERO_DISPARITY. See @ref stereoRectify. CALIB_USE_LU = (1 << 17), //!< use LU instead of SVD decomposition for solving. much faster but potentially less precise + CALIB_DISABLE_SCHUR_COMPLEMENT = (1 << 18), //!< disable Schur complement (use Bouguet calibration engine) CALIB_USE_EXTRINSIC_GUESS = (1 << 22), //!< For stereo and multi-view calibration. Use user provided extrinsics (R, T) as initial point for optimization // fisheye only flags CALIB_RECOMPUTE_EXTRINSIC = (1 << 23), //!< For fisheye model only. Recompute board position on each calibration iteration @@ -875,6 +876,7 @@ fx, fy, cx, cy that are optimized further. Otherwise, (cx, cy) is initially set center ( imageSize is used), and focal distances are computed in a least-squares fashion. Note, that if intrinsic parameters are known, there is no need to use this function just to estimate extrinsic parameters. Use @ref solvePnP instead. +- @ref CALIB_DISABLE_SCHUR_COMPLEMENT Disable Schur complement and use the Bouguet calibration engine (@cite Zhang2000, @cite BouguetMCT). - @ref CALIB_FIX_PRINCIPAL_POINT The principal point is not changed during the global optimization. It stays at the center or at a different location specified when @ref CALIB_USE_INTRINSIC_GUESS is set too. @@ -909,7 +911,9 @@ supplied distCoeffs matrix is used. Otherwise, it is set to 0. @return the overall RMS re-projection error. The function estimates the intrinsic camera parameters and extrinsic parameters for each of the -views. The algorithm is based on @cite Zhang2000 and @cite BouguetMCT . The coordinates of 3D object +views. By default, the optimization follows a sparse bundle adjustment formulation with Schur +complement; see @cite Triggs2000_bundle_adjustment and @cite Lourakis2009_sba for background. Use +@ref CALIB_DISABLE_SCHUR_COMPLEMENT to switch to the Bouguet calibration engine. The coordinates of 3D object points and their corresponding 2D projections in each view must be specified. That may be achieved by using an object with known geometry and easily detectable feature points. Such an object is called a calibration rig or calibration pattern, and OpenCV has built-in support for a chessboard as @@ -932,6 +936,10 @@ The algorithm performs the following steps: the projected (using the current estimates for camera parameters and the poses) object points objectPoints. See @ref projectPoints for details. +- In practice, robust acquisition is essential for stable results: use multiple board poses with + significant tilt, avoid collecting all views at a single working distance, span the expected + working-distance range (a larger board with larger squares can help for longer distances). + @note If you use a non-square (i.e. non-N-by-N) grid and @ref findChessboardCorners for calibration, and @ref calibrateCamera returns bad values (zero distortion coefficients, \f$c_x\f$ and @@ -1017,8 +1025,8 @@ less precise and less stable in some rare cases. @return the overall RMS re-projection error. The function estimates the intrinsic camera parameters and extrinsic parameters for each of the -views. The algorithm is based on @cite Zhang2000, @cite BouguetMCT and @cite strobl2011iccv. See -#calibrateCamera for other detailed explanations. +views. The object-releasing extension follows @cite strobl2011iccv and uses the same optimization +core as #calibrateCamera. See #calibrateCamera for other detailed explanations. @sa calibrateCamera, findChessboardCorners, solvePnP, initCameraMatrix2D, stereoCalibrate, undistort */ diff --git a/modules/calib/perf/perf_calibratecamera.cpp b/modules/calib/perf/perf_calibratecamera.cpp new file mode 100644 index 0000000000..95f7f07042 --- /dev/null +++ b/modules/calib/perf/perf_calibratecamera.cpp @@ -0,0 +1,164 @@ +// This file is part of OpenCV project. +// It is subject to the license terms in the LICENSE file found in the top-level directory +// of this distribution and at http://opencv.org/license.html +#include "perf_precomp.hpp" + +#include "opencv2/core/utils/filesystem.hpp" + +//#define SAVE_IMAGE_POINTS + +namespace opencv_test { + +#ifdef SAVE_IMAGE_POINTS + +static std::vector loadBulkImages(size_t max_images) +{ + const std::string data_dir = findDataDirectory("perf/calib3d/bulk_n500", false); + std::vector image_paths; + cv::utils::fs::glob(data_dir, "*.png", image_paths, false, false); + if (image_paths.empty()) + cv::utils::fs::glob(data_dir, "*.jpg", image_paths, false, false); + if (image_paths.empty()) + throw SkipTestException("No images found in perf/calib3d/bulk_n500"); + + std::sort(image_paths.begin(), image_paths.end()); + if (image_paths.size() > max_images) + image_paths.resize(max_images); + return image_paths; +} + +static std::vector> buildImagePoints(const std::vector& image_paths, const cv::Size pattern_size) +{ + std::vector> image_points; + image_points.reserve(image_paths.size()); + + for (const auto& path : image_paths) + { + Mat gray = imread(path, IMREAD_GRAYSCALE); + if (gray.empty()) + { + printf("Can't read image: %s\n", path.c_str()); + return std::vector>(); + } + + std::vector corners; + bool found = findChessboardCorners( + gray, pattern_size, corners, + CALIB_CB_ADAPTIVE_THRESH | CALIB_CB_NORMALIZE_IMAGE); + + if (found) + { + cornerSubPix(gray, corners, Size(11, 11), Size(-1, -1), + TermCriteria(TermCriteria::EPS + TermCriteria::COUNT, 30, 0.1)); + image_points.push_back(corners); + } + } + + return image_points; +} + +static void saveImagePoints(const std::vector>& image_points) +{ + const std::string points_file = "bulk_n500.yaml"; + cv::FileStorage fs(points_file, cv::FileStorage::WRITE | cv::FileStorage::FORMAT_YAML); + if (!fs.isOpened()) + { + printf("Cannot open yaml config \"%s\" for output\n", points_file.c_str()); + } + fs << "count" << (int)image_points.size(); + for (int i = 0; i < (int)image_points.size(); i++) + { + fs << cv::format("frame_%d", i) << image_points[i]; + } + + fs.release(); +} + +#else + +static std::vector> loadImagePoints() +{ + const std::string points_file = findDataFile("perf/calib3d/bulk_n500.yaml"); + cv::FileStorage fs(points_file, cv::FileStorage::READ); + if (!fs.isOpened()) + { + printf("Cannot open yaml config \"%s\" for output\n", points_file.c_str()); + return std::vector>(); + } + int count = fs["count"]; + std::vector> image_points(count); + + for (int i = 0; i < (int)image_points.size(); i++) + { + fs[cv::format("frame_%d", i)] >> image_points[i]; + } + + fs.release(); + + return image_points; +} + +#endif + +static std::vector> buildObjectPoints(const Size& pattern_size, float square_size, size_t count) +{ + std::vector board; + board.reserve(pattern_size.area()); + + for (int y = 0; y < pattern_size.height; ++y) + for (int x = 0; x < pattern_size.width; ++x) + board.push_back(Point3f(x * square_size, y * square_size, 0.f)); + + std::vector > object_points; + object_points.reserve(count); + for (size_t i = 0; i < count; i++) + object_points.push_back(board); + + return object_points; +} + +PERF_TEST(CalibrateCamera, BulkImages_N500) +{ + // NOTE: The images archive is published at https://dl.opencv.org/data/bulk_n500.zip + applyTestTag(CV_TEST_TAG_LONG); + + const cv::Size pattern_size(6, 8); + const cv::Size image_size(1280, 720); + +#ifdef SAVE_IMAGE_POINTS + std::vector image_paths = loadBulkImages(500); + std::vector> image_points = buildImagePoints(image_paths, pattern_size); + ASSERT_FALSE(image_points.empty()); + saveImagePoints(image_points); +#else + std::vector> image_points = loadImagePoints(); + ASSERT_FALSE(image_points.empty()); +#endif + + std::vector > object_points = buildObjectPoints(pattern_size, 1.0f, image_points.size()); + + Mat camera_matrix = Mat::eye(3, 3, CV_64F); + Mat dist_coeffs = Mat::zeros(8, 1, CV_64F); + std::vector rvecs; + std::vector tvecs; + double rms = 0.0; + + declare.in(image_points, object_points); + declare.out(camera_matrix, dist_coeffs); + declare.iterations(1); + + TEST_CYCLE() + { + camera_matrix = Mat::eye(3, 3, CV_64F); + dist_coeffs = Mat::zeros(8, 1, CV_64F); + rvecs.clear(); + tvecs.clear(); + rms = calibrateCamera(object_points, image_points, image_size, + camera_matrix, dist_coeffs, rvecs, tvecs, 0); + } + + EXPECT_NEAR(rms, 1.768263, 1e-4); + SANITY_CHECK_NOTHING(); +} + +} // namespace opencv_test diff --git a/modules/calib/src/calibration.cpp b/modules/calib/src/calibration.cpp index ede596d8cb..299f0c589a 100644 --- a/modules/calib/src/calibration.cpp +++ b/modules/calib/src/calibration.cpp @@ -46,6 +46,419 @@ #include #include +/* + This program implements Levenberg-Marquardt optimization using Schur complement to calibrate cameras. + It implements many speedups over the original implemntation. + The techniques used here are primarily based upon the following papers: + [1] "Bundle Adjustment -- A Modern Synthesis", Triggs et al., 1999 + [2] "SBA: A Software Package for Generic Sparse Bundle Adjustment", Lourakis and Argyros, 2000 +*/ + +namespace { + +class SchurLMSolver +{ +public: + enum State { STARTED, CALC_J, CHECK_ERR, DONE }; + + SchurLMSolver() + : currentError(0) + , state(STARTED) + , NINTRINSIC(18) + , nimages(0) + , n_global(0) + , lambdaLg10(-3) + , iters(0) + , prevErrNorm(DBL_MAX) + , solveMethod(cv::DECOMP_CHOLESKY) + { + } + + void init(int _nimages, int _NINTRINSIC, int nObjPoints, + cv::TermCriteria _criteria) + { + nimages = _nimages; + NINTRINSIC = _NINTRINSIC; + n_global = NINTRINSIC + nObjPoints * 3; + criteria = _criteria; + + U.create(n_global, n_global, CV_64F); + ea.create(n_global, 1, CV_64F); + deltaGlobal.create(n_global, 1, CV_64F); + + V.resize(nimages); + V_inv.resize(nimages); + W.resize(nimages); + eb.resize(nimages); + deltaLocal.resize(nimages); + + for (int i = 0; i < nimages; i++) + { + V[i].create(6, 6, CV_64F); + V_inv[i].create(6, 6, CV_64F); + W[i].create(n_global, 6, CV_64F); + eb[i].create(6, 1, CV_64F); + deltaLocal[i].create(6, 1, CV_64F); + } + + state = STARTED; + iters = 0; + lambdaLg10 = -3; + prevErrNorm = DBL_MAX; + solveMethod = cv::DECOMP_CHOLESKY; + } + + void reset() + { + U.setTo(0); + ea.setTo(0); + currentError = 0; + for (int i = 0; i < nimages; i++) + { + V[i].setTo(0); + W[i].setTo(0); + eb[i].setTo(0); + } + } + + void step(const cv::Mat& mask) + { + const double LOG10 = std::log(10.0); + double lambda = std::exp(lambdaLg10 * LOG10); + + // Invert V[i] blocks with regularization + for (int i = 0; i < nimages; i++) + { + cv::Mat V_reg = V[i].clone(); + V_reg.diag() *= (1.0 + lambda); + if (solveMethod == cv::DECOMP_QR) + { + cv::Mat I = cv::Mat::eye(6, 6, CV_64F); + cv::solve(V_reg, I, V_inv[i], solveMethod); + } + else + { + cv::invert(V_reg, V_inv[i], solveMethod); + } + } + + // Form Schur complement S and reduced RHS e + cv::Mat S = U.clone(); + S.diag() *= (1.0 + lambda); + cv::Mat e = ea.clone(); + + for (int i = 0; i < nimages; i++) + { + cv::Mat WtVinv = W[i] * V_inv[i]; + S -= WtVinv * W[i].t(); + e -= WtVinv * eb[i]; + } + + // Extract active submatrix using mask + std::vector mask_global(n_global); + const uchar* mask_ptr = mask.ptr(); + std::copy(mask_ptr, mask_ptr + NINTRINSIC, mask_global.begin()); + if (n_global > NINTRINSIC) + { + // Object points parameters follow intrinsic parameters in the global update vector + std::copy(mask_ptr + NINTRINSIC + nimages * 6, + mask_ptr + NINTRINSIC + nimages * 6 + (n_global - NINTRINSIC), + mask_global.begin() + NINTRINSIC); + } + + int nactive = 0; + for (int i = 0; i < n_global; i++) + if (mask_global[i]) nactive++; + + if (nactive == 0) + { + deltaGlobal.setTo(0); + for (int i = 0; i < nimages; i++) + deltaLocal[i].setTo(0); + return; + } + + cv::Mat S_sub(nactive, nactive, CV_64F); + cv::Mat e_sub(nactive, 1, CV_64F); + + int ii = 0; + for (int i = 0; i < n_global; i++) + { + if (!mask_global[i]) continue; + e_sub.at(ii) = e.at(i); + int jj = 0; + for (int j = 0; j < n_global; j++) + { + if (!mask_global[j]) continue; + S_sub.at(ii, jj) = S.at(i, j); + jj++; + } + ii++; + } + + // Solve the reduced system + cv::Mat delta_sub; + cv::solve(S_sub, e_sub, delta_sub, solveMethod); + + // Distribute the solution into the full update vector + deltaGlobal.setTo(0); + int j = 0; + for (int i = 0; i < n_global; i++) + { + if (mask_global[i]) + deltaGlobal.at(i) = delta_sub.at(j++); + } + + // Back-substitute to find the update for extrinsic parameters + for (int i = 0; i < nimages; i++) + { + deltaLocal[i] = V_inv[i] * (eb[i] - W[i].t() * deltaGlobal); + } + } + + bool iterate(double errNorm, bool& needsJacobian) + { + needsJacobian = false; + + if (state == DONE) + { + return false; + } + + if (state == STARTED) + { + state = CALC_J; + needsJacobian = true; + return true; + } + + if (state == CALC_J) + { + prevErrNorm = errNorm; + state = CHECK_ERR; + return true; + } + + // state == CHECK_ERR + if (errNorm > prevErrNorm) + { + // Reject step + if (++lambdaLg10 <= 16) + { + // Re-solve with higher lambda, DON'T recompute Jacobian + return true; + } + // Lambda exhausted + } + + // Accept step + lambdaLg10 = std::max(lambdaLg10 - 1, -16); + prevErrNorm = errNorm; + + if (++iters >= criteria.maxCount) + { + state = DONE; + return false; + } + + state = CALC_J; + needsJacobian = true; + return true; + } + + cv::Mat getGlobalUpdate() const { return deltaGlobal; } + cv::Mat getLocalUpdate(int i) const { return deltaLocal[i]; } + void setSolveMethod(int method) { solveMethod = method; } + + cv::Mat U; + std::vector V; + std::vector W; + cv::Mat ea; + std::vector eb; + double currentError; + + State state; + int NINTRINSIC; + int nimages; + int n_global; + int lambdaLg10; + +private: + int iters; + double prevErrNorm; + cv::TermCriteria criteria; + int solveMethod; + + std::vector V_inv; + cv::Mat deltaGlobal; + std::vector deltaLocal; +}; + +class JAccumulator : public cv::ParallelLoopBody +{ +public: + JAccumulator(SchurLMSolver& _solver, + const cv::Mat& _matM, const cv::Mat& _m, const cv::Mat& _npoints, + const cv::Mat& _param, int _flags, double _aspectRatio, + int _NINTRINSIC, bool _releaseObject, int _maxPoints, + cv::Mutex& _globalMutex) + : solver(_solver), matM(_matM), m(_m), npoints(_npoints), + param(_param), flags(_flags), aspectRatio(_aspectRatio), + NINTRINSIC(_NINTRINSIC), releaseObject(_releaseObject), + maxPoints(_maxPoints), globalMutex(_globalMutex) {} + + void operator()(const cv::Range& range) const CV_OVERRIDE + { + int nimages_total = npoints.checkVector(1, CV_32S); + int so = NINTRINSIC + nimages_total * 6; + + double fx = param.at(0), fy = param.at(1); + double cx = param.at(2), cy = param.at(3); + cv::Matx33d intrin(fx, 0, cx, 0, fy, cy, 0, 0, 1); + cv::Mat dist = param.rowRange(4, 4 + 14); + + cv::Mat JiBuf(maxPoints * 2, NINTRINSIC, CV_64F); + cv::Mat JeBuf(maxPoints * 2, 6, CV_64F); + cv::Mat JoBuf; + if (releaseObject) + JoBuf.create(maxPoints * 2, maxPoints * 3, CV_64F); + cv::Mat errBuf(maxPoints * 2, 1, CV_64F); + + int pos = 0; + for (int i = 0; i < range.start; i++) + pos += npoints.at(i); + + double localErr = 0; + + for (int i = range.start; i < range.end; i++) + { + int ni = npoints.at(i); + int si = NINTRINSIC + i * 6; + + cv::Mat _ri = param.rowRange(si, si + 3); + cv::Mat _ti = param.rowRange(si + 3, si + 6); + + cv::Mat _Mi = matM.colRange(pos, pos + ni); + if (releaseObject) + { + _Mi = param.rowRange(so, so + ni * 3); + _Mi = _Mi.reshape(3, 1); + } + cv::Mat _mi = m.colRange(pos, pos + ni); + + // Buffers for derivatives + cv::Mat Ji = JiBuf.rowRange(0, ni * 2); + Ji.setTo(0); + cv::Mat Je = JeBuf.rowRange(0, ni * 2); + Je.setTo(0); + cv::Mat Jo; + if (releaseObject) + { + Jo = JoBuf.rowRange(0, ni * 2).colRange(0, ni * 3); + Jo.setTo(0); + } + cv::Mat err = errBuf.rowRange(0, ni * 2); + cv::Mat _mp = err.reshape(2, 1); + + cv::Mat _dpdr, _dpdt, _dpdf, _dpdc, _dpdk, _dpdo; + + // Allocate contiguous buffers for Jacobians + _dpdr.create(ni * 2, 3, CV_64F); + _dpdt.create(ni * 2, 3, CV_64F); + _dpdk.create(ni * 2, NINTRINSIC - 4, CV_64F); + + if (!(flags & cv::CALIB_FIX_FOCAL_LENGTH)) + _dpdf.create(ni * 2, 2, CV_64F); + if (!(flags & cv::CALIB_FIX_PRINCIPAL_POINT)) + _dpdc.create(ni * 2, 2, CV_64F); + if (releaseObject) + _dpdo.create(ni * 2, ni * 3, CV_64F); + + cv::projectPoints(_Mi, _ri, _ti, intrin, dist, _mp, + _dpdr, _dpdt, + _dpdf.empty() ? cv::noArray() : _dpdf, + _dpdc.empty() ? cv::noArray() : _dpdc, + _dpdk, + _dpdo.empty() ? cv::noArray() : _dpdo, + (flags & cv::CALIB_FIX_ASPECT_RATIO) ? aspectRatio : 0.); + + // Copy Jacobians from projectPoints buffers to solver blocks + _dpdr.copyTo(Je.colRange(0, 3)); + _dpdt.copyTo(Je.colRange(3, 6)); + + if (!_dpdf.empty()) + { + if (flags & cv::CALIB_FIX_ASPECT_RATIO) + { + cv::Mat dpdf_fx = _dpdf.col(0); + cv::Mat dpdf_fy = _dpdf.col(1); + cv::scaleAdd(dpdf_fx, aspectRatio, dpdf_fy, dpdf_fy); + } + _dpdf.copyTo(Ji.colRange(0, 2)); + } + if (!_dpdc.empty()) + _dpdc.copyTo(Ji.colRange(2, 4)); + + _dpdk.copyTo(Ji.colRange(4, NINTRINSIC)); + + if (!_dpdo.empty()) + _dpdo.copyTo(Jo); + + cv::subtract(_mp, _mi, _mp); + localErr += cv::norm(err, cv::NORM_L2SQR); + + // Accumulate V and eb blocks (per-image) + solver.V[i] = Je.t() * Je; + solver.eb[i] = Je.t() * err; + cv::Mat JitJe = Ji.t() * Je; + JitJe.copyTo(solver.W[i].rowRange(0, NINTRINSIC)); + + if (releaseObject) + { + cv::Mat JotJe = Jo.t() * Je; + JotJe.copyTo(solver.W[i].rowRange(NINTRINSIC, NINTRINSIC + Jo.cols)); + } + + // Accumulate U and ea blocks (shared parameters) + { + cv::AutoLock lock(globalMutex); + solver.U(cv::Rect(0, 0, NINTRINSIC, NINTRINSIC)) += Ji.t() * Ji; + solver.ea.rowRange(0, NINTRINSIC) += Ji.t() * err; + + if (releaseObject) + { + cv::Mat JitJo = Ji.t() * Jo; + solver.U(cv::Rect(NINTRINSIC, 0, Jo.cols, NINTRINSIC)) += JitJo; + solver.U(cv::Rect(0, NINTRINSIC, NINTRINSIC, Jo.cols)) += JitJo.t(); + solver.U(cv::Rect(NINTRINSIC, NINTRINSIC, Jo.cols, Jo.cols)) += Jo.t() * Jo; + solver.ea.rowRange(NINTRINSIC, NINTRINSIC + Jo.cols) += Jo.t() * err; + } + + solver.currentError += localErr; + } + + pos += ni; + localErr = 0; + } + } + +private: + SchurLMSolver& solver; + const cv::Mat& matM; + const cv::Mat& m; + const cv::Mat& npoints; + const cv::Mat& param; + int flags; + double aspectRatio; + int NINTRINSIC; + bool releaseObject; + int maxPoints; + cv::Mutex& globalMutex; +}; + +} // anonymous namespace + + /* This is straight-forward port v3 of Matlab calibration engine by Jean-Yves Bouguet that is (in a large extent) based on the paper: @@ -144,7 +557,6 @@ static void subMatrix(const Mat& src, Mat& dst, CV_Assert(src.type() == CV_64F && dst.type() == CV_64F); int m = (int)rows.size(), n = (int)cols.size(); int i1 = 0, j1 = 0; - const uchar* colsdata = cols.empty() ? 0 : &cols[0]; for(int i = 0; i < m; i++) { if(rows[i]) @@ -152,20 +564,31 @@ static void subMatrix(const Mat& src, Mat& dst, const double* srcptr = src.ptr(i); double* dstptr = dst.ptr(i1++); - for(int j = j1 = 0; j < n; j++) + if(n > 0) { - if(colsdata[j]) - dstptr[j1++] = srcptr[j]; + const uchar* colsdata = cols.data(); + for(int j = j1 = 0; j < n; j++) + { + if(colsdata[j]) + dstptr[j1++] = srcptr[j]; + } } } } } -static double calibrateCameraInternal( const Mat& objectPoints, - const Mat& imagePoints, const Mat& npoints, - Size imageSize, int iFixedPoint, Mat& cameraMatrix, Mat& distCoeffs, - Mat rvecs, Mat tvecs, Mat newObjPoints, Mat stdDevs, - Mat perViewErr, int flags, const TermCriteria& termCrit ) +/* + T his is straight-forward port v3 of Matlab calibration engine by Jean-Yves Bouguet * + that is (in a large extent) based on the paper: + Z. Zhang. "A flexible new technique for camera calibration". + IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(11):1330-1334, 2000. + The 1st initial port was done by Valery Mosyagin. + */ +static double calibrateCameraInternalBouguet( const Mat& objectPoints, + const Mat& imagePoints, const Mat& npoints, + Size imageSize, int iFixedPoint, Mat& cameraMatrix, Mat& distCoeffs, + Mat rvecs, Mat tvecs, Mat newObjPoints, Mat stdDevs, + Mat perViewErr, int flags, const TermCriteria& termCrit ) { int NINTRINSIC = CALIB_NINTRINSIC; @@ -636,6 +1059,678 @@ static double calibrateCameraInternal( const Mat& objectPoints, return std::sqrt(reprojErr/total); } +static double calibrateCameraInternalSchur( const Mat& objectPoints, + const Mat& imagePoints, const Mat& npoints, + Size imageSize, int iFixedPoint, Mat& cameraMatrix, Mat& distCoeffs, + Mat rvecs, Mat tvecs, Mat newObjPoints, Mat stdDevs, + Mat perViewErr, int flags, const TermCriteria& termCrit ) +{ + int NINTRINSIC = CALIB_NINTRINSIC; + double reprojErr = 0; + + double aspectRatio = 0.; + int nimages = npoints.checkVector(1, CV_32S); + CV_Assert(nimages >= 1); + int ndistCoeffs = (int)distCoeffs.total(); + bool releaseObject = iFixedPoint > 0 && iFixedPoint < npoints.at(0) - 1; + + // 0. check the parameters & allocate buffers + if (imageSize.width <= 0 || imageSize.height <= 0) + { + CV_Error(cv::Error::StsOutOfRange, "image width and height must be positive"); + } + + if (flags & CALIB_TILTED_MODEL) + { + // when the tilted sensor model is used the distortion coefficients matrix must have 14 parameters + if (ndistCoeffs != 14) + { + CV_Error(cv::Error::StsBadArg, "The tilted sensor model must have 14 parameters"); + } + } + else + { + // when the thin prism model is used the distortion coefficients matrix must have 12 parameters + if (flags & CALIB_THIN_PRISM_MODEL) + { + if (ndistCoeffs != 12) + { + CV_Error(cv::Error::StsBadArg, "Thin prism model must have 12 parameters"); + } + } + } + + if (!rvecs.empty()) + { + int cn = rvecs.channels(); + CV_Assert(rvecs.depth() == CV_32F || rvecs.depth() == CV_64F); + CV_Assert((rvecs.rows == nimages && (rvecs.cols * cn == 3)) || + (rvecs.rows == 1 && rvecs.cols == nimages && cn == 3)); + } + + if (!tvecs.empty()) + { + int cn = tvecs.channels(); + CV_Assert(tvecs.depth() == CV_32F || tvecs.depth() == CV_64F); + CV_Assert((tvecs.rows == nimages && tvecs.cols * cn == 3) || + (tvecs.rows == 1 && tvecs.cols == nimages && cn == 3)); + } + + CV_Assert(cameraMatrix.type() == CV_32F || cameraMatrix.type() == CV_64F); + CV_Assert(cameraMatrix.rows == 3 && cameraMatrix.cols == 3); + + CV_Assert(distCoeffs.type() == CV_32F || distCoeffs.type() == CV_64F); + CV_Assert(distCoeffs.rows == 1 || distCoeffs.cols == 1); + CV_Assert(ndistCoeffs == 4 || ndistCoeffs == 5 || ndistCoeffs == 8 || + ndistCoeffs == 12 || ndistCoeffs == 14); + + int total = 0, maxPoints = 0; + for(int i = 0; i < nimages; i++ ) + { + int ni = npoints.at(i); + if( ni < 4 ) + { + CV_Error_( cv::Error::StsOutOfRange, ("The number of points in the view #%d is < 4", i)); + } + maxPoints = MAX( maxPoints, ni ); + total += ni; + } + + if( !newObjPoints.empty() ) + { + int cn = newObjPoints.channels(); + CV_Assert(newObjPoints.depth() == CV_32F || newObjPoints.depth() == CV_64F); + CV_Assert((newObjPoints.rows == maxPoints && newObjPoints.cols*cn == 3) || + (newObjPoints.rows == 1 && newObjPoints.cols == maxPoints && cn == 3)); + } + + if( !stdDevs.empty() ) + { + int cn = stdDevs.channels(); + CV_Assert(stdDevs.depth() == CV_32F || stdDevs.depth() == CV_64F); + int nstddev = nimages*6 + NINTRINSIC + (releaseObject ? maxPoints*3 : 0); + + CV_Assert((stdDevs.rows == nstddev && stdDevs.cols*cn == 1) || + (stdDevs.rows == 1 && stdDevs.cols == nstddev && cn == 1)); + } + + Mat matM( 1, total, CV_64FC3 ); + Mat _m( 1, total, CV_64FC2 ); + Mat allErrors(1, total, CV_64FC2); + + if (objectPoints.channels() == 3) + { + objectPoints.convertTo(matM, CV_64F); + }else + { + convertPointsToHomogeneous(objectPoints, matM); + matM.convertTo(matM, CV_64F); + } + + if (imagePoints.channels() == 2) + { + imagePoints.convertTo(_m, CV_64F); + }else + { + convertPointsFromHomogeneous(imagePoints, _m); + _m.convertTo(_m, CV_64F); + } + + int nparams = NINTRINSIC + nimages*6; + if( releaseObject ) + nparams += maxPoints * 3; + + std::vector k(14, 0.0); + Mat _k( distCoeffs.rows, distCoeffs.cols, CV_64F, k.data()); + if( distCoeffs.total() < 8 ) + { + if( distCoeffs.total() < 5 ) + flags |= CALIB_FIX_K3; + flags |= CALIB_FIX_K4 | CALIB_FIX_K5 | CALIB_FIX_K6; + } + const double minValidAspectRatio = 0.01; + const double maxValidAspectRatio = 100.0; + + Matx33d A; + cameraMatrix.convertTo(A, CV_64F); + + // 1. initialize intrinsic parameters & LM solver + if( flags & CALIB_USE_INTRINSIC_GUESS ) + { + if( A(0, 0) <= 0 || A(1, 1) <= 0 ) + CV_Error( cv::Error::StsOutOfRange, "Focal length (fx and fy) must be positive" ); + if( A(0, 2) < 0 || A(0, 2) >= imageSize.width || + A(1, 2) < 0 || A(1, 2) >= imageSize.height ) + CV_Error( cv::Error::StsOutOfRange, "Principal point must be within the image" ); + if( fabs(A(0, 1)) > 1e-5 ) + CV_Error( cv::Error::StsOutOfRange, "Non-zero skew is not supported by the function" ); + if( fabs(A(1, 0)) > 1e-5 || fabs(A(2, 0)) > 1e-5 || + fabs(A(2, 1)) > 1e-5 || fabs(A(2,2)-1) > 1e-5 ) + CV_Error( cv::Error::StsOutOfRange, + "The intrinsic matrix must have [fx 0 cx; 0 fy cy; 0 0 1] shape" ); + A(0, 1) = A(1, 0) = A(2, 0) = A(2, 1) = 0.; + A(2, 2) = 1.; + + if (flags & CALIB_FIX_ASPECT_RATIO) + { + aspectRatio = A(0, 0) / A(1, 1); + if (aspectRatio < minValidAspectRatio || aspectRatio > maxValidAspectRatio) + { + CV_Error(cv::Error::StsOutOfRange, + "The specified aspect ratio (= cameraMatrix[0][0] / cameraMatrix[1][1]) is incorrect"); + } + } + distCoeffs.convertTo(_k, CV_64F); + }else + { + Scalar mean, sdv; + meanStdDev(matM, mean, sdv); + if (fabs(mean[2]) > 1e-5 || fabs(sdv[2]) > 1e-5) + { + CV_Error(cv::Error::StsBadArg, + "For non-planar calibration rigs the initial intrinsic matrix must be specified"); + } + for (int i = 0; i < total; i++) + { + matM.at(i).z = 0.; + } + + if (flags & CALIB_FIX_ASPECT_RATIO) + { + aspectRatio = A(0, 0); + aspectRatio /= A(1, 1); + if (aspectRatio < minValidAspectRatio || aspectRatio > maxValidAspectRatio) + { + CV_Error(cv::Error::StsOutOfRange, + "The specified aspect ratio (= cameraMatrix[0][0] / cameraMatrix[1][1]) is incorrect"); + } + } + initIntrinsicParams2D(matM, _m, npoints, imageSize, A, aspectRatio); + } + + SchurLMSolver solver; + solver.init(nimages, NINTRINSIC, releaseObject ? maxPoints : 0, termCrit); + if (flags & CALIB_USE_LU) + solver.setSolveMethod(cv::DECOMP_LU); + else if (flags & CALIB_USE_QR) + solver.setSolveMethod(cv::DECOMP_QR); + + Mat_ param_m(nparams, 1); + std::vector mask_vec(nparams, 1); + + param_m(0) = A(0, 0); param_m(1) = A(1, 1); param_m(2) = A(0, 2); param_m(3) = A(1, 2); + for (int i = 0; i < 14; i++) + param_m(4 + i) = k[i]; + + if (flags & CALIB_FIX_ASPECT_RATIO) + mask_vec[0] = 0; + if (flags & CALIB_FIX_FOCAL_LENGTH) + mask_vec[0] = mask_vec[1] = 0; + if (flags & CALIB_FIX_PRINCIPAL_POINT) + mask_vec[2] = mask_vec[3] = 0; + if (flags & CALIB_ZERO_TANGENT_DIST) + { + param_m(6) = param_m(7) = 0; + mask_vec[6] = mask_vec[7] = 0; + } + if (!(flags & CALIB_RATIONAL_MODEL)) + flags |= CALIB_FIX_K4 + CALIB_FIX_K5 + CALIB_FIX_K6; + if (!(flags & CALIB_THIN_PRISM_MODEL)) + flags |= CALIB_FIX_S1_S2_S3_S4; + if (!(flags & CALIB_TILTED_MODEL)) + flags |= CALIB_FIX_TAUX_TAUY; + + mask_vec[4] = !(flags & CALIB_FIX_K1); + mask_vec[5] = !(flags & CALIB_FIX_K2); + if (flags & CALIB_FIX_TANGENT_DIST) + mask_vec[6] = mask_vec[7] = 0; + mask_vec[8] = !(flags & CALIB_FIX_K3); + mask_vec[9] = !(flags & CALIB_FIX_K4); + mask_vec[10] = !(flags & CALIB_FIX_K5); + mask_vec[11] = !(flags & CALIB_FIX_K6); + + if (flags & CALIB_FIX_S1_S2_S3_S4) + { + mask_vec[12] = 0; + mask_vec[13] = 0; + mask_vec[14] = 0; + mask_vec[15] = 0; + } + + if (flags & CALIB_USE_INTRINSIC_GUESS) + { + CV_Assert(nparams == NINTRINSIC + nimages * 6); + } + + if (flags & CALIB_FIX_ASPECT_RATIO) + { + mask_vec[0] = 0; + if (aspectRatio <= 0) + CV_Error(Error::StsBadArg, "CALIB_FIX_ASPECT_RATIO is set but aspectRatio is not positive"); + } + + if (flags & CALIB_FIX_TAUX_TAUY) + { + mask_vec[16] = mask_vec[17] = 0; + } + + if (releaseObject) + { + int s = NINTRINSIC + nimages * 6; + std::copy(matM.ptr(), matM.ptr(0, maxPoints - 1) + 3, + param_m.ptr() + s); + mask_vec[s + 0] = 0; + mask_vec[s + 1] = 0; + mask_vec[s + 2] = 0; + mask_vec[s + iFixedPoint * 3 + 0] = 0; + mask_vec[s + iFixedPoint * 3 + 1] = 0; + mask_vec[s + iFixedPoint * 3 + 2] = 0; + mask_vec[nparams - 1] = 0; + } + + Mat mask(mask_vec); + int nparams_nz = countNonZero(mask); + + if (nparams_nz >= 2 * total) + { + CV_Error_(Error::StsBadArg, + ("There should be less vars to optimize (having %d) than the " + "number of residuals (%d = 2 per point)", nparams_nz, 2 * total)); + } + + // 2. initialize extrinsic parameters + for(int i = 0, pos = 0; i < nimages; i++ ) + { + int ni = npoints.at(i); + + int s = NINTRINSIC + i*6; + Mat _ri = param_m.rowRange(s, s + 3); + Mat _ti = param_m.rowRange(s + 3, s + 6); + + Mat _Mi = matM.colRange(pos, pos + ni); + Mat _mi = _m.colRange(pos, pos + ni); + + solvePnP(_Mi, _mi, A, _k, _ri, _ti, false); + + pos += ni; + } + + // 3. run the optimization + Mat_ prev_param = param_m.clone(); + Mutex globalMutex; + + // Compute initial error + Mat allErrorsBuf(1, total, CV_64FC2); + reprojErr = 0; + int so = NINTRINSIC + nimages * 6; + int pos = 0; + for (int i = 0; i < nimages; i++) + { + int ni = npoints.at(i); + int si = NINTRINSIC + i * 6; + + Mat _ri = param_m.rowRange(si, si + 3); + Mat _ti = param_m.rowRange(si + 3, si + 6); + + Mat _Mi = matM.colRange(pos, pos + ni); + if (releaseObject) + { + _Mi = param_m.rowRange(so, so + ni * 3); + _Mi = _Mi.reshape(3, 1); + } + Mat _mi = _m.colRange(pos, pos + ni); + Mat _me = allErrorsBuf.colRange(pos, pos + ni); + + double fx = param_m(0), fy = param_m(1), cx = param_m(2), cy = param_m(3); + if (flags & CALIB_FIX_ASPECT_RATIO) + fx = fy * aspectRatio; + Matx33d intrin(fx, 0, cx, 0, fy, cy, 0, 0, 1); + Mat dist = param_m.rowRange(4, 4 + 14); + + projectPoints(_Mi, _ri, _ti, intrin, dist, _me); + subtract(_me, _mi, _me); + + reprojErr += norm(_me, NORM_L2SQR); + pos += ni; + } + + bool recomputeFinalErrors = (termCrit.maxCount == 0); + bool jacobianAtCurrentParams = false; + if (termCrit.maxCount > 0) + { + solver.reset(); + + parallel_for_(Range(0, nimages), + JAccumulator(solver, matM, _m, npoints, param_m, + flags, aspectRatio, NINTRINSIC, + releaseObject, maxPoints, globalMutex)); + jacobianAtCurrentParams = true; + // JAccumulator acc(solver, matM, _m, npoints, param_m, + // flags, aspectRatio, NINTRINSIC, + // releaseObject, maxPoints, globalMutex); + // acc(Range(0, nimages)); + solver.step(mask); + + double prevErr = reprojErr; + // Apply the initial step + param_m.copyTo(prev_param); + // 1. Global parameters (Intrinsics + Objects) + for (int kk = 0; kk < NINTRINSIC; kk++) + param_m(kk) -= solver.getGlobalUpdate().at(kk); + + if (releaseObject) + { + int param_obj_start = NINTRINSIC + nimages * 6; + int num_obj_params = solver.n_global - NINTRINSIC; + for (int kk = 0; kk < num_obj_params; kk++) + param_m(param_obj_start + kk) -= solver.getGlobalUpdate().at(NINTRINSIC + kk); + } + for (int i = 0; i < nimages; i++) + { + int si = NINTRINSIC + i * 6; + for (int kk = 0; kk < 6; kk++) + param_m(si + kk) -= solver.getLocalUpdate(i).at(kk); + } + jacobianAtCurrentParams = false; + for (int iter = 0; iter < termCrit.maxCount; iter++) + { + if (flags & CALIB_FIX_ASPECT_RATIO) + { + param_m(0) = param_m(1) * aspectRatio; + } + + reprojErr = 0; + int pos_iter = 0; + for (int i = 0; i < nimages; i++) + { + int ni = npoints.at(i); + int si = NINTRINSIC + i * 6; + + Mat _ri = param_m.rowRange(si, si + 3); + Mat _ti = param_m.rowRange(si + 3, si + 6); + + Mat _Mi = matM.colRange(pos_iter, pos_iter + ni); + if (releaseObject) + { + _Mi = param_m.rowRange(so, so + ni * 3).reshape(3, 1); + } + Mat _mi = _m.colRange(pos_iter, pos_iter + ni); + Mat _me = allErrorsBuf.colRange(pos_iter, pos_iter + ni); + + double fx = param_m(0), fy = param_m(1), cx = param_m(2), cy = param_m(3); + Matx33d intrin(fx, 0, cx, 0, fy, cy, 0, 0, 1); + Mat dist = param_m.rowRange(4, 4 + 14); + + projectPoints(_Mi, _ri, _ti, intrin, dist, _me); + subtract(_me, _mi, _me); + + double viewErr = norm(_me, NORM_L2SQR); + if (!perViewErr.empty()) + { + perViewErr.at(i) = std::sqrt(viewErr / ni); + } + reprojErr += viewErr; + pos_iter += ni; + } + + // Accept or reject the step + if (reprojErr < prevErr) + { + // Step accepted + solver.lambdaLg10 = std::max(solver.lambdaLg10 - 1, -16); + prevErr = reprojErr; + + // Check convergence BEFORE saving new prev_param + double paramChange = norm(param_m, prev_param, NORM_L2) / + (norm(prev_param, NORM_L2) + DBL_EPSILON); + param_m.copyTo(prev_param); + if (paramChange < termCrit.epsilon) + { + break; + } + + // Compute new Jacobian at new position + solver.reset(); + parallel_for_(Range(0, nimages), + JAccumulator(solver, matM, _m, npoints, param_m, + flags, aspectRatio, NINTRINSIC, + releaseObject, maxPoints, globalMutex)); + jacobianAtCurrentParams = true; + }else + { + // Step rejected, increase lambda and recompute step + if (++solver.lambdaLg10 > 16) + { + // Lambda exhausted, restore best params + prev_param.copyTo(param_m); + recomputeFinalErrors = true; + break; + } + + // Revert params + prev_param.copyTo(param_m); + } + + // Compute step (possibly with new lambda) + solver.step(mask); + + // Apply step. Update layouts: + // deltaGlobal: [Intrinsics (0..NINTRINSIC-1) | Object points (NINTRINSIC..n_global-1)] + // param_m: [Intrinsics | Extrinsic parameters | Object points] + + // Apply Intrinsics update + for (int kk = 0; kk < NINTRINSIC; kk++) + { + param_m(kk) -= solver.getGlobalUpdate().at(kk); + } + + // Apply Objects update (if any) + if (releaseObject) + { + int param_obj_start = NINTRINSIC + nimages * 6; + int num_obj_params = solver.n_global - NINTRINSIC; + for (int kk = 0; kk < num_obj_params; kk++) + { + param_m(param_obj_start + kk) -= solver.getGlobalUpdate().at(NINTRINSIC + kk); + } + } + + // Apply Local parameters (Extrinsics) + for (int i = 0; i < nimages; i++) + { + int si = NINTRINSIC + i * 6; + for (int kk = 0; kk < 6; kk++) + param_m(si + kk) -= solver.getLocalUpdate(i).at(kk); + } + jacobianAtCurrentParams = false; + } + } + + if (flags & CALIB_FIX_ASPECT_RATIO) + { + param_m(0) = param_m(1) * aspectRatio; + } + + if (recomputeFinalErrors) + { + Mat finalErrorsBuf; + Mat* errorsPtr = &allErrorsBuf; + if (!stdDevs.empty()) + { + finalErrorsBuf.create(1, total, CV_64FC2); + errorsPtr = &finalErrorsBuf; + } + + Mat& errors = *errorsPtr; + reprojErr = 0; + pos = 0; + for (int i = 0; i < nimages; i++) + { + int ni = npoints.at(i); + int si = NINTRINSIC + i * 6; + + Mat _ri = param_m.rowRange(si, si + 3); + Mat _ti = param_m.rowRange(si + 3, si + 6); + + Mat _Mi = matM.colRange(pos, pos + ni); + if (releaseObject) + { + _Mi = param_m.rowRange(so, so + ni * 3).reshape(3, 1); + } + Mat _mi = _m.colRange(pos, pos + ni); + Mat _me = errors.colRange(pos, pos + ni); + + double fx = param_m(0), fy = param_m(1), cx = param_m(2), cy = param_m(3); + Matx33d intrin(fx, 0, cx, 0, fy, cy, 0, 0, 1); + Mat dist = param_m.rowRange(4, 4 + 14); + + projectPoints(_Mi, _ri, _ti, intrin, dist, _me); + subtract(_me, _mi, _me); + + double viewErr = norm(_me, NORM_L2SQR); + if (!perViewErr.empty()) + { + perViewErr.at(i) = std::sqrt(viewErr / ni); + } + reprojErr += viewErr; + pos += ni; + } + } + + if (!stdDevs.empty()) + { + if (!jacobianAtCurrentParams) + { + solver.reset(); + parallel_for_(Range(0, nimages), + JAccumulator(solver, matM, _m, npoints, param_m, + flags, aspectRatio, NINTRINSIC, + releaseObject, maxPoints, globalMutex)); + jacobianAtCurrentParams = true; + } + Mat JtJ = Mat::zeros(nparams, nparams, CV_64F); + + // Manual copy for U to avoid assertion + Mat srcU = solver.U(Rect(0, 0, NINTRINSIC, NINTRINSIC)); + Mat dstU = JtJ(Rect(0, 0, NINTRINSIC, NINTRINSIC)); + for(int r=0; r(r); + double* dptr = dstU.ptr(r); + for(int c=0; c(r); + double* dptr = dstV.ptr(r); + for(int c=0; c mask_copy(mask.ptr(), mask.ptr() + nparams); + subMatrix(JtJ, JtJN, mask_copy, mask_copy); + completeSymm(JtJN, false); + cv::invert(JtJN, JtJinv, DECOMP_SVD); + + int nErrors = 2 * total - nparams_nz; + double sigma2 = norm(allErrorsBuf, NORM_L2SQR) / nErrors; + + int j = 0; + for (int s = 0; s < nparams; s++) + { + stdDevs.at(s) = mask.at(s) ? + std::sqrt(JtJinv.at(j, j) * sigma2) : 0.0; + if (mask.at(s)) + { + j++; + } + } + } + + // Store optimization results + A = Matx33d(param_m(0), 0, param_m(2), 0, param_m(1), param_m(3), 0, 0, 1); + Mat(A).convertTo(cameraMatrix, cameraMatrix.type()); + + _k = Mat(distCoeffs.size(), CV_64F, param_m.ptr() + 4); + _k.convertTo(distCoeffs, distCoeffs.type()); + + + if (!newObjPoints.empty() && releaseObject) + { + int s = NINTRINSIC + nimages * 6; + Mat _Mi = param_m.rowRange(s, s + maxPoints * 3); + _Mi.reshape(3, 1).convertTo(newObjPoints, newObjPoints.type()); + + } + for (int i = 0; i < nimages; i++) + { + if (!rvecs.empty()) + { + Mat src = param_m.rowRange(NINTRINSIC + i * 6, NINTRINSIC + i * 6 + 3); + if (rvecs.rows == nimages && rvecs.cols * rvecs.channels() == 9) + { + Mat dst(3, 3, rvecs.depth(), rvecs.ptr(i)); + Rodrigues(src, A); + Mat(A).convertTo(dst, dst.type()); + } + else + { + Mat dst(3, 1, rvecs.depth(), rvecs.rows == 1 ? + rvecs.data + i * rvecs.elemSize() : rvecs.ptr(i)); + src.convertTo(dst, dst.type()); + } + } + if (!tvecs.empty()) + { + Mat src = param_m.rowRange(NINTRINSIC + i * 6 + 3, NINTRINSIC + i * 6 + 6); + Mat dst(3, 1, tvecs.depth(), tvecs.rows == 1 ? + tvecs.data + i * tvecs.elemSize() : tvecs.ptr(i)); + src.convertTo(dst, dst.type()); + } + } + return std::sqrt(reprojErr/total); +} + //////////////////////////////// Stereo Calibration /////////////////////////////////// static double stereoCalibrateImpl( @@ -731,9 +1826,14 @@ static double stereoCalibrateImpl( { Mat mIntr(A[k], /* copyData = */ false); Mat mDist(distInitial[k], /* copyData = */ false); - calibrateCameraInternal(objectPoints, imagePoints[k], - _npoints, imageSize, 0, mIntr, mDist, - Mat(), Mat(), Mat(), Mat(), Mat(), flags, termCrit); + if( flags & CALIB_DISABLE_SCHUR_COMPLEMENT ) + calibrateCameraInternalBouguet(objectPoints, imagePoints[k], + _npoints, imageSize, 0, mIntr, mDist, + Mat(), Mat(), Mat(), Mat(), Mat(), flags, termCrit); + else + calibrateCameraInternalSchur(objectPoints, imagePoints[k], + _npoints, imageSize, 0, mIntr, mDist, + Mat(), Mat(), Mat(), Mat(), Mat(), flags, termCrit); } } @@ -1852,13 +2952,27 @@ double calibrateCameraRO(InputArrayOfArrays _objectPoints, errorsM = _perViewErrors.getMat(); } - double reprojErr = calibrateCameraInternal( + double reprojErr = 0; + if( flags & CALIB_DISABLE_SCHUR_COMPLEMENT ) + { + reprojErr = calibrateCameraInternalBouguet( objPt, imgPt, npoints, imageSize, iFixedPoint, cameraMatrix, distCoeffs, rvecM, tvecM, newObjPt, stdDeviationsM, errorsM, flags, criteria); + } + else + { + reprojErr = calibrateCameraInternalSchur( + objPt, imgPt, npoints, imageSize, iFixedPoint, + cameraMatrix, distCoeffs, + rvecM, tvecM, + newObjPt, + stdDeviationsM, + errorsM, flags, criteria); + } if( stddev_needed ) { @@ -1884,13 +2998,13 @@ double calibrateCameraRO(InputArrayOfArrays _objectPoints, { _rvecs.create(3, 1, CV_64F, i, true); Mat rv = _rvecs.getMat(i); - memcpy(rv.ptr(), rvecM.ptr(i), 3*sizeof(double)); + Mat(rvecM.row(i).t()).copyTo(rv); } if( tvecs_needed && tvecs_mat_vec) { _tvecs.create(3, 1, CV_64F, i, true); Mat tv = _tvecs.getMat(i); - memcpy(tv.ptr(), tvecM.ptr(i), 3*sizeof(double)); + Mat(tvecM.row(i).t()).copyTo(tv); } } diff --git a/modules/calib/test/test_multiview_calib.cpp b/modules/calib/test/test_multiview_calib.cpp index 0df18e12f1..b8c94b9783 100644 --- a/modules/calib/test/test_multiview_calib.cpp +++ b/modules/calib/test/test_multiview_calib.cpp @@ -312,7 +312,7 @@ TEST_F(MultiViewTest, OneLine) std::vector board_pattern = genAsymmetricObjectPoints(); std::vector> objPoints(num_frames, board_pattern); - std::vector flagsForIntrinsics(3, CALIB_RATIONAL_MODEL); + std::vector flagsForIntrinsics(3, CALIB_RATIONAL_MODEL | CALIB_DISABLE_SCHUR_COMPLEMENT); std::vector Ks, distortions, Rs, Rs_rvec, Ts; double rms = calibrateMultiview(objPoints, image_points_all, image_sizes, visibility, models, @@ -334,7 +334,6 @@ TEST_F(MultiViewTest, OneLine) TEST_F(MultiViewTest, OneLineInitialGuess) { - applyTestTag(CV_TEST_TAG_VERYLONG); const string root = cvtest::TS::ptr()->get_data_path() + "cv/cameracalibration/multiview/3cams-one-line/"; const std::vector cam_names = {"cam_0", "cam_1", "cam_3"}; const std::vector image_sizes = {{1920, 1080}, {1920, 1080}, {1920, 1080} }; @@ -389,7 +388,7 @@ TEST_F(MultiViewTest, OneLineInitialGuess) { Mat K, dist; double mono_rms = calibrateMono(board_pattern, image_points_all[c], image_sizes[c], - cv::CALIB_MODEL_PINHOLE, cv::CALIB_RATIONAL_MODEL, + cv::CALIB_MODEL_PINHOLE, cv::CALIB_RATIONAL_MODEL | CALIB_DISABLE_SCHUR_COMPLEMENT, K, dist); CV_LOG_INFO(NULL, "K:" << K); @@ -440,7 +439,6 @@ TEST_F(MultiViewTest, OneLineInitialGuess) TEST_F(MultiViewTest, CamsToFloor) { - applyTestTag(CV_TEST_TAG_VERYLONG); const string root = cvtest::TS::ptr()->get_data_path() + "cv/cameracalibration/multiview/3cams-to-floor/"; const std::vector cam_names = {"cam_0", "cam_1", "cam_2"}; std::vector image_sizes = {{1920, 1080}, {1920, 1080}, {1280, 720}}; @@ -508,7 +506,6 @@ TEST_F(MultiViewTest, CamsToFloor) TEST_F(MultiViewTest, Hetero) { - applyTestTag(CV_TEST_TAG_VERYLONG); const string root = cvtest::TS::ptr()->get_data_path() + "cv/cameracalibration/multiview/3cams-hetero/"; const std::vector cam_names = {"cam_7", "cam_4", "cam_8"}; std::vector image_sizes = {{1920, 1080}, {1920, 1080}, {2048, 2048}}; diff --git a/modules/core/include/opencv2/core/cv_cpu_dispatch.h b/modules/core/include/opencv2/core/cv_cpu_dispatch.h index 345a5fa18c..0136146edb 100644 --- a/modules/core/include/opencv2/core/cv_cpu_dispatch.h +++ b/modules/core/include/opencv2/core/cv_cpu_dispatch.h @@ -140,7 +140,7 @@ #endif #if defined _WIN32 && (defined(_M_ARM) || defined(_M_ARM64) || defined(_M_ARM64EC)) && (defined(CV_CPU_COMPILE_NEON) || !defined(_MSC_VER)) -# include +# include # include # define CV_NEON 1 #elif defined(__ARM_NEON) @@ -233,7 +233,7 @@ struct VZeroUpperGuard { # define CV_SSE 1 # define CV_SSE2 1 #elif defined _WIN32 && (defined(_M_ARM) || defined(_M_ARM64) || defined(_M_ARM64EC)) && (defined(CV_CPU_COMPILE_NEON) || !defined(_MSC_VER)) -# include +# include # include # define CV_NEON 1 #elif defined(__ARM_NEON) diff --git a/modules/core/include/opencv2/core/hal/intrin_wasm.hpp b/modules/core/include/opencv2/core/hal/intrin_wasm.hpp index 7c4d8e05df..f21cca7aea 100644 --- a/modules/core/include/opencv2/core/hal/intrin_wasm.hpp +++ b/modules/core/include/opencv2/core/hal/intrin_wasm.hpp @@ -1859,9 +1859,12 @@ inline _Tpwvec v_expand_high(const _Tpvec& a) \ { return _Tpwvec(__CV_CAT(intrin, _high)(a.val)); } \ inline _Tpwvec v_load_expand(const _Tp* ptr) \ { \ - v128_t a = wasm_v128_load(ptr); \ - return _Tpwvec(intrin(a)); \ -} + using lane_t = typename _Tpwvec::lane_type; \ + alignas(16) lane_t tmp[_Tpwvec::nlanes]; \ + for(int i = 0; i < _Tpwvec::nlanes; i++) \ + tmp[i] = static_cast(ptr[i]); \ + return _Tpwvec(wasm_v128_load(tmp)); \ +} \ OPENCV_HAL_IMPL_WASM_EXPAND(v_uint8x16, v_uint16x8, uchar, v128_cvtu8x16_i16x8) OPENCV_HAL_IMPL_WASM_EXPAND(v_int8x16, v_int16x8, schar, v128_cvti8x16_i16x8) @@ -1873,9 +1876,12 @@ OPENCV_HAL_IMPL_WASM_EXPAND(v_int32x4, v_int64x2, int, v128_cvti32x4_i64x2) #define OPENCV_HAL_IMPL_WASM_EXPAND_Q(_Tpvec, _Tp, intrin) \ inline _Tpvec v_load_expand_q(const _Tp* ptr) \ { \ - v128_t a = wasm_v128_load(ptr); \ - return _Tpvec(intrin(a)); \ -} + using lane_t = typename _Tpvec::lane_type; \ + alignas(16) lane_t tmp[_Tpvec::nlanes]; \ + for(int i = 0; i < _Tpvec::nlanes; i++) \ + tmp[i] = static_cast(ptr[i]); \ + return _Tpvec(wasm_v128_load(tmp)); \ +} \ OPENCV_HAL_IMPL_WASM_EXPAND_Q(v_uint32x4, uchar, v128_cvtu8x16_i32x4) OPENCV_HAL_IMPL_WASM_EXPAND_Q(v_int32x4, schar, v128_cvti8x16_i32x4) diff --git a/modules/core/perf/perf_stat.cpp b/modules/core/perf/perf_stat.cpp index 025700c989..8846475844 100644 --- a/modules/core/perf/perf_stat.cpp +++ b/modules/core/perf/perf_stat.cpp @@ -29,7 +29,7 @@ PERF_TEST_P(Size_MatType, mean, TYPICAL_MATS) declare.in(src, WARMUP_RNG).out(s); - TEST_CYCLE() s = mean(src); + TEST_CYCLE() s = cv::mean(src); SANITY_CHECK(s, 1e-5); } @@ -45,7 +45,7 @@ PERF_TEST_P(Size_MatType, mean_mask, TYPICAL_MATS) declare.in(src, WARMUP_RNG).in(mask).out(s); - TEST_CYCLE() s = mean(src, mask); + TEST_CYCLE() s = cv::mean(src, mask); SANITY_CHECK(s, 5e-5); } diff --git a/modules/core/src/matrix_transform.cpp b/modules/core/src/matrix_transform.cpp index 3de1024ff1..cf7eddc3c5 100644 --- a/modules/core/src/matrix_transform.cpp +++ b/modules/core/src/matrix_transform.cpp @@ -13,11 +13,330 @@ namespace cv { ////////////////////////////////////// transpose ///////////////////////////////////////// +#if CV_SIMD128 +static void transpose_8bit_simd(const uchar* src, size_t sstep, uchar* dst, size_t dstep, Size sz) +{ + const int m = sz.width, n = sz.height; + int i = 0; + for (; i <= m - 16; i += 16) + { + int j = 0; + for (; j <= n - 16; j += 16) + { + v_uint8x16 r0 = v_load(src + i + sstep*(j+ 0)); + v_uint8x16 r1 = v_load(src + i + sstep*(j+ 1)); + v_uint8x16 r2 = v_load(src + i + sstep*(j+ 2)); + v_uint8x16 r3 = v_load(src + i + sstep*(j+ 3)); + v_uint8x16 r4 = v_load(src + i + sstep*(j+ 4)); + v_uint8x16 r5 = v_load(src + i + sstep*(j+ 5)); + v_uint8x16 r6 = v_load(src + i + sstep*(j+ 6)); + v_uint8x16 r7 = v_load(src + i + sstep*(j+ 7)); + v_uint8x16 r8 = v_load(src + i + sstep*(j+ 8)); + v_uint8x16 r9 = v_load(src + i + sstep*(j+ 9)); + v_uint8x16 r10 = v_load(src + i + sstep*(j+10)); + v_uint8x16 r11 = v_load(src + i + sstep*(j+11)); + v_uint8x16 r12 = v_load(src + i + sstep*(j+12)); + v_uint8x16 r13 = v_load(src + i + sstep*(j+13)); + v_uint8x16 r14 = v_load(src + i + sstep*(j+14)); + v_uint8x16 r15 = v_load(src + i + sstep*(j+15)); + + v_uint8x16 t0, t1, t2, t3, t4, t5, t6, t7, + t8, t9, t10, t11, t12, t13, t14, t15; + + v_zip(r0, r1, t0, t1); + v_zip(r2, r3, t2, t3); + v_zip(r4, r5, t4, t5); + v_zip(r6, r7, t6, t7); + v_zip(r8, r9, t8, t9); + v_zip(r10, r11, t10, t11); + v_zip(r12, r13, t12, t13); + v_zip(r14, r15, t14, t15); + + v_uint16x8 s0, s1, s2, s3, s4, s5, s6, s7, + s8, s9, s10, s11, s12, s13, s14, s15; + v_zip(v_reinterpret_as_u16(t0), v_reinterpret_as_u16(t2), s0, s1); + v_zip(v_reinterpret_as_u16(t1), v_reinterpret_as_u16(t3), s2, s3); + v_zip(v_reinterpret_as_u16(t4), v_reinterpret_as_u16(t6), s4, s5); + v_zip(v_reinterpret_as_u16(t5), v_reinterpret_as_u16(t7), s6, s7); + v_zip(v_reinterpret_as_u16(t8), v_reinterpret_as_u16(t10), s8, s9); + v_zip(v_reinterpret_as_u16(t9), v_reinterpret_as_u16(t11), s10, s11); + v_zip(v_reinterpret_as_u16(t12), v_reinterpret_as_u16(t14), s12, s13); + v_zip(v_reinterpret_as_u16(t13), v_reinterpret_as_u16(t15), s14, s15); + + v_uint32x4 u0, u1, u2, u3, u4, u5, u6, u7, + u8, u9, u10, u11, u12, u13, u14, u15; + + v_zip(v_reinterpret_as_u32(s0), v_reinterpret_as_u32(s4), u0, u1); + v_zip(v_reinterpret_as_u32(s1), v_reinterpret_as_u32(s5), u2, u3); + v_zip(v_reinterpret_as_u32(s2), v_reinterpret_as_u32(s6), u4, u5); + v_zip(v_reinterpret_as_u32(s3), v_reinterpret_as_u32(s7), u6, u7); + v_zip(v_reinterpret_as_u32(s8), v_reinterpret_as_u32(s12), u8, u9); + v_zip(v_reinterpret_as_u32(s9), v_reinterpret_as_u32(s13), u10, u11); + v_zip(v_reinterpret_as_u32(s10), v_reinterpret_as_u32(s14), u12, u13); + v_zip(v_reinterpret_as_u32(s11), v_reinterpret_as_u32(s15), u14, u15); + + v_uint32x4 v0 = v_combine_low (u0, u8); + v_uint32x4 v1 = v_combine_high(u0, u8); + v_uint32x4 v2 = v_combine_low (u1, u9); + v_uint32x4 v3 = v_combine_high(u1, u9); + v_uint32x4 v4 = v_combine_low (u2, u10); + v_uint32x4 v5 = v_combine_high(u2, u10); + v_uint32x4 v6 = v_combine_low (u3, u11); + v_uint32x4 v7 = v_combine_high(u3, u11); + v_uint32x4 v8 = v_combine_low (u4, u12); + v_uint32x4 v9 = v_combine_high(u4, u12); + v_uint32x4 v10 = v_combine_low (u5, u13); + v_uint32x4 v11 = v_combine_high(u5, u13); + v_uint32x4 v12 = v_combine_low (u6, u14); + v_uint32x4 v13 = v_combine_high(u6, u14); + v_uint32x4 v14 = v_combine_low (u7, u15); + v_uint32x4 v15 = v_combine_high(u7, u15); + + v_store(dst + dstep*(i+ 0) + j, v_reinterpret_as_u8(v0)); + v_store(dst + dstep*(i+ 1) + j, v_reinterpret_as_u8(v1)); + v_store(dst + dstep*(i+ 2) + j, v_reinterpret_as_u8(v2)); + v_store(dst + dstep*(i+ 3) + j, v_reinterpret_as_u8(v3)); + v_store(dst + dstep*(i+ 4) + j, v_reinterpret_as_u8(v4)); + v_store(dst + dstep*(i+ 5) + j, v_reinterpret_as_u8(v5)); + v_store(dst + dstep*(i+ 6) + j, v_reinterpret_as_u8(v6)); + v_store(dst + dstep*(i+ 7) + j, v_reinterpret_as_u8(v7)); + v_store(dst + dstep*(i+ 8) + j, v_reinterpret_as_u8(v8)); + v_store(dst + dstep*(i+ 9) + j, v_reinterpret_as_u8(v9)); + v_store(dst + dstep*(i+10) + j, v_reinterpret_as_u8(v10)); + v_store(dst + dstep*(i+11) + j, v_reinterpret_as_u8(v11)); + v_store(dst + dstep*(i+12) + j, v_reinterpret_as_u8(v12)); + v_store(dst + dstep*(i+13) + j, v_reinterpret_as_u8(v13)); + v_store(dst + dstep*(i+14) + j, v_reinterpret_as_u8(v14)); + v_store(dst + dstep*(i+15) + j, v_reinterpret_as_u8(v15)); + } + for (; j < n; j++) + for (int k = 0; k < 16; k++) + dst[dstep*(i+k) + j] = src[i + sstep*j + k]; + } + for (; i < m; i++) + for (int j = 0; j < n; j++) + dst[dstep*i + j] = src[i + sstep*j]; +} + +static void transpose_16bit_simd(const uchar* src, size_t sstep, uchar* dst, size_t dstep, Size sz) +{ + const ushort* src16 = reinterpret_cast(src); + ushort* dst16 = reinterpret_cast(dst); + + const size_t sstep_e = sstep / sizeof(ushort); + const size_t dstep_e = dstep / sizeof(ushort); + + const int m = sz.width, n = sz.height; + int i = 0; + + for (; i <= m - 8; i += 8) + { + int j = 0; + for (; j <= n - 8; j += 8) + { + v_uint16x8 r0 = v_load(src16 + i + sstep_e*(j+0)); + v_uint16x8 r1 = v_load(src16 + i + sstep_e*(j+1)); + v_uint16x8 r2 = v_load(src16 + i + sstep_e*(j+2)); + v_uint16x8 r3 = v_load(src16 + i + sstep_e*(j+3)); + v_uint16x8 r4 = v_load(src16 + i + sstep_e*(j+4)); + v_uint16x8 r5 = v_load(src16 + i + sstep_e*(j+5)); + v_uint16x8 r6 = v_load(src16 + i + sstep_e*(j+6)); + v_uint16x8 r7 = v_load(src16 + i + sstep_e*(j+7)); + + v_uint16x8 t0, t1, t2, t3, t4, t5, t6, t7; + v_zip(r0, r1, t0, t1); + v_zip(r2, r3, t2, t3); + v_zip(r4, r5, t4, t5); + v_zip(r6, r7, t6, t7); + v_uint32x4 u0, u1, u2, u3, u4, u5, u6, u7; + v_zip(v_reinterpret_as_u32(t0), v_reinterpret_as_u32(t4), u0, u1); + v_zip(v_reinterpret_as_u32(t1), v_reinterpret_as_u32(t5), u2, u3); + v_zip(v_reinterpret_as_u32(t2), v_reinterpret_as_u32(t6), u4, u5); + v_zip(v_reinterpret_as_u32(t3), v_reinterpret_as_u32(t7), u6, u7); + v_uint32x4 v0, v1, v2, v3, v4, v5, v6, v7; + v_zip(u0, u4, v0, v1); + v_zip(u1, u5, v2, v3); + v_zip(u2, u6, v4, v5); + v_zip(u3, u7, v6, v7); + + v_store(dst16 + dstep_e*(i+0) + j, v_reinterpret_as_u16(v0)); + v_store(dst16 + dstep_e*(i+1) + j, v_reinterpret_as_u16(v1)); + v_store(dst16 + dstep_e*(i+2) + j, v_reinterpret_as_u16(v2)); + v_store(dst16 + dstep_e*(i+3) + j, v_reinterpret_as_u16(v3)); + v_store(dst16 + dstep_e*(i+4) + j, v_reinterpret_as_u16(v4)); + v_store(dst16 + dstep_e*(i+5) + j, v_reinterpret_as_u16(v5)); + v_store(dst16 + dstep_e*(i+6) + j, v_reinterpret_as_u16(v6)); + v_store(dst16 + dstep_e*(i+7) + j, v_reinterpret_as_u16(v7)); + } + for (; j < n; j++) + for (int k = 0; k < 8; k++) + dst16[dstep_e*(i+k) + j] = src16[i + sstep_e*j + k]; + } + for (; i < m; i++) + for (int j = 0; j < n; j++) + dst16[dstep_e*i + j] = src16[i + sstep_e*j]; +} + +static void transpose_32bit_simd(const uchar* src, size_t sstep, uchar* dst, size_t dstep, Size sz) +{ + const uint32_t* src32 = reinterpret_cast(src); + uint32_t* dst32 = reinterpret_cast(dst); + + const size_t sstep_e = sstep / sizeof(uint32_t); + const size_t dstep_e = dstep / sizeof(uint32_t); + + const int m = sz.width, n = sz.height; + int i = 0; + for (; i <= m - 4; i += 4) + { + int j = 0; + for (; j <= n - 4; j += 4) + { + v_uint32x4 r0 = v_load(src32 + i + sstep_e*(j+0)); + v_uint32x4 r1 = v_load(src32 + i + sstep_e*(j+1)); + v_uint32x4 r2 = v_load(src32 + i + sstep_e*(j+2)); + v_uint32x4 r3 = v_load(src32 + i + sstep_e*(j+3)); + v_uint32x4 o0, o1, o2, o3; + v_transpose4x4(r0, r1, r2, r3, o0, o1, o2, o3); + + v_store(dst32 + dstep_e*(i+0) + j, o0); + v_store(dst32 + dstep_e*(i+1) + j, o1); + v_store(dst32 + dstep_e*(i+2) + j, o2); + v_store(dst32 + dstep_e*(i+3) + j, o3); + } + for (; j < n; j++) + for (int k = 0; k < 4; k++) + dst32[dstep_e*(i+k) + j] = src32[i + sstep_e*j + k]; + } + for (; i < m; i++) + for (int j = 0; j < n; j++) + dst32[dstep_e*i + j] = src32[i + sstep_e*j]; +} + +static void transpose_48bit_simd(const uchar* src, size_t sstep, uchar* dst, size_t dstep, Size sz) +{ + const short* src16 = reinterpret_cast(src); + short* dst16 = reinterpret_cast(dst); + + const size_t sstep_e = sstep / sizeof(short); + const size_t dstep_e = dstep / sizeof(short); + + const int m = sz.width, n = sz.height; + int i = 0; + + for (; i <= m - 8; i += 8) + { + int j = 0; + for (; j <= n - 8; j += 8) + { + v_int16x8 C0_0, C1_0, C2_0; + v_int16x8 C0_1, C1_1, C2_1; + v_int16x8 C0_2, C1_2, C2_2; + v_int16x8 C0_3, C1_3, C2_3; + v_int16x8 C0_4, C1_4, C2_4; + v_int16x8 C0_5, C1_5, C2_5; + v_int16x8 C0_6, C1_6, C2_6; + v_int16x8 C0_7, C1_7, C2_7; + + v_load_deinterleave(src16 + sstep_e*(j+0) + i*3, C0_0, C1_0, C2_0); + v_load_deinterleave(src16 + sstep_e*(j+1) + i*3, C0_1, C1_1, C2_1); + v_load_deinterleave(src16 + sstep_e*(j+2) + i*3, C0_2, C1_2, C2_2); + v_load_deinterleave(src16 + sstep_e*(j+3) + i*3, C0_3, C1_3, C2_3); + v_load_deinterleave(src16 + sstep_e*(j+4) + i*3, C0_4, C1_4, C2_4); + v_load_deinterleave(src16 + sstep_e*(j+5) + i*3, C0_5, C1_5, C2_5); + v_load_deinterleave(src16 + sstep_e*(j+6) + i*3, C0_6, C1_6, C2_6); + v_load_deinterleave(src16 + sstep_e*(j+7) + i*3, C0_7, C1_7, C2_7); + + v_uint16x8 t0, t1, t2, t3, t4, t5, t6, t7; + v_zip(v_reinterpret_as_u16(C0_0), v_reinterpret_as_u16(C0_1), t0, t1); + v_zip(v_reinterpret_as_u16(C0_2), v_reinterpret_as_u16(C0_3), t2, t3); + v_zip(v_reinterpret_as_u16(C0_4), v_reinterpret_as_u16(C0_5), t4, t5); + v_zip(v_reinterpret_as_u16(C0_6), v_reinterpret_as_u16(C0_7), t6, t7); + v_uint32x4 u0, u1, u2, u3, u4, u5, u6, u7; + v_zip(v_reinterpret_as_u32(t0), v_reinterpret_as_u32(t4), u0, u1); + v_zip(v_reinterpret_as_u32(t1), v_reinterpret_as_u32(t5), u2, u3); + v_zip(v_reinterpret_as_u32(t2), v_reinterpret_as_u32(t6), u4, u5); + v_zip(v_reinterpret_as_u32(t3), v_reinterpret_as_u32(t7), u6, u7); + v_uint32x4 s0, s1, s2, s3, s4, s5, s6, s7; + v_zip(u0, u4, s0, s1); v_zip(u1, u5, s2, s3); + v_zip(u2, u6, s4, s5); v_zip(u3, u7, s6, s7); + v_int16x8 r0_0 = v_reinterpret_as_s16(s0), r0_1 = v_reinterpret_as_s16(s1); + v_int16x8 r0_2 = v_reinterpret_as_s16(s2), r0_3 = v_reinterpret_as_s16(s3); + v_int16x8 r0_4 = v_reinterpret_as_s16(s4), r0_5 = v_reinterpret_as_s16(s5); + v_int16x8 r0_6 = v_reinterpret_as_s16(s6), r0_7 = v_reinterpret_as_s16(s7); + + v_zip(v_reinterpret_as_u16(C1_0), v_reinterpret_as_u16(C1_1), t0, t1); + v_zip(v_reinterpret_as_u16(C1_2), v_reinterpret_as_u16(C1_3), t2, t3); + v_zip(v_reinterpret_as_u16(C1_4), v_reinterpret_as_u16(C1_5), t4, t5); + v_zip(v_reinterpret_as_u16(C1_6), v_reinterpret_as_u16(C1_7), t6, t7); + v_zip(v_reinterpret_as_u32(t0), v_reinterpret_as_u32(t4), u0, u1); + v_zip(v_reinterpret_as_u32(t1), v_reinterpret_as_u32(t5), u2, u3); + v_zip(v_reinterpret_as_u32(t2), v_reinterpret_as_u32(t6), u4, u5); + v_zip(v_reinterpret_as_u32(t3), v_reinterpret_as_u32(t7), u6, u7); + v_zip(u0, u4, s0, s1); v_zip(u1, u5, s2, s3); + v_zip(u2, u6, s4, s5); v_zip(u3, u7, s6, s7); + v_int16x8 r1_0 = v_reinterpret_as_s16(s0), r1_1 = v_reinterpret_as_s16(s1); + v_int16x8 r1_2 = v_reinterpret_as_s16(s2), r1_3 = v_reinterpret_as_s16(s3); + v_int16x8 r1_4 = v_reinterpret_as_s16(s4), r1_5 = v_reinterpret_as_s16(s5); + v_int16x8 r1_6 = v_reinterpret_as_s16(s6), r1_7 = v_reinterpret_as_s16(s7); + + v_zip(v_reinterpret_as_u16(C2_0), v_reinterpret_as_u16(C2_1), t0, t1); + v_zip(v_reinterpret_as_u16(C2_2), v_reinterpret_as_u16(C2_3), t2, t3); + v_zip(v_reinterpret_as_u16(C2_4), v_reinterpret_as_u16(C2_5), t4, t5); + v_zip(v_reinterpret_as_u16(C2_6), v_reinterpret_as_u16(C2_7), t6, t7); + v_zip(v_reinterpret_as_u32(t0), v_reinterpret_as_u32(t4), u0, u1); + v_zip(v_reinterpret_as_u32(t1), v_reinterpret_as_u32(t5), u2, u3); + v_zip(v_reinterpret_as_u32(t2), v_reinterpret_as_u32(t6), u4, u5); + v_zip(v_reinterpret_as_u32(t3), v_reinterpret_as_u32(t7), u6, u7); + v_zip(u0, u4, s0, s1); v_zip(u1, u5, s2, s3); + v_zip(u2, u6, s4, s5); v_zip(u3, u7, s6, s7); + v_int16x8 r2_0 = v_reinterpret_as_s16(s0), r2_1 = v_reinterpret_as_s16(s1); + v_int16x8 r2_2 = v_reinterpret_as_s16(s2), r2_3 = v_reinterpret_as_s16(s3); + v_int16x8 r2_4 = v_reinterpret_as_s16(s4), r2_5 = v_reinterpret_as_s16(s5); + v_int16x8 r2_6 = v_reinterpret_as_s16(s6), r2_7 = v_reinterpret_as_s16(s7); + + v_store_interleave(dst16 + dstep_e*(i+0) + j*3, r0_0, r1_0, r2_0); + v_store_interleave(dst16 + dstep_e*(i+1) + j*3, r0_1, r1_1, r2_1); + v_store_interleave(dst16 + dstep_e*(i+2) + j*3, r0_2, r1_2, r2_2); + v_store_interleave(dst16 + dstep_e*(i+3) + j*3, r0_3, r1_3, r2_3); + v_store_interleave(dst16 + dstep_e*(i+4) + j*3, r0_4, r1_4, r2_4); + v_store_interleave(dst16 + dstep_e*(i+5) + j*3, r0_5, r1_5, r2_5); + v_store_interleave(dst16 + dstep_e*(i+6) + j*3, r0_6, r1_6, r2_6); + v_store_interleave(dst16 + dstep_e*(i+7) + j*3, r0_7, r1_7, r2_7); + } + for (; j < n; j++) + for (int k = 0; k < 8; k++) + { + dst16[dstep_e*(i+k) + j*3 + 0] = src16[sstep_e*j + (i+k)*3 + 0]; + dst16[dstep_e*(i+k) + j*3 + 1] = src16[sstep_e*j + (i+k)*3 + 1]; + dst16[dstep_e*(i+k) + j*3 + 2] = src16[sstep_e*j + (i+k)*3 + 2]; + } + } + for (; i < m; i++) + for (int j = 0; j < n; j++) + { + dst16[dstep_e*i + j*3 + 0] = src16[sstep_e*j + i*3 + 0]; + dst16[dstep_e*i + j*3 + 1] = src16[sstep_e*j + i*3 + 1]; + dst16[dstep_e*i + j*3 + 2] = src16[sstep_e*j + i*3 + 2]; + } +} +#endif template static void transpose_( const uchar* src, size_t sstep, uchar* dst, size_t dstep, Size sz ) { - int i=0, j, m = sz.width, n = sz.height; +#if CV_SIMD128 + switch (sizeof(T)) + { + case 1: transpose_8bit_simd(src, sstep, dst, dstep, sz); return; + case 2: transpose_16bit_simd(src, sstep, dst, dstep, sz); return; + case 4: transpose_32bit_simd(src, sstep, dst, dstep, sz); return; + case 6: transpose_48bit_simd(src, sstep, dst, dstep, sz); return; + default: break; + } +#endif + + int i = 0, j, m = sz.width, n = sz.height; #if CV_ENABLE_UNROLLED for(; i <= m - 4; i += 4 ) @@ -146,10 +465,16 @@ static bool ocl_transpose( InputArray _src, OutputArray _dst ) return false; } + String deviceMacro; + if (dev.isIntel()) + deviceMacro = " -D INTEL_GPU"; + else + deviceMacro = ""; + ocl::Kernel k(kernelName.c_str(), ocl::core::transpose_oclsrc, - format("-D T=%s -D T1=%s -D cn=%d -D TILE_DIM=%d -D BLOCK_ROWS=%d -D rowsPerWI=%d%s", + format("-D T=%s -D T1=%s -D cn=%d -D TILE_DIM=%d -D BLOCK_ROWS=%d -D rowsPerWI=%d%s%s", ocl::memopTypeToStr(type), ocl::memopTypeToStr(depth), - cn, TILE_DIM, BLOCK_ROWS, rowsPerWI, inplace ? " -D INPLACE" : "")); + cn, TILE_DIM, BLOCK_ROWS, rowsPerWI, inplace ? " -D INPLACE" : "", deviceMacro.c_str())); if (k.empty()) return false; diff --git a/modules/core/src/mean.simd.hpp b/modules/core/src/mean.simd.hpp index b7ca5ba52e..db7b411f39 100644 --- a/modules/core/src/mean.simd.hpp +++ b/modules/core/src/mean.simd.hpp @@ -224,9 +224,10 @@ static int sumsqr_(const T* src0, const uchar* mask, ST* sum, SQT* sqsum, int le v0 = (ST)src[0], v1 = (ST)src[1]; s0 += v0; sq0 += (SQT)v0*v0; s1 += v1; sq1 += (SQT)v1*v1; - v0 = (ST)src[2], v1 = (ST)src[3]; - s2 += v0; sq2 += (SQT)v0*v0; - s3 += v1; sq3 += (SQT)v1*v1; + ST v2, v3; + v2 = (ST)src[2], v3 = (ST)src[3]; + s2 += v2; sq2 += (SQT)v2*v2; + s3 += v3; sq3 += (SQT)v3*v3; } sum[k] = s0; sum[k+1] = s1; sum[k+2] = s2; sum[k+3] = s3; diff --git a/modules/core/src/opencl/transpose.cl b/modules/core/src/opencl/transpose.cl index 01ea7dd9d5..ae4d7f08d3 100644 --- a/modules/core/src/opencl/transpose.cl +++ b/modules/core/src/opencl/transpose.cl @@ -53,10 +53,10 @@ #define TSIZE ((int)sizeof(T1)*3) #endif -#ifndef INPLACE - #define LDS_STEP (TILE_DIM + 1) +#ifndef INPLACE + __kernel void transpose(__global const uchar * srcptr, int src_step, int src_offset, int src_rows, int src_cols, __global uchar * dstptr, int dst_step, int dst_offset) { @@ -120,6 +120,7 @@ __kernel void transpose(__global const uchar * srcptr, int src_step, int src_off __kernel void transpose_inplace(__global uchar * srcptr, int src_step, int src_offset, int src_rows) { +#ifdef INTEL_GPU int x = get_global_id(0); int y = get_global_id(1) * rowsPerWI; @@ -141,6 +142,70 @@ __kernel void transpose_inplace(__global uchar * srcptr, int src_step, int src_o storepix(tmp, src); } } +#else + int gp_x = get_group_id(0); + int gp_y = get_group_id(1); + int lx = get_local_id(0); + int ly = get_local_id(1); + + __local T tile_a[TILE_DIM * LDS_STEP]; + __local T tile_b[TILE_DIM * LDS_STEP]; + + if (gp_x > gp_y) + { + int x_a = gp_x * TILE_DIM + lx; + int y_a = gp_y * TILE_DIM + ly; + int x_b = gp_y * TILE_DIM + lx; + int y_b = gp_x * TILE_DIM + ly; + + // Load + #pragma unroll + for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) + { + if (y_a + i < src_rows && x_a < src_rows) + tile_a[mad24(ly + i, LDS_STEP, lx)] = + loadpix(srcptr + mad24(y_a + i, src_step, mad24(x_a, TSIZE, src_offset))); + + if (y_b + i < src_rows && x_b < src_rows) + tile_b[mad24(ly + i, LDS_STEP, lx)] = + loadpix(srcptr + mad24(y_b + i, src_step, mad24(x_b, TSIZE, src_offset))); + } + + barrier(CLK_LOCAL_MEM_FENCE); + + // Store (transposed) + #pragma unroll + for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) + { + if (y_b + i < src_rows && x_b < src_rows) + storepix(tile_a[mad24(lx, LDS_STEP, ly + i)], + srcptr + mad24(y_b + i, src_step, mad24(x_b, TSIZE, src_offset))); + + if (y_a + i < src_rows && x_a < src_rows) + storepix(tile_b[mad24(lx, LDS_STEP, ly + i)], + srcptr + mad24(y_a + i, src_step, mad24(x_a, TSIZE, src_offset))); + } + } + else if (gp_x == gp_y) + { + int x = gp_x * TILE_DIM + lx; + int y = gp_y * TILE_DIM + ly; + + #pragma unroll + for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) + if (y + i < src_rows && x < src_rows) + tile_a[mad24(ly + i, LDS_STEP, lx)] = + loadpix(srcptr + mad24(y + i, src_step, mad24(x, TSIZE, src_offset))); + + barrier(CLK_LOCAL_MEM_FENCE); + + #pragma unroll + for (int i = 0; i < TILE_DIM; i += BLOCK_ROWS) + if (y + i < src_rows && x < src_rows) + storepix(tile_a[mad24(lx, LDS_STEP, ly + i)], + srcptr + mad24(y + i, src_step, mad24(x, TSIZE, src_offset))); + } +#endif } #endif // INPLACE diff --git a/modules/dnn/src/layers/batch_norm_layer.cpp b/modules/dnn/src/layers/batch_norm_layer.cpp index 1167fbf488..0c8a4659ec 100644 --- a/modules/dnn/src/layers/batch_norm_layer.cpp +++ b/modules/dnn/src/layers/batch_norm_layer.cpp @@ -83,7 +83,7 @@ public: if( hasBias ) { CV_Assert((size_t)biasBlobIndex < blobs.size()); - const Mat& b = blobs[weightsBlobIndex]; + const Mat& b = blobs[biasBlobIndex]; CV_Assert(b.isContinuous() && b.type() == CV_32F && b.total() == (size_t)n); } diff --git a/modules/dnn/src/layers/resize_layer.cpp b/modules/dnn/src/layers/resize_layer.cpp index 4cc9723d85..c42f1de26d 100644 --- a/modules/dnn/src/layers/resize_layer.cpp +++ b/modules/dnn/src/layers/resize_layer.cpp @@ -514,17 +514,30 @@ public: attrs.nearest_mode = ov::op::v4::Interpolate::NearestMode::ROUND_PREFER_FLOOR; - - std::vector shape = {outHeight, outWidth}; - auto out_shape = std::make_shared(ov::element::i64, ov::Shape{2}, shape.data()); - - auto& input_shape = ieInpNode.get_shape(); - CV_Assert_N(input_shape[2] != 0, input_shape[3] != 0); - std::vector scales = {static_cast(outHeight) / input_shape[2], static_cast(outWidth) / input_shape[3]}; - auto scales_shape = std::make_shared(ov::element::f32, ov::Shape{2}, scales.data()); + std::shared_ptr out_shape_node; + std::shared_ptr scales_node; + if (nodes.size() == 2) + { + auto& ieRefNode = nodes[1].dynamicCast()->node; + auto ref_shape = std::make_shared(ieRefNode, ov::element::i64); + auto hw_indices = std::make_shared(ov::element::i64, ov::Shape{2}, std::vector{2, 3}); + auto gather_axis = std::make_shared(ov::element::i64, ov::Shape{}, 0LL); + out_shape_node = std::make_shared(ref_shape, hw_indices, gather_axis); + std::vector dummy_scales = {1.0f, 1.0f}; + scales_node = std::make_shared(ov::element::f32, ov::Shape{2}, dummy_scales.data()); + } + else + { + std::vector shape = {outHeight, outWidth}; + out_shape_node = std::make_shared(ov::element::i64, ov::Shape{2}, shape.data()); + auto& input_shape = ieInpNode.get_shape(); + CV_Assert_N(input_shape[2] != 0, input_shape[3] != 0); + std::vector scales = {static_cast(outHeight) / input_shape[2],static_cast(outWidth) / input_shape[3]}; + scales_node = std::make_shared(ov::element::f32, ov::Shape{2}, scales.data()); + } auto axes = std::make_shared(ov::element::i64, ov::Shape{2}, std::vector{2, 3}); - auto interp = std::make_shared(ieInpNode, out_shape, scales_shape, axes, attrs); + auto interp = std::make_shared(ieInpNode, out_shape_node, scales_node, axes, attrs); return Ptr(new InfEngineNgraphNode(interp)); } #endif // HAVE_DNN_NGRAPH diff --git a/modules/highgui/misc/java/src/java/highgui+HighGui.java b/modules/highgui/misc/java/src/java/highgui+HighGui.java index 7a25095166..e986618fb4 100644 --- a/modules/highgui/misc/java/src/java/highgui+HighGui.java +++ b/modules/highgui/misc/java/src/java/highgui+HighGui.java @@ -43,8 +43,7 @@ public final class HighGui { public static void imshow(String winname, Mat img) { if (img.empty()) { - System.err.println("Error: Empty image in imshow"); - System.exit(-1); + throw new IllegalArgumentException("Image is empty"); } else { ImageWindow tmpWindow = windows.get(winname); if (tmpWindow == null) { @@ -113,8 +112,7 @@ public final class HighGui { // If there are no windows to be shown return if (windows.isEmpty()) { - System.err.println("Error: waitKey must be used after an imshow"); - System.exit(-1); + throw new IllegalStateException("No windows created. Call imshow() first"); } // Remove the unused windows @@ -145,8 +143,7 @@ public final class HighGui { win.lbl.setIcon(icon); } } else { - System.err.println("Error: no imshow associated with" + " namedWindow: \"" + win.name + "\""); - System.exit(-1); + throw new IllegalStateException("No image set for window: \"" + win.name + "\". Call imshow() first"); } } @@ -158,6 +155,7 @@ public final class HighGui { } } catch (InterruptedException e) { e.printStackTrace(); + Thread.currentThread().interrupt(); } // Set all windows as already used diff --git a/modules/imgcodecs/src/grfmt_pam.cpp b/modules/imgcodecs/src/grfmt_pam.cpp index 12371e1f5d..efae09f033 100644 --- a/modules/imgcodecs/src/grfmt_pam.cpp +++ b/modules/imgcodecs/src/grfmt_pam.cpp @@ -481,6 +481,23 @@ bool PAMDecoder::readHeader() } } while (fieldtype != PAM_HEADER_ENDHDR); + if (selected_fmt != IMWRITE_PAM_FORMAT_NULL && flds_depth) { + if (selected_fmt == IMWRITE_PAM_FORMAT_BLACKANDWHITE && m_channels != 1) { + CV_Error(Error::StsError, "fmt is IMWRITE_PAM_FORMAT_BLACKANDWHITE but number of channels is not 1"); + } + if (selected_fmt == IMWRITE_PAM_FORMAT_GRAYSCALE && m_channels != 1) { + CV_Error(Error::StsError, "fmt is IMWRITE_PAM_FORMAT_GRAYSCALE but number of channels is not 1"); + } + if (selected_fmt == IMWRITE_PAM_FORMAT_GRAYSCALE_ALPHA && m_channels != 2) { + CV_Error(Error::StsError, "fmt is IMWRITE_PAM_FORMAT_GRAYSCALE_ALPHA but number of channels is not 2"); + } + if (selected_fmt == IMWRITE_PAM_FORMAT_RGB && m_channels != 3) { + CV_Error(Error::StsError, "fmt is IMWRITE_PAM_FORMAT_RGB but number of channels is not 3"); + } + if (selected_fmt == IMWRITE_PAM_FORMAT_RGB_ALPHA && m_channels != 4) { + CV_Error(Error::StsError, "fmt is IMWRITE_PAM_FORMAT_RGB_ALPHA but number of channels is not 4"); + } + } if (flds_endhdr && flds_height && flds_width && flds_depth && flds_maxval) { if (selected_fmt == IMWRITE_PAM_FORMAT_NULL) diff --git a/modules/java/generator/android/java/org/opencv/android/Utils.java b/modules/java/generator/android/java/org/opencv/android/Utils.java index eef4c45622..80e85f5f88 100644 --- a/modules/java/generator/android/java/org/opencv/android/Utils.java +++ b/modules/java/generator/android/java/org/opencv/android/Utils.java @@ -23,24 +23,20 @@ public class Utils { public static String exportResource(Context context, int resourceId, String dirname) { String fullname = context.getResources().getString(resourceId); String resName = fullname.substring(fullname.lastIndexOf("/") + 1); - try { - InputStream is = context.getResources().openRawResource(resourceId); + try (InputStream is = context.getResources().openRawResource(resourceId)) { File resDir = context.getDir(dirname, Context.MODE_PRIVATE); File resFile = new File(resDir, resName); - FileOutputStream os = new FileOutputStream(resFile); - - byte[] buffer = new byte[4096]; - int bytesRead; - while ((bytesRead = is.read(buffer)) != -1) { - os.write(buffer, 0, bytesRead); + try (FileOutputStream os = new FileOutputStream(resFile)) { + byte[] buffer = new byte[4096]; + int bytesRead; + while ((bytesRead = is.read(buffer)) != -1) { + os.write(buffer, 0, bytesRead); + } } - is.close(); - os.close(); return resFile.getAbsolutePath(); } catch (IOException e) { - e.printStackTrace(); throw new CvException("Failed to export resource " + resName + ". Exception thrown: " + e); } @@ -53,19 +49,22 @@ public class Utils { public static Mat loadResource(Context context, int resourceId, int flags) throws IOException { - InputStream is = context.getResources().openRawResource(resourceId); - ByteArrayOutputStream os = new ByteArrayOutputStream(is.available()); + byte[] data; + int dataLength; - byte[] buffer = new byte[4096]; - int bytesRead; - while ((bytesRead = is.read(buffer)) != -1) { - os.write(buffer, 0, bytesRead); + try (InputStream is = context.getResources().openRawResource(resourceId); + ByteArrayOutputStream os = new ByteArrayOutputStream(is.available())) { + byte[] buffer = new byte[4096]; + int bytesRead; + while ((bytesRead = is.read(buffer)) != -1) { + os.write(buffer, 0, bytesRead); + } + data = os.toByteArray(); + dataLength = os.size(); } - is.close(); - Mat encoded = new Mat(1, os.size(), CvType.CV_8U); - encoded.put(0, 0, os.toByteArray()); - os.close(); + Mat encoded = new Mat(1, dataLength, CvType.CV_8U); + encoded.put(0, 0, data); Mat decoded = Imgcodecs.imdecode(encoded, flags); encoded.release(); diff --git a/modules/python/src2/cv2_convert.cpp b/modules/python/src2/cv2_convert.cpp index 0e99e0726f..c02a83a15e 100644 --- a/modules/python/src2/cv2_convert.cpp +++ b/modules/python/src2/cv2_convert.cpp @@ -209,9 +209,9 @@ bool pyopencv_to(PyObject* o, Mat& m, const ArgInfo& info) if (ismultichannel) { int channels = ndims >= 1 ? (int)_sizes[ndims - 1] : 1; - if (channels > CV_CN_MAX) + if (channels < 1 || channels > CV_CN_MAX) { - failmsg("%s unable to wrap channels, too high (%d > CV_CN_MAX=%d)", info.name, (int)channels, (int)CV_CN_MAX); + failmsg("%s unable to wrap channels, invalid count (%d, must be in [1, %d])", info.name, (int)channels, (int)CV_CN_MAX); return false; } ndims--; diff --git a/modules/python/test/test_mat.py b/modules/python/test/test_mat.py index 72614fda36..73efdf3fb8 100644 --- a/modules/python/test/test_mat.py +++ b/modules/python/test/test_mat.py @@ -89,6 +89,23 @@ try: print(res1) + def test_mat_wrap_channels_zero(self): + # Passing a 0-channel array must raise cv.error, not segfault. + data = np.zeros((100, 100, 0), dtype=np.uint8) + + with self.assertRaises(cv.error): + cv.resize(data, (200, 200)) # channels=0 -> invalid, must not segfault + + with self.assertRaises(cv.error): + mat_data = cv.Mat(data, wrap_channels=True) # unable to wrap channels, invalid count (0) + cv.utils.dumpInputArray(mat_data) + + # Verify that channels=1 (lower bound) still works correctly + data_1ch = np.zeros((100, 100, 1), dtype=np.uint8) + mat_1ch = cv.Mat(data_1ch, wrap_channels=True) + res = cv.utils.dumpInputArray(mat_1ch) + self.assertIn("CV_8UC1", res) + def test_ufuncs(self): data = np.arange(10) mat_data = cv.Mat(data) diff --git a/modules/videoio/include/opencv2/videoio.hpp b/modules/videoio/include/opencv2/videoio.hpp index 5560fedb0d..652e232e7e 100644 --- a/modules/videoio/include/opencv2/videoio.hpp +++ b/modules/videoio/include/opencv2/videoio.hpp @@ -138,7 +138,7 @@ enum VideoCaptureProperties { CAP_PROP_FOURCC =6, //!< 4-character code of codec. see VideoWriter::fourcc . CAP_PROP_FRAME_COUNT =7, //!< Number of frames in the video file. CAP_PROP_FORMAT =8, //!< Format of the %Mat objects (see Mat::type()) returned by VideoCapture::retrieve(). - //!< Set value -1 to fetch undecoded RAW video streams (as Mat 8UC1). + //!< Set value -1 to fetch undecoded RAW video streams (as Mat 8UC1). Default is 8UC3. FFmpeg backend supports 8UC4 with alpha, if it's available. CAP_PROP_MODE =9, //!< Backend-specific value indicating the current capture mode. CAP_PROP_BRIGHTNESS =10, //!< Brightness of the image (only for those cameras that support). CAP_PROP_CONTRAST =11, //!< Contrast of the image (only for cameras). @@ -226,6 +226,8 @@ enum VideoWriterProperties { VIDEOWRITER_PROP_KEY_FLAG = 11, //!< Set to non-zero to signal that the following frames are key frames or zero if not, when encapsulating raw video (\ref VIDEOWRITER_PROP_RAW_VIDEO != 0). FFmpeg back-end only. VIDEOWRITER_PROP_PTS = 12, //!< Specifies the frame presentation timestamp for each frame using the FPS time base. This property is **only** necessary when encapsulating **externally** encoded video where the decoding order differs from the presentation order, such as in GOP patterns with bi-directional B-frames. The value should be provided by your external encoder and for video sources with fixed frame rates it is equivalent to dividing the current frame's presentation time (\ref CAP_PROP_POS_MSEC) by the frame duration (1000.0 / VideoCapture::get(\ref CAP_PROP_FPS)). It can be queried from the resulting encapsulated video file using VideoCapture::get(\ref CAP_PROP_PTS). FFmpeg back-end only. VIDEOWRITER_PROP_DTS_DELAY = 13, //!< Specifies the maximum difference between presentation (pts) and decompression timestamps (dts) using the FPS time base. This property is necessary **only** when encapsulating **externally** encoded video where the decoding order differs from the presentation order, such as in GOP patterns with bi-directional B-frames. The value should be calculated based on the specific GOP pattern used during encoding. For example, in a GOP with presentation order IBP and decoding order IPB, this value would be 1, as the B-frame is the second frame presented but the third to be decoded. It can be queried from the resulting encapsulated video file using VideoCapture::get(\ref CAP_PROP_DTS_DELAY). Non-zero values usually imply the stream is encoded using B-frames. FFmpeg back-end only. + VIDEOWRITER_PROP_COLOR_SPACE = 14, //!< (**open-only**) GStreamer backend only. Pixel format for the encoding profile. Default is "I420". Other values: "NV12", "BGRx". See GStreamer raw video formats for more options. + VIDEOWRITER_PROP_ENABLE_ALPHA = 15, //!< (**open-only**) FFmpeg backend only. Defines that input frames contain alpha channel. #ifndef CV_DOXYGEN CV__VIDEOWRITER_PROP_LATEST #endif diff --git a/modules/videoio/src/cap_ffmpeg_impl.hpp b/modules/videoio/src/cap_ffmpeg_impl.hpp index afd8fd4098..33f6eecc70 100644 --- a/modules/videoio/src/cap_ffmpeg_impl.hpp +++ b/modules/videoio/src/cap_ffmpeg_impl.hpp @@ -612,6 +612,7 @@ struct CvCapture_FFMPEG bool rawModeInitialized; bool rawSeek; bool convertRGB; + bool enableAlpha; AVPacket packet_filtered; #if LIBAVFORMAT_BUILD >= CALC_FFMPEG_VERSION(58, 20, 100) AVBSFContext* bsfc; @@ -669,6 +670,7 @@ void CvCapture_FFMPEG::init() rawModeInitialized = false; rawSeek = false; convertRGB = true; + enableAlpha = false; memset(&packet_filtered, 0, sizeof(packet_filtered)); av_init_packet(&packet_filtered); bsfc = NULL; @@ -1085,9 +1087,10 @@ bool CvCapture_FFMPEG::open(const char* _filename, int index, const Ptr(CAP_PROP_FORMAT); @@ -1098,8 +1101,19 @@ bool CvCapture_FFMPEG::open(const char* _filename, int index, const Ptrcolor_primaries) << ", transfer: " << av_color_transfer_name(sw_picture->color_trc) ); - const AVPixelFormat result_format = convertRGB ? AV_PIX_FMT_BGR24 : (AVPixelFormat)sw_picture->format; + + const AVPixelFormat color_format = enableAlpha ? AV_PIX_FMT_BGRA : AV_PIX_FMT_BGR24; + const AVPixelFormat result_format = convertRGB ? color_format : (AVPixelFormat)sw_picture->format; switch (result_format) { + case AV_PIX_FMT_BGRA: *depth = CV_8U; *cn = 4; break; case AV_PIX_FMT_BGR24: *depth = CV_8U; *cn = 3; break; case AV_PIX_FMT_GRAY8: *depth = CV_8U; *cn = 1; break; case AV_PIX_FMT_GRAY16LE: *depth = CV_16U; *cn = 1; break; @@ -2122,6 +2139,12 @@ double CvCapture_FFMPEG::getProperty( int property_id ) const case CAP_PROP_FORMAT: if (rawMode) return -1; + else if (!convertRGB) + return CV_8UC1; + else if (enableAlpha) + return CV_8UC4; + else + return CV_8UC3; break; case CAP_PROP_CONVERT_RGB: return convertRGB; @@ -2365,8 +2388,20 @@ bool CvCapture_FFMPEG::setProperty( int property_id, double value ) seek((int64_t)(value*ic->duration)); return true; case CAP_PROP_FORMAT: + if (!convertRGB) + return false; if (value == -1) return setRaw(); + else if (value == CV_8UC3) + { + enableAlpha = false; + return true; + } + else if (value == CV_8UC4) + { + enableAlpha = true; + return true; + } return false; case CAP_PROP_CONVERT_RGB: convertRGB = (value != 0); @@ -2775,6 +2810,12 @@ bool CvVideoWriter_FFMPEG::writeFrame( const unsigned char* data, int step, int return false; } } + else if (input_pix_fmt == AV_PIX_FMT_BGRA) { + if (cn != 4) { + CV_LOG_WARNING(NULL, "write frame skipped - expected 4 channels but got " << cn); + return false; + } + } else if (input_pix_fmt == AV_PIX_FMT_GRAY8 || input_pix_fmt == AV_PIX_FMT_GRAY16LE) { if (cn != 1) { CV_LOG_WARNING(NULL, "write frame skipped - expected 1 channel but got " << cn); @@ -3151,6 +3192,8 @@ bool CvVideoWriter_FFMPEG::open( const char * filename, int fourcc, use_opencl = params.get(VIDEOWRITER_PROP_HW_ACCELERATION_USE_OPENCL); } + bool enable_alpha = params.get(VIDEOWRITER_PROP_ENABLE_ALPHA, false); + if (params.warnUnusedParameters()) { CV_LOG_ERROR(NULL, "VIDEOIO/FFMPEG: unsupported parameters in VideoWriter, see logger INFO channel for details"); @@ -3185,7 +3228,7 @@ bool CvVideoWriter_FFMPEG::open( const char * filename, int fourcc, { switch (depth) { - case CV_8U: input_pix_fmt = AV_PIX_FMT_BGR24; break; + case CV_8U: input_pix_fmt = enable_alpha ? AV_PIX_FMT_BGRA : AV_PIX_FMT_BGR24; break; default: CV_LOG_WARNING(NULL, "Unsupported input depth for color image: " << depth); return false; @@ -3387,7 +3430,7 @@ bool CvVideoWriter_FFMPEG::open( const char * filename, int fourcc, break; default: // good for lossy formats, MPEG, etc. - codec_pix_fmt = AV_PIX_FMT_YUV420P; + codec_pix_fmt = enable_alpha ? AV_PIX_FMT_YUVA420P : AV_PIX_FMT_YUV420P; break; } diff --git a/modules/videoio/src/cap_gstreamer.cpp b/modules/videoio/src/cap_gstreamer.cpp index 24e414c5f9..042cce3c83 100644 --- a/modules/videoio/src/cap_gstreamer.cpp +++ b/modules/videoio/src/cap_gstreamer.cpp @@ -2512,25 +2512,43 @@ bool CvVideoWriter_GStreamer::open( const std::string &filename, int fourcc, CV_WARN("OpenCV backend does not support this file type (extension): " << filename); return false; } - - //create pipeline elements encodebin.reset(gst_element_factory_make("encodebin", NULL)); - + if (!encodebin) + { + CV_WARN("GStreamer: cannot create encodebin element"); + return false; + } GSafePtr containercaps; GSafePtr containerprofile; GSafePtr videoprofile; containercaps.attach(gst_caps_from_string(mime)); - - //create encodebin profile containerprofile.attach(gst_encoding_container_profile_new("container", "container", containercaps.get(), NULL)); - videoprofile.reset(gst_encoding_video_profile_new(videocaps.get(), NULL, NULL, 1)); - gst_encoding_container_profile_add_profile(containerprofile.get(), (GstEncodingProfile*)videoprofile.get()); - + unsigned int colorspace_fourcc = (unsigned int)params.get(VIDEOWRITER_PROP_COLOR_SPACE, CV_FOURCC('I', '4', '2', '0')); + const char* colorspace = gst_video_format_to_string(gst_video_format_from_fourcc(colorspace_fourcc)); + GSafePtr prof_caps; + std::string caps_str = std::string("video/x-raw, format=") + std::string(colorspace); + prof_caps.attach(gst_caps_from_string(caps_str.c_str())); + videoprofile.attach(gst_encoding_video_profile_new(prof_caps.get(), NULL, NULL, 1)); + // Transfer ownership to the container profile + gst_encoding_container_profile_add_profile( + containerprofile.get(), + (GstEncodingProfile*)videoprofile.detach() + ); g_object_set(G_OBJECT(encodebin.get()), "profile", containerprofile.get(), NULL); source.reset(gst_element_factory_make("appsrc", NULL)); + if (!source) + { + CV_WARN("GStreamer: cannot create appsrc element"); + return false; + } file.reset(gst_element_factory_make("filesink", NULL)); + if (!file) + { + CV_WARN("GStreamer: cannot create filesink element"); + return false; + } g_object_set(G_OBJECT(file.get()), "location", (const char*)filename.c_str(), NULL); } @@ -2709,21 +2727,22 @@ void CvVideoWriter_GStreamer::write(InputArray image) Mat imageMat = image.getMat(); const size_t buf_size = imageMat.total() * imageMat.elemSize(); - duration = ((double)1/framerate) * GST_SECOND; - timestamp = num_frames * duration; + duration = gst_util_uint64_scale_int(GST_SECOND, 1, framerate); + timestamp = gst_util_uint64_scale_int(num_frames, GST_SECOND, framerate); //gst_app_src_push_buffer takes ownership of the buffer, so we need to supply it a copy GstBuffer *buffer = gst_buffer_new_allocate(NULL, buf_size, NULL); GstMapInfo info; - gst_buffer_map(buffer, &info, (GstMapFlags)GST_MAP_READ); - memcpy(info.data, (guint8*)imageMat.data, buf_size); - gst_buffer_unmap(buffer, &info); + if (gst_buffer_map(buffer, &info, (GstMapFlags)GST_MAP_WRITE)) { + memcpy(info.data, (guint8*)imageMat.data, buf_size); + gst_buffer_unmap(buffer, &info); + } GST_BUFFER_DURATION(buffer) = duration; GST_BUFFER_PTS(buffer) = timestamp; GST_BUFFER_DTS(buffer) = timestamp; //set the current number in the frame GST_BUFFER_OFFSET(buffer) = num_frames; - + GST_BUFFER_OFFSET_END(buffer) = num_frames + 1; ret = gst_app_src_push_buffer(GST_APP_SRC(source.get()), buffer); if (ret != GST_FLOW_OK) { diff --git a/modules/videoio/test/test_ffmpeg.cpp b/modules/videoio/test/test_ffmpeg.cpp index abf28acbc1..6b1c482fe4 100644 --- a/modules/videoio/test/test_ffmpeg.cpp +++ b/modules/videoio/test/test_ffmpeg.cpp @@ -1030,6 +1030,75 @@ inline static std::string videoio_ffmpeg_mismatch_name_printer(const testing::Te INSTANTIATE_TEST_CASE_P(/**/, videoio_ffmpeg_channel_mismatch, testing::ValuesIn(mismatch_cases), videoio_ffmpeg_mismatch_name_printer); +#ifndef _WIN32 + +typedef tuple AlphaChannelParams; +typedef testing::TestWithParam< AlphaChannelParams > videoio_ffmpeg_alpha_channel; + +// New feature in https://github.com/opencv/opencv/pull/28751 requires FFmpeg wrapper rebuild on Windows +TEST_P(videoio_ffmpeg_alpha_channel, write_read) +{ + if (!videoio_registry::hasBackend(CAP_FFMPEG)) + throw SkipTestException("FFmpeg backend was not found"); + + const int fourcc = fourccFromString(get<0>(GetParam())); + const string filename = "video_with_alpha_channel." + get<1>(GetParam()); + cv::VideoWriter writer(filename, cv::CAP_FFMPEG, fourcc, 1, Size(320, 240), + {VIDEOWRITER_PROP_IS_COLOR, 1, + VIDEOWRITER_PROP_ENABLE_ALPHA, 1}); + + ASSERT_TRUE(writer.isOpened()); + + for (int i = 0; i < 10; i ++) + { + cv::Mat frame; + cv::Mat gray_frame(240, 320, CV_8UC1, cv::Scalar::all(0)); + gray_frame(Rect(i*10, i*10, i*10, i*10)).setTo(255); + cv::Mat channels[4] = {gray_frame, gray_frame, gray_frame, gray_frame}; + cv::merge(channels, 4, frame); + writer.write(frame); + } + + writer.release(); + + cv::VideoCapture cap(filename, cv::CAP_FFMPEG, {cv::CAP_PROP_FORMAT, CV_8UC4}); + ASSERT_TRUE(cap.isOpened()); + ASSERT_EQ(10, cap.get(cv::CAP_PROP_FRAME_COUNT)); + + for (int i = 0; i < 10; i++) + { + cv::Mat frame; + cap >> frame; + EXPECT_EQ(4, frame.channels()); + EXPECT_EQ(320, frame.cols); + EXPECT_EQ(240, frame.rows); + EXPECT_EQ(0, frame.data[0]); + EXPECT_EQ(0, frame.data[1]); + EXPECT_EQ(0, frame.data[2]); + EXPECT_EQ(0, frame.data[3]); + + cv::Mat channels[4]; + cv::split(frame, channels); + int g_non_zero = cv::countNonZero(channels[1]); + int alpha_non_zero = cv::countNonZero(channels[3]); + + EXPECT_EQ(g_non_zero, alpha_non_zero); + } + remove(filename.c_str()); +} + +AlphaChannelParams alpha_params[] = +{ + make_tuple("FFV1", "mkv"), + make_tuple("FFV1", "avi") + // webm and hevc formats are disable as require fresh FFmpeg + //make_tuple("VP90", "webm") + //make_tuple("hevc", "mp4") +}; + +INSTANTIATE_TEST_CASE_P(/**/, videoio_ffmpeg_alpha_channel, testing::ValuesIn(alpha_params)); +#endif + // PR: https://github.com/opencv/opencv/pull/27523 // TODO: Enable the tests back on Windows after FFmpeg plugin rebuild #ifndef _WIN32