diff options
Diffstat (limited to 'tracker-pt/module/point_extractor.cpp')
-rw-r--r-- | tracker-pt/module/point_extractor.cpp | 307 |
1 files changed, 194 insertions, 113 deletions
diff --git a/tracker-pt/module/point_extractor.cpp b/tracker-pt/module/point_extractor.cpp index b67c4036..3329fafc 100644 --- a/tracker-pt/module/point_extractor.cpp +++ b/tracker-pt/module/point_extractor.cpp @@ -9,11 +9,11 @@ #include "point_extractor.h" #include "point_tracker.h" #include "frame.hpp" - #include "cv/numeric.hpp" #include "compat/math.hpp" +#include "compat/math-imports.hpp" -#include <opencv2/videoio.hpp> +#include <opencv2/imgproc.hpp> #undef PREVIEW //#define PREVIEW @@ -29,12 +29,13 @@ #include <QDebug> -using namespace types; -using namespace pt_module; +using namespace numeric_types; + +// meanshift code written by Michael Welter /* http://en.wikipedia.org/wiki/Mean-shift -In this application the idea, is to eliminate any bias of the point estimate +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 @@ -43,31 +44,28 @@ 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 +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 +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 vec2 ¤t_center, f filter_width) +static vec2 MeanShiftIteration(const cv::Mat1b &frame_gray, const vec2 ¤t_center, f filter_width) { - // Most amazingly this function runs faster with doubles than with floats. - const f s = 1.0 / filter_width; + const f s = 1 / filter_width; f m = 0; - vec2 com { 0, 0 }; + vec2 com { 0, 0 }; for (int i = 0; i < frame_gray.rows; i++) { - auto frame_ptr = (uint8_t const* restrict_ptr)frame_gray.ptr(i); + uint8_t const* const __restrict frame_ptr = frame_gray.ptr(i); for (int j = 0; j < frame_gray.cols; j++) { f val = frame_ptr[j]; - val = val * val; // taking the square wights brighter parts of the image stronger. - { - f dx = (j - current_center[0])*s; - f dy = (i - current_center[1])*s; - f f = std::fmax(0, 1 - dx*dx - dy*dy); - val *= f; - } + val = val * val; // taking the square weighs brighter parts of the image stronger. + f dx = (j - current_center[0])*s; + f dy = (i - current_center[1])*s; + f max = std::fmax(f(0), 1 - dx*dx - dy*dy); + val *= max; m += val; com[0] += j * val; com[1] += i * val; @@ -75,13 +73,15 @@ static cv::Vec2d MeanShiftIteration(const cv::Mat &frame_gray, const vec2 &curre } if (m > f(.1)) { - com *= f(1) / m; + com *= 1 / m; return com; } else return current_center; } +namespace pt_module { + PointExtractor::PointExtractor(const QString& module_name) : s(module_name) { blobs.reserve(max_blobs); @@ -89,24 +89,19 @@ PointExtractor::PointExtractor(const QString& module_name) : s(module_name) void PointExtractor::ensure_channel_buffers(const cv::Mat& orig_frame) { - if (ch[0].rows != orig_frame.rows || ch[0].cols != orig_frame.cols) - for (unsigned k = 0; k < 3; k++) - ch[k] = cv::Mat1b(orig_frame.rows, orig_frame.cols); + for (cv::Mat1b& x : ch) + x.create(orig_frame.rows, orig_frame.cols); } void PointExtractor::ensure_buffers(const cv::Mat& frame) { const int W = frame.cols, H = frame.rows; - if (frame_gray.rows != W || frame_gray.cols != H) - { - frame_gray = cv::Mat1b(H, W); - frame_bin = cv::Mat1b(H, W); - frame_blobs = cv::Mat1b(H, W); - } + frame_gray.create(H, W); + frame_bin.create(H, W); } -void PointExtractor::extract_single_channel(const cv::Mat& orig_frame, int idx, cv::Mat& dest) +void PointExtractor::extract_single_channel(const cv::Mat& orig_frame, int idx, cv::Mat1b& dest) { ensure_channel_buffers(orig_frame); @@ -117,17 +112,50 @@ void PointExtractor::extract_single_channel(const cv::Mat& orig_frame, int idx, cv::mixChannels(&orig_frame, 1, &dest, 1, from_to, 1); } -void PointExtractor::extract_channels(const cv::Mat& orig_frame, const int* order, int order_npairs) +void PointExtractor::filter_single_channel(const cv::Mat& orig_frame, float r, float g, float b, bool overexp, cv::Mat1b& dest) { ensure_channel_buffers(orig_frame); - cv::mixChannels(&orig_frame, 1, (cv::Mat*) ch, order_npairs, order, order_npairs); + // just filter for colour or also include overexposed regions? + if (!overexp) + cv::transform(orig_frame, dest, cv::Mat(cv::Matx13f(b, g, r))); + else + { + for (int i = 0; i < orig_frame.rows; i++) + { + cv::Vec3b const* const __restrict orig_ptr = orig_frame.ptr<cv::Vec3b>(i); + uint8_t* const __restrict dest_ptr = dest.ptr(i); + for (int j = 0; j < orig_frame.cols; j++) + { + // get the intensity of the key color (i.e. +ve coefficients) + uchar blue = orig_ptr[j][0], green = orig_ptr[j][1], red = orig_ptr[j][2]; + float key = std::max(b, 0.0f) * blue + std::max(g, 0.0f) * green + std::max(r, 0.0f) * red; + // get the intensity of the non-key color (i.e. -ve coefficients) + float nonkey = std::max(-b, 0.0f) * blue + std::max(-g, 0.0f) * green + std::max(-r, 0.0f) * red; + // the result is key color minus non-key color inversely weighted by key colour intensity + dest_ptr[j] = std::max(0.0f, std::min(255.0f, key - (255.0f - key) / 255.0f * nonkey)); + } + } + } } void PointExtractor::color_to_grayscale(const cv::Mat& frame, cv::Mat1b& output) { + if (frame.channels() == 1) + { + output.create(frame.rows, frame.cols); + frame.copyTo(output); + return; + } + + const float half_chr_key_str = *s.chroma_key_strength * 0.5; switch (s.blob_color) { + case pt_color_green_only: + { + extract_single_channel(frame, 1, output); + break; + } case pt_color_blue_only: { extract_single_channel(frame, 0, output); @@ -138,18 +166,44 @@ void PointExtractor::color_to_grayscale(const cv::Mat& frame, cv::Mat1b& output) extract_single_channel(frame, 2, output); break; } - case pt_color_average: + case pt_color_red_chromakey: + { + filter_single_channel(frame, 1, -half_chr_key_str, -half_chr_key_str, s.chroma_key_overexposed, output); + break; + } + case pt_color_green_chromakey: { - const int W = frame.cols, H = frame.rows; - const cv::Mat tmp = frame.reshape(1, W * H); - cv::Mat output_ = output.reshape(1, W * H); - cv::reduce(tmp, output_, 1, cv::REDUCE_AVG); + filter_single_channel(frame, -half_chr_key_str, 1, -half_chr_key_str, s.chroma_key_overexposed, output); break; } + case pt_color_blue_chromakey: + { + filter_single_channel(frame, -half_chr_key_str, -half_chr_key_str, 1, s.chroma_key_overexposed, output); + break; + } + case pt_color_cyan_chromakey: + { + filter_single_channel(frame, -*s.chroma_key_strength, 0.5, 0.5, s.chroma_key_overexposed, output); + break; + } + case pt_color_yellow_chromakey: + { + filter_single_channel(frame, 0.5, 0.5, -*s.chroma_key_strength, s.chroma_key_overexposed, output); + break; + } + case pt_color_magenta_chromakey: + { + filter_single_channel(frame, 0.5, -*s.chroma_key_strength, 0.5, s.chroma_key_overexposed, output); + break; + } + case pt_color_hardware: + eval_once(qDebug() << "camera driver doesn't support grayscale"); + goto do_grayscale; default: - once_only(qDebug() << "wrong pt_color_type enum value" << int(s.blob_color)); - /*FALLTHROUGH*/ - case pt_color_natural: + eval_once(qDebug() << "wrong pt_color_type enum value" << int(s.blob_color)); + [[fallthrough]]; + case pt_color_bt709: +do_grayscale: cv::cvtColor(frame, output, cv::COLOR_BGR2GRAY); break; } @@ -175,31 +229,97 @@ void PointExtractor::threshold_image(const cv::Mat& frame_gray, cv::Mat1b& outpu cv::noArray(), hist, 1, - (int const*) &hist_size, + &hist_size, &ranges); - const f radius = (f) threshold_radius_value(frame_gray.cols, frame_gray.rows, threshold_slider_value); + const f radius = threshold_radius_value(frame_gray.cols, frame_gray.rows, threshold_slider_value); - auto ptr = (float const* const restrict_ptr) hist.ptr(0); - const unsigned area = uround(3 * M_PI * radius*radius); + float const* const __restrict ptr = hist.ptr<float>(0); + const unsigned area = unsigned(iround(3 * pi * radius*radius)); const unsigned sz = unsigned(hist.cols * hist.rows); - unsigned thres = 32; - for (unsigned i = sz-1, cnt = 0; i > 32; i--) + unsigned thres = 1; + for (unsigned i = sz-1, cnt = 0; i > 1; i--) { - cnt += ptr[i]; + cnt += (unsigned)ptr[i]; if (cnt >= area) break; thres = i; } - cv::threshold(frame_gray, output, thres, 255, CV_THRESH_BINARY); + cv::threshold(frame_gray, output, thres, 255, cv::THRESH_BINARY); } } -void PointExtractor::extract_points(const pt_frame& frame_, pt_preview& preview_frame_, std::vector<vec2>& points) +static void draw_blobs(cv::Mat& preview_frame, const blob* blobs, unsigned nblobs, const cv::Size& size) +{ + for (unsigned k = 0; k < nblobs; k++) + { + const blob& b = blobs[k]; + + if (b.radius < 0) + continue; + + const f dpi = preview_frame.cols / f(320); + const f offx = 10 * dpi, offy = f(7.5) * dpi; + + const f cx = preview_frame.cols / f(size.width), + cy = preview_frame.rows / f(size.height), + c = std::fmax(f(1), cx+cy)/2; + + cv::Point p(iround(b.pos[0] * cx), iround(b.pos[1] * cy)); + + auto outline_color = k >= PointModel::N_POINTS + ? cv::Scalar(192, 192, 192) + : cv::Scalar(255, 255, 0); + + cv::ellipse(preview_frame, p, + {iround(b.rect.width/(f)2+2*c), iround(b.rect.height/(f)2+2*c)}, + 0, 0, 360, outline_color, iround(dpi), cv::LINE_AA); + + char buf[16]; + std::snprintf(buf, sizeof(buf), "%.2fpx", (double)b.radius); + + auto text_color = k >= PointModel::N_POINTS + ? cv::Scalar(160, 160, 160) + : cv::Scalar(0, 0, 255); + + cv::Point pos(iround(b.pos[0]*cx+offx), iround(b.pos[1]*cy+offy)); + cv::putText(preview_frame, buf, pos, + cv::FONT_HERSHEY_PLAIN, iround(dpi), text_color, + 1); + } +} + +static vec2 meanshift_initial_guess(const cv::Rect rect, cv::Mat& frame_roi) +{ + vec2 ret = {rect.width/(f)2, rect.height/(f)2}; + + // compute center initial guess + double ynorm = 0, xnorm = 0, y = 0, x = 0; + for (int j = 0; j < rect.height; j++) + { + const unsigned char* __restrict ptr = frame_roi.ptr<unsigned char>(j); + for (int i = 0; i < rect.width; i++) + { + double val = ptr[i] * 1./255; + x += i * val; + y += j * val; + xnorm += val; + ynorm += val; + } + } + constexpr double eps = 1e-4; + if (xnorm > eps && ynorm > eps) + ret = { (f)(x / xnorm), (f)(y / ynorm) }; + return ret; +} + +void PointExtractor::extract_points(const pt_frame& frame_, + pt_preview& preview_frame_, + bool preview_visible, + std::vector<vec2>& points) { const cv::Mat& frame = frame_.as_const<Frame>()->mat; - cv::Mat& preview_frame = *preview_frame_.as<Preview>(); ensure_buffers(frame); color_to_grayscale(frame, frame_gray); @@ -211,24 +331,24 @@ void PointExtractor::extract_points(const pt_frame& frame_, pt_preview& preview_ threshold_image(frame_gray, frame_bin); - blobs.clear(); - frame_bin.copyTo(frame_blobs); - - const f region_size_min = s.min_point_size; - const f region_size_max = s.max_point_size; + const f region_size_min = (f)s.min_point_size; + const f region_size_max = (f)s.max_point_size; unsigned idx = 0; - for (int y=0; y < frame_blobs.rows; y++) + + blobs.clear(); + + for (int y=0; y < frame_bin.rows; y++) { - const unsigned char* ptr_bin = frame_blobs.ptr(y); - for (int x=0; x < frame_blobs.cols; x++) + const unsigned char* __restrict ptr_bin = frame_bin.ptr(y); + for (int x=0; x < frame_bin.cols; x++) { if (ptr_bin[x] != 255) continue; idx = blobs.size() + 1; cv::Rect rect; - cv::floodFill(frame_blobs, + cv::floodFill(frame_bin, cv::Point(x,y), cv::Scalar(idx), &rect, @@ -244,8 +364,8 @@ void PointExtractor::extract_points(const pt_frame& frame_, pt_preview& preview_ for (int i=rect.y; i < ymax; i++) { - unsigned char const* const restrict_ptr ptr_blobs = frame_blobs.ptr(i); - unsigned char const* const restrict_ptr ptr_gray = frame_gray.ptr(i); + unsigned char const* const __restrict ptr_blobs = frame_bin.ptr(i); + unsigned char const* const __restrict ptr_gray = frame_gray.ptr(i); for (int j=rect.x; j < xmax; j++) { if (ptr_blobs[j] != idx) @@ -257,12 +377,12 @@ void PointExtractor::extract_points(const pt_frame& frame_, pt_preview& preview_ } } - const double radius = std::sqrt(cnt / M_PI); + const f radius = std::sqrt((f)cnt) / std::sqrt(pi); if (radius > region_size_max || radius < region_size_min) continue; blobs.emplace_back(radius, - vec2(rect.width/2., rect.height/2.), + vec2(rect.width/f(2), rect.height/f(2)), std::pow(f(norm), f(1.1))/cnt, rect); @@ -272,9 +392,7 @@ void PointExtractor::extract_points(const pt_frame& frame_, pt_preview& preview_ // XXX we could go to the next scanline unless the points are really small. // i'd expect each point being present on at least one unique scanline // but it turns out some people are using 2px points -sh 20180110 -#if BROKEN && 0 - break; -#endif + //break; } } end: @@ -288,29 +406,22 @@ end: for (idx = 0; idx < sz; ++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 - + blob& b = blobs[idx]; + cv::Rect rect = b.rect & cv::Rect(0, 0, W, H); // crop at frame boundaries cv::Mat frame_roi = frame_gray(rect); // smaller values mean more changes. 1 makes too many changes while 1.5 makes about .1 static constexpr f radius_c = f(1.75); const f kernel_radius = b.radius * radius_c; - vec2 pos(rect.width/2., rect.height/2.); // position relative to ROI. + vec2 pos = meanshift_initial_guess(rect, frame_roi); // position relative to ROI. for (int iter = 0; iter < 10; ++iter) { vec2 com_new = MeanShiftIteration(frame_roi, pos, kernel_radius); vec2 delta = com_new - pos; pos = com_new; - if (delta.dot(delta) < 1e-2) + if (delta.dot(delta) < f(1e-3)) break; } @@ -318,43 +429,11 @@ end: b.pos[1] = pos[1] + rect.y; } - for (unsigned k = 0; k < blobs.size(); k++) - { - blob& b = blobs[k]; - - const f dpi = preview_frame.cols / f(320); - const f offx = 10 * dpi, offy = 7.5 * dpi; - - const f cx = preview_frame.cols / f(frame.cols), - cy = preview_frame.rows / f(frame.rows), - c_ = (cx+cy)/2; - - static constexpr unsigned fract_bits = 16; - static constexpr double c_fract(1 << fract_bits); - - cv::Point p(iround(b.pos[0] * cx * c_fract), iround(b.pos[1] * cy * c_fract)); + if (preview_visible) + draw_blobs(preview_frame_.as<Frame>()->mat, + blobs.data(), blobs.size(), + frame_gray.size()); - auto circle_color = k >= PointModel::N_POINTS - ? cv::Scalar(192, 192, 192) - : cv::Scalar(255, 255, 0); - - const f overlay_size = dpi > 1.5 ? 2 : 1; - - cv::circle(preview_frame, p, iround((b.radius + 3.3) * c_ * c_fract), circle_color, overlay_size, cv::LINE_AA, fract_bits); - - char buf[16]; - buf[sizeof(buf)-1] = '\0'; - std::snprintf(buf, sizeof(buf) - 1, "%.2fpx", b.radius); - - auto text_color = k >= PointModel::N_POINTS - ? cv::Scalar(160, 160, 160) - : cv::Scalar(0, 0, 255); - - cv::Point pos(iround(b.pos[0]*cx+offx), iround(b.pos[1]*cy+offy)); - cv::putText(preview_frame, buf, pos, - cv::FONT_HERSHEY_PLAIN, overlay_size, text_color, - 1); - } // End of mean shift code. At this point, blob positions are updated with hopefully less noisy less biased values. points.reserve(max_blobs); @@ -375,3 +454,5 @@ blob::blob(f radius, const vec2& pos, f brightness, const cv::Rect& rect) : { //qDebug() << "radius" << radius << "pos" << pos[0] << pos[1]; } + +} // ns pt_module |