diff options
| -rw-r--r-- | tracker-pt/point_extractor.cpp | 94 | ||||
| -rw-r--r-- | tracker-pt/point_extractor.h | 7 |
2 files changed, 94 insertions, 7 deletions
diff --git a/tracker-pt/point_extractor.cpp b/tracker-pt/point_extractor.cpp index a9f431d9..e09af056 100644 --- a/tracker-pt/point_extractor.cpp +++ b/tracker-pt/point_extractor.cpp @@ -8,6 +8,8 @@ #include "point_extractor.h" #include "compat/util.hpp" +#include "point_tracker.h" +#include <QDebug> #include <opencv2/videoio.hpp> @@ -19,6 +21,59 @@ using namespace types; +/* +http://en.wikipedia.org/wiki/Mean-shift +In this application the idea, is to eliminate any bias of the point estimate +which is introduced by the rather arbitrary thresholded area. One must recognize +that the thresholded area can only move in one pixel increments since it is +binary. Thus, its center of mass might make "jumps" as pixels are added/removed +from the thresholded area. +With mean-shift, a moving "window" or kernel is multiplied with the gray-scale +image, and the COM is calculated of the result. This is iterated where the +kernel center is set the previously computed COM. Thus, peaks in the image intensity +distribution "pull" the kernel towards themselves. Eventually it stops moving, i.e. +then the computed COM coincides with the kernel center. We hope that the +corresponding location is a good candidate for the extracted point. +The idea similar to the window scaling suggested in Berglund et al. "Fast, bias-free +algorithm for tracking single particles with variable size and shape." (2008). +*/ +static cv::Vec2d MeanShiftIteration(const cv::Mat &frame_gray, const cv::Vec2d ¤t_center, double filter_width) +{ + // Most amazingling this function runs faster with doubles than with floats. + const double s = 1.0 / filter_width; + + auto EvalFilter = [¤t_center, s](int row, int col) + { + double dx = (col - current_center[0])*s; + double dy = (row - current_center[1])*s; + double f = std::max(0.0, 1.0 - dx*dx - dy*dy); + return f; + }; + + double m = 0; + cv::Vec2d com(0.0, 0.0); + for (int i = 0; i < frame_gray.rows; i++) + { + auto frame_ptr = (uint8_t *)frame_gray.ptr(i); + for (int j = 0; j < frame_gray.cols; j++) + { + double val = frame_ptr[j]; + val = val * val; // taking the square wights brighter parts of the image stronger. + val *= EvalFilter(i, j); + m += val; + com[0] += j * val; + com[1] += i * val; + } + } + if (m > 0.1) + { + com *= 1.0 / m; + return com; + } + else + return current_center; +} + PointExtractor::PointExtractor() { blobs.reserve(max_blobs); @@ -135,12 +190,13 @@ void PointExtractor::extract_points(cv::Mat& frame, std::vector<vec2>& points) cnt++; } } - if (norm > 1e0) + if (norm > 0) { - const double radius = sqrt(cnt / M_PI); + const double radius = sqrt(cnt / M_PI), N = double(norm); if (radius > region_size_max || radius < region_size_min) continue; - blob b(radius, cv::Vec2d(m10 / norm, m01 / norm), norm/sqrt(double(cnt))); + + blob b(radius, cv::Vec2d(m10 / N, m01 / N), N/sqrt(double(cnt)), rect); blobs.push_back(b); { char buf[64]; @@ -162,6 +218,33 @@ end: sort(blobs.begin(), blobs.end(), [](const blob& b1, const blob& b2) -> bool { return b2.brightness < b1.brightness; }); + for (idx = 0; idx < std::min(PointModel::N_POINTS, unsigned(blobs.size())); ++idx) + { + blob &b = blobs[idx]; + cv::Rect rect = b.rect; + rect.x -= rect.width / 2; + rect.y -= rect.height / 2; + rect.width *= 2; + rect.height *= 2; + rect &= cv::Rect(0, 0, W, H); // crop at frame boundaries + + cv::Mat frame_roi = frame_gray(rect); + + const double kernel_radius = b.radius * 1.5; + cv::Vec2d pos(b.pos[0] - rect.x, b.pos[1] - rect.y); // position relative to ROI. + + for (int iter = 0; iter < 10; ++iter) + { + cv::Vec2d com_new = MeanShiftIteration(frame_roi, pos, kernel_radius); + cv::Vec2d delta = com_new - b.pos; + pos = com_new; + if (delta.dot(delta) < 1.0e-2) break; + } + + b.pos[0] = pos[0] + rect.x; + b.pos[1] = pos[1] + rect.y; + } + // End of mean shift code. At this point, blob positions are updated which hopefully less noisy less biased values. points.reserve(max_blobs); points.clear(); @@ -171,3 +254,8 @@ end: points.push_back(p); } } + +PointExtractor::blob::blob(double radius, const cv::Vec2d& pos, double brightness, cv::Rect& rect) : radius(radius), brightness(brightness), pos(pos), rect(rect) +{ + //qDebug() << "radius" << radius << "pos" << pos[0] << pos[1]; +} diff --git a/tracker-pt/point_extractor.h b/tracker-pt/point_extractor.h index 775e4292..ad350344 100644 --- a/tracker-pt/point_extractor.h +++ b/tracker-pt/point_extractor.h @@ -42,10 +42,9 @@ private: { double radius, brightness; vec2 pos; - blob(double radius, const cv::Vec2d& pos, double brightness) : radius(radius), brightness(brightness), pos(pos) - { - //qDebug() << "radius" << radius << "pos" << pos[0] << pos[1]; - } + cv::Rect rect; + + blob(double radius, const cv::Vec2d& pos, double brightness, cv::Rect &rect); }; std::vector<blob> blobs; |