2019-10-25 23:21:22 qq_43294951 阅读数 24
• ###### matlab 基于阈值分割的区域生长的医学图像分割3讲

matlab 基于阈值分割的区域生长的医学图像分割3讲

186 人正在学习 去看看 刘昱显

``````#include <iostream>
#include <opencv2/opencv.hpp>
#include <Eigen/Dense>
#include <math.h>
#include <stack>
#define _MATH_DEFINES_DEFINED
using namespace std;
using namespace cv;
using namespace Eigen;
float Otsu(const Mat& image,int T)
{
int nr = image.rows;
int nl = image.cols;
Mat hist = Mat::zeros(1, 256, CV_32F);
float* hi = hist.ptr<float>(0);
float junzhi=0;//定义总均值，类一均值，类二均值
float junzhi0 = 0;
float junzhi1 = 0;
float lei0 = 0;//定义类一，类二的概率
float lei1 = 0;
float fangcha0=0;//定义类一，类二的方差
float fangcha1 = 0;
float neifangcha = 0;//定义类内方差
float jianfangcha = 0;//定义类间方差
float zongfangcha = 0;//定义总方差，类内加类间
float ppppp = 0;
for (int i = 0; i < nr; i++)//这里是统计直方图
{
const float* im = image.ptr<float>(i);
for (int j = 0; j < nl; j++)
{
hi[int(im[j])] = hi[int(im[j])] +1;
}
}
//cout << hist << endl;
for (int i = 0; i < 256; i++)//计算均值
{
junzhi = i * hi[i] / nr / nl;
}
for (int i = 0; i <= T; i++)//计算
{
lei0 = lei0 + hi[i] / nr / nl;
}
for (int i = T+1; i <256; i++)
{
lei1 = lei1 + hi[i] / nr / nl;
}
for (int i = 0; i <= T; i++)
{
junzhi0 = junzhi0 + i * hi[i] / nr / nl / lei0;
}
for (int i = T + 1; i < 256; i++)
{
junzhi1 = junzhi1 + i * hi[i] / nr / nl / lei1;
}
for (int i = 0; i <= T; i++)
{
fangcha0 = fangcha0 + (i -junzhi0)* (i - junzhi0) * hi[i] / nr / nl / lei0;
}
for (int i = T + 1; i < 256; i++)
{
fangcha1 = fangcha1 + (i - junzhi1) * (i - junzhi1) * hi[i] / nr / nl / lei1;
}
neifangcha = lei0 * fangcha0 + lei1 * fangcha1;
jianfangcha = lei0 * lei1 * (junzhi1 - junzhi0) * (junzhi1 - junzhi0);
zongfangcha = neifangcha + jianfangcha;
ppppp= jianfangcha / zongfangcha;
return ppppp;
}
int getT(const Mat& image, Mat& new_image)
{
int temp = 0;
float max = 0;
max = Otsu(image, 20);
for (int i =  0; i < 255; i++)
{
float gg = Otsu(image, i);
if (max < gg)
{
temp = i;
max =gg;
}
}//遍历得到最大的分割点temp
for (int i = 0; i < image.rows; i++)
{
const float* im =image.ptr<float>(i);
float* p = new_image.ptr<float>(i);
for (int j = 0; j < image.cols; j++)
{
float ge = im[j];
if (ge <= temp)
{
p[j] = 0;
}
else
{
p[j] = 255;
}

}
}//图像进行二值化
return temp;
}
int main()
{
imshow("原图", image);
Mat IM;
image.convertTo(IM, CV_32F);
//imshow("fsSSd", IM);
Mat new_image = Mat::zeros(IM.rows, IM.cols, CV_32F);
int T = getT(IM, new_image);
cout << T << endl;
imshow("分割以后", new_image);
cv::waitKey(0);
}
``````

ok，得到的阈值分割图如下，感觉海星

2019-04-22 16:24:12 Lemon_jay 阅读数 529
• ###### matlab 基于阈值分割的区域生长的医学图像分割3讲

matlab 基于阈值分割的区域生长的医学图像分割3讲

186 人正在学习 去看看 刘昱显

# 使用C++、opencv对图像进行简单的阈值分割

0----RGB        R、G、B值大于阈值的点设为黑色
1----|G-R|        G-R值的绝对值小于threshod_value的点设为黑色
2----|2G-R-B|        2G-R-B值的绝对值大于threshod_value的点设为黑色
3----G/R        G/R值小于threshod_value的点设为黑色

int createTrackbar(const string& trackbarname, const string& winname, int* value, int count, TrackbarCallback onChange=0,                                    void* userdata=0)

``````#include "stdafx.h"
#include <opencv2\opencv.hpp>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <iostream>
using namespace cv;

Mat src,dst;

//滑动条相关全局变量
int mode = 0;
int r = 255, g = 255, b = 255;
int threshod_value = 70;

//滑动条回调函数声明
void on_mode(int, void*);
void on_threshold(int, void*);

//阈值处理函数声明
void RGB_threshold();
void G_B_threshold();
void G_B_R_threshold();
void G_R_threshold();

//帮助文本
void show_help()
{
std::cout << "---------help---------" << std::endl;
std::cout << "mode:" << std::endl;
std::cout << "0----RGB		R、G、B值大于阈值的点设为黑色" << std::endl;
std::cout << "1----|G-R|		G-R值的绝对值小于threshod_value的点设为黑色" << std::endl;
std::cout << "2----|2G-R-B|		2G-R-B值的绝对值大于threshod_value的点设为黑色" << std::endl;
std::cout << "3----G/R		G/R值小于threshod_value的点设为黑色" << std::endl;
}
int main()
{
system("color 02");

show_help();

//读取图像
namedWindow("src", WINDOW_NORMAL);
imshow("src", src);
namedWindow("dst", WINDOW_NORMAL);

//创建滑动条
createTrackbar("mode", "src", &mode, 4, on_mode);
on_mode(0, 0);

createTrackbar("R", "src", &r, 255, on_threshold);
on_threshold(0, 0);

createTrackbar("G", "src", &g, 255, on_threshold);
on_threshold(0, 0);

createTrackbar("B", "src", &b, 255, on_threshold);
on_threshold(0, 0);

createTrackbar("thre_val", "src", &threshod_value, 255, on_threshold);
on_threshold(0, 0);

waitKey();
return 0;
}

void on_mode(int, void*)
{
}
void on_threshold(int, void*)
{
switch (mode)
{
case 0:
RGB_threshold();
break;
case 1:
G_B_threshold();
break;
case 2:
G_B_R_threshold();
break;
case 3:
G_R_threshold();
break;
}
}

void RGB_threshold()
{
//复制源图到目标图像
dst = src.clone();
//遍历图像每个像素点
for(int i=0;i<src.rows;i++)
for (int j = 0; j < src.cols; j++)
{
//大于阈值的点设为黑色
if (src.at<Vec3b>(i, j)[2] > r)
{
dst.at<Vec3b>(i, j)[0] = 0;
dst.at<Vec3b>(i, j)[1] = 0;
dst.at<Vec3b>(i, j)[2] = 0;
}
if (src.at<Vec3b>(i, j)[1] > g)
{
dst.at<Vec3b>(i, j)[0] = 0;
dst.at<Vec3b>(i, j)[1] = 0;
dst.at<Vec3b>(i, j)[2] = 0;
}
if (src.at<Vec3b>(i, j)[0] > b)
{
dst.at<Vec3b>(i, j)[0] = 0;
dst.at<Vec3b>(i, j)[1] = 0;
dst.at<Vec3b>(i, j)[2] = 0;
}
}
imshow("dst", dst);
}
void G_B_threshold()
{
dst = src.clone();
for (int i = 0; i<src.rows; i++)
for (int j = 0; j < src.cols; j++)
{
if ((    abs(src.at<Vec3b>(i, j)[1] - src.at<Vec3b>(i, j)[2])	<	threshod_value	))
{
dst.at<Vec3b>(i, j)[0] = 0;
dst.at<Vec3b>(i, j)[1] = 0;
dst.at<Vec3b>(i, j)[2] = 0;
}
}
//开操作
//Mat element = getStructuringElement(MORPH_RECT, Size(3, 3));
//morphologyEx(dst, dst, MORPH_OPEN, element);
imshow("dst", dst);
}
void G_B_R_threshold()
{
dst = src.clone();
for (int i = 0; i<src.rows; i++)
for (int j = 0; j < src.cols; j++)
{
if ((abs(2*src.at<Vec3b>(i, j)[1] - src.at<Vec3b>(i, j)[2] - src.at<Vec3b>(i, j)[0])	>	threshod_value))
{
dst.at<Vec3b>(i, j)[0] = 0;
dst.at<Vec3b>(i, j)[1] = 0;
dst.at<Vec3b>(i, j)[2] = 0;
}
}
//腐蚀操作
//Mat element = getStructuringElement(MORPH_RECT, Size(3, 3));
//morphologyEx(dst, dst, MORPH_ERODE, element);
imshow("dst", dst);
}
void G_R_threshold()
{
dst = src.clone();
int height = dst.rows;
int width = dst.cols;
for (int row = 0; row < height; row++)
{
for (int col = 0; col < width; col++)
{
double a = (double)dst.at<Vec3b>(row, col)[1] / dst.at<Vec3b>(row, col)[2];
a = a * 100;
if (a < threshod_value)
{
dst.at<Vec3b>(row, col)[0] = 0;
dst.at<Vec3b>(row, col)[1] = 0;
dst.at<Vec3b>(row, col)[2] = 0;
}
}
}
imshow("dst", dst);
}``````

2013-01-24 14:14:54 ily6418031hwm 阅读数 13028
• ###### matlab 基于阈值分割的区域生长的医学图像分割3讲

matlab 基于阈值分割的区域生长的医学图像分割3讲

186 人正在学习 去看看 刘昱显

由于阈值处理直观、实现简单且计算速度快，因此图像阈值处理在图像分割中处于核心地位，下面我会为大家介绍阈值处理的方法，并用OpenCV给出实现的代码。

1.  处理流程：

1.为全局阈值选择一个初始估计值T(图像的平均灰度)。
2.用T分割图像。产生两组像素：G1有灰度值大于T的像素组成，G2有小于等于T像素组成。
3.计算G1和G2像素的平均灰度值m1和m2；
4.计算一个新的阈值:T = (m1 + m2) / 2;
5.重复步骤2和4,直到连续迭代中的T值间的差为零。

此种方法主要适用于图像直方图有明显波谷。

2.源代码：

```/*******************************************************************************

2.用T分割图像。产生两组像素：G1有灰度值大于T的像素组成，G2有小于等于T像素组成。
3.计算G1和G2像素的平均灰度值m1和m2；
4.计算一个新的阈值:T = (m1 + m2) / 2;
5.重复步骤2和4,直到连续迭代中的T值间的差为零。
*******************************************************************************/
int CThreasholdProcess::BasicGlobalThreshold(int*pg,int start,int end)
{
int  i,t,t1,t2,k1,k2;
double u,u1,u2;
t=0;
u=0;

for (i=start;i<end;i++)
{
t+=pg[i];
u+=i*pg[i];
}

k2=(int) (u/t);                          //  计算此范围灰度的平均值
do
{
k1=k2;
t1=0;
u1=0;
for (i=start;i<=k1;i++)
{             //  计算低灰度组的累加和
t1+=pg[i];
u1+=i*pg[i];
}

t2=t-t1;
u2=u-u1;
if (t1)
u1=u1/t1;                     //  计算低灰度组的平均值
else
u1=0;
if (t2)
u2=u2/t2;                     //  计算高灰度组的平均值
else
u2=0;
k2=(int) ((u1+u2)/2);                 //  得到新的阈值估计值
}while(k1!=k2);                           //  数据未稳定，继续

return k1;                              //  返回阈值
}```
3.演示结果

a).被处理的源图像                                                                                                                         b).基本的全局阈值处理后图像                                                                      c）源图直方图

2012-02-11 16:18:32 renshengrumenglibing 阅读数 48377
• ###### matlab 基于阈值分割的区域生长的医学图像分割3讲

matlab 基于阈值分割的区域生长的医学图像分割3讲

186 人正在学习 去看看 刘昱显

1．引言

2．阈值分割的基本概念

（原始图像）    （阈值分割后的二值化图像）

T(x，y，N(x，y)，f(x，y))

点相关的全局阈值T＝T(f(x，y))
(只与点的灰度值有关)

(与点的灰度值和该点的局部邻域特征有关)
局部阈值或动态阈值T＝T(x，y，N(x，y)，f(x，y))
(与点的位置、该点的灰度值和该点邻域特征有关)

1) 基于点的全局阈值方法；
2) 基于区域的全局阈值方法
3) 局部阈值方法和多阈值方法

3．基于点的全局阈值选取方法
3.1  p-分位数法
1962年Doyle[10]提出的p-分位数法(也称p-tile法)可以说是最古老的一种阈值选取方法。该方法使目标或背景的像素比例等于其先验概率来设定阈值，简单高效，但是对于先验概率难于估计的图像却无能为力。

3.2  迭代方法选取阈值[11]

3.3  直方图凹面分析法

Rosenfeld和Torre[12]提出可以构造一个包含直方图 的最小凸多边形 ，由集差 确定 的凹面。若 和 分别表示 与 在灰度级之处的高度，则 取局部极大值时所对应的灰度级可以作为阈值。也有人使用低通滤波的方法平滑直方图，但是滤波尺度的选择并不容易[13]。

3.4 最大类间方差法

，  ，  .                 (3)

在实际运用中，往往使用以下简化计算公式：
（T） ＝ WA(μa-μ)2  + Wb(μb-μ)2

3.5 熵方法

(4)

(5)

3.6 最小误差阈值

Eb(T)是目标类错分到背景类的概率，Eo(T)是背景类错分到目标类的概率

(9)

(10)

3.7 矩量保持法

3.8 模糊集方法

Pal等[21]首先，他们把一幅具有 个灰度级的 图像看作一个模糊集 ，其中隶属函数 定义如下：
(11)

= min{μ(l),1-μ(l)}

3.9 小结

4 基于区域的全局阈值选取方法

4.1 二维熵阈值分割方法[25]

（二维灰度直方图： 灰度－领域平均灰度）

Abutaleb[26]，和Pal]结合Kapur]和Kirby的方法，分别提出了各自的二维熵阈值化方法，其准则函数都是使目标熵和背景熵之和最大化。Brink[27]的方法则是使这两者中的较小者最大化，该方法的计算复杂度为 ，后来有人改进为递推快速算法将时间复杂度降为 (其中 为最大灰度级数)。

4.2  简单统计法
Kittler等人[28],[29]提出一种基于简单的图像统计的阈值选取方法。使用这种方法，阈值可以直接计算得到，从而避免了分析灰度直方图，也不涉及准则函数的优化。该方法的计算公式为
(19)

4.3  直方图变化法

例如，由于目标区的象素具有一定的一致性和相关性，因此梯度值应该较小，背景区也类似。而边界区域或者噪声，就具有较大的梯度值。最简单的直方图变换方法，就是根据梯度值加权，梯度值小的象素权加大，梯度值大的象素权减小。这样，就可以使直方图的双峰更加突起，谷底更加凹陷。

4.4 其它基于区域的全局阈值法

5 局部阈值法和多阈值法

5.1 局部阈值（动态阈值）

（阈值低，对亮区效果好，则暗区差）          （阈值高，对暗区效果好，则亮区差）

（按两个区域取局部阈值的分割结果）

```/************************************************************************/
/* 全局阈值分割  自动求取阈值        */
/************************************************************************/
//自动求取阈值，增加对场景的适应性
//只需求取一次，之后就可以一直使用
#include<cv.h>
#include <highgui.h>
#include <iostream>
#include <math.h>
using namespace std;
int main(){
IplImage * image,* image2;
cvNamedWindow("image",1);
cvShowImage("image",image);
image2 = cvCreateImage(cvSize(image->width,image->height),image->depth,1);
double T = 0;
double dT0 = 1.0;//阈值求取结束标志
double dT = 255.0;

//求取平均灰度，作为阈值T的初始值T0
int i, j;
double T0 = 0,T1 = 0,T2 = 0;//初始阈值
int count1,count2;
unsigned char * ptr,*dst;
for (i = 0 ; i< image->height ; i++)
{
for (j =0 ; j < image->width;j++)
{
ptr = (unsigned char *)image->imageData + i*image->widthStep + j;
T0 += ((double)(*ptr))/image->width/image->height;
}
}
cout<<"T0:     "<<T0<<endl;
T = (int)(T0 + 0.5);
//计算T两侧的灰度平均值，然后把二者的均值赋值给T
while (dT > dT0)
{

T1 = 0;
T2 = 0;
count1 = 0;
count2 = 0;
for (i = 0 ; i< image->height ; i++)
{
for (j =0 ; j < image->width;j++)
{
ptr = (unsigned char *)image->imageData + i*image->widthStep + j;
if (*ptr > T)
{
T1 += ((double)(*ptr))/image->width/image->height;
count1++;
}
else if(*ptr < T)
{
T2 +=  ((double)(*ptr))/image->width/image->height;
count2++;
}
}
}

T1 = T1*image->width*image->height/count1;
T2 = T2*image->width*image->height/count2;
dT = fabs(T - (T1 + T2)/2);

cout<<"T1"<<T1<<endl;
cout<<"T2"<<T2<<endl;
cout<<"dT  " << dT<<endl;
T = (T1 + T2)/2;
cout<<"T:     "<<T<<endl;
}

//根据求取的阈值进行分割
for (i = 0 ; i< image2->height ; i++)
{
for (j =0 ; j < image2->width;j++)
{
ptr = (unsigned char *)image->imageData + i*image->widthStep + j;
dst = (unsigned char *)image2->imageData+i*image2->widthStep+j;
if (*ptr > T)
{
*dst = 255;
}
else
{
*dst =0;
}
}
}

cvNamedWindow("image2",1);
cvShowImage("image2",image2);
cvSaveImage("E:\\image\\dowels2.tif",image2);
cvWaitKey(0);
return 0;
}```

5.1.1 阈值插值法
首先将图像分解成系列子图，由于子图相对原图很小，因此受阴影或对比度空间变化等带来的问题的影响会比较小。然后对每个子图计算一个局部阈值（此时的阈值可用任何一种固定阈值选取方法）。通过对这些子图所得到的阈值进行插值，就可以得到对原图中每个象素进行分割所需要的合理阈值。这里对应每个象素的阈值合起来构成的一个曲面，叫做阈值曲面。

5.1.2 水线阈值算法

5.1.3 其它的局部阈值法

文献[31] [32]中提出一种基于局部梯度最大值的插值方法。首先平滑图像，并求得具有局部梯度最大值的像素点，然后利用这些像素点的位置和灰度在图像上内插，得到灰度级阈值表面。

5.2 多阈值法

5.2.1 基于小波的多域值方法。

5.2.2 基于边界点的递归多域值方法。

5.2.3 均衡对比度递归多域值方法。

6 阈值化算法评价简介

然而，如同所有的图像分割方法一样，阈值化结果的评价是一个比较困难的问题。事实上对图像分割本身还缺乏比较系统的精确的研究，因此对其评价则更差一些。人们先后已经提出了几十个评价准则。这些准则中又有定性的，也有定量的；有分析算法的，也有检测实验结果的，文献[37]将它们大致分为13类。

2019-07-06 21:20:18 weixin_43124455 阅读数 348
• ###### matlab 基于阈值分割的区域生长的医学图像分割3讲

matlab 基于阈值分割的区域生长的医学图像分割3讲

186 人正在学习 去看看 刘昱显

# 理论

otsu理论的话， 看冈萨雷斯的《数字图像处理》 第10.3.6节 多阈值处理

# MATLAB代码

otsu 双阈值函数：

``````

function [t1,t2]=DoubleOtsuThresh(img)
%
%  Otsu 双阈值求解
%  输入 图像img，输出 最优阈值t1和t2(归一化,范围在[0,1])
%
%

BinsNum = 256;
hist = imhist(img,BinsNum);
p = hist / sum(hist);          % 直方图的概率密度函数
mG= sum(p .* (1:BinsNum)');     % 全局均值

P1 = cumsum(p);                 % 概率分布
m1 = cumsum(p .* (1:BinsNum)')./P1; % 256*1  每个阈值的前景平均灰度

% 根据算法理论，从k2+1累加到L-1，可以先倒着累加再翻转回来
P3= cumsum(flip(p));
m3 = cumsum(flip(p) .* flip(1:BinsNum)')./P3;
P3=flip(P3);
P3=[P3(2:end) ;0];             %P3的索引用k2，则P3(1)实际上是从p(2)加到p(256)，所以要移动一个
%移动后也用不到P3(1)这个值，因为k2>k1>1
m3=flip(m3);
m3=[m3(2:end) ;0];

% m2 P2为 k1*k2的索引矩阵，为 256*256 的上三角矩阵 因为k1<k2
m2=zeros(BinsNum,BinsNum);
P2=zeros(BinsNum,BinsNum);

for k1=1:BinsNum-2
for k2=k1+1:BinsNum-1
P2(k1,k2)=1-P1(k1)-P3(k2);
m2(k1,k2) = sum( (k1+1:k2)' .* p(k1+1:k2) )./ P2(k1,k2) ;
end
end

% 遍历k1,k2各种组合  求方差
variance=zeros(BinsNum,BinsNum);
for k1=1:BinsNum-2
for k2=k1+1:BinsNum-1

% variance 为256*256 的上三角矩阵 因为l1<k2
variance(k1,k2) = P1(k1)*(m1(k1)-mG)^2 + ...
P2(k1,k2)*(m2(k1,k2)-mG)^2 + ...
P3(k2)*(m3(k2)-mG)^2;

end
end

% 最大方差点即为最优阈值
[~,index] = max(variance(:));
[index_row,index_col] = ind2sub(size(variance),index);

%灰度值为[0,255] 因此需要-1
t1= (index_row-1)./(BinsNum-1);
t2= (index_col-1)./(BinsNum-1);

end
``````

``````
function out=img2gray(img,t1,t2)
[M,N]=size(img);
out=zeros(M,N);

for m=1:M
for n=1:N
if img(m,n)<t1
out(m,n)=0;
else
if( img(m,n)>t1 && img(m,n)<t2)
out(m,n)=127;  %灰度级可调
else
out(m,n)=255;	%灰度级可调
end
end
end
end
end

``````

``````%%  图B 彩色图
imgB=rgb2gray(imgB);
imgB=im2double(imgB);

%% 先高斯滤波平滑图B，去高频，再用阈值处理
gaussH=fspecial('gaussian',[5 5],10);
smoothB=imfilter(imgB,gaussH);

[t1B,t2B] =DoubleOtsuThresh(smoothB);

% 如果imgB是unit8类型,则灰度范围在[0,255] 要把算法得到的两个阈值乘255
% 如果imgB是[0,1]范围的double类型则不用乘
%  T1B=t1B*255;
%  T2B=t2B*255;

T1B=t1B;
T2B=t2B;

J9=img2gray(smoothB,T1B,T2B); %转为灰度图

figure
imshow(J9,[]);title('平滑+双阈值分割'); % imshow里面J9后面的[]相当于直方图均衡效果，让灰度值均匀分布,可以去掉对比一下
``````

``````
%%
% —— 图D —— %

%   imgD=rgb2gray(imgD);
[t1D,t2D] =DoubleOtsuThresh(imgD);

T1D=t1D*255;
T2D=t2D*255;
% 分为 0 127 255 显示

J8=img2gray(imgD,T1D,T2D);

figure
subplot 121
imshow(imgD);title('原图D');

subplot 122
imshow(J8,[]);title('双阈值分割');

``````

``````
%% 用书上的图测试算法,输出的两个阈值和书上结果一样 80和177

[t1,t2] =DoubleOtsuThresh(imgC);
T1=t1*255;
T2=t2*255;
J7=img2gray(imgC,T1,T2);
figure
subplot 121
imshow(imgC);title('原图C');
subplot 122
imshow(J7,[]);title('双阈值');
disp(['     otsu双阈值为：   ' ,'t1=',num2str(T1),'     t2=',num2str(T2)])

``````

The end .
( markdown如何支持matlab的代码片啊 ,为什么matlab这个key没有效果，只能沦落到用js这个key来高亮代码片)