Canny Konturenfilter
Zum Code
Dies ist meine experimentelle Implementation des Canny-Algorithmus. Dieser berechnet die Konturen in Bildern, welche dann später in Ketten umgewandelt und durch Polygone annähernd beschrieben werden könnten. Das Canny-Verfahren ist eins der besten und exaktesten, es liefert sehr gute Ergebnisse.

Für weitere Informationen sollte man mal nach Canny googeln. Hier ein kleiner Abriss des Ablaufs:

Zuerst berechnet man zwei Gauss-Filter, einmal einen horizontalen und einen vertikalen. Diese werden auf das Eingabebild (welches vorher in Graustufen umgewandelt wurde) angewandet, wodurch wir die Ableitung des Bildes in X- und in Y-Richtung erhalten. Daraus wird ein Norm- Bild errechnet, wobei gilt Norm^2 = X'^2 + Y'^2. Dann setzt das Beispiel alle Werte auf ein gewisses Minimum, Abschnitt "Thresholding", um im letzten Schritte die erkannten Konturen auf die lokalen Maxima zu reduzieren. Dafür werden für eventuelle Konturpunkte die theoretischen Werte vor und nach ihm auf der Kante interpoliert (bilinear) und überprüft, ob der aktuelle Wert größer ist als seine interpolierten Nachbar. Beim Filtern der Maxima wird das Thresholding wieder rückgängig gemacht.

History
15.05.2003 Hinzugefügt

Autor: Dominik Auras <Dominik_auf_vbinside.de> 

Code aus Form1.frm
Option Explicit
Const pi = 3.14159265
Private Sub Command1_Click()
  Dim filterx() As Double, filtery() As Double, img() As Double
  Dim img2() As Double, img3() As Double, norm_img() As Double
  Dim min As Double, max As Double, level As Double, alfa As Double
  Dim n1 As Long, n2 As Long, sigma1 As Double, sigma2 As Double
  Dim thres_img() As Double, thin_img() As Double
    
    alfa = 0.1
    n1 = 10: n2 = 10
    sigma1 = 1
    sigma2 = 1
    
    ReadGreyImg img
    Debug.Print "Ausgelesen"
    
    d2dgauss n1, sigma1, n2, sigma2, 0.5 * pi, filterx
    Debug.Print "Convolution Kernel x erstellt"
    
    d2dgauss n1, sigma1, n2, sigma2, 1 * pi, filtery
    Debug.Print "Convolution Kernel y erstellt"
    
    ApplyConvKernel img, filterx, img2
    Debug.Print "Kernel x angewendet"
    
    ApplyConvKernel img, filtery, img3
    Debug.Print "Kernel y angewendet"
    
    CalcNorm img2, img3, norm_img
    Debug.Print "Normbild berechnet"
    
    max = FindMax(norm_img)
    min = FindMin(norm_img)
    level = alfa * (max - min) + min
    
    Threshold norm_img, level, thres_img
    Debug.Print "Thresholding abgeschlossen"
    
    Thinning thres_img, img2, img3, norm_img, min, max, level, thin_img
    Debug.Print "Thinning abgeschlossen"
    
    ShowGreyImg_Forced thin_img, 255, 0
End Sub
Private Sub Form_Load()
    Picture1.Picture = LoadPicture(App.Path & "\001.gif")
End Sub
Sub Thinning(img() As Double, img2() As Double, img3() As Double, _
    norm_img() As Double, min As Double, max As Double, level As _
    Double, out() As Double)
  Dim x As Long, y As Long, x1 As Double, x2 As Double
  Dim y1 As Double, y2 As Double, n3 As Double
  Dim z3 As Double
    
    ReDim out(UBound(img, 1), UBound(img, 2))
    
    For x = 1 To UBound(img, 1) - 1
      For y = 1 To UBound(img, 2) - 1
        If img(x, y) > level Then
          x1 = img2(x, y) / norm_img(x, y)
          y1 = img3(x, y) / norm_img(x, y)
          x2 = -x1: y2 = -y1
          
          Debug.Assert x1 <= 1 And x1 >= -1 And y1 <= 1 And y1 >= -1
          Debug.Assert x2 <= 1 And x2 >= -1 And y2 <= 1 And y2 >= -1
          
          n3 = interpol(img, x, y, x1, y1)
          z3 = interpol(img, x, y, x2, y2)
          
          'Die Punkte (x1,y1) = n3 und (x2,y2) = z3 müssen _
              interpoliert werden
          ' Richtung: (x2,y2) -> (x,y) -> (x1,y1)
          
          If img(x, y) >= n3 And img(x, y) >= z3 Then
            out(x, y) = max
          Else
            out(x, y) = min
          End If
        Else
          out(x, y) = min
        End If
      Next y
    Next x
    
    'interp2
    'berechnet die farbwerte an den punkten
    ' img2(i,j)/norm_img(i,j), img3(i,j)/norm_img(i,j)
    ' -img2(i,j)/norm_img(i,j), -img3(i,j)/norm_img(i,j)
    
    ' norm_img(i,j) = sqr(img2(i,j)^2 + img3(i,j)^2)
    ' betrag der änderung
    
    ' img2(i,j)/norm_img(i,j) == änderung x / änderung gesamt
    ' -> normalisierung auf 0 bis 1
    
    ' X =
    '  -1 0 1
    '  -1 0 1
    '  -1 0 1
    '---------
    ' Y =
    '  -1 -1 -1
    '   0  0  0
    '   1  1  1
    '---------
    ' Z = entsprechender Ausschnitt in img(x-1 bis x+1, y-1 bis y+1)
    
    ' (-1|-1|Z)  (0|-1|Z)  (1|-1|Z)
    '                  (ZI1)
    ' (-1| 0|Z)  (0| 0|Z)  (1| 0|Z)
    '        (ZI2)
    ' (-1| 1|Z)  (0| 1|Z)  (1| 1|Z)
    '
    ' ZI1 = (änderung x / änderung gesamt, änderung y / änderung gesamt)
    ' ZI1 = (-änderung x / änderung gesamt, -änderung y / änderung _
        gesamt)
    '
    ' Wir haben nur eine Kante, wenn laut Interpolation die Mitte _
        das Maximum ist
    ' wir betrachten praktisch die Punkte auf der Kante, die vor _
        und nach dem
    ' aktuellen Pixel liegen
    '
    ' bilineare Interpolation
    '  1. Quadrant bestimmen
    '  2. die 4 Randwerte auslesen
    '  3. Abstand berechnen
    '     Koords = Abstand von (0|0|Z), muss Abstand von links oben sein
    '     1. Quadrant (links oben): Abstand von (-1|-1|Z)
    '     2. Quadrant (rechts oben): Abstand von (0|-1|Z)
    '     3. Quadrant (rechts unten): Abstand von (0|0|Z)
    '     4. Quadrant (links unten): Abstand von (-1|0|Z)
    '  4. Interpolieren, lineare Interpolation für beide Zeilen, _
        dann lineare
    '     Interpolation für die interpolierten Werte der Zeilen
    '
    'Gesucht sind die interpolierten Werte ZI1 und ZI2
    '
    ' Wenn der originale Wert >= dem ersten und zweiten interpolierten
    ' Wert ist ==> max
    ' sonst ==> min
End Sub
Function interpol(img() As Double, x As Long, y As Long, x1 As _
    Double, y1 As Double) As Double
  Dim n1 As Double, n2 As Double, n3 As Double
    
    If x1 < 0 And y1 > 0 Then
      '1. Quadrant
      n1 = Bilinear(1 - x1, img(x - 1, y - 1), x1, img(x, y - 1))
      n2 = Bilinear(1 - x1, img(x - 1, y), x1, img(x, y))
      n3 = Bilinear(1 - y1, n1, y1, n2)
      'Punkt(x1,y1) = n3
    End If
    
    If x1 > 0 And y1 > 0 Then
      '2. Quadrant
      n1 = Bilinear(x1, img(x, y - 1), 1 - x1, img(x + 1, y - 1))
      n2 = Bilinear(x1, img(x, y), 1 - x1, img(x + 1, y))
      n3 = Bilinear(1 - y1, n1, y1, n2)
      'Punkt(x1,y1) = n3
    End If
    
    If x1 > 0 And y1 < 0 Then
      '3. Quadrant
      n1 = Bilinear(x1, img(x, y), 1 - x1, img(x + 1, y))
      n2 = Bilinear(x1, img(x, y + 1), 1 - x1, img(x + 1, y + 1))
      n3 = Bilinear(y1, n1, 1 - y1, n2)
      'Punkt(x1,y1) = n3
    End If
    
    If x1 < 0 And y1 < 0 Then
      '4. Quadrant
      n1 = Bilinear(1 - x1, img(x - 1, y), x1, img(x, y))
      n2 = Bilinear(1 - x1, img(x - 1, y + 1), x1, img(x, y + 1))
      n3 = Bilinear(y1, n1, 1 - y1, n2)
      'Punkt(x1,y1) = n3
    End If
    
    interpol = n3
End Function
Function Bilinear(ByVal d0 As Double, ByVal f0 As Double, ByVal d1 _
    As Double, ByVal f1 As Double) As Double
    Bilinear = f0 * d0 + f1 * d1
End Function
Sub Threshold(img() As Double, max As Double, out() As Double)
  Dim x As Long, y As Long
    
    ReDim out(UBound(img, 1), UBound(img, 2))
    
    For x = 0 To UBound(img, 1)
      For y = 0 To UBound(img, 2)
        If img(x, y) > max Then
          out(x, y) = img(x, y)
        Else
          out(x, y) = max
        End If
      Next y
    Next x
End Sub
Function FindMax(inp() As Double) As Double
  Dim x As Long, y As Long, max As Double
    
    max = inp(0, 0)
    
    For x = 0 To UBound(inp, 1)
      For y = 0 To UBound(inp, 2)
        If inp(x, y) > max Then
          max = inp(x, y)
        End If
      Next y
    Next x
    
    FindMax = max
End Function
Function FindMin(inp() As Double) As Double
  Dim x As Long, y As Long, min As Double
    
    min = inp(0, 0)
    
    For x = 0 To UBound(inp, 1)
      For y = 0 To UBound(inp, 2)
        If inp(x, y) < min Then
          min = inp(x, y)
        End If
      Next y
    Next x
    
    FindMin = min
End Function
Sub CalcNorm(in1() As Double, in2() As Double, out() As Double)
  Dim x As Long, y As Long
    
    ReDim out(UBound(in1, 1), UBound(in1, 2))
    
    For x = 0 To UBound(in1, 1)
      For y = 0 To UBound(in1, 2)
        out(x, y) = Sqr(in1(x, y) ^ 2 + in2(x, y) ^ 2)
      Next y
    Next x
End Sub
Sub ShowGreyImg_Enforced(img() As Double)
  Dim x As Long, y As Long, color As Long
    
    Picture1.Cls
    
    For x = 0 To UBound(img, 1)
      For y = 0 To UBound(img, 2)
        color = 128 + Int(img(x, y) * 128 / 400)
        If color > 255 Then
          color = 255: Stop
        End If
        If color < 0 Then
          color = 0: Stop
        End If
        
        Picture1.PSet (x, y), RGB(color, color, color)
      Next y
    Next x
End Sub
Sub ShowGreyImg(img() As Double)
  Dim x As Long, y As Long, color As Long
    
    Picture1.Cls
    
    For x = 0 To UBound(img, 1)
      For y = 0 To UBound(img, 2)
        color = Int(img(x, y))
        If color > 255 Then
          color = 255
        End If
        If color < 0 Then
          color = 0
        End If
        
        Picture1.PSet (x, y), RGB(color, color, color)
      Next y
    Next x
End Sub
Sub ShowGreyImg_Forced(img() As Double, rate As Double, min As Double)
  Dim x As Long, y As Long, color As Long
    
    Picture1.Cls
    
    For x = 0 To UBound(img, 1)
      For y = 0 To UBound(img, 2)
        color = Int(img(x, y) * rate) + min
        If color > 255 Then
          color = 255
        End If
        If color < 0 Then
          color = 0
        End If
        
        Picture1.PSet (x, y), RGB(color, color, color)
      Next y
    Next x
End Sub
Sub ApplyConvKernel(img() As Double, k() As Double, out() As Double)
  Dim midx As Long, midy As Long, x As Long, y As Long
  Dim kx As Long, ky As Long, cx As Long, cy As Long
  Dim c As Double
    
    midx = Int((UBound(k, 1) + 1) / 2)
    midy = Int((UBound(k, 2) + 1) / 2)
    
    ReDim out(UBound(img, 1), UBound(img, 2))
    
    'der erste punkt der berechnet werden kann ist midx,midy
    'alles links und rechts davon ist schwarz
    
    For x = 0 To UBound(img, 1)
      For y = 0 To UBound(img, 2)
        For kx = 0 To UBound(k, 1)
          For ky = 0 To UBound(k, 2)
            cx = x - midx + kx
            cy = y - midy + ky
            
            If cx < 0 Or cy < 0 Or cx > UBound(img, 1) Or cy > _
                UBound(img, 2) Then
              c = 0
            Else
              c = img(cx, cy)
            End If
            
            out(x, y) = out(x, y) + k(kx, ky) * c
          Next ky
        Next kx
      Next y
    Next x
End Sub
Sub ReadGreyImg(img() As Double)
  Dim x As Long, y As Long, R As Long, G As Long, B As Long
  Dim color As Long, lum As Double
    
    ReDim img(Picture1.ScaleWidth - 1, Picture1.ScaleHeight - 1)
    For x = 0 To Picture1.ScaleWidth - 1
      For y = 0 To Picture1.ScaleHeight - 1
        color = Picture1.Point(x, y)
        R = color And &HFF
        G = (color \ 2 ^ 8) And &HFF
        B = (color \ 2 ^ 16) And &HFF
        
        lum = 0.299 * R + 0.587 * G + 0.114 * B
        
        img(x, y) = lum / 255
      Next y
    Next x
End Sub
Sub d2dgauss(n1 As Long, sigma1 As Double, n2 As Long, sigma2 As _
    Double, theta As Double, h() As Double)
  Dim i As Long, j As Long, u1 As Double, u2 As Double, x() As Double
  Dim n As Double, cos_theta As Double, sin_theta As Double
    
    If theta = pi / 2 Then
      cos_theta = 0
      sin_theta = 1
    Else
      cos_theta = Cos(theta)
      sin_theta = Sin(theta)
    End If
    
    ReDim h(n2 - 1, n1 - 1)
    
    For i = 1 To n2
      For j = 1 To n1
        
        u1 = cos_theta * (j - (n1 + 1) / 2) - sin_theta * (i - (n2 _
            + 1) / 2)
        u2 = sin_theta * (j - (n1 + 1) / 2) + cos_theta * (i - (n2 _
            + 1) / 2)
        
        h(i - 1, j - 1) = gauss(u1, sigma1) * dgauss(u2, sigma2)
      Next j
    Next i
    
    AbsMultArray h, h, x
    n = Sqr(SumXY(x))
    DivideArray1 h, n, x
    
    h = x
End Sub
Function dgauss(x As Double, std As Double) As Double
    dgauss = -x * gauss(x, std) / std ^ 2
End Function
Function gauss(x As Double, std As Double) As Double
    gauss = Exp(-x ^ 2 / (2 * std ^ 2)) / (std * Sqr(2 * pi))
End Function
Function SumXY(arr() As Double) As Double
  Dim x As Long, y As Long, sum As Double
    
    For y = LBound(arr, 1) To UBound(arr, 2)
      For x = LBound(arr, 2) To UBound(arr, 1)
        sum = sum + arr(x, y)
      Next x
    Next y
    
    SumXY = sum
End Function
Sub AbsMultArray(inp() As Double, mult() As Double, out() As Double)
  Dim x As Long, y As Long
    
    ReDim out(UBound(inp, 1), UBound(inp, 2))
    
    For x = LBound(inp, 1) To UBound(inp, 1)
      For y = LBound(inp, 2) To UBound(inp, 2)
        out(x, y) = Abs(inp(x, y)) * Abs(mult(x, y))
      Next y
    Next x
End Sub
Sub DivideArray1(inp() As Double, divi As Double, out() As Double)
  Dim x As Long, y As Long
    
    ReDim out(UBound(inp, 1), UBound(inp, 2))
    
    For x = LBound(inp, 1) To UBound(inp, 1)
      For y = LBound(inp, 2) To UBound(inp, 2)
        out(x, y) = inp(x, y) / divi
      Next y
    Next x
End Sub