根据模糊阈值分割的原理，本文在“模糊阈值分割代码（1）"的基础上进行了拓展，主要是利用图像直方图，搜索最佳波谷阈值，实现多阈值分割。
其过程主要如下：
1. 获取图像直方图
2. 进行直方图平滑，以去毛刺
3. 进行波峰搜索。这部分比较复杂，需要考虑的问题也多。包括波峰的判定准则、伪峰值的去除、半峰值的判断等等。
一般来说，一个峰值包括起始位置（上坡段）、终止位置（下坡段）和峰值位置。本文上、下坡段判断是通过连续增长（下降）的数量来判断的，比如连续超过10个增长的就认为是上坡段，而连续10个下降的就是下坡段。当然了，有些峰值非常小，这个就需要缩小判定阈值。峰值的位置是在上下坡段的区间内取重心的方式得到的（见函数GetPeaks）。
4. 在波峰之间采用S函数的模糊化（MemshipFunc），获得模糊曲线。我们知道，模糊化的过程窗口的宽度很重要，一般窗宽取峰值距离的0.3~0.8。由于上一步把峰值取到了，所以这一步就可以根据峰值间的距离设定窗宽，达到自适应的目的。
5. 获得模糊率曲线后，搜索不同峰值间的最小值就可以获得分割阈值。
struct Peack
{
int nStat;
int nEnd;
int nValue;
Peack():nStat(0),nEnd(0),nValue(0)
{}
};
/*
名称：BlurCure
参数：lpBmp - 图像数据指针
lSrcWidth - 图像宽度
lSrcHeight - 图像高度
dLineBites - 图像单行数据量
说明：获取模糊阈值分割的模糊曲线dBlurData
*/
void BlurCurve(LPSTR lpBmp, long lSrcWidth, long lSrcHeight, DWORD dLineBites)
{
if(lpBmp == NULL)
return;
int i,j;
LPSTR lpSrcPtr = NULL;
BYTE btTempValue = 0;
int nTistogram[256] = {0};
for (i = 0;i < lSrcHeight; ++i)
{
lpSrcPtr = lpBmp + i * dLineBites;
for (j = 0;j <lSrcWidth; ++j)
{
btTempValue = (BYTE)*(lpSrcPtr + j);
nTistogram[btTempValue]++;
}
}
//-----------------直方图平均-----------
int nSum = 0, nMaxTtg = 0;
const int c_nScop = 5;
int arrValue[c_nScop] = {0};
for (j = 0;j <256; ++j)
{
nSum += nTistogram[j];
if(j < c_nScop)
{
arrValue[j] = nTistogram[j];
nTistogram[j] = nSum/(j+1);
}
else
{
i = j%c_nScop;
nSum -= arrValue[i];
arrValue[i] = nTistogram[j];
nTistogram[j] = nSum/c_nScop;
}
if(nMaxTtg < nTistogram[j])
nMaxTtg = nTistogram[j];
outFile<<nTistogram[j]<<endl;
}//end
//--------------------------------------- vector <Peack> vctPeak;
GetPeaks(nTistogram,nMaxTtg/500,vctPeak);
int nNum = vctPeak.size();
if(nNum > 1)
{
int arrThesh[255] = {0};
int nSetScop = 255/(nNum -1);
for (i=0; i<nNum-1; i++)
{
int nWinWidth = int(0.3*(vctPeak[i+1].nValue-vctPeak[i].nValue));
double dGetMemship = 0,dMinData = 0xFFFFFFFF;
for (int nGrayPos = vctPeak[i].nValue;nGrayPos <= vctPeak[i+1].nValue; ++nGrayPos)
{
double dOutData = 0;
for (j = 0;j <256; ++j)
{
dGetMemship = MemshipFunc(j,nGrayPos,nWinWidth);
if(dGetMemship > 0.5)
dGetMemship = 1 - dGetMemship;
if(dGetMemship < 0)
dGetMemship = 0;//-dGetMemship;
dOutData += nTistogram[j]*dGetMemship;
}
//dOutData *= 2.0/(lSrcHeight*lSrcWidth);
if(dMinData > dOutData)
{
dMinData = dOutData;
arrThesh[i] = nGrayPos;
}
}
}
//---------
for (i = 0;i < lSrcHeight; ++i)
{
lpSrcPtr = lpSrcData + i * dLineBites;
lpDstPtr = lpDstData + i * dLineBites;
for (j = 0;j <lSrcWidth; ++j)
{
btTempValue = (BYTE)*(lpSrcPtr + j);
int k;
for ( k=0; k<nNum-1; k++)
{
if(btTempValue < arrThesh[k])
{
*(lpDstPtr + j)= (BYTE)(k*nSetScop);
break;
}
}
if(k == nNum-1)
*(lpDstPtr + j) = 255;
}
}//END FOR
}
else if(nNum == 1)//单峰，用峰值左侧点为分割点
{
for (i = 0;i < lSrcHeight; ++i)
{
lpSrcPtr = lpSrcData + i * dLineBites;
lpDstPtr = lpDstData + i * dLineBites;
for (j = 0;j <lSrcWidth; ++j)
{
btTempValue = (BYTE)*(lpSrcPtr + j);
if(btTempValue < vctPeak[0].nEnd)
*(lpDstPtr + j)= 0;
else
*(lpDstPtr + j) = 255;
}
}
}
else
{
MessageBox(_T("分割失败"),_T("提示"),MB_ICONWARNING);
}
}
/*
函数:GetPeaks
参数: pTistogram - 图像直方图
nLowValue - 谷底最小值
vctPeak - 返回的峰值列表
说明：函数主要对处理过的直方图查找谷峰
*/
void GetPeaks(const int pTistogram[256], const int nLowValue,vector<Peack> &vctPeak)
{
int arrPos[2] = {0};
int nSratMin = 0xFFFFFF,nContinueNum=0;
int nPreValue = pTistogram[0];
int nSetThres = 8;
for (int i=1; i<256; i++)
{
if(nSratMin > pTistogram[i])
{
nSratMin = pTistogram[i];
if(nSratMin <= nLowValue)
nSetThres = 5;
}
//-------------------------------
if(arrPos[0] != 0)
{
if(pTistogram[i] < nPreValue)
{
nContinueNum ++;
if(nContinueNum >= nSetThres)
arrPos[1] = i;
}
else// if(nTistogram[i] > nPreValue)
{
if(arrPos[1] != 0)//data out
{
Peack pk;
pk.nStat = arrPos[0];
pk.nEnd = arrPos[1];
vctPeak.push_back(pk);
//--------
arrPos[0] = arrPos[1] =0;
nSratMin = 0xFFFFFF;
nSetThres = 10;
}
nContinueNum = 0;
}
}
else
{
if(pTistogram[i] > nPreValue)
{
nContinueNum ++;
if(nContinueNum >= nSetThres)
arrPos[0] = i-nSetThres;
}
else //if(nTistogram[i] < nPreValue)
{
static int s_nLowNum = 0;
if(pTistogram[i] < nPreValue)
{
s_nLowNum ++;
if(s_nLowNum == i && s_nLowNum > 10)
arrPos[0] = 1;
}
else if(nContinueNum > 2)
nContinueNum ++;
else
nContinueNum = 0;
}
}
/////////////////////////////////
nPreValue = pTistogram[i];
}
//============峰值计算==============
int nSum = 0,nTistoSum = 0;
vector<Peack>::iterator it=vctPeak.begin();
while (it != vctPeak.end())
{
nSum = nTistoSum = 0;
for (int i=(*it).nStat; i<=(*it).nEnd; i++)
{
nSum += i*pTistogram[i];
nTistoSum += pTistogram[i];
}
(*it).nValue = nSum/nTistoSum;
it++;
}
}
下面是“模糊阈值分割代码（1）”中的瓶底图像获得的峰值信息：
三个阈值信息：
最终的分割结果：
从图上可以看出，该种方法还是比较准确的得到了分割的阈值，结果还是比较理想的。
当然，这是针对的多峰值的情况，对于单峰值的情况也是适用的。
其实这种算法在上世纪80年代就已经推出了，之所以没有推广开，问题是在求解的过程中有很多不确定性。比如峰值的确定，本身峰值怎么求的准确就是一个困难的问题，特别是在复杂的情况下，比如弱灰度、有噪音的情况，峰值往往并不明显。而要获得自私应的模糊阈值的，其窗宽又依赖于峰值选取的结果，结果可想而知，通用性就大大降低了！要解决这类算法的通用性问题，首先要减小不确定性，这也是现在依旧在探讨的问题。