前言

时间过得非常快已经快要11月了，有点混所以进度算有点小落后。这次主要介绍meanShift原理从直方图循序了解，藉由本篇文章介绍主要能达到对于原理的理解并解实作。

画直线与方块

在介绍直方图之前，首先要先建立一个标头档专门在绘製直线与方块。而直线主要是有点到点所以先建立Point，矩形则是参考OpenCV方式给予点和长宽所以建立Rect。

新增Point.h

#pragma once#ifndef POINT_H#define POINT_H#include "general.h"#include "Image.h"class Point {public:Point(C_UINT32 x, C_UINT32 y) :_x(x), _y(y) {};void X(C_UINT32 x);UINT32 X() const;void Y(C_UINT32 y);UINT32 Y() const;private:UINT32 _x;UINT32 _y;};#endif // !POINT_H

新增Point.cpp

#include "Point.h"void Point::X(C_UINT32 x) {_x = x;}UINT32 Point::X() const {return _x;}void Point::Y(C_UINT32 y) {_y = y;}UINT32 Point::Y() const {return _y;}

新增Rect.h

EndX：取得实际终点X大小。EndY：取得实际终点Y大小。

#pragma once#ifndef RECT_H#define RECT_H#include "general.h"#include "Image.h"class Rect {public:Rect(C_UINT32 x, C_UINT32 y, C_UINT32 width, C_UINT32 height) :_x(x), _y(y), _width(width), _height(height) {  };void X(C_UINT32 x);UINT32 X() const;void Y(C_UINT32 y);UINT32 Y() const;void Width(C_UINT32 width);UINT32 Width() const;void Height(C_UINT32 height);UINT32 Height() const;UINT32 EndX() const;UINT32 EndY() const;private:UINT32 _x;UINT32 _y;UINT32 _width;UINT32 _height;};#endif // !RECT_H

新增Rect.cpp

#include "Rect.h"void Rect::X(C_UINT32 x){_x = x;}UINT32 Rect::X() const{return _x;}void Rect::Y(C_UINT32 y){_y = y;}UINT32 Rect::Y() const {return _y;}void Rect::Width(C_UINT32 width){_width = width;}UINT32 Rect::Width() const{return _width;}void Rect::Height(C_UINT32 height){_height = height;}UINT32 Rect::Height() const{return _height;}UINT32 Rect::EndX() const{return _x + _width;}UINT32 Rect::EndY() const{return _y + _height;}

新增draw.h

#pragma once#ifndef DRAW_H#define DRAW_H#include "general.h"#include "Image.h"#include "Point.h"#include "Rect.h"namespace MNDT {inline void DrawLine8bit(const Image& image, const Point& p1, const Point& p2, C_UCHAE& color){assert(image.Width() > p1.X() && image.Width() > p2.X());assert(image.Height() > p1.Y() && image.Height() > p2.Y());C_UINT32 diffX = abs(static_cast<int32_t>(p1.X()) - static_cast<int32_t>(p2.X()));C_UINT32 diffY = abs(static_cast<int32_t>(p1.Y()) - static_cast<int32_t>(p2.Y()));C_UINT32 base = diffX > diffY ? diffX : diffY;C_FLOAT baseX = static_cast<float>(diffX) / static_cast<float>(base);C_FLOAT baseY = static_cast<float>(diffY) / static_cast<float>(base);float x = static_cast<float>(p1.X());float y = static_cast<float>(p1.Y());for (UINT32 index = 0; index < base; index++){image.image[static_cast<UINT32>(y)][static_cast<UINT32>(x)] = color;x += baseX;y += baseY;}}inline void DrawLine8bit(const Image& image, const Point& p1, const Point& p2){DrawLine8bit(image, p1, p2, 255);}inline void DrawRect8bit(const Image& image, const Rect& rect, C_UCHAE& color){C_UINT32 endX = rect.EndX();C_UINT32 endY = rect.EndY();for (UINT32 x = rect.X(); x < endX; x++){image.image[rect.Y()][x] = color;image.image[endY][x] = color;}for (UINT32 y = rect.Y(); y < endY; y++){image.image[y][rect.X()] = color;image.image[y][endX] = color;}}inline void DrawRect8bit(const Image& image, const Rect& rect){DrawRect8bit(image, rect, 255);}}#endif // !DRAW_H

八位元直方图

在八位元影像当中，是由0到255像素所组成的，而直方图顾名思义就是计算0到255的数量。

步骤
1.走访图片每个像素并记录。

Library.h加入

SetHistogram8bit：取得直方图资料，有一重载因后面会使用到已区间计算直方图。

/*SetHistogram8bit Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's heightminRange= min pixelmaxRange= max pixelbin= histogran split size*/void SetHistogram8bit(C_UCHAE* src, int32_t* histogram, C_UINT32 width, C_UINT32 height);void SetHistogram8bit(C_UCHAE* src, int32_t* histogram, C_UINT32 width, C_UINT32 height, C_UCHAE minRange, C_UCHAE maxRange, C_UCHAE bin);

Libary.cpp

void Library::SetHistogram8bit(C_UCHAE* src, int32_t* histogram, C_UINT32 width, C_UINT32 height){C_UINT32 size = width * height;for (UINT32 count = 0; count < size; count++){histogram[*src]++;src++;}}void Library::SetHistogram8bit(C_UCHAE* src, int32_t* histogram, C_UINT32 width, C_UINT32 height, C_UCHAE minRange, C_UCHAE maxRange, C_UCHAE bin){assert(maxRange > minRange);C_UCHAE diffRange = maxRange - minRange;C_UCHAE interval = diffRange / bin;C_UINT32 size = width * height;for (UINT32 count = 0; count < size; count++){histogram[*(src + count) / interval]++;}}

规一化

[1]使用OpenCV的规一化有多种方法，这里使用最大值规一化。

步骤
1.取得size大小中最大值。

2.将size个都规一化。

Library.h加入

SetNormalizedHistogram8bit：规一化直方图数据。

/*SetNormalizedHistogram8bit Parameter:histogram= histogram datasize= histogram sizebase= max pixel*/void SetNormalizedHistogram8bit(int32_t* histogram, C_UINT32 size, C_UCHAE base);

Library.cpp加入

void Library::SetNormalizedHistogram8bit(int32_t* histogram, C_UINT32 size, C_UCHAE base){int32_t max = 0;for (UINT32 index = 0; index < size; index++){max = max < histogram[index] ? histogram[index] : max;}float normalizedBase = static_cast<float>(base) / max;for (UINT32 index = 0; index < size; index++){histogram[index] = static_cast<int32_t>(histogram[index] * normalizedBase);}}

八位元直方图可视化

使用上述取得值方图资料后再规一化，并画在256 * 256大小图片上。

步骤
1.计算像素直方图数量。

2.规一化直方图数据。

3.画出0到255直方图。

结果图。

Library.h加入

Histogram8bit：取得直方图图像。

/*Histogram8bit Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's height*/void Histogram8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height);

Library.cpp加入

void Library::Histogram8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height){int32_t* histogram = new int32_t[256]{ 0 };int32_t max = 0;SetHistogram8bit(src, histogram, width, height);SetNormalizedHistogram8bit(histogram, 256, 255);Image purImage(pur, 256, 256, MNDT::GRAY_8BIT);for (UINT32 index = 0; index < 256; index++){Point start(index, (255 - histogram[index]));Point end(index, 255);MNDT::DrawLine8bit(purImage, start, end);}delete[] histogram;}

直方图等化

直方图等化主要是使用累积分布函数来去改变图片性质，图一为直方图数量，图二为依序累积后的数量，然而最后依照原数量除上累积总数量则是累积分布的结果。
结果测试图片来源[2]，公式相关参考之前上课老师所提供的讲义。

图一来源：[2]。

图二来源：[2]。

结果图。
步骤
1.计算像素直方图数量。

2.计算每个累积分布机率在乘上255做规一化。

3.计算后的直方图带入原始图像。

Library.h加入

/*HistogramEqualization8bit Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's height*/void HistogramEqualization8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height);

Library.h加入

void Library::HistogramEqualization8bit(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height){C_UINT32 size = width * height;int32_t* histogram = new int32_t[256]{ 0 };SetHistogram8bit(src, histogram, width, height);float base = 0;for (UINT32 index = 0; index < 256; index++){if (histogram[index] > 0){base += static_cast<float>(histogram[index]) / static_cast<float>(size);histogram[index] = static_cast<int32_t>(base * 255);}}for (UINT32 index = 0; index < size; index++){*(pur + index) = histogram[*(src + index)];}delete[] histogram;}

直方图比较

[1]有介绍OpenCV有多种比较方法，这边使用直方图相减去计算误差，设定为越接近1代表误差越小。
步骤
1.取得两张图片的直方图。

2.两个直方图相减计算累加。

3.一扣掉累积除上像素数量。

Library.h加入

/*CompareHistogram Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's height*/double CompareHistogram(C_UCHAE* src, C_UCHAE* pur, C_UINT32 width, C_UINT32 height);

Library.cpp加入

double Library::CompareHistogram(C_UCHAE* src, C_UCHAE* pur, C_UINT32 width, C_UINT32 height){int32_t* srcHistogram = new int32_t[256]{ 0 };int32_t* purHistogram = new int32_t[256]{ 0 };SetHistogram8bit(src, srcHistogram, width, height);SetHistogram8bit(pur, purHistogram, width, height);int32_t diff = 0;for (UINT32 index = 0; index < 256; index++){diff += abs(srcHistogram[index] - purHistogram[index]);}return 1.0 - static_cast<double>(diff) / static_cast<double>(width * height);}

影像撷取通道

这边介绍取出单一通道影像，因后面meanShift会使用到。主要用途为将BGR或其他色彩通道指定取得单一通道，例如BGR，通道1为只取出B色彩图像，通道2为只取出G色彩图像，通道3为只取出R色彩图像。

取出B通道结果图。

Library.h加入

/*Channel Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's heightchannel= color channel*/void Channel(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32 channel);

Library.cpp加入

void Library::Channel(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32 channel){assert(src != nullptr && pur != nullptr);assert(width > 0 && height > 0);assert(channel > 0 && channel < 4);C_UCHAE* purEnd = pur + width * height + 1;src += (channel - 1);while (pur < purEnd){*pur = *src;pur++;src += 3;}}

直方图反投影

主要将输入图片依照输入的分布直方图去做输出，简单来说先计算一张图片的直方图资料，再将另一张图片像素使用计算好的直方图资料做投影，也就是套上直方图资料在输入图像的机率分布，而这边使用HSV的H，因为H代表色度比较能代表图片的特徵。

结果图。

取得直方图步骤-使用C#
1.将一张图片转为HSV。

2.使用Channel取得H。

3.计算直方图资料。

4.直方图资料规一化。

步骤
1.将每个像素转为bin区间，并转为输入直方图的值。

C#取得直方图资料
hisData与之前输入图片原理相同，hisPtr为H通道图像的指标。

IntPtr hisPtr = hisData.Scan0;int[] histogram = new int[bin];mndtSetHistogram8bit(hisPtr, histogram                   , hisImage.Width, hisImage.Height                   , 0, bin                   , bin);mndtSetNormalizedHistogram8bit(histogram, bin, bin);

Library.h加入

/*BackProjection Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's heighthistogram= histogram dataminRange= min pixelmaxRange= max pixelbin= histogran split size*/void BackProjection(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32* histogram, C_UCHAE minRange, C_UCHAE maxRange, C_UCHAE bin);

Library.cpp加入

void Library::BackProjection(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32* histogram, C_UCHAE minRange, C_UCHAE maxRange, C_UCHAE bin){assert(src != nullptr && pur != nullptr);assert(width > 0 && height > 0);C_UCHAE* purEnd = pur + width * height + 1;C_UCHAE diffRange = maxRange - minRange;C_UCHAE interval = diffRange / bin;while (pur < purEnd){*pur = histogram[(*src) / interval];pur++;src++;}}

均值漂移

这边均值漂移作法叫简单，如下图一，只需要已中心计算範围内的值即可求出移动量和方向，而图一的红点在图片上即是像素，像素也则代表对移动和方向的影响量。上面使用反向投影求出分布图后，在使用均值漂移则可以求出新图片的目标位置。
注：实际上均值漂移也是一种机器学习，通常会使用高斯核来去计算分群目标[4]。

图一来源：[4]。

结果图，蓝色为原始位置，红色为搜寻到的位置。

公式简介
参考[3]和OpenCV的meanShift原始码。假设矩形大小为120　*　60，则中心位置为(60, 30)。矩形内每个座标已中心为原点(二维座标系)去累加计算x和y乘上像素(影像度)，最后在除上总像素计算x和y的位移量。

原点(60, 30)。
二维座标为，G(x) = x - 60，G(y) = y - 30。

步骤
1.计算矩形内的x像素、y像素和像素总和。

2.取得x和y位移量并且修正数值。

3.偏移达到threshold或time则结束。

Library.h加入

/*MeanShift Parameter:src= source of imagepur= purpose of imagewidth= Image's widthheight= Image's heightrectPoint= original rect pointtime= trans timethreshold= diff threshold*/void MeanShift(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32* rectPoint, C_UINT32 times, C_DOUBLE threshold);

Library.cpp加入

void Library::MeanShift(C_UCHAE* src, UCHAE* pur, C_UINT32 width, C_UINT32 height, C_UINT32* rectPoint, C_UINT32 times, C_DOUBLE threshold){assert(src != nullptr && pur != nullptr);assert(width > 0 && height > 0);memcpy(pur, src, width * height * sizeof(UCHAE));Image srcImage(const_cast<UCHAE*>(src), width, height, ColerType::BGR2GRAY_8BIT);Image purImage(const_cast<UCHAE*>(pur), width, height, ColerType::BGR2GRAY_8BIT);Rect rect(rectPoint[0], rectPoint[1], rectPoint[2], rectPoint[3]);MNDT::DrawRect8bit(purImage, rect);for (UINT32 time = 0; time < times; time++){C_UINT32 endX = rect.EndX();C_UINT32 endY = rect.EndY();int32_t centerX = (rect.Width() >> 1) + rect.X();int32_t centerY = (rect.Height() >> 1) + rect.Y();int32_t sumX = 0;int32_t sumY = 0;int32_t sumBase = 0;for (UINT32 y = rect.Y(); y < endY; y++){for (UINT32 x = rect.X(); x < endX; x++){if (srcImage.image[y][x] > 0){C_INT32 xDiff = x - centerX;C_INT32 yDiff = y - centerY;sumX += (xDiff * srcImage.image[y][x]);sumY += (yDiff * srcImage.image[y][x]);sumBase += srcImage.image[y][x];}}}C_DOUBLE offsetX = (static_cast<double>(sumX) / static_cast<double>(sumBase));C_DOUBLE offsetY = (static_cast<double>(sumY) / static_cast<double>(sumBase));int32_t updateX = rect.X() + static_cast<int32_t>(offsetX);int32_t updateY = rect.Y() + static_cast<int32_t>(offsetY);updateX = std::min(std::max(updateX, 0), static_cast<int32_t>(width - 1));updateY = std::min(std::max(updateY, 0), static_cast<int32_t>(height - 1));rect.X(updateX);rect.Y(updateY);// lossif ((offsetX * offsetX + offsetY * offsetY) <= threshold) break;}MNDT::DrawRect8bit(purImage, rect);}

结论

C#视窗原始码
C++函数原始码
感觉时间过得很快，今年又快过完了，但还有很多东西都还没学到一个段落，看来要努力加把劲了，希望大家今年也能完成自己的目标準备朝向新的一年前进。文章若有文题欢迎留言提问谢谢。

参考文献

[1]阿洲(2015). OpenCV教学 | 阿洲的程式教学 from: http://monkeycoding.com/?page_id=12 (2018.10.19).
[2]维基百科(2018). 直方图均衡化 from:https://zh.wikipedia.org/wiki/%E7%9B%B4%E6%96%B9%E5%9B%BE%E5%9D%87%E8%A1%A1%E5%8C%96 (2018.10.10).
[3]CamShift算法--Mean Shift算法 from:http://www.voidcn.com/article/p-rekabebj-zs.html
[4]简单易学的机器学习算法——Mean Shift聚类算法 from:https://blog.csdn.net/google19890102/article/details/51030884