程序?qū)崿F(xiàn)：

    上面的公式看似復(fù)雜，其實(shí)我們關(guān)心的就是最后的5個(gè)計(jì)算步驟,，這里說(shuō)明一下,，有的書上以隸屬度矩陣的某一范數(shù)小于一定值作為收斂的條件，這也可,，不過計(jì)算量稍微要大一點(diǎn)了。

        程序采用VB6.0編制,，完全按照以上的步驟進(jìn)行,。

‘程序?qū)崿F(xiàn)功能：模糊聚類和硬聚類
‘作    者： laviewpbt
‘聯(lián)系方式： [email protected]
‘QQ:33184777
‘版本：Version 2.3.1
‘說(shuō)明：復(fù)制請(qǐng)保留源作者信息，轉(zhuǎn)載請(qǐng)說(shuō)明,，歡迎大家提出意見和建議

Private Declare Function GetTickCount Lib "kernel32" () As Long

Private Enum IniCenterMethod    ‘初始中心的方法
    CreateRandom                ‘隨機(jī)的中心點(diǎn)
    CreateByHcm                 ‘由HCM創(chuàng)建的中心點(diǎn)
    CreateByRandomZadeh         ‘先隨機(jī)創(chuàng)建隸屬矩陣,，然后計(jì)算得到的中心點(diǎn)

    CreateByHand                ‘手工確定初始中心點(diǎn)

End Enum

Private Enum AntiFuzzyMethod    ‘反模糊的方法
    Max                         ‘最大隸屬度法
    Middle                      ‘中位數(shù)法
    Mean                        ‘加權(quán)均值法
End Enum

Private Type FcmInfo
     Center() As Double         ‘聚類中心
     Degree() As Double         ‘隸屬度，為Double類型
     Class() As Byte            ‘記錄數(shù)據(jù)屬于那一類
     TimeUse As Long            ‘所用時(shí)間
     Iterations  As Long        ‘迭帶次數(shù)
     ErrMsg As String           ‘錯(cuò)誤信息
End Type

Private Type HcmInfo
    Center() As Double          ‘聚類中心
    Class() As Byte             ‘記錄數(shù)據(jù)屬于那一類
    TimeUse As Long             ‘所用時(shí)間
    Iterations  As Long         ‘迭帶次數(shù)
    ErrMsg As String            ‘錯(cuò)誤信息
End Type

‘*************************************************************************************
‘*    作    者 ：    laviewpbt
‘*    函 數(shù) 名 ：    Fcm
‘*    參    數(shù) ：    Data     -   待分類的樣本,，第一維的大小表示樣本的個(gè)數(shù),，
‘*                                第二維的大小表示樣本的維數(shù)
‘*                   Cluster  -   分類數(shù)
‘*                   CreateIniCenter - 初始聚類中心的創(chuàng)建方法
‘*                   AntiFuzzy -  反模糊化的方法
‘*                   Exponent  -  一個(gè)控制聚類效果的參數(shù)，一般取2
‘*                   Maxiterations  - 最大的迭代次數(shù)
‘*                   MinImprovement - 最小的改進(jìn)參數(shù)（兩次迭代間聚類中心的距離）
‘*    返回值 ：      FcmInfo結(jié)構(gòu),，記錄了相關(guān)的參數(shù)
‘*    功能描述 ：    利用模糊理論的聚類方法把數(shù)據(jù)分類
‘*    日    期 ：    2004-10-27 10.25.32
‘*    修 改 人 ：    laviewpbt
‘*    日    期 ：    2006-11-7 19.25.31
‘*    版    本 ：    Version 2.3.1
‘**************************************************************************************

Private Function Fcm(ByRef Data() As Double, ByVal Cluster As Long, Optional ByVal CreateIniCenter As IniCenterMethod = IniCenterMethod.CreateByHcm, Optional AntiFuzzy As AntiFuzzyMethod = Max, Optional Exponent As Byte = 2, Optional Maxiterations As Long = 400, Optional MinImprovement As Double = 0.01, Optional ByRef CenterByHandle As Variant) As FcmInfo
    If ArrayRange(Data) <> 2 Then
        Fcm.ErrMsg = "數(shù)據(jù)只能為二維數(shù)組"
        Exit Function
    End If
    Dim i As Long, j As Long, k As Long, l As Long, m As Long
    Dim DataNumber As Long, DataSize As Long
    Dim Temp As Double, Sum1 As Double, Sum2 As Double, Sum3 As Double, Index As Integer
    Dim OldCenter() As Double
    Fcm.TimeUse = GetTickCount
    DataNumber = UBound(Data, 1): DataSize = UBound(Data, 2)
    ReDim Fcm.center(1 To Cluster, 1 To DataSize) As Double
    ReDim Fcm.Degree(1 To Cluster, 1 To DataNumber) As Double
    ReDim Fcm.Class(1 To DataNumber) As Byte
    ReDim OldCenter(1 To Cluster, 1 To DataSize) As Double
    On Error GoTo ErrHandle:
    Randomize
    If CreateIniCenter = CreateRandom Then
        For i = 1 To Cluster
            For j = 1 To DataSize
                OldCenter(i, j) = Data(Rnd * DataNumber, j)    ‘產(chǎn)生隨機(jī)初始中心點(diǎn)
            Next
        Next
    ElseIf CreateIniCenter = CreateByHcm Then
        Dim HcmCenter As HcmInfo
        HcmCenter = Hcm(Data, Cluster)
        For i = 1 To Cluster
            For j = 1 To DataSize
                OldCenter(i, j) = HcmCenter.center(i, j)   ‘產(chǎn)生HCM初始中心點(diǎn)
            Next
        Next
    ElseIf CreateIniCenter = CreateByRandomZadeh Then
        ReDim RndDegree(1 To Cluster, 1 To DataNumber) As Double
        Dim RndSum As Double
        For i = 1 To Cluster
            For j = 1 To DataNumber
                RndDegree(i, j) = Rnd           ‘創(chuàng)建隨機(jī)的隸屬度
            Next
        Next
        For j = 1 To DataNumber
            RndSum = 0
            For i = 1 To Cluster
                RndSum = RndSum + RndDegree(i, j)
            Next
            For i = 1 To Cluster
                RndDegree(i, j) = RndDegree(i, j) / RndSum   ‘隸屬度矩陣每列之后必須為1
            Next
        Next
       
        For i = 1 To Cluster
            For j = 1 To DataSize
                Sum1 = 0: Sum2 = 0
                For k = 1 To DataNumber
                    Temp = Exp(Log(RndDegree(i, k)) * Exponent)  ‘其實(shí)就是RndDegree(i, k)^Exponent
                    Sum1 = Sum1 + Temp * Data(k, j)           ‘隸屬度的平方乘以數(shù)值
                    Sum2 = Sum2 + Temp                        ‘隸屬度的和
                Next
                OldCenter(i, j) = Sum1 / Sum2                 ‘得到聚類中心
            Next
        Next
    ElseIf CreateIniCenter = CreateByHand Then
        If IsMissing(CenterByHandle) Then
            Fcm.ErrMsg = "請(qǐng)?zhí)峁┏跏季垲愔行摹?"
            Exit Function
        ElseIf UBound(CenterByHandle, 1) <> Cluster Or UBound(CenterByHandle, 2) <> DataSize Then
            Fcm.ErrMsg = "手工提供的初始聚類中心維數(shù)有錯(cuò)誤."
            Exit Function
        End If
        For i = 1 To Cluster
            For j = 1 To DataSize
                OldCenter(i, j) = CenterByHandle(i, j)
            Next
        Next
    End If

    
    Do
        Fcm.Iterations = Fcm.Iterations + 1
        For i = 1 To Cluster
            For j = 1 To DataNumber
                Sum1 = 0: Sum3 = 1
                For k = 1 To DataSize
                    Temp = Data(j, k) - OldCenter(i, k)
                    Sum1 = Sum1 + Temp * Temp             ‘計(jì)算第j點(diǎn)到第i個(gè)聚類中心的距離
                Next
                If Sum1 = 0 Then
                    Fcm.Degree(i, j) = 1                      ‘如果j點(diǎn)與第i個(gè)聚類中心重合,，則完全屬于該類
                Else
                    For k = 1 To Cluster
                        Sum2 = 0
                        If k <> i Then
                            For l = 1 To DataSize
                                Temp = Data(j, l) - OldCenter(k, l)
                                Sum2 = Sum2 + Temp * Temp  ‘計(jì)算第j點(diǎn)到其他聚類中心的距離
                            Next
                            Sum3 = Sum3 + Exp(Log(Sum1 / Sum2) * (2 / (Exponent - 1)))      ‘計(jì)算累加值,
                        End If
                    Next
                    Fcm.Degree(i, j) = 1 / Sum3    ‘計(jì)算新的隸屬度
                End If
            Next
        Next
       
        For i = 1 To Cluster
            For j = 1 To DataSize
                Sum1 = 0: Sum2 = 0
                For k = 1 To DataNumber
                    Temp = Exp(Log(Fcm.Degree(i, k)) * Exponent)
                    Sum1 = Sum1 + Temp * Data(k, j)           ‘隸屬度的平方乘以數(shù)值
                    Sum2 = Sum2 + Temp                        ‘隸屬度的和
                Next
                Fcm.Center(i, j) = Sum1 / Sum2                    ‘得到新的聚類中心
            Next
        Next
       
        Temp = 0
        For i = 1 To Cluster
            For j = 1 To DataSize
                Temp = Temp + (OldCenter(i, j) - Fcm.Center(i, j)) ^ 2      ‘ 計(jì)算兩次迭代之間的聚類中心的距離
                OldCenter(i, j) = Fcm.Center(i, j)                          ‘ 保留上一次的聚類中心
            Next
        Next

    Loop While Fcm.Iterations < Maxiterations And Temp > MinImprovement
   
    If AntiFuzzy = Max Then
        For i = 1 To DataNumber
            Temp = -1
            For k = 1 To Cluster
                If Temp < Fcm.Degree(k, i) Then    ‘得到列方向的最大值
                    Temp = Fcm.Degree(k, i)
                    Index = k
                End If
            Next
            Fcm.Class(i) = Index                  ‘Index記錄了列方向最大隸屬度的類
        Next
    ElseIf AntiFuzzy = Mean Then
         For i = 1 To DataNumber
             Temp = 0
             For j = 1 To Cluster
                Temp = Temp + Fcm.Degree(j, i) * j   ‘去隸書乘以對(duì)應(yīng)的類別數(shù)之和
             Next
             Fcm.Class(i) = CInt(Temp)
      Next
    ElseIf AntiFuzzy = Middle Then
        For i = 1 To DataNumber
            Temp = 0
            For j = 1 To Cluster
                If Temp <= 0.5 And Temp + Fcm.Degree(j, i) >= 0.5 Then
                    Index = j
                    Exit For
                Else
                    Temp = Temp + Fcm.Degree(j, i)   ‘取面積的一半處
                End If
            Next
            Fcm.Class(i) = Index
        Next
    End If
    Fcm.TimeUse = GetTickCount - Fcm.TimeUse
    Exit Function
ErrHandle:
    Fcm.ErrMsg = Err.Description
    Fcm.TimeUse = GetTickCount - Fcm.TimeUse
End Function

‘*************************************************************************************
‘*    作    者 ：    laviewpbt
‘*    函 數(shù) 名 ：    Hcm
‘*    參    數(shù) ：    Data     -   待分類的樣本，第一維的大小表示樣本的個(gè)數(shù),，
‘*                                第二維的大小表示樣本的維數(shù)
‘*                   Cluster  -   分類數(shù)
‘*                   Maxiterations  - 最大的迭代次數(shù)
‘                    MinImprovement - 最小的改進(jìn)參數(shù)（兩次迭代間聚類中心的距離）
‘*    返回值 ：      HcmInfo結(jié)構(gòu),，記錄了相關(guān)的參數(shù)
‘*    功能描述 ：    直接利用硬聚類方法把數(shù)據(jù)分類
‘*    日    期 ：    2004-10-24 20.10.56
‘*    修 改 人 ：    laviewpbt
‘*    日    期 ：    2006-11-7 20.11.23
‘*    版    本 ：    Version 2.3.1
‘**************************************************************************************

Private Function Hcm(ByRef Data() As Double, ByVal Cluster As Byte, Optional Maxiterations As Long = 400, Optional MinImprovement As Double = 0.01) As HcmInfo
    If ArrayRange(Data) <> 2 Then
        Hcm.ErrMsg = "數(shù)據(jù)只能為二維數(shù)組"
        Exit Function
    End If
    Dim i As Long, j As Long, k As Long, l As Long, m As Long
    Dim DataNumber As Long, DataSize As Long
    Dim Temp As Double, DX As Double, DY As Double, SumValue() As Double, SumNumber() As Long
    Dim OldCenter() As Double, Distance As Double, Dist As Double, Index As Long
    On Error GoTo ErrHandle:
    Hcm.TimeUse = GetTickCount
    DataNumber = UBound(Data, 1): DataSize = UBound(Data, 2)
    ReDim Hcm.Center(1 To Cluster, 1 To DataSize) As Double
    ReDim Hcm.Class(1 To DataNumber) As Byte
    ReDim OldCenter(1 To Cluster, 1 To DataSize) As Double
    For i = 1 To Cluster
        For j = 1 To DataSize
            OldCenter(i, j) = Data(i * DataNumber / Cluster, j) ‘產(chǎn)生初始中心點(diǎn)
        Next
    Next
    Do
        Hcm.Iterations = Hcm.Iterations + 1
        ReDim SumNumber(Cluster) As Long
        ReDim SumValue(Cluster, DataSize) As Double
        For i = 1 To DataNumber
            Distance = 40000000000#
            For j = 1 To Cluster
                Dist = 0
                For k = 1 To DataSize
                    Temp = Data(i, k) - OldCenter(j, k)
                    Dist = Dist + Temp * Temp             ‘計(jì)算第j點(diǎn)到第i個(gè)聚類中心的距離
                Next
                If Distance > Dist Then
                    Distance = Dist
                    Index = j                         ‘把i點(diǎn)歸于距離該點(diǎn)最近的中心點(diǎn)所在的類
                End If
            Next
            Hcm.Class(i) = Index
            For j = 1 To DataSize
                SumValue(Index, j) = SumValue(Index, j) + Data(i, j)
            Next
            SumNumber(Index) = SumNumber(Index) + 1
        Next
       
        For i = 1 To Cluster
            For k = 1 To DataSize
                If SumNumber(i) = 0 Then
                    Hcm.Center(i, k) = Data(Rnd * DataNumber, k)
                Else
                    Hcm.Center(i, k) = SumValue(i, k) / SumNumber(i)         ‘求新的中心
                End If
            Next
        Next
        Temp = 0
        For i = 1 To Cluster
            For j = 1 To DataSize
                Temp = Temp + (OldCenter(i, j) - Hcm.Center(i, j)) ^ 2      ‘ 計(jì)算兩次迭代之間的聚類中心的距離
                OldCenter(i, j) = Hcm.Center(i, j)                          ‘ 保留上一次的聚類中心
            Next
        Next
    Loop While Hcm.Iterations < Maxiterations And Temp > MinImprovement
    Hcm.TimeUse = GetTickCount - Hcm.TimeUse
    Exit Function
ErrHandle:
    Hcm.ErrMsg = Err.Description
    Hcm.TimeUse = GetTickCount - Hcm.TimeUse
    End Function

‘*************************************************************************************
‘*    作    者 ：    網(wǎng)絡(luò)
‘*    函 數(shù) 名 ：    ArrayRange
‘*    參    數(shù) ：    Data     -   待測(cè)試的數(shù)據(jù)
‘*    返回值 ：      返回?cái)?shù)組的維數(shù)
‘*    日    期 ：    2006-7-10 13.20.40
‘*    修 改 人 ：    laviewpbt
‘*    日    期 ：    2006-11-7 10。10,。45
‘*    版    本 ：    Version 1.2.1
‘**************************************************************************************
Public Function ArrayRange(Data() As Double) As Integer
    Dim i As Integer, ret As Integer
    Dim ErrF As Boolean
    ErrF = False
    On Error GoTo ErrHandle
    For i = 1 To 60               ‘VB中數(shù)組最大為60
        ret = UBound(mArray, i)   ‘用UBound函數(shù)判斷某一維的上界,，如果大數(shù)組的實(shí)際維數(shù)時(shí)產(chǎn)生超出范圍錯(cuò)誤，此時(shí)我們通過Resume Next 來(lái)捕捉錯(cuò)這個(gè)錯(cuò)誤
        ret = ret + 1
        If ErrF Then Exit For
    Next
    ArrayRange = ret
    Exit Function
ErrHandle:
    ret = i
    ErrF = True
    Resume Next
End Function

 測(cè)試情況：

1,、簡(jiǎn)單數(shù)據(jù)的聚類

原始數(shù)據(jù)為:
1    2   
2    3   
1.5    2.5   
1.5    2   
5.1    1   
4.1    1   
5    3   
6    2   
聚類中心為:
1.500    2.374   
5.062    1.750   
隸屬矩陣為:
1.00 1.00 1.00 1.00 0.00 0.03 0.02 0.00
0.00 0.00 0.00 0.00 1.00 0.97 0.98 1.00
迭代次數(shù)為:2

以下是一組氣候的統(tǒng)計(jì)數(shù)據(jù),，根據(jù)要求我們把氣候分為三類：

原始數(shù)據(jù)為:
1    3.5    1    0   
2    2.5    2    2   
2    3.5    1    1   
3    3    3    1   
3    3    1    1   
5    .5    5    2   
6    1.5    4    0   
6    1.5    4    1   
5    3    2    2   
4    3    1    2   
聚類中心為:
4.434    2.992    1.600    1.948   
5.740    1.251    4.250    0.931   
1.902    3.280    1.204    0.845   
隸屬矩陣為:
0.09 0.26 0.02 0.40 0.26 0.12 0.08 0.02 0.92 0.88
0.04 0.08 0.00 0.16 0.04 0.82 0.89 0.97 0.04 0.03
0.87 0.67 0.98 0.44 0.70 0.06 0.04 0.01 0.03 0.09
迭代次數(shù)為:4
每個(gè)樣本的類別為：
3 3 3 3 3 2 2 2 1 1

這個(gè)分析的結(jié)果是可以接受的。

2 ,、二維數(shù)據(jù)的聚類

隨機(jī)生成二維的數(shù)據(jù),，然后比較分類結(jié)果

 上圖中第一排第二個(gè)式HCM的結(jié)果,，其他的圖是不同的參數(shù)相結(jié)合的FCM的結(jié)果,，至于結(jié)果的好壞對(duì)于這些點(diǎn)群不好說(shuō)。

  這副圖選擇的聚類類別為4,，中間那4個(gè)點(diǎn)被聚集到左上角的一類，F(xiàn)CM和HCM的結(jié)果是一樣的,。

      上面這副圖的聚類類別為4,，但他顯示出無(wú)論HCM還是FCM都沒有給出合理的分類,，產(chǎn)生這一現(xiàn)象的原因與其說(shuō)是算法的問題,，不如說(shuō)是對(duì)距離的定義，在FCM/HCM中我們使用的是歐式距離,，這使得中間那一部分到左側(cè)的距離要近一些,，因此被聚集到這個(gè)類中,。通過適當(dāng)改變對(duì)距離的定義，我們也可以得到合理的結(jié)果,。

3 ,、 圖像分割中的應(yīng)用,。

圖像分割也是數(shù)據(jù)聚類的過程，我們選擇圖像的原始象素矩陣為聚類樣本,。效果如下：

        第一行第二列是HCM的結(jié)果,，其他是不同參數(shù)FCM的結(jié)果,，其中第四副是采用加權(quán)平均反模糊化的，看到這副效果是比較差的,，邊緣有很多毛刺,，這是因?yàn)榧訖?quán)平均的過程有個(gè)取整的操作所造成的,。

      這副圖片是在工業(yè)現(xiàn)場(chǎng)拍攝的一副火焰圖像，我以前做一個(gè)項(xiàng)目用到的,，分割后的圖片根據(jù)一定的原理可以得到一定溫度 ,，并且能大概的不同燃料和其他物份的分布狀況,。

    在常用的二值化及其vb.net實(shí)現(xiàn)一文中,，我曾提到用模糊聚類的方法二值化圖像，應(yīng)該說(shuō)這種方法有其合理之處,，他不同于其他的閥值分割,，二帶有一定模糊信息,，并且能方便實(shí)現(xiàn)彩色二值分割（即樣本的維數(shù)為3）。

       如上圖,，程序完美的把天空綠草和樹木分開,。

      當(dāng)然FCM在其他方面也有著廣泛的應(yīng)用,。

      以上是最基本的模糊聚類算法,，針對(duì)FCM的計(jì)算速度慢的特點(diǎn)，人們提出了快速FCM的步驟,，不過那玩意我沒看懂是什么意思,。 

      前面我說(shuō)過遺傳算法，遺傳算法在優(yōu)化方面有很大的優(yōu)勢(shì),，而FCM算法實(shí)際上就是找到目標(biāo)函數(shù)的最小值，因此可以把兩者結(jié)合起來(lái)從而獲得更好的效果,。

最后談一點(diǎn)速度優(yōu)化方面的問題：

1,、FCM算法的計(jì)算量是比較大的，在算法模型不變的情況下我們可以通過以下方法減小計(jì)算量.   a.如果對(duì)某一個(gè)特定的問題,，我們知道聚類中心的大概位置,，則通過程序提供的CreateByHand方法運(yùn)行。 b,、在不確定聚類中心的情況下,，選擇聚類中心由HCM產(chǎn)生，HCM算法的速度是相當(dāng)塊的,。  c.針對(duì)不同對(duì)象選擇不同數(shù)據(jù)類型,，這點(diǎn)在下面要講到。

2,、如果處理對(duì)象是圖像,，則數(shù)據(jù)量一般很大，速度就是關(guān)鍵了,?？紤]到圖像數(shù)據(jù)是byte類型的，則可以把FCM算法的Data()參數(shù)數(shù)據(jù)類型改為byte,我們知道浮點(diǎn)數(shù)的運(yùn)算總是很慢的,。并且考慮到象素值沒有小數(shù)部分,，程序中有些/可以改為\，整除總比一般除法塊,，還有既然確定了樣本的第二維,，則把程序中所有的 DataSize改為3，并且對(duì)于所有的有關(guān)DataSize的小循環(huán),，全部改為手寫。

3,、實(shí)踐證明，RndDegree(i, k)^Exponent 的計(jì)算速度比 Exp(Log(RndDegree(i, k)) * Exponent)要慢,。X^2 比X*X要慢,，我是指大數(shù)據(jù)量的。