2. 程式人生 > >GMM混合高斯背景建模C++結合Opencv實現（內附Matlab實現）

GMM混合高斯背景建模C++結合Opencv實現（內附Matlab實現）

阿新 • • 發佈：2018-11-11

最近在做視訊流檢測方面的工作，一般情況下這種視訊流檢測是較複雜的場景，比如交通監控，或者各種監控攝像頭，場景比較複雜，因此需要構建背景影象，然後去檢測場景中每一幀動態變化的前景部分，GMM高斯模型是建模的一種方法，關於高斯建模的介紹有很多部落格了，大家可以去找一找，本篇部落格主要依賴於上一個老兄，他用matlab實現了GMM模型，我在其基礎上利用C++和OpenCV進行了重寫，下面會給出C++程式碼，希望能給大家一點幫助，本文能力有限，如有問題可以一起交流，一起改進。

先展示結果吧，不要問我的影象是神馬影象，哈哈哈 感覺很low b。 其實都是視訊中的某一幀影象。方法可以用就行啦！！

背景影象：
在這裡插入圖片描述
待檢測影象：
在這裡插入圖片描述
檢測結果：
在這裡插入圖片描述
我其實是想檢測出圖片中的缺口，以下是檢測結果，我改變一下形式吧 ，這樣能更清楚的看出結果
在這裡插入圖片描述
其實還是很精準的嘛，來來來，下面上程式碼~

C++結合OpenCV程式碼


//Writen by 蘇丶 2018//11/04
//轉載請附轉載連結
#include<opencv2/opencv.hpp>
#include<opencv2/core.hpp>
#include<math.h>
#include<iostream>
#include<stdio.h>

#define K 3
#define
 
 initial_sd 36
#define  D 2.5
#define WIDTH 640
#define HEIGHT 32
#define thresh 0.25
#define learning_rate 0.01

using namespace cv;
using namespace std;


double minIndex(double a, double b, double c)
{
	if (a < b && a < c) return 1;
	if (b < a && b < c) return 2;
	if (c <
 
 a && c < b)return 3;
	
}

struct Idx {
	double data;
	int index;
};

bool sortRule(Idx a, Idx b)
{
	return a.data > b.data;
}



int main()
{
	//輸入背景影象
	Mat srcImage = imread("1.tif", 0);

	double *weight,*pixelValueMean,*sd,*udiff;
	int *matchCount;
	weight = (double*)malloc(sizeof(double)*HEIGHT*WIDTH*K);//建立三個高斯模型的權重
	pixelValueMean = (double*)malloc(sizeof(double)*HEIGHT*WIDTH*K);
	sd = (double*)malloc(sizeof(double)*HEIGHT*WIDTH*K);
	matchCount = (int*)malloc(sizeof(int)*HEIGHT*WIDTH*K);
	udiff = (double*)malloc(sizeof(double)*HEIGHT*WIDTH*K);

	//為權重，畫素均值，畫素標準差，匹配次數，圖片與高斯均值的差  進行賦初始值
	for (int i = 0; i < HEIGHT; i++)
	{
		uchar* pixelValuePtr = (uchar*)srcImage.data + i * WIDTH;
		for (int j = 0; j < WIDTH; j++)
		{
			for (int k = 0; k < K; k++)
			{
				if (k == 0)
				{
					weight[i*WIDTH*K + j * K + k] = 1;
					pixelValueMean[i*WIDTH*K + j * K + k] = pixelValuePtr[j];
					sd[i*WIDTH*K + j * K + k] = initial_sd;
					matchCount[i*WIDTH*K + j * K + k] = 1;
					udiff[i*WIDTH*K + j * K + k] = 0;
				}
				else
				{
					matchCount[i*WIDTH*K + j * K + k] = 1;
					pixelValueMean[i*WIDTH*K + j * K + k] =0;
					weight[i*WIDTH*K + j * K + k] =pow(2.0,-10);
					sd[i*WIDTH*K + j * K + k] = initial_sd;
					udiff[i*WIDTH*K + j * K + k] = 0;
				}
			}
		}
	}


	Mat dstImage = imread("2.tif",0);//輸入待檢測影象

	//計算待檢測影象與各個高斯均值的差
	for (int i = 0; i < HEIGHT; i++)
	{
		uchar* pixelValuePtr1 = (uchar*)dstImage.data + i * WIDTH;
		for (int j = 0; j < WIDTH; j++)
		{
			for (int m = 0; m < K; m++)
			{
				//計算待檢測影象與每一層高斯均值圖的差值
				udiff[i*WIDTH*K + j * K + m] = abs(pixelValuePtr1[j] - pixelValueMean[i*WIDTH*K + j * K + m]);
			}
		}
	}

	//更新每一個畫素的背景模型
	double p;
	int match,match_ind;
	int kValue1, kValue2, kValue3;
	int min_w_index;
	double sum_weight;

	//建立前景影象
	Mat foreGroundImage = Mat::zeros(dstImage.size(), dstImage.type());
	Mat backGroundImage= Mat::zeros(dstImage.size(), dstImage.type());

	double rank1, rank2, rank3;
	vector<Idx>Rank_index;
	Idx d1, d2, d3;

	for (int i = 0; i < HEIGHT; i++)
	{
		uchar* pixelValuePtr2 = (uchar*)dstImage.data + i * WIDTH;
		uchar* pixelValuePtr3 = (uchar*)foreGroundImage.data + i * WIDTH;
		uchar* pixelValuePtr4 = (uchar*)backGroundImage.data + i * WIDTH;
		for (int j = 0; j < WIDTH; j++)
		{
			match = 0;//畫素值與高斯模型匹配的標識
			match_ind = 0;//為該畫素最匹配的高斯模型的標號

			for (int k = 0; k < K; k++)
			{
				if (abs(udiff[i*WIDTH*K + j * K + k]) <= D * sd[i*WIDTH*K + j * K + k] && match == 0)
				{
					match = 1;
					match_ind = k;
					p = learning_rate / weight[i * WIDTH * K + j * K + k];
					weight[i * WIDTH * K + j * K + k] = (1 - learning_rate) * weight[i * WIDTH * K + j * K + k] + learning_rate;
					//更新用的都是待檢測影象的畫素值
					pixelValueMean[i * WIDTH *K + j * K + k] = (1 - p) * pixelValueMean[i * WIDTH * K + j * K + k] + p * pixelValuePtr2[j];
					sd[i * WIDTH * K + j * K + k] = sqrt((1 - p) * pow(sd[i * WIDTH * K + j * K + k], 2) + p * pow((pixelValuePtr2[j] - pixelValueMean[i*WIDTH*K + j * K + k]), 2));
					if (matchCount[i*WIDTH*K + j * K + k] != 255)
					{
						matchCount[i*WIDTH*K + j * K + k]++;
					}
				}
				else
				{
					weight[i*WIDTH*K + j * K + k] = (1 - learning_rate)*weight[i*WIDTH*K + j * K + k];
				}
			}
			if (match == 0)
			{
				kValue1 = weight[i * WIDTH * K + j * K];
				kValue2 = weight[i*WIDTH*K + j * K + 1];
				kValue3 = weight[i*WIDTH*K + j * K + 2];
				
				min_w_index = minIndex(kValue1, kValue2, kValue3);//找到權重最小值

				matchCount[i * WIDTH * K + j * K + min_w_index] = 1;
				weight[i*WIDTH*K + j * K + min_w_index] = 1 / (matchCount[i*WIDTH*K + j * K] + matchCount[i*WIDTH*K + j * K + 1] + matchCount[i*WIDTH*K + j * K + 2] - 1);
				pixelValueMean[i*WIDTH*K + j * K + min_w_index] = pixelValuePtr2[j];
				sd[i*WIDTH*K + j * K + min_w_index] = initial_sd;
			}

			sum_weight = weight[i*WIDTH*K + j * K] + weight[i*WIDTH*K + j * K + 1] + weight[i*WIDTH*K + j * K + 2];

			for (int m = 0; m < 3; m++)
			{
				weight[i*WIDTH*K + j * K + m] /= sum_weight;
			}

			rank1= weight[i*WIDTH*K + j * K];
			rank2= weight[i*WIDTH*K + j * K + 1];
			rank3= weight[i*WIDTH*K + j * K + 2];
			
			d1.data = rank1;
			d1.index = 0;
			d2.data = rank2;
			d2.index = 1;
			d3.data = rank3;
			d3.index = 2;
			Rank_index.push_back(d1);
			Rank_index.push_back(d2);
			Rank_index.push_back(d3);
			//然後對rank進行降序排序，主要取其索引值
			sort(Rank_index.begin(), Rank_index.end(), sortRule);
		    
			
			//將前景的初始值設定為255
			pixelValuePtr3[j] = pixelValuePtr2[j];

			if (match == 1)
			{
				if (match_ind == Rank_index[0].index)
				{
					pixelValuePtr3[j] = 0;
				}
				else if (match_ind == Rank_index[1].index)
				{
					if (weight[i*WIDTH*K + j * K + Rank_index[0].index] < thresh)
					{
						pixelValuePtr3[j] = 0;
					}
				}
				else if(match_ind == Rank_index[3].index)
				{
					if (weight[i*WIDTH*K + j * K + Rank_index[3].index] > 1 - thresh)
					{
						pixelValuePtr3[j] = 0;
					}
				}
			}

			for (int kk = 0; kk < K; kk++)
			{
			    pixelValuePtr4[j] = pixelValuePtr4[j] + pixelValueMean[i*WIDTH*K + j * K + kk] * weight[i*WIDTH*K + j * K + kk];
			}
		}
	}

	imshow("back", backGroundImage);
	imshow("fore", foreGroundImage);
	imshow("original", srcImage);
	imshow("daijiance", dstImage);
	imwrite("fore.tif", foreGroundImage);

	waitKey();
	return 0;
}

Matlab程式碼：
我主要是針對以為老哥的Matlab程式碼改出來的C++版本，所以把老哥的部落格連結附在下面，然後把程式碼也展示一下吧
Matlab GMM實現

麻雞，再次吐槽這個配色，我不知道怎麼改啊 ，有沒有人會改 請給我留言，我看著這個配色好不爽

clear;
fr=imread('first.jpg');
% 讀取該影象作為背景 
fr_bw1 = rgb2gray(fr);   
% 將背景轉換為灰度影象 
fr_size = size(fr);             
%取幀大小 
width = fr_size(2); 
height = fr_size(1); %獲取原影象的尺寸
fg = zeros(height, width);   %前景,讀取的第二張圖片獲得
bg_bw = zeros(height, width);%背景,讀取的第一張圖片獲得
fr_bw1 = double(fr_bw1);
% --------------------- mog variables ----------------------------------- 
C = 3;           % 組成混合高斯的單高斯數目 (一般3-5) 
D = 2.5;         % 閾值（一般2.5個標準差） 
alpha = 0.01;    % learning rate 學習率決定更新速度(between 0 and 1) (from paper 0.01) 
thresh = 0.25;   % foreground threshold 前景閾值(0.25 or 0.75 in paper) 
sd_init = 36;    % initial standard deviation 初始化標準差(for new components) var = 36 in paper 
w = zeros(height,width,C);          % initialize weights array 初始化權值陣列 
w(:,:,1) = 1;
w(:,:,2:C) = 2^-10;                  % 第一個高斯分佈的初始權重為1，其餘分佈的權重為0
mean = zeros(height,width,C);        % pixel means 畫素均值 
mean(:,:,1) = fr_bw1;                % 第一個高斯分佈的初始均值為參考幀的值，其餘分佈的均值為0s
sd = sd_init*ones(height,width,C);   % pixel standard deviations 畫素標準差 
matchcnt = ones(height, width,C);    % 匹配的次數，初始值都設為1
u_diff = zeros(height,width,C);      % difference of each pixel from mean 圖片與高斯均值的差 
 
fr = imread('second.jpg');       % read in frame 讀取幀 
fr_bw = rgb2gray(fr);  	   % convert frame to grayscale 轉換為灰度影象
fr_bw = double(fr_bw);     % 將灰度圖值設定為雙精度
 %求匯入進來的圖片與各個高斯均值的差
for m=1:C  
	u_diff(:,:,m) = abs(fr_bw - double(mean(:,:,m)));     
end       
% update gaussian components for each pixel 更新每個畫素的背景模型 
%rank_ind = zeros(C,1);
for i=1:height
	for j=1:width              
        match = 0; %畫素與高斯模型匹配的標識
        match_ind = 0;%為該畫素最匹配的高斯模型的標號
        for k=1:C %與第k個高斯模型進行比對，然後更新引數                 
            if (abs(u_diff(i,j,k)) <= D*sd(i,j,k) && (match == 0))       % pixel matches component畫素匹配了高斯中的第k個模型 			
                match = 1;                               
                % variable to signal component match 設定匹配標記 				
                match_ind = k;
                % update weights, mean, sd, p  更新權值，均值，標準差和引數學習率                     
  				p = alpha/w(i,j,k);           %理應使用p = alpha/gaussian才對，這裡勉強      
                w(i,j,k) = (1-alpha)*w(i,j,k) + alpha;                     
				%p = alpha/w(i,j,k);           %理應使用p = alpha/gaussian才對，這裡勉強       
                mean(i,j,k) = (1-p)*mean(i,j,k) + p*double(fr_bw(i,j)); 
                sd(i,j,k) =   sqrt((1-p)*(sd(i,j,k)^2) + p*((double(fr_bw(i,j)) - mean(i,j,k)))^2);
                if matchcnt(i, j, k) ~= 255                            % 匹配次數達到255就不加了。在實時視訊序列中，上限是必須的，否則可能溢位。
                    matchcnt(i, j, k) = matchcnt(i, j, k) + 1;
                end 
            else                                    
			%當與第k個模型沒有匹配的話，則第k個模型所佔的比重自然而然地下降
                w(i,j,k) = (1-alpha)*w(i,j,k);      % weight slighly decreases 權值減小                                      
            end
        end
   		% if no components match, create new component 如果沒有匹配的模型則建立新模型
		if(match==0)  % 沒有匹配的高斯，建立新的高斯取代：排序後排在最後面的那個
            [min_w,min_w_index]=min(w(i,j,:));
            matchcnt(i,j,min_w_index) = 1; % 匹配次數設為1，一個小值
            w(i,j,min_w_index) = 1 / ( sum( matchcnt(i, j, :) ) - 1 );% 權值為其它高斯分佈匹配次數之和的倒數
			mean(i,j,min_w_index)=double(fr_bw(i,j));
			sd(i,j,min_w_index)=sd_init;
        end
        %無論匹配是否成功，都要將該畫素在不同模型上的權重標準歸一化
        w_sum = sum(w(i, j, :));
        w(i, j, :) = w(i, j, :) / w_sum;
        
        %針對該畫素，計算多個模型的優先順序（依據權重）
		rank = w(i,j,:)./sd(i,j,:);             
		[sorted_rank, rank_ind] = sort(rank, 'descend');
        
        %將前景的初始值設定為255，即為白色；
        fg(i,j) = fr_bw(i,j);
        %當該畫素匹配成功的時候，利用高斯混合模型，將該畫素值重新設定
        if(match == 1)
            switch match_ind
                case rank_ind(1)% 與最優的高斯匹配，肯定是歸為背景點
                    fg(i,j) = 0;
                case rank_ind(2)% 與中間的高斯匹配，如果最上面一個高斯的權值小於thresh，則這點歸為背景點
                    if w(i, j, rank_ind(1)) < thresh
                        fg(i, j) = 0;
                    end
                case rank_ind(3)% 與最下面的高斯匹配，如果最下面的高斯權值大於1-thresh（或者前兩個高斯權值和小於thresh），則這點歸為背景點
                    if w(i, j, rank_ind(3)) > 1 - thresh
                        fg(i, j) = 0;
                    end
            end
        end
        for k=1:C
            bg_bw(i,j) = bg_bw(i,j)+ mean(i,j,k)*w(i,j,k);%更新背景
        end
%         % 根據rank的排序結果調整引數的順序
%         tmp_T = [mean(i, j, :); sd(i, j, :); w(i, j, :); matchcnt(i, j, :)]; % 為了排序時，幾個引數同步調整，所以組合在一起
%         mean(i, j, :) = tmp_T(1, rank_ind);    %即同時利用rank，將mean進行了排序
%         sd(i, j, :) = tmp_T(2, rank_ind);      %同理
%         w(i, j, :) = tmp_T(3, rank_ind);       %同理
%         matchcnt(i, j, :) = tmp_T(4, rank_ind);%同理
%         if w(i, j, 1) > thresh                 %使用大於閾值的進行背景構造即可
%             bg_bw(i, j) = w(i, j, 1) * mean(i, j, 1); 
%         else
%             if w(i, j, 1) + w(i, j, 2) > thresh
%                 bg_bw(i, j) = w(i, j, 1) * mean(i, j, 1) + w(i, j, 2) * mean(i, j, 2);
%             else
%                 bg_bw(i, j) = sum(w(i, j, :) .* mean(i, j, :));
%             end
%         end
	end
end
figure(1),subplot(3,1,1),imshow(fr);    %顯示輸入影象
subplot(3,1,2),imshow(uint8(bg_bw));  %顯示背景影象
subplot(3,1,3),imshow(uint8(fg));     %顯示前景影象

轉載請附本部落格連結地址：
https://blog.csdn.net/weixin_38285131/article/details/83721069

相關推薦

GMM混合高斯背景建模C++結合Opencv實現（內附Matlab實現）

最近在做視訊流檢測方面的工作，一般情況下這種視訊流檢測是較複雜的場景，比如交通監控，或者各種監控攝像頭，場景比較複雜，因此需要構建背景影象，然後去檢測場景中每一幀動態變化的前景部分，GMM高斯模型是建模的一種方法，關於高斯建模的介紹有很多部落格了，大家可以去找一找，本篇部落格主要依賴於上

混合高斯背景建模原理及實現(C# )

原文http://blog.csdn.net/jinshengtao/article/details/26278725 前些日子一直在忙答辯的事情，畢業後去了華為，影象處理什麼的都派不上用場了。打算分3-4篇文章，把我研究生階段學過的常用演算法為大家和4107的師弟師妹

背景建模演算法--平均背景建模、單高斯背景建模C實現

1、背景建模   視訊影象中運動目標的檢測的一種基本思想是對影象的背景進行建模，既是將所有畫素分為背景和運動前景兩類。 背景建模方法可以分為兩類 – 顏色背景建模和紋理背景建模。此文只涉及顏色背景建模。紋理背景建模可參考 顏色背景建模的基本原理：對影象中每個畫素的顏色值（灰度或彩色）進

【影象處理基礎知識】-混合高斯背景建模

0、高斯模糊原理-------------------------------------------------------分割線----------------------------------------------------------------1、一維高斯函式

GMM混合高斯模型演算法詳解

使用概率模型的原因 k均值等價於假設了球對稱形狀的聚類。使用帶權歐式距離，仍然假設了軸對齊的橢球。沒有考慮聚類的形狀。 促使概率模型的原因：混合模型 提供觀測點到聚類的軟分配soft assignment（分配包含不確定性） 考慮了聚類的形狀而不僅僅是中心

PCA,LDA的C++實現（內附程式碼~）

原始碼獲取通道：C++實現LDA， 實現PCA pca_z執行結果： LDA_Z執行結果： 原理簡單易懂，用於降維。 Visual C++ 無法解析的外部符號問題 在VS環境下用Eigen進行SVD分解 先輸入陣列，二維陣列，然後轉換為M

綜合實踐報告（7）C程式原始碼操作（內附五套程式碼）

#include <iostream> #include <cstring> #include <cstdio> using namespace std; char str[9999]; // str 刪除多餘空格後的 int cint[25],cmain[25

C++實現控制檯版2048（內附彩蛋）

前言　　之前做過一個JavaScript版本的2048遊戲，最近在學習C++，昨天晚上突然心血來潮，想用C++來實現，因為核心演算法已十分理解，所以兩個小時擼出來一個C++的簡易版本。簡介　　二維陣列遍歷，C++基本資料型別，string類，控制結構，函式。　　本方法不涉及指

高斯混合模型視訊背景建模的EM演算法與Matlab 實現

1.問題描述 影像的背景前景分離. 輸⼊為影像監控的1000 幀 (如下⽅圖中左邊所⽰), 要求輸出是背景和前景 (如下⽅圖中右邊所⽰). 2.背景知識 觀察待處理的監控影像，可以發現，前景主要是來來往往的行人，背景始終是攝像頭對準的固定區域，

機器學習筆記（十）EM演算法及實踐（以混合高斯模型（GMM）為例來次完整的EM）

今天要來討論的是EM演算法。第一眼看到EM我就想到了我大楓哥，EM Master，千里馬，RUA!!!不知道看這個部落格的人有沒有懂這個梗的。好的，言歸正傳，今天要講的EM演算法，全稱是Expectation maximization，期望最大化。怎麼個意思呢，就是給你一

sklearn聚類模型：基於密度的DBSCAN；基於混合高斯模型的GMM

type 依次 imp style wid -m ict ria htm 1 sklearn聚類方法詳解 2 對比不同聚類算法在不同數據集上的表現 3 用scikit-learn學習K-Means聚類 4 用scikit-learn學習DB

EM算法與混合高斯模型

exp n) white http tlab real 概率 mat strong EM算法與混合高斯模型 close all; clear; clc; %% Sample Generate N=5000; a_real =[3/10,5/10,2/10]; mu_re

數值分析中的高斯消元 c語言實現附帶註釋

#include <iostream> #include <cstdio> #include <cmath> using namespace std; double

機器學習之混合高斯模型(Gaussian Mixture Model)聚類演算法+程式碼

機器學習之混合高斯模型聚類演算法1 演算法原理2 演算法例項3 典型應用參考資料 機器學習分為監督學習、無監督學習和半監督學習(強化學習)。無監督學習最常應用的場景是聚類(clustering)和降維(dimension reduction)。聚類演算法包括：

SVM的高斯kernel 引數c和gamma說明

The SVM with a Gaussian kernel function has two such training parameters: C which controls overfitting of the model, and gamma (γ) which

Matlab產生二維混合高斯分佈及隨機數

參考：https://ww2.mathworks.cn/help/stats/gmdistribution.random.html 修改如下： Create a gmdistribution object and generate random variates. De

高斯模糊函式 c 程式碼

注：程式碼來源於 http://hi.baidu.com/iceboy_/blog/item/729b79cae4744a18be09e6e9.html ，本人只是稍走修改，測試過，能使用。模糊後的資料保存於 傳入的data中。  int gaussBlur:(int *d

資料學習(9)·最大期望演算法·混合高斯模型(上)

作者課堂筆記摘錄，有問題請聯絡 [email protected] Preview 高斯混合模型（Mixture of Gaussians） 最大期望演算法（EM）

安卓 高斯背景模糊popupwindow，彈簧彈出 rebounds

最終類似效果圖，沒有截動畫，gridelayout彈上來的時候是波紋的。 直接上程式碼。 public class MyBlurPopWin extends BlurPopupWindow { Context mContext; GridLayout

robots_estimation and learning之混合高斯模型

    近來slam技術在無人汽車領域比較火熱，所以現在打算把cousera上的robots_estimation and learning的4周課程刷一遍。做作業做了，但是沒有提交(之前買的課到期了)，所以也就沒有驗證正確性。 一、課程背景介紹 我們知道SLAM問題