opencv实现 Matlab 寻找峰值函数 findpeaks

findpeaks 寻找峰值函数

[pks,locs] = findpeaks(coutss,'minpeakdistance',30,'minpeakheight',mean_value-5);

OpenCV实现findpeaks(coutss,'minpeakdistance',30,'minpeakheight',mean_value-5); 

 

[cpp] view plain copy

print ?

1. /* 
2. 输入参数： 
3. Matdata: 输入的数据矩阵 
4. minpeakdistance：设定两峰值间的最小间隔数 
5. minpeakheight  ：设定峰值的最小高度 
6. 输出函数： 
7. mapPeaks_max        ：按照key  值（灰度级0-255）排序的 
8. Vector_map_Peaks_max：按照value值（各灰度级的像素值）排序的 
9. 备注：此函数只完成了波峰的提取，波谷自行完成 
10. - zhang 
11. */  
12. void findPeaks(const Mat Matdata, int minpeakdistance, int minpeakheight, std::map <int, double> &mapPeaks_max, vector<pair<int, double>>&Vector_map_Peaks_max)  
13. {  
14.     int row = Matdata.rows - 1;  
15.     int col = Matdata.cols;  
16.     vector<int> sign;  
17.     vector<double> markdata;  
18.     vector<double> markdata1;  
19.     //minMaxLoc寻找矩阵(一维数组当作向量,用Mat定义) 中最小值和最大值的位置.   
20.     //类型转换注意：data为unchar*   
21.     double* pMatdata = (double*)Matdata.data;  
22.     for (int i = 1; i < col; i++)  
23.     {  
24.         /*相邻值做差： 
25.         *小于0，赋-1 
26.         *大于0，赋1 
27.         *等于0，赋0 
28.         */  
29.         markdata1.push_back(*(pMatdata + i-1));  
30.         int diff = *(pMatdata + i) - *(pMatdata + i-1);  
31.         if (diff > 0)  
32.         {  
33.             sign.push_back(1);  
34.         }  
35.         else if (diff < 0)  
36.         {  
37.             sign.push_back(-1);  
38.         }  
39.         else  
40.         {  
41.             sign.push_back(0);  
42.         }  
43.     }  
44.     //再对sign相邻位做差    
45.     //保存极大值和极小值的位置    
46.     vector<int> indMax;  
47.     vector<int> indMin;  
48.     std::map <int, double> mapPeaks_min;  
49.     for (int j = 1; j < sign.size(); j++)  
50.     {  
51.         int diff = sign[j] - sign[j - 1];  
52.         if (diff < 0)  
53.         {  
54.             indMax.push_back(j);  
55.             //根据峰值最小高度进行筛选  
56.             if ( *(pMatdata +  indMax[indMax.size() - 1]) > minpeakheight){  
57.                 mapPeaks_max.insert(pair<int, double>(j, *(pMatdata + indMax[indMax.size() - 1])));  
58.             }  
59.               
60.         }  
61.         else if (diff>0)  
62.         {  
63.             indMin.push_back(j);  
64.             mapPeaks_min.insert(pair<int, double>(j, *(pMatdata + indMin[indMin.size() - 1])));  
65.               
66.         }  
67.     }  
68.   
69.     if (mapPeaks_max.size() >= 2){  
70.   
71.         //找到像素值最大的对应的容器指针  
72.         std::map<int, double>::iterator iter_high_max1_pos = std::max_element(std::begin(mapPeaks_max), std::end(mapPeaks_max), cmpDescending);  
73.         int high_max1_pos = iter_high_max1_pos->first;//最高点的位置  
74.         map<int, double>::iterator  iter;  
75.         map<int, double>::iterator  temp;  
76.         //以位置为序  
77.         map<int, double>::iterator  iter_high_max1_pos;  
78.         map<int, double>::iterator  iter_high_max1_pos_temp1;  
79.         map<int, double>::iterator  iter_high_max1_pos_temp2;  
80.         //知道峰值最高的对应容器指针  
81.         //我们以最高峰值处进行作业搜索所以在此声明两个局部变量保存该指针  
82.         iter_high_max1_pos_temp1= iter_high_max1_pos;  
83.         iter_high_max1_pos_temp2 = iter_high_max1_pos;  
84.          
85.         //rate为宽度，以iter_high_max1_pos为起点,根据minpeakdistance的距离倍数搜索  
86.         //minpeakdistance宽度内取一个最大的峰值max_value_temp  
87.         int rate = 1;  
88.         double max_value_temp = 0;  
89.         map<int, double>::iterator max_value_temp_iter;  
90.         //每个minpeakdistance内的第一个指针防止被删除  
91.         bool is_first_max = true;  
92.         //挨着最大峰值的  
93.         if (mapPeaks_max.size() >= 2){  
94.             if (iter_high_max1_pos_temp1 != mapPeaks_max.begin()){  
95.                   
96.                 for (iter = --iter_high_max1_pos_temp1; iter != mapPeaks_max.begin(); iter--){  
97.                     //小于minpeakdistance  
98.                     if (labs((*iter).first - high_max1_pos) < (rate)*minpeakdistance){  
99.                             temp = iter;  
100.                             iter++;  
101.                             mapPeaks_max.erase(temp);  
102.                     }  
103.                     if (rate*minpeakdistance<=labs((*iter).first - high_max1_pos) && labs((*iter).first - high_max1_pos)<=(rate+1)*minpeakdistance){  
104.                         if ((*iter).second > max_value_temp){  
105.                             max_value_temp = (*iter).second;  
106.                             if (is_first_max){  
107.                                 max_value_temp_iter = iter;  
108.                                 is_first_max = false;  
109.                             }  
110.                             else{  
111.                                 //删除原先最大的  
112.                                 mapPeaks_max.erase(max_value_temp_iter);  
113.                                 temp = iter;  
114.                                 iter++;  
115.                                 max_value_temp_iter = temp;   
116.                             }  
117.                         }  
118.                         else{  
119.                             is_first_max = true;  
120.                             max_value_temp = 0;  
121.                             temp = iter;  
122.                             iter++;  
123.                             mapPeaks_max.erase(temp);  
124.                         }  
125.                     }  
126.                     //大于（rate + 1）*minpeakdistance  
127.                     else if (labs((*iter).first - high_max1_pos)>(rate + 1)*minpeakdistance){  
128.                            
129.                             rate++;  
130.                             iter++;  
131.                     }  
132.                 }  
133.             }  
134.             is_first_max = true;  
135.             if (iter_high_max1_pos_temp2 != mapPeaks_max.end() && mapPeaks_max.size()>=2){  
136.               
137.                 for (iter = ++iter_high_max1_pos_temp2; iter != mapPeaks_max.end(); iter++){  
138.                     //TRACE("2Vector:%d \n", debug++);  
139.                     //小于minpeakdistance  
140.                     if (labs((*iter).first - high_max1_pos) < (rate)*minpeakdistance){  
141.                         temp = iter;  
142.                         iter--;  
143.                         mapPeaks_max.erase(temp);  
144.                     }  
145.                     if (rate*minpeakdistance <= labs((*iter).first - high_max1_pos) && labs((*iter).first - high_max1_pos) <= (rate + 1)*minpeakdistance){  
146.   
147.                         if ((*iter).second > max_value_temp){  
148.                             max_value_temp = (*iter).second;  
149.                             if (is_first_max){  
150.                                 max_value_temp_iter = iter;  
151.                                 is_first_max = false;  
152.                             }  
153.                             else{  
154.                                 //删除原先最大的  
155.                                 mapPeaks_max.erase(max_value_temp_iter);  
156.                                   
157.                                 temp = iter;  
158.                                 iter--;  
159.                                 max_value_temp_iter = temp;  
160.                             }  
161.                         }  
162.   
163.                     }  
164.                     //大于（rate + 1）*minpeakdistance  
165.                     else if (labs((*iter).first - high_max1_pos)>(rate + 1)*minpeakdistance){  
166.                         //初始化  
167.                         is_first_max = true;  
168.                         max_value_temp = 0;  
169.                         rate++;  
170.                         iter--;  
171.   
172.                     }  
173.                 }  
174.                 }  
175.             //因为容器指针的begin()指针在上面的循环中未进行判断，所以这里单独处理  
176.             //end()指的是最后一个元素的下一个位置所以不需要淡定进行判断  
177.             map<int, double>::iterator max_value_temp_iter_more;  
178.             if (mapPeaks_max.size() >= 2){  
179.                 max_value_temp_iter = mapPeaks_max.begin();  
180.                 max_value_temp_iter_more = max_value_temp_iter++;  
181.                 if (labs(max_value_temp_iter->first - max_value_temp_iter_more->first) < minpeakdistance){  
182.                     if (max_value_temp_iter->second > max_value_temp_iter_more->second){  
183.                         mapPeaks_max.erase(max_value_temp_iter_more);  
184.                     }  
185.                     else{  
186.                         mapPeaks_max.erase(max_value_temp_iter);  
187.                     }  
188.   
189.                 }  
190.             }  
191.         }  
192.         Vector_map_Peaks_max.clear();  
193.         //二维向量排序（根据函数选择以key或者value排序）  
194.         sortMapByValue(mapPeaks_max, Vector_map_Peaks_max);  
195.     }  
196.   
197. }  

/*
输入参数：
Matdata: 输入的数据矩阵
minpeakdistance：设定两峰值间的最小间隔数
minpeakheight  ：设定峰值的最小高度
输出函数：
mapPeaks_max        ：按照key  值（灰度级0-255）排序的
Vector_map_Peaks_max：按照value值（各灰度级的像素值）排序的
备注：此函数只完成了波峰的提取，波谷自行完成
- zhang
*/
void findPeaks(const Mat Matdata, int minpeakdistance, int minpeakheight, std::map <int, double> &mapPeaks_max, vector<pair<int, double>>&Vector_map_Peaks_max)
{
	int row = Matdata.rows - 1;
	int col = Matdata.cols;
	vector<int> sign;
	vector<double> markdata;
	vector<double> markdata1;
	//minMaxLoc寻找矩阵(一维数组当作向量,用Mat定义) 中最小值和最大值的位置. 
	//类型转换注意：data为unchar* 
	double* pMatdata = (double*)Matdata.data;
	for (int i = 1; i < col; i++)
	{
		/*相邻值做差：
		*小于0，赋-1
		*大于0，赋1
		*等于0，赋0
		*/
		markdata1.push_back(*(pMatdata + i-1));
		int diff = *(pMatdata + i) - *(pMatdata + i-1);
		if (diff > 0)
		{
			sign.push_back(1);
		}
		else if (diff < 0)
		{
			sign.push_back(-1);
		}
		else
		{
			sign.push_back(0);
		}
	}
	//再对sign相邻位做差  
	//保存极大值和极小值的位置  
	vector<int> indMax;
	vector<int> indMin;
	std::map <int, double> mapPeaks_min;
	for (int j = 1; j < sign.size(); j++)
	{
		int diff = sign[j] - sign[j - 1];
		if (diff < 0)
		{
			indMax.push_back(j);
			//根据峰值最小高度进行筛选
			if ( *(pMatdata +  indMax[indMax.size() - 1]) > minpeakheight){
				mapPeaks_max.insert(pair<int, double>(j, *(pMatdata + indMax[indMax.size() - 1])));
            }
			
		}
		else if (diff>0)
		{
			indMin.push_back(j);
			mapPeaks_min.insert(pair<int, double>(j, *(pMatdata + indMin[indMin.size() - 1])));
			
		}
	}

	if (mapPeaks_max.size() >= 2){

		//找到像素值最大的对应的容器指针
		std::map<int, double>::iterator iter_high_max1_pos = std::max_element(std::begin(mapPeaks_max), std::end(mapPeaks_max), cmpDescending);
		int high_max1_pos = iter_high_max1_pos->first;//最高点的位置
		map<int, double>::iterator  iter;
		map<int, double>::iterator  temp;
		//以位置为序
		map<int, double>::iterator  iter_high_max1_pos;
		map<int, double>::iterator  iter_high_max1_pos_temp1;
		map<int, double>::iterator  iter_high_max1_pos_temp2;
		//知道峰值最高的对应容器指针
        //我们以最高峰值处进行作业搜索所以在此声明两个局部变量保存该指针
		iter_high_max1_pos_temp1= iter_high_max1_pos;
		iter_high_max1_pos_temp2 = iter_high_max1_pos;
	   
	    //rate为宽度，以iter_high_max1_pos为起点,根据minpeakdistance的距离倍数搜索
		//minpeakdistance宽度内取一个最大的峰值max_value_temp
		int rate = 1;
		double max_value_temp = 0;
		map<int, double>::iterator max_value_temp_iter;
		//每个minpeakdistance内的第一个指针防止被删除
		bool is_first_max = true;
		//挨着最大峰值的
		if (mapPeaks_max.size() >= 2){
			if (iter_high_max1_pos_temp1 != mapPeaks_max.begin()){
				
				for (iter = --iter_high_max1_pos_temp1; iter != mapPeaks_max.begin(); iter--){
					//小于minpeakdistance
					if (labs((*iter).first - high_max1_pos) < (rate)*minpeakdistance){
							temp = iter;
							iter++;
							mapPeaks_max.erase(temp);
					}
					if (rate*minpeakdistance<=labs((*iter).first - high_max1_pos) && labs((*iter).first - high_max1_pos)<=(rate+1)*minpeakdistance){
						if ((*iter).second > max_value_temp){
							max_value_temp = (*iter).second;
							if (is_first_max){
								max_value_temp_iter = iter;
								is_first_max = false;
							}
							else{
								//删除原先最大的
								mapPeaks_max.erase(max_value_temp_iter);
								temp = iter;
								iter++;
								max_value_temp_iter = temp; 
							}
						}
						else{
							is_first_max = true;
							max_value_temp = 0;
							temp = iter;
							iter++;
							mapPeaks_max.erase(temp);
						}
					}
					//大于（rate + 1）*minpeakdistance
					else if (labs((*iter).first - high_max1_pos)>(rate + 1)*minpeakdistance){
						 
							rate++;
							iter++;
					}
				}
			}
			is_first_max = true;
			if (iter_high_max1_pos_temp2 != mapPeaks_max.end() && mapPeaks_max.size()>=2){
			
				for (iter = ++iter_high_max1_pos_temp2; iter != mapPeaks_max.end(); iter++){
					//TRACE("2Vector:%d \n", debug++);
					//小于minpeakdistance
					if (labs((*iter).first - high_max1_pos) < (rate)*minpeakdistance){
						temp = iter;
						iter--;
						mapPeaks_max.erase(temp);
					}
					if (rate*minpeakdistance <= labs((*iter).first - high_max1_pos) && labs((*iter).first - high_max1_pos) <= (rate + 1)*minpeakdistance){

						if ((*iter).second > max_value_temp){
							max_value_temp = (*iter).second;
							if (is_first_max){
								max_value_temp_iter = iter;
								is_first_max = false;
							}
							else{
								//删除原先最大的
								mapPeaks_max.erase(max_value_temp_iter);
								
								temp = iter;
								iter--;
								max_value_temp_iter = temp;
							}
						}

					}
					//大于（rate + 1）*minpeakdistance
					else if (labs((*iter).first - high_max1_pos)>(rate + 1)*minpeakdistance){
						//初始化
						is_first_max = true;
						max_value_temp = 0;
						rate++;
						iter--;

					}
				}
				}
			//因为容器指针的begin()指针在上面的循环中未进行判断，所以这里单独处理
			//end()指的是最后一个元素的下一个位置所以不需要淡定进行判断
			map<int, double>::iterator max_value_temp_iter_more;
			if (mapPeaks_max.size() >= 2){
				max_value_temp_iter = mapPeaks_max.begin();
				max_value_temp_iter_more = max_value_temp_iter++;
				if (labs(max_value_temp_iter->first - max_value_temp_iter_more->first) < minpeakdistance){
					if (max_value_temp_iter->second > max_value_temp_iter_more->second){
						mapPeaks_max.erase(max_value_temp_iter_more);
					}
					else{
						mapPeaks_max.erase(max_value_temp_iter);
					}

				}
			}
		}
		Vector_map_Peaks_max.clear();
		//二维向量排序（根据函数选择以key或者value排序）
		sortMapByValue(mapPeaks_max, Vector_map_Peaks_max);
	}

}

[cpp] view plain copy

print ?

1. //辅助函数  
2. //降序  
3. int cmpDescending (const pair<int, double>& x, const pair<int, double>&y)  
4. {  
5.     return  x.second<y.second;  
6. }  
7. //升序  
8. int cmpAscending(const pair<int, double>& x, const pair<int, double>&y)  
9. {  
10.     return  x.second>y.second;  
11. }  
12. //完成MAP到vector 的转化并且完成排序  
13. void sortMapByValue(map<int, double>& tMap, vector<pair<int, double>>& tVector)  
14. {  
15.     for (map<int, double>::iterator curr = tMap.begin(); curr != tMap.end(); curr++)  
16.     {  
17.         tVector.push_back(make_pair(curr->first, curr->second));  
18.     }  
19.       
20.     sort(tVector.begin(), tVector.end(), cmpAscending);  
21. }