static check

CMakeLists.txt

@@ -1,7 +1,19 @@
 cmake_minimum_required(VERSION 3.12)
 
 project(match)
-find_package(OpenCV REQUIRED)
+find_package(OpenCV 4.8 REQUIRED)
+
+option(ENABLE_OPENMP "enable openmp" OFF)
+if(ENABLE_OPENMP)
+        # find OpenMP
+        find_package(OpenMP REQUIRED)
+
+        if(OPENMP_FOUND)
+            set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")
+            set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+            set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
+        endif(OPENMP_FOUND)
+endif(ENABLE_OPENMP)
 
 #==============================================================
 #library
@@ -12,16 +24,24 @@ add_library(algo SHARED
     serialize.cpp
     privateType.h
     apiExport.h
+    integral.h
+    integral.cpp
 )
 
 target_include_directories(algo PRIVATE ${OpenCV_INCLUDE_DIRS})
-target_link_libraries(algo ${OpenCV_LIBRARIES})
+target_link_libraries(algo ${OpenCV_LIBRARIES} $<$<BOOL:${OPENMP_FOUND}>:OpenMP::OpenMP_CXX>)
 target_compile_options(algo PRIVATE
         $<$<CXX_COMPILER_ID:MSVC>:/W4 /WX /external:W0>
-        $<$<STREQUAL:${CMAKE_SYSTEM_NAME},Linux>: -mlsx -fPIC -fvisibility=hidden -Wl,--exclude-libs,ALL -Wall -Wextra -Wpedantic -Wmisleading-indentation -Wunused -Wuninitialized -Wshadow -Wconversion -Werror>
+        $<$<STREQUAL:${CMAKE_SYSTEM_NAME},Linux>: -fPIC -fvisibility=hidden -Wl,--exclude-libs,ALL -Wall -Wextra -Wpedantic -Wmisleading-indentation -Wunused -Wuninitialized -Wshadow -Wconversion -Werror>
         $<$<AND:$<CXX_COMPILER_ID:Clang>,$<STREQUAL:${CMAKE_SYSTEM_NAME},Windows>>:/W4 /WX /external:W0>
+        #$<$<AND:$<CXX_COMPILER_ID:Clang>,$<STREQUAL:${CMAKE_SYSTEM_NAME},Darwin>>: -fPIC -fvisibility=hidden -Wall -Wextra -Wpedantic -Wmisleading-indentation -Wunused -Wuninitialized -Wshadow -Wconversion -Werror>
+        $<$<STREQUAL:${CMAKE_SYSTEM_PROCESSOR},loongarch64>:-mlsx>
 )
-target_compile_definitions(algo PUBLIC API_EXPORTS)
+target_compile_definitions(algo PUBLIC API_EXPORTS
+        $<$<STREQUAL:${CMAKE_SYSTEM_PROCESSOR},loongarch64>:CV_LSX>
+)
+
+message("Arch:${CMAKE_SYSTEM_PROCESSOR}")
 
 #==============================================================
 #exe
@@ -37,3 +57,4 @@ target_compile_options(${PROJECT_NAME} PRIVATE
         $<$<STREQUAL:${CMAKE_SYSTEM_NAME},Linux>: -fPIC -fvisibility=hidden -Wall -Wextra -Wpedantic -Wmisleading-indentation -Wunused -Wuninitialized -Wshadow -Wconversion -Werror>
         $<$<AND:$<CXX_COMPILER_ID:Clang>,$<STREQUAL:${CMAKE_SYSTEM_NAME},Windows>>:/W4 /WX /external:W0>
 )
+target_compile_definitions(${PROJECT_NAME} PRIVATE IMG_DIR="${CMAKE_CURRENT_SOURCE_DIR}/img")

LICENSE

@@ -0,0 +1,24 @@
+BSD 2-Clause License
+
+Copyright (c) 2024, SurfaceMan
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

README.md

@@ -0,0 +1,20 @@
+# Template match with gray model(ncc)
+
+## Note: branch [feature-rotate-model](https://github.com/SurfaceMan/gray_match/tree/feature-rotate-model) method 2x faster at matching!
+
+## highlights:
+ 1. original code based [Fastest_Image_Pattern_Matching](https://github.com/DennisLiu1993/Fastest_Image_Pattern_Matching), you can check out tag [v1.0](https://github.com/SurfaceMan/gray_match/releases/tag/v1.0) for more details.
+ 2. refactor simd match process with opencv [Universal intrinsics](https://docs.opencv.org/4.x/df/d91/group__core__hal__intrin.html), have be tested on x86_64(sse),arm(neon),LoongArch(lsx).
+ 3. support model save/load as binary file
+ 4. provide pure c interface
+ 5. support openmp
+ 6. position with 3x3 subpixel interpolation
+ 
+## usage:
+all you need can be found in [main.cpp](main.cpp)
+
+## gallery:
+![sample](img/result.png)
+
+
+ 

TODO.md

@@ -0,0 +1 @@
+1. opencv intrinsics api [changed](https://github.com/opencv/opencv/pull/24371) since 4.9, need refactor 

grayMatch.cpp

@@ -141,6 +141,15 @@ int computeLayers(const int width, const int height, const int minArea) {
     return layer;
 }
 
+inline cv::Point2d transform(const cv::Point2d &point, const cv::Mat &rotate) {
+    const auto ptr = rotate.ptr<double>();
+
+    auto x = point.x * ptr[ 0 ] + point.y * ptr[ 1 ] + ptr[ 2 ];
+    auto y = point.x * ptr[ 3 ] + point.y * ptr[ 4 ] + ptr[ 5 ];
+
+    return {x, y};
+}
+
 cv::Size computeRotationSize(const cv::Size &dstSize, const cv::Size &templateSize, double angle,
                              const cv::Mat &rotate) {
     if (angle > 360) {
@@ -157,25 +166,20 @@ cv::Size computeRotationSize(const cv::Size &dstSize, const cv::Size &templateSi
         return dstSize;
     }
 
-    const std::vector<cv::Point2d> points{
-        {0, 0},
-        {static_cast<double>(dstSize.width) - 1, 0},
-        {0, static_cast<double>(dstSize.height) - 1},
-        {static_cast<double>(dstSize.width) - 1, static_cast<double>(dstSize.height) - 1}};
-    std::vector<cv::Point2d> trans;
-    cv::transform(points, trans, rotate);
+    const std::vector<cv::Point2d> pt{
+        transform({0, 0}, rotate), transform({static_cast<double>(dstSize.width) - 1, 0}, rotate),
+        transform({0, static_cast<double>(dstSize.height) - 1}, rotate),
+        transform({static_cast<double>(dstSize.width) - 1, static_cast<double>(dstSize.height) - 1},
+                  rotate)};
 
-    cv::Point2d min = trans[ 0 ];
-    cv::Point2d max = trans[ 0 ];
-    for (const auto &point : trans) {
-        min.x = std::min(min.x, point.x);
-        min.y = std::min(min.y, point.y);
-        max.x = std::max(max.x, point.x);
-        max.y = std::max(max.y, point.y);
-    }
+    cv::Point2d min;
+    cv::Point2d max;
+    min.x = std::min(std::min(std::min(pt[ 0 ].x, pt[ 1 ].x), pt[ 2 ].x), pt[ 3 ].x);
+    min.y = std::min(std::min(std::min(pt[ 0 ].y, pt[ 1 ].y), pt[ 2 ].y), pt[ 3 ].y);
+    max.x = std::max(std::max(std::max(pt[ 0 ].x, pt[ 1 ].x), pt[ 2 ].x), pt[ 3 ].x);
+    max.y = std::max(std::max(std::max(pt[ 0 ].y, pt[ 1 ].y), pt[ 2 ].y), pt[ 3 ].y);
 
     if (angle > 0 && angle < 90) {
-        ;
     } else if (angle > 90 && angle < 180) {
         angle -= 90;
     } else if (angle > 180 && angle < 270) {
@@ -199,15 +203,74 @@ cv::Size computeRotationSize(const cv::Size &dstSize, const cv::Size &templateSi
                            (templateSize.width > size.width && templateSize.height < size.height) ||
                            templateSize.area() > size.area();
     if (wrongSize) {
-        size = {static_cast<int>(max.x - min.x + 0.5), static_cast<int>(max.y - min.y + 0.5)};
+        size = {static_cast<int>(lround(max.x - min.x + 0.5)),
+                static_cast<int>(lround(max.y - min.y + 0.5))};
     }
 
     return size;
 }
 
-void coeffDenominator(const cv::Mat &src, const cv::Size &templateSize, cv::Mat &result,
-                      const double mean, const double normal, const double invArea,
-                      const bool equal1) {
+#ifdef CV_SIMD
+#include "integral.h"
+
+void matchTemplateSimd(const cv::Mat &src, const cv::Mat &templateImg, cv::Mat &result) {
+    result.create(src.size() - templateImg.size() + cv::Size(1, 1), CV_32FC1);
+
+    auto       *resultStart = result.ptr<float>();
+    const auto  resultStep  = result.step[ 0 ] / result.step[ 1 ];
+    const auto *srcStart    = src.ptr<uchar>();
+    const auto  srcStep     = src.step[ 0 ];
+    const auto *temStart    = templateImg.ptr<uchar>();
+    const auto  temStep     = templateImg.step[ 0 ];
+
+    std::vector<cv::v_uint32> tmp(result.cols, cv::v_setzero_u32());
+    for (int y = 0; y < result.rows; y++) {
+        auto *resultPtr = resultStart + resultStep * y;
+
+        for (int templateRow = 0; templateRow < templateImg.rows; templateRow++) {
+            auto *temPtr = temStart + temStep * templateRow;
+
+            for (int i = 0; i < templateImg.cols; i += simdSize(cv::v_uint8)) {
+                auto  vTem   = cv::v_load_aligned(temPtr + i);
+                auto *srcPtr = srcStart + srcStep * (y + templateRow) + i;
+
+                int x = 0;
+                for (int n = 0; n < result.cols / 4; n++) {
+                    auto vSrc = cv::v_load(srcPtr + x);
+                    tmp[ x ]  = cv::v_add(tmp[ x ], cv::v_dotprod_expand_fast(vSrc, vTem));
+                    x++;
+
+                    vSrc     = cv::v_load(srcPtr + x);
+                    tmp[ x ] = cv::v_add(tmp[ x ], cv::v_dotprod_expand_fast(vSrc, vTem));
+                    x++;
+
+                    vSrc     = cv::v_load(srcPtr + x);
+                    tmp[ x ] = cv::v_add(tmp[ x ], cv::v_dotprod_expand_fast(vSrc, vTem));
+                    x++;
+
+                    vSrc     = cv::v_load(srcPtr + x);
+                    tmp[ x ] = cv::v_add(tmp[ x ], cv::v_dotprod_expand_fast(vSrc, vTem));
+                    x++;
+                }
+                for (; x < result.cols; x++) {
+                    auto vSrc = cv::v_load(srcPtr + x);
+                    tmp[ x ]  = cv::v_add(tmp[ x ], cv::v_dotprod_expand_fast(vSrc, vTem));
+                }
+            }
+        }
+
+        for (int x = 0; x < result.cols; x++) {
+            const auto sum = cv::v_reduce_sum(tmp[ x ]);
+            tmp[ x ]       = cv::v_setzero_u32();
+            resultPtr[ x ] = static_cast<float>(sum);
+        }
+    }
+}
+#endif
+
+void ccoeffDenominator(const cv::Mat &src, const cv::Size &templateSize, cv::Mat &result,
+                       const double mean, const double normal, const double invArea,
+                       const bool equal1) {
     if (equal1) {
         result = cv::Scalar::all(1);
         return;
@@ -215,43 +278,44 @@ void coeffDenominator(const cv::Mat &src, const cv::Size &templateSize, cv::Mat
 
     cv::Mat sum;
     cv::Mat sqSum;
+#ifdef CV_SIMD
+    integralSimd(src, sum, sqSum);
+#else
     cv::integral(src, sum, sqSum, CV_64F);
+#endif
 
-    const auto *q0 = sqSum.ptr<double>(0);
+    const auto *q0 = reinterpret_cast<double *>(sqSum.data);
     const auto *q1 = q0 + templateSize.width;
-    const auto *q2 = sqSum.ptr<double>(templateSize.height);
+    const auto *q2 = reinterpret_cast<double *>(sqSum.data) + sqSum.step1() * templateSize.height;
     const auto *q3 = q2 + templateSize.width;
 
-    const auto *p0 = sum.ptr<double>(0);
+    const auto *p0 = reinterpret_cast<double *>(sum.data);
     const auto *p1 = p0 + templateSize.width;
-    const auto *p2 = sum.ptr<double>(templateSize.height);
+    const auto *p2 = reinterpret_cast<double *>(sum.data) + sum.step1() * templateSize.height;
     const auto *p3 = p2 + templateSize.width;
 
-    const auto step = sum.step / sizeof(double);
+    const auto step           = sum.step / sizeof(double);
+    auto      *resultStartPtr = reinterpret_cast<float *>(result.data);
 
     for (int y = 0; y < result.rows; y++) {
-        auto *scorePtr = result.ptr<float>(y);
+        auto *scorePtr = resultStartPtr + y * result.step1();
         auto  idx      = y * step;
         for (int x = 0; x < result.cols; x++, idx++) {
             auto      &score     = scorePtr[ x ];
             const auto partSum   = p0[ idx ] - p1[ idx ] - p2[ idx ] + p3[ idx ];
-            auto       partMean  = partSum * partSum;
-            const auto num       = score - partSum * mean;
-            partMean            *= invArea;
+            const auto numerator = score - partSum * mean;
 
-            auto partSum2 = q0[ idx ] - q1[ idx ] - q2[ idx ] + q3[ idx ];
+            const auto partSqSum    = q0[ idx ] - q1[ idx ] - q2[ idx ] + q3[ idx ];
+            auto       partSqNormal = partSqSum - partSum * partSum * invArea;
 
-            const auto diff = std::max(partSum2 - partMean, 0.);
-            if (diff <= std::min(0.5, 10 * FLT_EPSILON * partSum2)) {
-                partSum2 = 0;
-            } else {
-                partSum2 = sqrt(diff) * normal;
-            }
+            const auto   diff = std::max(partSqNormal, 0.);
+            const double denominator =
+                diff <= std::min(0.5, 10 * FLT_EPSILON * partSqSum) ? 0 : sqrt(diff) * normal;
 
-            if (abs(num) < partSum2) {
-                score = static_cast<float>(num / partSum2);
-            } else if (abs(num) < partSum2 * 1.125) {
-                score = num > 0.f ? 1.f : -1.f;
+            if (abs(numerator) < denominator) {
+                score = static_cast<float>(numerator / denominator);
+            } else if (abs(numerator) < denominator * 1.125) {
+                score = numerator > 0.f ? 1.f : -1.f;
             } else {
                 score = 0;
             }
@@ -259,48 +323,14 @@ void coeffDenominator(const cv::Mat &src, const cv::Size &templateSize, cv::Mat
     }
 }
 
-#ifdef CV_SIMD
-float convSimd(const uchar *kernel, const uchar *src, const int kernelWidth) {
-    const auto blockSize = cv::v_uint8::nlanes;
-    auto       vSum      = cv::vx_setall_u32(0);
-    int        i         = 0;
-    for (; i < kernelWidth - blockSize; i += blockSize) {
-        vSum += cv::v_dotprod_expand(cv::v_load(kernel + i), cv::v_load(src + i));
-    }
-    auto sum = cv::v_reduce_sum(vSum);
-
-    for (; i < kernelWidth; i++) {
-        sum += kernel[ i ] * src[ i ];
-    }
-
-    return static_cast<float>(sum);
-}
-
-void matchTemplateSimd(const cv::Mat &src, const cv::Mat &templateImg, cv::Mat &result) {
-    result = cv::Mat::zeros(src.size() - templateImg.size() + cv::Size(1, 1), CV_32FC1);
-
-    for (int y = 0; y < result.rows; y++) {
-        auto *resultPtr = result.ptr<float>(y);
-        for (int x = 0; x < result.cols; x++) {
-            auto &score = resultPtr[ x ];
-            for (int templateRow = 0; templateRow < templateImg.rows; templateRow++) {
-                auto *srcPtr  = src.ptr<uchar>(y + templateRow) + x;
-                auto *temPtr  = templateImg.ptr<uchar>(templateRow);
-                score        += convSimd(temPtr, srcPtr, templateImg.cols);
-            }
-        }
-    }
-}
-#endif
-
 void matchTemplate(cv::Mat &src, cv::Mat &result, const Model *model, int level) {
 #ifdef CV_SIMD
     matchTemplateSimd(src, model->pyramids[ level ], result);
 #else
     cv::matchTemplate(src, model->pyramids[ level ], result, cv::TM_CCORR);
 #endif
-    coeffDenominator(src, model->pyramids[ level ].size(), result, model->mean[ level ][ 0 ],
-                     model->normal[ level ], model->invArea[ level ], model->equal1[ level ] == 1);
+    ccoeffDenominator(src, model->pyramids[ level ].size(), result, model->mean[ level ][ 0 ],
+                      model->normal[ level ], model->invArea[ level ], model->equal1[ level ] == 1);
 }
 
 void nextMaxLoc(const cv::Point &pos, const cv::Size templateSize, const double maxOverlap,
@@ -334,17 +364,26 @@ void nextMaxLoc(cv::Mat &score, const cv::Point &pos, const cv::Size templateSiz
     cv::minMaxLoc(score, nullptr, &maxScore, nullptr, &maxPos);
 }
 
-inline cv::Point2d transform(const cv::Point2d &point, const cv::Mat &rotate) {
+inline cv::Mat getRotationMatrix2D(const cv::Point2f &center, double angle) {
+    angle              *= CV_PI / 180;
+    const double alpha  = std::cos(angle);
+    const double beta   = std::sin(angle);
+
+    cv::Mat    rotate(2, 3, CV_64FC1);
     const auto ptr = rotate.ptr<double>();
+    ptr[ 0 ]       = alpha;
+    ptr[ 1 ]       = beta;
+    ptr[ 2 ]       = (1 - alpha) * center.x - beta * center.y;
+    ptr[ 3 ]       = -beta;
+    ptr[ 4 ]       = alpha;
+    ptr[ 5 ]       = beta * center.x + (1 - alpha) * center.y;
 
-    auto x = point.x * ptr[ 0 ] + point.y * ptr[ 1 ] + ptr[ 2 ];
-    auto y = point.x * ptr[ 3 ] + point.y * ptr[ 4 ] + ptr[ 5 ];
-
-    return {x, y};
+    return rotate;
 }
 
-inline cv::Point2d transform(const cv::Point2d &point, const cv::Point &center, double angle) {
-    const auto rotate = cv::getRotationMatrix2D(center, angle, 1.);
+inline cv::Point2d transform(const cv::Point2d &point, const cv::Point &center,
+                             const double angle) {
+    const auto rotate = getRotationMatrix2D(center, angle);
 
     return transform(point, rotate);
 }
@@ -354,14 +393,13 @@ inline cv::Point2d sizeCenter(const cv::Size &size) {
 }
 
 void cropRotatedRoi(const cv::Mat &src, const cv::Size &templateSize, const cv::Point2d topLeft,
-                    const cv::Mat &rotate, cv::Mat &roi) {
+                    cv::Mat &rotate, cv::Mat &roi) {
     const auto     point = transform(topLeft, rotate);
     const cv::Size paddingSize(templateSize.width + 6, templateSize.height + 6);
-    auto           rt    = rotate;
-    rt.at<double>(0, 2) -= point.x - 3;
-    rt.at<double>(1, 2) -= point.y - 3;
+    rotate.at<double>(0, 2) -= point.x - 3;
+    rotate.at<double>(1, 2) -= point.y - 3;
 
-    cv::warpAffine(src, roi, rt, paddingSize);
+    cv::warpAffine(src, roi, rotate, paddingSize);
 }
 
 void filterOverlap(std::vector<Candidate> &candidates, const std::vector<cv::RotatedRect> &rects,
@@ -431,11 +469,24 @@ Model *trainModel(const cv::Mat &src, int level) {
 
     auto  *result = new Model;
     Model &model  = *result;
-    cv::buildPyramid(src, model.pyramids, level);
-    model.borderColor = cv::mean(src).val[ 0 ] < 128 ? 255 : 0;
-    model.reserve(model.pyramids.size());
 
-    for (const auto &pyramid : model.pyramids) {
+    std::vector<cv::Mat> pyramids;
+    cv::buildPyramid(src, pyramids, level);
+    model.borderColor = cv::mean(src).val[ 0 ] < 128 ? 255 : 0;
+    model.reserve(pyramids.size());
+
+    for (const auto &pyramid : pyramids) {
+        int     alignedWidth = static_cast<int>(cv::alignSize(pyramid.cols, simdSize(cv::v_uint8)));
+        auto    img          = cv::Mat::zeros(pyramid.rows, alignedWidth, CV_8UC1);
+        cv::Mat sub          = img(cv::Rect(0, 0, pyramid.cols, pyramid.rows));
+        for (int y = 0; y < pyramid.rows; y++) {
+            auto *dstPtr = sub.ptr<uchar>(y);
+            auto *srcPtr = pyramid.ptr<uchar>(y);
+            memcpy(dstPtr, srcPtr, pyramid.cols);
+        }
+
+        model.pyramids.push_back(sub);
+
         auto invArea = 1. / pyramid.size().area();
 
         cv::Scalar mean;
@@ -460,6 +511,9 @@ inline double sizeAngleStep(const cv::Size &size) {
     return atan(2. / std::max(size.width, size.height)) * 180. / CV_PI;
 }
 
+#pragma omp declare reduction(combine : std::vector<Candidate> : omp_out.insert(                   \
+        omp_out.end(), omp_in.begin(), omp_in.end()))
+
 std::vector<Candidate> matchTopLevel(const cv::Mat &dstTop, double startAngle, double spanAngle,
                                      double maxOverlap, double minScore, int maxCount,
                                      const Model *model, int level) {
@@ -469,12 +523,13 @@ std::vector<Candidate> matchTopLevel(const cv::Mat &dstTop, double startAngle, d
     auto        angleStep         = sizeAngleStep(templateTop.size());
     cv::Point2d center            = sizeCenter(dstTop.size());
     const auto  topScoreThreshold = minScore * pow(0.9, level);
-    bool calMaxByBlock = (dstTop.size().area() / templateTop.size().area() > 500) && maxCount > 10;
+    bool calMaxByBlock = dstTop.size().area() / templateTop.size().area() > 500 && maxCount > 10;
 
     const auto count = static_cast<int>(spanAngle / angleStep) + 1;
+#pragma omp parallel for reduction(combine : candidates)
     for (int i = 0; i < count; i++) {
         const auto angle  = startAngle + angleStep * i;
-        auto       rotate = cv::getRotationMatrix2D(center, angle, 1.);
+        auto       rotate = getRotationMatrix2D(center, angle);
         auto       size   = computeRotationSize(dstTop.size(), templateTop.size(), angle, rotate);
 
         auto tx                  = (size.width - 1) / 2. - center.x;
@@ -532,13 +587,57 @@ std::vector<Candidate> matchTopLevel(const cv::Mat &dstTop, double startAngle, d
     return candidates;
 }
 
+cv::Point2f computeSubpixel(const cv::Mat &score) {
+    cv::Point2f result(0, 0);
+    auto       *mag = score.ptr<float>();
+
+    const auto gx  = (-mag[ 0 ] + mag[ 2 ] - mag[ 3 ] + mag[ 5 ] - mag[ 6 ] + mag[ 8 ]) / 3.0f;
+    const auto gy  = (mag[ 6 ] + mag[ 7 ] + mag[ 8 ] - mag[ 0 ] - mag[ 1 ] - mag[ 2 ]) / 3.0f;
+    const auto gxx = (mag[ 0 ] - 2.0f * mag[ 1 ] + mag[ 2 ] + mag[ 3 ] - 2.0f * mag[ 4 ] +
+                      mag[ 5 ] + mag[ 6 ] - 2.0f * mag[ 7 ] + mag[ 8 ]) /
+                     6.0f;
+    const auto gxy = (-mag[ 0 ] + mag[ 2 ] + mag[ 6 ] - mag[ 8 ]) / 2.0f;
+    const auto gyy = (mag[ 0 ] + mag[ 1 ] + mag[ 2 ] - 2.0f * (mag[ 3 ] + mag[ 4 ] + mag[ 5 ]) +
+                      mag[ 6 ] + mag[ 7 ] + mag[ 8 ]) /
+                     6.0f;
+
+    cv::Mat hessian(2, 2, CV_32FC1);
+    hessian.at<float>(0, 0) = gxx;
+    hessian.at<float>(0, 1) = gxy;
+    hessian.at<float>(1, 0) = gxy;
+    hessian.at<float>(1, 1) = gyy;
+
+    cv::Mat eigenvalue;
+    cv::Mat eigenvector;
+    cv::eigen(hessian, eigenvalue, eigenvector);
+    float nx, ny;
+    if (fabs(eigenvalue.at<float>(0, 0)) >= fabs(eigenvalue.at<float>(1, 0))) {
+        nx = eigenvector.at<float>(0, 0);
+        ny = eigenvector.at<float>(0, 1);
+    } else {
+        nx = eigenvector.at<float>(1, 0);
+        ny = eigenvector.at<float>(1, 1);
+    }
+
+    const auto denominator = gxx * nx * nx + 2 * gxy * nx * ny + gyy * ny * ny;
+    if (denominator != 0.0) {
+        const auto t = -(gx * nx + gy * ny) / denominator;
+        result.x     = t * nx;
+        result.y     = t * ny;
+    }
+
+    return result;
+}
+
 std::vector<Candidate> matchDownLevel(const std::vector<cv::Mat>   &pyramids,
                                       const std::vector<Candidate> &candidates, double minScore,
                                       int subpixel, const Model *model, int level) {
     std::vector<Candidate> levelMatched;
+    auto                   count = static_cast<int>(candidates.size());
 
-    for (const auto &candidate : candidates) {
-        auto pose    = candidate;
+#pragma omp parallel for reduction(combine : levelMatched) schedule(dynamic)
+    for (int index = 0; index < count; index++) {
+        auto pose    = candidates[ index ];
         bool matched = true;
         for (int currentLevel = level - 1; currentLevel >= 0; currentLevel--) {
             const auto &currentTemplateImg = model->pyramids[ currentLevel ];
@@ -560,7 +659,7 @@ std::vector<Candidate> matchDownLevel(const std::vector<cv::Mat>   &pyramids,
             cv::Mat   newScoreRect;
             for (int i = -1; i <= 1; i++) {
                 auto angle  = pose.angle + i * currentAngleStep;
-                auto rotate = cv::getRotationMatrix2D(center, angle, 1.);
+                auto rotate = getRotationMatrix2D(center, angle);
 
                 cv::Mat roi;
                 cropRotatedRoi(currentDstImg, tmpSize, topLeft, rotate, roi);
@@ -592,7 +691,9 @@ std::vector<Candidate> matchDownLevel(const std::vector<cv::Mat>   &pyramids,
             }
 
             if (!newScoreRect.empty()) {
-                // TODO subpixel
+                auto offset         = computeSubpixel(newScoreRect);
+                newCandidate.pos.x += offset.x;
+                newCandidate.pos.y += offset.y;
             }
 
             // back to
@@ -663,7 +764,7 @@ std::vector<Pose> matchModel(const cv::Mat &dst, const Model *model, int level,
         for (const auto &candidate : levelMatched) {
             std::vector<cv::Point2f> points{topRight + cv::Point2f(candidate.pos),
                                             bottomRight + cv::Point2f(candidate.pos)};
-            auto rotate = cv::getRotationMatrix2D(candidate.pos, -candidate.angle, 1.);
+            auto                     rotate = getRotationMatrix2D(candidate.pos, -candidate.angle);
             std::vector<cv::Point2f> rotatedPoints;
             cv::transform(points, rotatedPoints, rotate);
 
@@ -695,14 +796,15 @@ std::vector<Pose> matchModel(const cv::Mat &dst, const Model *model, int level,
     return result;
 }
 
-Model_t trainModel(const unsigned char *data, int width, int height, int channels, int bytesPerline,
+Model_t trainModel(const unsigned char *data, int width, int height, int channels, int bytesPerLine,
                    int roiLeft, int roiTop, int roiWidth, int roiHeight, int levelNum) {
-    if ((1 != channels && 3 == channels && 4 == channels) || nullptr == data) {
+    if ((1 != channels && 3 != channels && 4 != channels) || nullptr == data) {
         return nullptr;
     }
 
-    auto    type = channels == 1 ? CV_8UC1 : (channels == 3 ? CV_8UC3 : CV_8UC4);
-    cv::Mat img(cv::Size(width, height), type, (void *)data, bytesPerline);
+    const auto    type = channels == 1 ? CV_8UC1 : channels == 3 ? CV_8UC3 : CV_8UC4;
+    const cv::Mat img(cv::Size(width, height), type, const_cast<unsigned char *>(data),
+                      bytesPerLine);
 
     cv::Mat src;
     if (1 == channels) {
@@ -711,9 +813,9 @@ Model_t trainModel(const unsigned char *data, int width, int height, int channel
         cv::cvtColor(img, src, channels == 3 ? cv::COLOR_RGB2GRAY : cv::COLOR_RGBA2GRAY);
     }
 
-    cv::Rect rect(roiLeft, roiTop, roiWidth, roiHeight);
-    cv::Rect imageRect(0, 0, width, height);
-    auto     roi = rect & imageRect;
+    const cv::Rect rect(roiLeft, roiTop, roiWidth, roiHeight);
+    const cv::Rect imageRect(0, 0, width, height);
+    const auto     roi = rect & imageRect;
     if (roi.empty()) {
         return nullptr;
     }
@@ -721,8 +823,8 @@ Model_t trainModel(const unsigned char *data, int width, int height, int channel
     return trainModel(src(roi), levelNum);
 }
 
-void matchModel(const unsigned char *data, int width, int height, int channels, int bytesPerline,
-                int roiLeft, int roiTop, int roiWidth, int roiHeight, const Model_t model,
+void matchModel(const unsigned char *data, int width, int height, int channels, int bytesPerLine,
+                int roiLeft, int roiTop, int roiWidth, int roiHeight, Model *const model,
                 int *count, Pose *poses, int level, double startAngle, double spanAngle,
                 double maxOverlap, double minScore, int subpixel) {
     if (nullptr == count) {
@@ -739,8 +841,9 @@ void matchModel(const unsigned char *data, int width, int height, int channels,
         return;
     }
 
-    auto    type = channels == 1 ? CV_8UC1 : (channels == 3 ? CV_8UC3 : CV_8UC4);
-    cv::Mat img(cv::Size(width, height), type, (void *)data, bytesPerline);
+    const auto    type = channels == 1 ? CV_8UC1 : channels == 3 ? CV_8UC3 : CV_8UC4;
+    const cv::Mat img(cv::Size(width, height), type, const_cast<unsigned char *>(data),
+                      bytesPerLine);
 
     cv::Mat dst;
     if (1 == channels) {
@@ -749,18 +852,18 @@ void matchModel(const unsigned char *data, int width, int height, int channels,
         cv::cvtColor(img, dst, channels == 3 ? cv::COLOR_RGB2GRAY : cv::COLOR_RGBA2GRAY);
     }
 
-    cv::Rect rect(roiLeft, roiTop, roiWidth, roiHeight);
-    cv::Rect imageRect(0, 0, width, height);
-    auto     roi = rect & imageRect;
+    const cv::Rect rect(roiLeft, roiTop, roiWidth, roiHeight);
+    const cv::Rect imageRect(0, 0, width, height);
+    const auto     roi = rect & imageRect;
     if (roi.empty()) {
         *count = 0;
         return;
     }
 
-    auto result = matchModel(dst(roi), model, level, startAngle, spanAngle, maxOverlap, minScore,
-                             *count, subpixel);
+    const auto result = matchModel(dst(roi), model, level, startAngle, spanAngle, maxOverlap,
+                                   minScore, *count, subpixel);
 
-    auto size = std::min(*count, static_cast<int>(result.size()));
+    const auto size = std::min(*count, static_cast<int>(result.size()));
     for (int i = 0; i < size; i++) {
         const auto &pose = result[ i ];
         poses[ i ]       = {pose.x + static_cast<float>(roi.x), pose.y + static_cast<float>(roi.y),
@@ -779,7 +882,7 @@ void freeModel(Model_t *model) {
     *model = nullptr;
 }
 
-int modelLevel(const Model_t model) {
+int modelLevel(Model *const model) {
     if (nullptr == model) {
         return 0;
     }
@@ -787,7 +890,7 @@ int modelLevel(const Model_t model) {
     return static_cast<int>(model->pyramids.size());
 }
 
-void modelImage(const Model_t model, int level, unsigned char *data, int length, int *width,
+void modelImage(Model *const model, int level, unsigned char *data, int length, int *width,
                 int *height, int *channels) {
     if (nullptr == model) {
         return;

grayMatch.h

@@ -20,16 +20,16 @@ struct Pose {
  * @param width image width
  * @param height image height
  * @param channels image channels 1(gray)/3(rgb)/4(rgba)
- * @param bytesPerline bytes per line
+ * @param bytesPerLine bytes per line
  * @param roiLeft rectangle roi left
  * @param roiTop rectangle roi top
  * @param roiWidth rectangle roi width
  * @param roiHeight rectangle roi height
- * @param levelNum pyramid levels (> 0)
+ * @param levelNum pyramid levels (> 0:user setting,-1:auto)
  * @return
  */
 API_PUBLIC Model_t trainModel(const unsigned char *data, int width, int height, int channels,
-                              int bytesPerline, int roiLeft, int roiTop, int roiWidth,
+                              int bytesPerLine, int roiLeft, int roiTop, int roiWidth,
                               int roiHeight, int levelNum);
 /**
  * @brief match model
@@ -37,7 +37,7 @@ API_PUBLIC Model_t trainModel(const unsigned char *data, int width, int height,
  * @param width image width
  * @param height image height
  * @param channels image channels 1(gray)/3(rgb)/4(rgba)
- * @param bytesPerline bytes per line
+ * @param bytesPerLine bytes per line
  * @param roiLeft rectangle roi left
  * @param roiTop rectangle roi top
  * @param roiWidth rectangle roi width
@@ -45,7 +45,7 @@ API_PUBLIC Model_t trainModel(const unsigned char *data, int width, int height,
  * @param model trained model
  * @param count in(max detect count)/out(found count)
  * @param poses pose array inited with size not less than count
- * @param level match start at which level (level>=0 && level<modelLevel-1)
+ * @param level match start at which level (level>=0 && level<modelLevel-1,-1:auto)
  * @param startAngle rotation start angle
  * @param spanAngle rotation angle range
  * @param maxOverlap overlap threshold
@@ -54,17 +54,16 @@ API_PUBLIC Model_t trainModel(const unsigned char *data, int width, int height,
  * @return
  */
 API_PUBLIC void matchModel(const unsigned char *data, int width, int height, int channels,
-                           int bytesPerline, int roiLeft, int roiTop, int roiWidth, int roiHeight,
-                           const Model_t model, int *count, Pose *poses, int level,
-                           double startAngle, double spanAngle, double maxOverlap, double minScore,
-                           int subpixel);
+                           int bytesPerLine, int roiLeft, int roiTop, int roiWidth, int roiHeight,
+                           Model_t model, int *count, Pose *poses, int level, double startAngle,
+                           double spanAngle, double maxOverlap, double minScore, int subpixel);
 
 /**
  * @brief get trained model levels
  * @param model
  * @return pyramid level
  */
-API_PUBLIC int modelLevel(const Model_t model);
+API_PUBLIC int modelLevel(Model_t model);
 
 /**
  * @brief get trained model image
@@ -77,8 +76,8 @@ API_PUBLIC int modelLevel(const Model_t model);
  * @param channels image channels, can input nullptr
  * @return
  */
-API_PUBLIC void modelImage(const Model_t model, int level, unsigned char *data, int length,
-                           int *width, int *height, int *channels);
+API_PUBLIC void modelImage(Model_t model, int level, unsigned char *data, int length, int *width,
+                           int *height, int *channels);
 
 /**
  * @brief free model
@@ -94,10 +93,10 @@ API_PUBLIC void freeModel(Model_t *model);
  * @param size in(buffer size)/out(written size)
  * @return true(success)false(failed)
  */
-API_PUBLIC bool serialize(const Model_t model, unsigned char *buffer, int *size);
+API_PUBLIC bool serialize(Model_t model, unsigned char *buffer, int *size);
 
 /**
- * @brief desrialize model
+ * @brief deserialize model
  * @param buffer
  * @param size buffer size
  * @return model

integral.cpp

@@ -0,0 +1,180 @@
+#include "integral.h"
+#include "privateType.h"
+
+#include <opencv2/core/hal/intrin.hpp>
+
+inline void expand(const cv::v_int32 &src, cv::v_float64 &low, cv::v_float64 &high) {
+    low  = cv::v_cvt_f64(src);
+    high = cv::v_cvt_f64_high(src);
+}
+
+inline void integralSum(const cv::v_uint16 &src, double *dst, const double *prevDst,
+                        cv::v_uint32 &pre) {
+    auto sum = cv::v_add(src, cv::v_rotate_left<1>(src));
+    sum      = cv::v_add(sum, cv::v_rotate_left<2>(sum));
+    sum      = cv::v_add(sum, cv::v_rotate_left<4>(sum));
+
+    cv::v_uint32 v1;
+    cv::v_uint32 v2;
+    cv::v_expand(sum, v1, v2);
+    v1  = cv::v_add(v1, pre);
+    v2  = cv::v_add(v2, pre);
+    pre = cv::v_setall_u32(cv::v_extract_n<simdSize(cv::v_uint32) - 1>(v2));
+
+    cv::v_float64 v3;
+    cv::v_float64 v4;
+    expand(cv::v_reinterpret_as_s32(v1), v3, v4);
+    cv::v_store(dst, cv::v_add(v3, cv::v_load(prevDst)));
+    cv::v_store(dst + simdSize(cv::v_float64),
+                cv::v_add(v4, cv::v_load(prevDst + simdSize(cv::v_float64))));
+
+    expand(cv::v_reinterpret_as_s32(v2), v3, v4);
+    cv::v_store(dst + simdSize(cv::v_float64) * 2,
+                cv::v_add(v3, cv::v_load(prevDst + simdSize(cv::v_float64) * 2)));
+    cv::v_store(dst + simdSize(cv::v_float64) * 3,
+                cv::v_add(v4, cv::v_load(prevDst + simdSize(cv::v_float64) * 3)));
+}
+
+inline void integralSqSum(cv::v_uint16 &src, double *dst, double *prevDst, cv::v_uint32 &pre) {
+    cv::v_uint32 v1;
+    cv::v_uint32 v2;
+    cv::v_expand(src, v1, v2);
+
+    {
+        auto         shift1 = cv::v_rotate_left<1>(src);
+        cv::v_uint32 v3;
+        cv::v_uint32 v4;
+        cv::v_expand(shift1, v3, v4);
+
+        v1 = cv::v_add(v1, v3);
+        v2 = cv::v_add(v2, v4);
+
+        v4 = cv::v_extract<2>(v1, v2);
+        v2 = cv::v_add(v2, v4);
+
+        v3 = cv::v_rotate_left<2>(v1);
+        v1 = cv::v_add(v1, v3);
+
+        v1 = cv::v_add(v1, pre);
+        v2 = cv::v_add(v2, v1);
+
+        pre = cv::v_setall_u32(cv::v_extract_n<simdSize(cv::v_uint32) - 1>(v2));
+    }
+
+    cv::v_float64 v3;
+    cv::v_float64 v4;
+    expand(cv::v_reinterpret_as_s32(v1), v3, v4);
+    cv::v_store(dst, cv::v_add(v3, cv::v_load(prevDst)));
+    cv::v_store(dst + simdSize(cv::v_float64),
+                cv::v_add(v4, cv::v_load(prevDst + simdSize(cv::v_float64))));
+
+    expand(cv::v_reinterpret_as_s32(v2), v3, v4);
+    cv::v_store(dst + simdSize(cv::v_float64) * 2,
+                cv::v_add(v3, cv::v_load(prevDst + simdSize(cv::v_float64) * 2)));
+    cv::v_store(dst + simdSize(cv::v_float64) * 3,
+                cv::v_add(v4, cv::v_load(prevDst + simdSize(cv::v_float64) * 3)));
+}
+
+/*
+inline void integralSqSum(cv::v_uint32 &src, double *dst, double *prevDst, cv::v_uint32 &pre) {
+    src += cv::v_rotate_left<1>(src);
+    src += cv::v_rotate_left<2>(src);
+    src += pre;
+    pre  = cv::v_setall_u32(cv::v_extract_n<simdSize(cv::v_uint32) - 1>(src));
+
+    cv::v_float64 v1;
+    cv::v_float64 v2;
+    expand(cv::v_reinterpret_as_s32(src), v1, v2);
+
+    cv::v_store(dst, v1 + cv::v_load(prevDst));
+    cv::v_store(dst + simdSize(cv::v_float64), v2 + cv::v_load(prevDst +
+simdSize(cv::v_float64)));
+}
+
+inline void integralSqSum(cv::v_uint16 &src, double *dst, double *prevDst, cv::v_uint32 &pre) {
+    cv::v_uint32 v1;
+    cv::v_uint32 v2;
+    cv::v_expand(src, v1, v2);
+    integralSqSum(v1, dst, prevDst, pre);
+    integralSqSum(v2, dst + simdSize(cv::v_uint32), prevDst + simdSize(cv::v_uint32),
+pre);
+}
+*/
+
+inline void integralSum(const cv::v_uint16 &v1, const cv::v_uint16 &v2, double *dst,
+                        const double *prevDst, cv::v_uint32 &pre) {
+    integralSum(v1, dst, prevDst, pre);
+    integralSum(v2, dst + simdSize(cv::v_uint16), prevDst + simdSize(cv::v_uint16), pre);
+}
+
+inline void integralSqSum(cv::v_uint16 &v1, cv::v_uint16 &v2, double *dst, double *prevDst,
+                          cv::v_uint32 &pre) {
+    v1 = cv::v_mul_wrap(v1, v1);
+    v2 = cv::v_mul_wrap(v2, v2);
+
+    integralSqSum(v1, dst, prevDst, pre);
+    integralSqSum(v2, dst + simdSize(cv::v_uint16), prevDst + simdSize(cv::v_uint16), pre);
+}
+
+void integralSimd(const cv::Mat &src, cv::Mat &sum, cv::Mat &sqSum) {
+    const auto size = src.size() + cv::Size(1, 1);
+    sum.create(size, CV_64FC1);
+    sqSum.create(size, CV_64FC1);
+    memset(sum.data, 0, sum.step[ 0 ]);
+    memset(sqSum.data, 0, sqSum.step[ 0 ]);
+
+    const auto *srcStart   = src.data;
+    const auto  srcStep    = src.step[ 0 ];
+    auto       *sumStart   = reinterpret_cast<double *>(sum.data) + sum.step1() + 1;
+    const auto  sumStep    = sum.step[ 0 ] / sum.step[ 1 ];
+    auto       *sqSumStart = reinterpret_cast<double *>(sqSum.data) + sqSum.step1() + 1;
+    const auto  sqSumStep  = sqSum.step[ 0 ] / sqSum.step[ 1 ];
+    const auto  end        = size.width - simdSize(cv::v_uint8);
+    for (int y = 0; y < src.rows; y++) {
+        auto       *srcPtr    = srcStart + srcStep * y;
+        auto       *sumPtr    = sumStart + sumStep * y;
+        const auto *preSumPtr = sumStart + sumStep * (y - 1);
+        sumPtr[ -1 ]          = 0;
+
+        cv::v_uint32 prevSum = cv::vx_setzero_u32();
+        for (int x = 0; x < end; x += simdSize(cv::v_uint8)) {
+            cv::v_uint16 v1;
+            cv::v_uint16 v2;
+            cv::v_expand(cv::v_load(srcPtr + x), v1, v2);
+
+            integralSum(v1, v2, sumPtr + x, preSumPtr + x, prevSum);
+        }
+    }
+
+    for (int y = 0; y < src.rows; y++) {
+        auto *srcPtr      = srcStart + srcStep * y;
+        auto *sqSumPtr    = sqSumStart + sqSumStep * y;
+        auto *preSqSumPtr = sqSumStart + sqSumStep * (y - 1);
+        sqSumPtr[ -1 ]    = 0;
+
+        cv::v_uint32 prevSqSum = cv::vx_setzero_u32();
+        for (int x = 0; x < end; x += simdSize(cv::v_uint8)) {
+            cv::v_uint16 v1;
+            cv::v_uint16 v2;
+            cv::v_expand(cv::v_load(srcPtr + x), v1, v2);
+
+            integralSqSum(v1, v2, sqSumPtr + x, preSqSumPtr + x, prevSqSum);
+        }
+    }
+
+    const auto start = src.cols - src.cols % simdSize(cv::v_uint8);
+    for (int y = 0; y < src.rows; y++) {
+        auto       *srcPtr      = srcStart + srcStep * y;
+        auto       *sumPtr      = sumStart + sumStep * y;
+        auto       *sqSumPtr    = sqSumStart + sqSumStep * y;
+        const auto *preSumPtr   = sumStart + sumStep * (y - 1);
+        const auto *preSqSumPtr = sqSumStart + sqSumStep * (y - 1);
+        for (int x = start; x < src.cols; x++) {
+            const auto val   = srcPtr[ x ];
+            const auto sqVal = val * val;
+
+            sumPtr[ x ]   = sumPtr[ x - 1 ] + val + preSumPtr[ x ] - preSumPtr[ x - 1 ];
+            sqSumPtr[ x ] = sqSumPtr[ x - 1 ] + sqVal + preSqSumPtr[ x ] - preSqSumPtr[ x - 1 ];
+        }
+    }
+}

integral.h

@@ -0,0 +1,8 @@
+#ifndef INTEGRAL_H
+#define INTEGRAL_H
+
+#include <opencv2/opencv.hpp>
+
+void integralSimd(const cv::Mat &src, cv::Mat &sum, cv::Mat &sqSum);
+
+#endif // INTEGRAL_H

main.cpp

@@ -1,59 +1,141 @@
 #include "grayMatch.h"
 
+#include <fstream>
 #include <iostream>
+#include <opencv2/core/utility.hpp>
 #include <opencv2/opencv.hpp>
 
-int main() {
-    auto src =
-        cv::imread("C:/Users/qiuyong/Desktop/test/template/model3.bmp", cv::IMREAD_GRAYSCALE);
-    auto dst =
-        cv::imread("C:/Users/qiuyong/Desktop/test/template/model3_src1.bmp", cv::IMREAD_GRAYSCALE);
+int main(int argc, const char *argv[]) {
+    const std::string keys = "{model m || model image}"
+                             "{scene s || scene image}"
+                             "{view v || view result}"
+                             "{threshold t | 0.7 | match minium score}"
+                             "{bench b || match benchmark}"
+                             "{help h || print this help}";
 
-    auto t0    = cv::getTickCount();
-    auto model = trainModel(src.data, src.cols, src.rows, src.channels(), int(src.step), 0, 0,
-                            src.cols, src.rows, -1);
-    auto t1    = cv::getTickCount();
+    cv::CommandLineParser cmd(argc, argv, keys);
+    if (!cmd.check()) {
+        cmd.printErrors();
+        return -1;
+    }
 
-    int size;
-    serialize(model, nullptr, &size);
-    std::vector<uchar> buffer(size);
-    serialize(model, buffer.data(), &size);
+    if (cmd.has("help")) {
+        cmd.printMessage();
+        return 0;
+    }
 
-    freeModel(&model);
+    auto srcFile = std::string(IMG_DIR) + "/3.bmp";
+    auto dstFile = std::string(IMG_DIR) + "/h.bmp";
+    if (cmd.has("model"))
+        srcFile = cmd.get<std::string>("model");
+    if (cmd.has("scene"))
+        dstFile = cmd.get<std::string>("scene");
 
-    model = deserialize(buffer.data(), int(buffer.size()));
+    auto src = cv::imread(srcFile, cv::IMREAD_GRAYSCALE);
+    auto dst = cv::imread(dstFile, cv::IMREAD_GRAYSCALE);
+    if (src.empty() || dst.empty()) {
+        return -1;
+    }
 
-    int               num = 70;
-    std::vector<Pose> poses(num);
+    const std::string modelName("model.bin");
+    {
+        auto t0    = cv::getTickCount();
+        auto model = trainModel(src.data, src.cols, src.rows, src.channels(),
+                                static_cast<int>(src.step), 0, 0, src.cols, src.rows, -1);
+        auto t1    = cv::getTickCount();
 
-    auto t2 = cv::getTickCount();
-    matchModel(dst.data, dst.cols, dst.rows, dst.channels(), int(dst.step), 0, 0, dst.cols,
-               dst.rows, model, &num, poses.data(), -1, 0, 360, 0, 0.5, 1);
-    auto t3 = cv::getTickCount();
+        // get size
+        int size;
+        serialize(model, nullptr, &size);
 
-    auto trainCost = double(t1 - t0) / cv::getTickFrequency();
-    auto matchCost = double(t3 - t2) / cv::getTickFrequency();
-    std::cout << "train(s):" << trainCost << " match(s):" << matchCost << std::endl;
+        // serialize to buffer
+        std::vector<uchar> buffer(size);
+        serialize(model, buffer.data(), &size);
 
-    cv::Mat color;
-    cv::cvtColor(dst, color, cv::COLOR_GRAY2RGB);
-    for (int i = 0; i < num; i++) {
-        auto           &pose = poses[ i ];
-        cv::RotatedRect rect(cv::Point2f(pose.x, pose.y), src.size(), -pose.angle);
+        // write to file
+        std::ofstream ofs(modelName, std::ios::binary | std::ios::out);
+        if (!ofs.is_open()) {
+            return -1;
+        }
+        ofs.write(reinterpret_cast<const char *>(buffer.data()), size);
 
-        cv::Point2f pts[ 4 ];
-        rect.points(pts);
+        freeModel(&model);
 
-        cv::line(color, pts[ 0 ], pts[ 1 ], cv::Scalar(255, 0, 0), 1, cv::LINE_AA);
-        cv::line(color, pts[ 1 ], pts[ 2 ], cv::Scalar(255, 0, 0), 1, cv::LINE_AA);
-        cv::line(color, pts[ 2 ], pts[ 3 ], cv::Scalar(255, 0, 0), 1, cv::LINE_AA);
-        cv::line(color, pts[ 3 ], pts[ 0 ], cv::Scalar(255, 0, 0), 1, cv::LINE_AA);
+        auto trainCost = static_cast<double>(t1 - t0) / cv::getTickFrequency();
+        std::cout << "train(s):" << trainCost << std::endl;
+    }
 
+    int               count = 70;
+    std::vector<Pose> poses(count);
+    Model_t           model;
+    auto              score = cmd.get<float>("threshold");
+    {
+        // open file
+        std::ifstream ifs(modelName, std::ios::binary | std::ios::in);
+        if (!ifs.is_open()) {
+            return -2;
+        }
+
+        // get size
+        ifs.seekg(0, std::ios::end);
+        auto size = ifs.tellg();
+        ifs.seekg(0, std::ios::beg);
+
+        // read to buffer
+        std::vector<uchar> buffer(size);
+        ifs.read(reinterpret_cast<char *>(buffer.data()), size);
+
+        // deserialize from buffer
+        model = deserialize(buffer.data(), static_cast<int>(buffer.size()));
+
+        auto t2 = cv::getTickCount();
+        matchModel(dst.data, dst.cols, dst.rows, dst.channels(), static_cast<int>(dst.step), 0, 0,
+                   dst.cols, dst.rows, model, &count, poses.data(), -1, 0, 360, 0, score, 1);
+        auto t3 = cv::getTickCount();
+
+        auto matchCost = static_cast<double>(t3 - t2) / cv::getTickFrequency();
+        std::cout << "match(s):" << matchCost << std::endl;
+    }
+
+    for (int i = 0; i < count; i++) {
+        const auto &pose = poses[ i ];
         std::cout << pose.x << "," << pose.y << "," << pose.angle << "," << pose.score << std::endl;
     }
 
-    cv::imshow("img", color);
-    cv::waitKey();
+    if (cmd.has("bench")) {
+        constexpr int times = 100;
+
+        auto start = cv::getTickCount();
+        for (int i = 0; i < times; i++) {
+            matchModel(dst.data, dst.cols, dst.rows, dst.channels(), static_cast<int>(dst.step), 0,
+                       0, dst.cols, dst.rows, model, &count, poses.data(), -1, 0, 360, 0, score, 1);
+            count = 70;
+        }
+        auto end = cv::getTickCount();
+
+        const auto cost = static_cast<double>(end - start) / cv::getTickFrequency() / times;
+        std::cout << "match bench avg(s):" << cost << std::endl;
+    }
+
+    if (cmd.has("view")) {
+        cv::Mat color;
+        cv::cvtColor(dst, color, cv::COLOR_GRAY2RGB);
+        for (int i = 0; i < count; i++) {
+            const auto     &pose = poses[ i ];
+            cv::RotatedRect rect(cv::Point2f(pose.x, pose.y), src.size(), -pose.angle);
+
+            cv::Point2f pts[ 4 ];
+            rect.points(pts);
+
+            cv::line(color, pts[ 0 ], pts[ 1 ], cv::Scalar(255, 0, 0), 1, cv::LINE_AA);
+            cv::line(color, pts[ 1 ], pts[ 2 ], cv::Scalar(255, 0, 0), 1, cv::LINE_AA);
+            cv::line(color, pts[ 2 ], pts[ 3 ], cv::Scalar(255, 0, 0), 1, cv::LINE_AA);
+            cv::line(color, pts[ 3 ], pts[ 0 ], cv::Scalar(255, 0, 0), 1, cv::LINE_AA);
+        }
+
+        cv::imshow("img", color);
+        cv::waitKey();
+    }
 
     return 0;
 }

privateType.h

@@ -26,9 +26,16 @@ struct Model {
     }
 
     void reserve(const std::size_t size) {
+        pyramids.reserve(size);
         normal.reserve(size);
         invArea.reserve(size);
         mean.reserve(size);
         equal1.reserve(size);
     }
-};
+};
+
+#if CV_VERSION_MAJOR >= 4 && CV_VERSION_MINOR >= 8
+#define simdSize(type) cv::VTraits<type>::nlanes
+#else
+#define simdSize(type) type::nlanes
+#endif

serialize.cpp

@@ -1,12 +1,16 @@
 #include "grayMatch.h"
 #include "privateType.h"
 
+#include <opencv2/core/hal/intrin.hpp>
+
 class Buffer {
 public:
-    Buffer(int size_, unsigned char *data_)
+    Buffer(const int size_, unsigned char *data_)
         : m_size(size_)
         , m_data(data_) {}
 
+    virtual ~Buffer() = default;
+
     virtual void operator&(uchar &val)                   = 0;
     virtual void operator&(std::vector<cv::Mat> &val)    = 0;
     virtual void operator&(std::vector<cv::Scalar> &val) = 0;
@@ -22,7 +26,7 @@ public:
         this->operator&(val.borderColor);
     }
 
-    int count() const {
+    [[nodiscard]] int count() const {
         return m_size;
     }
 
@@ -31,34 +35,29 @@ protected:
     unsigned char *m_data = nullptr;
 };
 
-class WriteOperationBase {
-public:
-    virtual void write(void *dst, void *src, int size) = 0;
-};
-
-void binWrite(void *dst, void *src, int size) {
+void binWrite(void *const dst, const void *src, const int size) {
     memcpy(dst, src, size);
 }
 
-void fakeWrite(void *dst, void *src, int size) {
-    (void)(dst);
-    (void)(src);
-    (void)(size);
+void fakeWrite(void *const dst, const void *src, const int size) {
+    (void)dst;
+    (void)src;
+    (void)size;
 }
 
-using Write = void (*)(void *, void *, int);
+using Write = void (*)(void *, const void *, int);
 
-template <Write write> class OutBuffer : public Buffer {
+template <Write write> class OutBuffer final : public Buffer {
 public:
-    OutBuffer(unsigned char *data_)
+    explicit OutBuffer(unsigned char *const data_)
         : Buffer(0, data_) {}
 
-    void operator&(uchar &val) final {
+    void operator&(uchar &val) override {
         write(m_data + m_size, &val, sizeof(val));
         m_size += static_cast<int>(sizeof(val));
     }
-    void operator&(std::vector<cv::Mat> &val) final {
-        int size = static_cast<int>(val.size());
+    void operator&(std::vector<cv::Mat> &val) override {
+        const int size = static_cast<int>(val.size());
         write(m_data + m_size, &size, sizeof(size));
         m_size += static_cast<int>(sizeof(size));
 
@@ -78,8 +77,8 @@ public:
             m_size += val.cols;
         }
     }
-    void operator&(std::vector<cv::Scalar> &val) final {
-        int size = static_cast<int>(val.size());
+    void operator&(std::vector<cv::Scalar> &val) override {
+        const int size = static_cast<int>(val.size());
         write(m_data + m_size, &size, sizeof(size));
         m_size += static_cast<int>(sizeof(size));
 
@@ -87,20 +86,20 @@ public:
             writeElement(element);
         }
     }
-    void writeElement(cv::Scalar &val) {
+    void writeElement(const cv::Scalar &val) {
         write(m_data + m_size, val.val, sizeof(double) * 4);
         m_size += static_cast<int>(sizeof(double)) * 4;
     }
-    void operator&(std::vector<double> &val) final {
-        int size = static_cast<int>(val.size());
+    void operator&(std::vector<double> &val) override {
+        const int size = static_cast<int>(val.size());
         write(m_data + m_size, &size, sizeof(size));
         m_size += static_cast<int>(sizeof(size));
 
         write(m_data + m_size, val.data(), static_cast<int>(sizeof(double)) * size);
         m_size += static_cast<int>(sizeof(double)) * size;
     }
-    void operator&(std::vector<uchar> &val) final {
-        int size = static_cast<int>(val.size());
+    void operator&(std::vector<uchar> &val) override {
+        const int size = static_cast<int>(val.size());
         write(m_data + m_size, &size, sizeof(size));
         m_size += static_cast<int>(sizeof(size));
 
@@ -112,16 +111,16 @@ public:
 using SizeCountBuffer = OutBuffer<fakeWrite>;
 using WriteBuffer     = OutBuffer<binWrite>;
 
-class ReadBuffer : public Buffer {
+class ReadBuffer final : public Buffer {
 public:
-    ReadBuffer(unsigned char *data_)
+    explicit ReadBuffer(unsigned char *data_)
         : Buffer(0, data_) {}
 
-    void operator&(uchar &val) final {
+    void operator&(uchar &val) override {
         memcpy(&val, m_data + m_size, sizeof(uchar));
         m_size += static_cast<int>(sizeof(uchar));
     }
-    void operator&(std::vector<cv::Mat> &val) final {
+    void operator&(std::vector<cv::Mat> &val) override {
         int count = 0;
         memcpy(&count, m_data + m_size, sizeof(int));
         val.resize(count);
@@ -140,10 +139,17 @@ public:
         memcpy(&height, m_data + m_size, sizeof(int));
         m_size += static_cast<int>(sizeof(int));
 
-        val     = cv::Mat(cv::Size(width, height), CV_8UC1, m_data + m_size);
-        m_size += width * height;
+        const int  alignedWidth = static_cast<int>(cv::alignSize(width, simdSize(cv::v_uint8)));
+        const auto img          = cv::Mat::zeros(height, alignedWidth, CV_8UC1);
+        val                     = img(cv::Rect(0, 0, width, height));
+
+        for (int y = 0; y < height; y++) {
+            auto *ptr = val.ptr<uchar>(y);
+            memcpy(ptr, m_data + m_size, width);
+            m_size += width;
+        }
     }
-    void operator&(std::vector<cv::Scalar> &val) final {
+    void operator&(std::vector<cv::Scalar> &val) override {
         int count = 0;
         memcpy(&count, m_data + m_size, sizeof(int));
         val.resize(count);
@@ -157,7 +163,7 @@ public:
         memcpy(val.val, m_data + m_size, sizeof(double) * 4);
         m_size += static_cast<int>(sizeof(double)) * 4;
     }
-    void operator&(std::vector<double> &val) final {
+    void operator&(std::vector<double> &val) override {
         int count = 0;
         memcpy(&count, m_data + m_size, sizeof(int));
         val.resize(count);
@@ -166,7 +172,7 @@ public:
         memcpy(val.data(), m_data + m_size, sizeof(double) * count);
         m_size += static_cast<int>(sizeof(double)) * count;
     }
-    void operator&(std::vector<uchar> &val) final {
+    void operator&(std::vector<uchar> &val) override {
         int count = 0;
         memcpy(&count, m_data + m_size, sizeof(int));
         val.resize(count);
@@ -178,10 +184,10 @@ public:
 };
 
 void operation(Buffer *buf, Model &model) {
-    (*buf) & (model);
+    *buf &model;
 }
 
-bool serialize(const Model_t model, unsigned char *buffer, int *size) {
+bool serialize(Model *const model, unsigned char *buffer, int *size) {
     if (nullptr == size) {
         return false;
     }
@@ -191,15 +197,15 @@ bool serialize(const Model_t model, unsigned char *buffer, int *size) {
         return false;
     }
 
-    SizeCountBuffer countor(buffer);
-    operation(&countor, *model);
-    *size = countor.count();
+    SizeCountBuffer counter(buffer);
+    operation(&counter, *model);
 
     if (nullptr == buffer) {
+        *size = counter.count();
         return true;
     }
 
-    if (countor.count() > *size) {
+    if (counter.count() > *size) {
         *size = 0;
         return false;
     }
@@ -209,7 +215,7 @@ bool serialize(const Model_t model, unsigned char *buffer, int *size) {
     return true;
 }
 
-Model_t deserialize(unsigned char *buffer, int size) {
+Model_t deserialize(unsigned char *buffer, const int size) {
     if (size < 1 || nullptr == buffer) {
         return nullptr;
     }