Canny edge algorithmus und Sobel Matrix

[Gutelo]

Hallo,

Ich implementiere gerade ein Verfahren um Kannten in Bildern zu erkennen. Ablauf: Farbbild -> Graustufenbild -> Glaetten mit Gauss Matrix Kern -> Gradienten mittels Sobel Matrizen in x und y richtung > threshold -> non maximum reduction

Ein Problem habe ich bei den Sobel Matrizen. Die Faltung mit der Pixelmatrix mit den sobel matrizen liefert neben byte konformen werten auch negative werte. Im internet habe ich nirgends gefunden wie genau die beiden Einzelbilder aus den zwei Sobelmatrizen generiert werden. Diese Bilder werden oft separat gezeigt. Ausserdem zeigen meine errechneten Bilder keine klaren Kannten. Ich nehme an dass nach Anwendung der Sobelmatrizen der Farbraum neu skaliert wird und thresholds gesetzt werden. Ist das richtig? Wenn ja, wieso schreibt das kein Mensch. Auch steht ueberall dass die errechneten pixel normiert werden aber der normierungsfaktor fehlt generell bei sobel matrizen. Bei der gaussmatrix hingegen steht der triviale normierungsfaktor immer dabei.

Gutelo

[Gutelo]

Als Sobel Matrizen verwende ich:

Code:

           (1   2  1)
G_Y = N_Y *(0   0  0)
           (-1 -2 -1)

            (1  0  -1)
G_X = N_X * (2  0  -2)
            (1  0  -1)

Die Normierungsfaktoren N_Y und N_X werden selten angegeben. Zwei Werte N_Y und N_X habe ich fuer exakt diese Matrizen gesehen. Einmal N_Y=N_X=(1/4). Dies macht wenig Sinn da dadurch die Werte der Faltung zwischen -255 und +255 liegen und sich nicht als Graubild darstellen lassen. Mehr Sinn macht N_Y=N_X=(1/8). In diesem Fall wuerden die Werte der Faltung zwischen -128 und +128 liegen und durch ein positives Verschieben aller Werte um 128 liesse sich das Gradientenbild erstellen. Ob ein zusaetzlicher Threshold notwendig ist zur Darstellung der Gradientenplots konnte ich immer noch nicht rausfinden.

Wieso gibt es da nicht wenigstens irgendein vernuenftiges Tutorial fuer dieses Verfahren?

[Aphton]

Es ist klar, dass negative Werte rauskommen können.. Wenn man die Matrizen genauer ansieht, transformiert man von [0..255] auf [-1020..1020]
Dh du müsstest einen Normalisierungsschritt machen, indem du zuerst 1020 addierst, dann mit 255/2040 multiplizierst.

Bin mir momentan nicht ganz sicher.. Meine Überlegung war - wenn die "obere Pixellinie" die mit der ersten Zeile der G_Y Matrix multipliziert wird, aus 255'en bestehen, dann kommt dort als max. Wert 1*255 + 2*255 + 1*255 = 4*255 = 1020 raus.. Dasselbe gilt für die untere Pixellinie = -1020

Edit: Nen Sonderfall hat man bei Rändern. Ich hab ne Beispielsimplementierung, wo ich die Operatoren zähle, damit ich den neuen Intervall jeweils für X und Y kenne und somit dann normalisieren kann..

Delphi-Quellcode:

procedure gray(bmp: TBitmap);
var
  sl: PRGBTriple;
  y, x: Integer;
begin
  for y := 0 to bmp.Height - 1 do
  begin
    sl := bmp.ScanLine[y];
    for x := 0 to bmp.Width - 1 do
    begin
      sl.rgbtBlue   := (sl.rgbtBlue + sl.rgbtGreen + sl.rgbtRed) div 3;
      sl.rgbtGreen  := sl.rgbtBlue;
      sl.rgbtRed    := sl.rgbtBlue;
      Inc(sl);
    end;
  end;
end;

procedure sobel(bmp: TBitmap);
const
  SobelXOperator: Array[-1..1,-1..1] of ShortInt =
    ((-1,0,1),
     (-2,0,2),
     (-1,0,1));
  SobelYOperator: Array[-1..1,-1..1] of ShortInt =
    ((1,2,1),
     (0,0,0),
     (-1,-2,-1));
var
  Data: Array of Array of Byte;
  y, x: Integer;
  i, j: Integer;
  posXOp, negXOp: Integer;  // positive & negative Sobel X Operatoren
  posYOp, negYOp: Integer;  // -- || --                  Y Operatoren
  sumX, sumY: Integer;
  sl: PRGBTriple;
begin
  SetLength(Data, bmp.Height, bmp.Width);
  for y := 0 to bmp.Height - 1 do
  begin
    sl := bmp.ScanLine[y];
    for x := 0 to bmp.Width - 1 do
    begin
      Data[y,x] := sl.rgbtBlue;
      Inc(sl);
    end;
  end;

  for y := 0 to bmp.Height - 1 do
  begin
    sl := bmp.ScanLine[y];
    for x := 0 to bmp.Width - 1 do
    begin
      sumX    := 0;
      sumY    := 0;
      posXOp  := 0;
      negXOp  := 0;
      posYOp  := 0;
      negYOp  := 0;
      for i := -1 to 1 do
        if (x + i >= 0) and (x + i < bmp.Width) then
          for j := -1 to 1 do
            if (y + j >= 0) and (y + j < bmp.Height) then
            begin
              if SobelXOperator[j,i] < 0 then
                inc(negXOp, abs(SobelXOperator[j,i]))
              else
                inc(posXOp, abs(SobelXOperator[j,i]));

              if SobelYOperator[j,i] < 0 then
                inc(negYOp, abs(SobelYOperator[j,i]))
              else
                inc(posYOp, abs(SobelYOperator[j,i]));

              inc(sumX, Data[y+j,x+i] * SobelXOperator[j,i]);
              inc(sumY, Data[y+j,x+i] * SobelYOperator[j,i]);
            end;
      // Normalisierungsschritt
      inc(sumX, negXOp*255);
      sumX := sumX div (negXOp + posXOp);

      inc(sumY, negYOp*255);
      sumY := sumY div (negYOp + posYOp);

      sumX := (sumX + sumY) div 2;
      // .. done

      sl.rgbtBlue   := sumX;
      sl.rgbtGreen  := sumX;
      sl.rgbtRed    := sumX;

      inc(sl);
    end;
  end;
end;

[Medium]

Ich habe exakt genau diese Schritte als Basis meiner Bachelor Arbeit benutzt (Canny-Edge-Detector lässt grüßen.)

Die beiden Richtungen des Sobel-Operators werden dabei als Magnitude gerechnet, also die Länge der Vektoren, die sich aus den Grauwerten der beiden Ergbnisbilder an einer Koordinate ergeben: Sqrt(Sobel_X(x, y)² + Sobel_Y(x, y)²). Damit ist das Problem der negativen Werte gleich mit erschlagen.

Um einen Skalierungsfaktor habe ich mich in meiner Arbeit nicht geschert, da dies eh nur ein Zwischenschritt ist, und man die Non-Maxima-Unterdrückung locker ohne Normierung machen kann (bzw. sollte, weil schneller).
Wenn es denn aber sein muss, schaut man sich einfach mal den extremen Fall an:

Code:

 1  1 1|0             1 0 -1      1  2  1
 1  0  0  mit Kernlen 2 0 -2 und   0  0  0
1|0 0  0              1 0 -1     -1 -2 -1

Sobel_X und _Y sind hierbei = 3, deren Magnitude also rund 4,243. Das wäre der höchte Wert, auf den anzupassen wäre. Da dieser Fall aber ein starkes Extrem ist, zeigt ein so angepasstes als Graustufen angezeigtes Bild fast nur nahezu schwarze Stellen. Ich hatte daher um mir das Zwischenergebnis anzuzeigen eine nicht-lineare Skalierung genommen (Magnitude^K / 4,234 mit 0<K<1). Aber dies hat wie gesagt für die Non-Maxima-Unterdrückung keine Relevanz.

Wichtiger ist für diese, dass du neben der Magnitude der Sobel-Operatoren auch deren Richtungen beachtest, weil in diese Richtungen "gesucht" werden muss. Die Richtung ergibt sich aus Atan2(Sobel_Y(x, y), Sobel_X(x, y)), was einen Winkel ergibt. Den genauen Winkel braucht man allerdings nicht, es genügt in die bei Pixeln 8 möglichen Richtungen zu klassifizieren: 0, 45, 90, 135, 180, 225, 270 und 315 Grad.

Generell sind die Deutschen und Englischen Wikipedia-Artikel zum Canny-Edge-Detector sehr informativ, und diese sowie die darin angegeben Quellen (sowie eine Hand voll Google) haben nach etwas Arbeit ganz gut zum Ziel geführt: Listen von 4-Nachbarschafts-Koordinaten von zusammenhängenden Kanten, die ich darauf hin reduziert und zu Catmull-Rom-Splines verwertet habe. (Der Weg zu diesen war erheblich steiniger als der CED selbst, auwei das waren Zeiten )

[Gutelo]

Hallo Medium,

danke fuer die Antwort. Dass man aus den x und y gradienten das endgueltige bild mittels G=sqrt(Gx*Gx + Gy*Gy) berechnet ist mir klar. Es geht mir um eine vernuenftige Darstellung der einzelnen Bilder Gx und Gy und da muessen alle Werte positiv sein.

Ich habe immer noch das Problem, dass sowohl die einzelnen Bilder als auch das endgueltige Sobel Bild keinerlei Kannten zeigen. Wahrscheinlich habe ich irgendwo einen Fehler in der Implementierung. Hier mein Code:

Delphi-Quellcode:

// ##############################################################
// Prozedur um Graubild zu erzeugen
// ##############################################################

Procedure ConvertToGrayScale(ImgIn : TBitmap; var ImgOut : TBitmap);
type TRGB = Array[1..3] of Byte;
     PixArray = Array[1..MaxInt div SizeOf(TRGB)-1] of TRGB;
     PPixArray = ^PixArray;
     PixArray2D = Array of PPixArray;
var ImgIn2D : PixArray2D;     // Access to ImgIn Pixels
    ImgOut2D : PixArray2D;    // Access to ImgOut Pixels
    rows,cols : Integer;
    GrayShade : Byte;
begin

   // Build 2D Arrays via 'scanline' to speed up pixel operations
   // access pixels: Img[a]^[b][c] , a=line, b=col, c: 3=R,2=G,or 1=B
   // Input Image:
   SetLength(ImgIn2D,ImgIn.Height);
   For rows := 0 to ImgIn.Height-1 do ImgIn2D[rows] := ImgIn.ScanLine[rows];
   // Output Image:
   SetLength(ImgOut2D,ImgOut.Height);
   For rows := 0 to ImgOut.Height-1 do ImgOut2D[rows] := ImgOut.ScanLine[rows];

   // Convert to grayscale
   for rows := 0 to ImgIn.Height-1 do
   for cols := 0 to ImgIn.Width-1 do
    begin
      GrayShade := Round(  0.3 * ImgIn2D[rows]^[cols][3]
                         + 0.6 * ImgIn2D[rows]^[cols][2]
                         + 0.1 * ImgIn2D[rows]^[cols][1]);
      ImgOut2D[rows]^[cols][1] := GrayShade; // B-value
      ImgOut2D[rows]^[cols][2] := GrayShade; // G-value
      ImgOut2D[rows]^[cols][3] := GrayShade; // R-value
    end;
end; 


// ##############################################################
// Prozeduren um Filter Matrizen zu erzeugen
// ##############################################################

Procedure MakeSmoothenFilter(var Filter : TFilter);
var i : Integer;
begin
   // Set Dimensions of Filter Array
   SetLength(Filter,5);
   for i := 0 to High(Filter) do SetLength(Filter[i],5);
   // Define Gauss Filter Matrix
   //   2  4  5  4 2
   //   4  9 12  9 4
   //   5 12 15 12 5
   //   4  9 12  9 4
   //   2  4  5  4 2
   // first Line
   Filter[0,0] := 2; Filter[0,1] := 4; Filter[0,2] := 5;
   Filter[0,3] := 4; Filter[0,4] := 2;
   // second Line
   Filter[1,0] := 4; Filter[1,1] := 9; Filter[1,2] := 12;
   Filter[1,3] := 9; Filter[1,4] := 4;
   // third Line
   Filter[2,0] := 5;  Filter[2,1] := 12; Filter[2,2] := 15;
   Filter[2,3] := 12; Filter[2,4] := 5;
   // fourth Line
   Filter[3,0] := 4; Filter[3,1] := 9; Filter[3,2] := 12;
   Filter[3,3] := 9; Filter[3,4] := 4;
   // fifth Line
   Filter[4,0] := 2; Filter[4,1] := 4; Filter[4,2] := 5;
   Filter[4,3] := 4; Filter[4,4] := 2;
end;

Procedure MakeSobelFilterX(var Filter : TFilter);
var i : integer;
begin
   // Set Dimensions of Filter Array
   SetLength(Filter,3);
   for i := 0 to High(Filter) do SetLength(Filter[i],3);
   // Define Sobel-X Filter Matrix
   //   -1  0  1
   //   -2  0  2
   //   -1  0  1
   // first Line
   Filter[0,0] := -1; Filter[0,1] := 0; Filter[0,2] := 1;
   // second Line
   Filter[1,0] := -2; Filter[1,1] := 0; Filter[1,2] := 2;
   // third Line
   Filter[2,0] := -1; Filter[2,1] := 0; Filter[2,2] := 1;
end;

Procedure MakeSobelFilterY(var Filter : TFilter);
var i : integer;
begin
   // Set Dimensions of Filter Array
   SetLength(Filter,3);
   for i := 0 to High(Filter) do SetLength(Filter[i],3);
   // Define Sobel-Y Filter Matrix
   //    1   2   1
   //    0   0   0
   //   -1  -2  -1
   // first Line
   Filter[0,0] :=  1; Filter[0,1] :=  2; Filter[0,2] :=  1;
   // second Line
   Filter[1,0] :=  0; Filter[1,1] :=  0; Filter[1,2] :=  0;
   // third Line
   Filter[2,0] := -1; Filter[2,1] := -2; Filter[2,2] := -1;
end;


// ##############################################################
// Prozedur um Filter auf ein Bild anzuwenden (Faltung)
// ##############################################################

Procedure ApplyFilterOnImage(var ImgIn, ImgOut : TBitmap; Filter : TFilter);
type TRGB = Array[1..3] of Byte;
     PixArray = Array[1..MaxInt div SizeOf(TRGB)-1] of TRGB;
     PPixArray = ^PixArray;
     PixArray2D = Array of PPixArray;
var rows,cols : Integer;      // rows and cols in the image
    FLowX, FHighX : Integer;
    FLowY, FHighY : Integer;
    weight : Integer;         // Weight for Filter Matrix
    ImgIn2D : PixArray2D;     // Access to ImgIn Pixels
    ImgOut2D : PixArray2D;    // Access to ImgOut Pixels
    SumMinus : Integer;
    ScaleShift : Integer;
    PixelSum, x, y : Integer;
    Pixel : Integer;
begin
   // Filter matrix, element -> pixel position with respect to center pixel
   // Calculate extremes in x and y direction
   FLowX  := -Round(0.5 * Length(Filter)-1);
   FHighX := +Round(0.5 * Length(Filter)-1);
   FLowY  := -Round(0.5 * Length(Filter[0])-1);
   FHighY := +Round(0.5 * Length(Filter[0])-1);

   // Determine matrix weight and scaleshift required after convolution:
   ScaleShift := 0;
   SumMinus := 0;
   weight := 0;
   For y := 0 to High(Filter[0]) do
   For x := 0 to High(Filter) do
   begin
      if Filter[x,y] < 0 then SumMinus := SumMinus + Filter[x,y];
      weight := weight + Abs(Filter[x,y]);
   end;
   ScaleShift := Abs(Round((SumMinus/weight)*255));

   // Build 2D Arrays via 'scanline' to speed up pixel operations
   // access pixels: Img[a]^[b][c] , a=line, b=col, c: 3=R,2=G,or 1=B
   // Input Image:
   SetLength(ImgIn2D,ImgIn.Height);
   For rows := 0 to ImgIn.Height-1 do ImgIn2D[rows] := ImgIn.ScanLine[rows];
   // Output Image:
   SetLength(ImgOut2D,ImgOut.Height);
   For rows := 0 to ImgOut.Height-1 do ImgOut2D[rows] := ImgOut.ScanLine[rows];

   // Calculate filtered Image, i.e. convolution of Filter and ImgIn -> ImgOut
   For rows := 0-FLowY to ImgIn.Height-1-FHighY do
   For cols := 0-FLowX to ImgIn.Width-1-FHighX do
   begin
      // Perform Convolution
      PixelSum := 0;
      For y := FLowY to FHighY do
      For x := FLowX to FHighX do
      begin
        PixelSum := PixelSum +
            Filter[x-FLowX,y-FLowY]* ImgIn2D[rows+y]^[cols+x][1];
      end;
      Pixel := Round(PixelSum/weight) + ScaleShift;
      ImgOut2D[rows]^[cols][1] := Pixel;
      ImgOut2D[rows]^[cols][2] := Pixel;
      ImgOut2D[rows]^[cols][3] := Pixel;
   end;
end; 

// ##############################################################
// Prozedur um finales Sobelbild aus Einzelbildern zu erzeugen
// ##############################################################

Procedure GenerateFinalSobel(ImgSx, ImgSy : TBitmap; var ImgOut : TBitmap);
type TRGB = Array[1..3] of Byte;
     PixArray = Array[1..MaxInt div SizeOf(TRGB)-1] of TRGB;
     PPixArray = ^PixArray;
     PixArray2D = Array of PPixArray;
var  rows,cols : Integer;      // rows and cols in the image
     ImgSx2D  : PixArray2D;    // Access to ImgSx Pixels
     ImgSy2D  : PixArray2D;    // Access to ImgSy Pixels
     ImgOut2D : PixArray2D;    // Access to ImgOut Pixels
     Gx,Gy,Gf : Byte;
begin
   // Build 2D Arrays via 'scanline' to speed up pixel operations
   // access pixels: Img[a]^[b][c] , a=line, b=col, c: 3=R,2=G,or 1=B
   // SobelX Image:
   SetLength(ImgSx2D,ImgSx.Height);
   For rows := 0 to ImgSx.Height-1 do ImgSx2D[rows] := ImgSx.ScanLine[rows];
   // SobelY Image:
   SetLength(ImgSy2D,ImgSy.Height);
   For rows := 0 to ImgSy.Height-1 do ImgSy2D[rows] := ImgSy.ScanLine[rows];
   // Output Image:
   SetLength(ImgOut2D,ImgOut.Height);
   For rows := 0 to ImgOut.Height-1 do ImgOut2D[rows] := ImgOut.ScanLine[rows];

   // Perform G = sqrt(Gx*Gx + Gy*Gy);
   for rows := 0 to ImgSx.Height-1 do
   for cols := 0 to ImgSx.Width-1 do
   begin
     Gx := ImgSx2D[rows]^[cols][1];
     Gy := ImgSy2D[rows]^[cols][1];
     Gf := Round(sqrt(Abs(Gx*Gx)+Abs(Gy*Gy)));
     ImgOut2D[rows]^[cols][1] := Gf;
     ImgOut2D[rows]^[cols][2] := Gf;
     ImgOut2D[rows]^[cols][3] := Gf;
   end;
end;  

// ##############################################################
// Prozedur um geglaettetes Graubild zu erzeugen
// ##############################################################

Procedure Smoothen(var ImgIn, ImgOut : TBitmap);
var GaussFilter : TFilter;
    SobelXFilter : TFilter;
    ImgGray : TBitmap;
begin
   //convert to Grayscale
   ImgGray := TBitmap.Create;
   ImgGray.Width := ImgIn.Width;
   ImgGray.Height := ImgIn.Height;
   ConvertToGrayScale(ImgIn, ImgGray);
   //smoothen grayscale image
   MakeSmoothenFilter(GaussFilter);
   ApplyFilterOnImage(ImgGray, ImgOut, GaussFilter);
end;

// ##############################################################
// Prozedur fuer kompletten Sobel
// ##############################################################

Procedure Sobel(var ImgIn, ImgOut : TBitmap);
var SobelX, SobelY : TFilter;
    ImgSmooth, ImgSx, ImgSy : TBitmap;
begin
    // create ImgSmooth
   ImgSmooth := TBitmap.create;
   ImgSmooth.width := ImgIn.width;
   ImgSmooth.height := ImgIn.height;
   // Smoothen Image
   Smoothen(ImgIn, ImgSmooth);
   // create ImgSx
   ImgSx := TBitmap.create;
   ImgSx.width := ImgIn.width;
   ImgSx.height := ImgIn.height;
   // create ImgSy
   ImgSy := TBitmap.create;
   ImgSy.width := ImgIn.width;
   ImgSy.height := ImgIn.height;
   // Make SobelX and SobelY Filter
   MakeSobelFilterX(SobelX);
   MakeSobelFilterY(SobelY);
   // Apply Sobel Filter on ImgIn -> ImgSx and ImgSy
   ApplyFilterOnImage(ImgSmooth, ImgSx, SobelX);
   ApplyFilterOnImage(ImgSmooth, ImgSy, SobelY);
   // ImgSx and ImgSy -> ImgOut
   GenerateFinalSobel(ImgSx, ImgSy, ImgOut);
end;

Hier eine Prozedur zum Aufruf des Ganzen:

Delphi-Quellcode:

procedure TForm1.Button2Click(Sender: TObject);
var ImgIn, ImgOut : TBitmap;
begin
  ImgIn := TBitmap.Create;
  ImgOut := TBitmap.Create;
  ImgIn.LoadFromFile('Book.bmp');
  ImgOut.Width := ImgIn.Width;
  ImgOut.Height := ImgIn.Height;
  //Smoothen(ImgIn, ImgOut);
  Sobel(ImgIn, ImgOut);
  Image3.Picture.Bitmap := ImgOut;
end;

Die Umwandlung in ein Graubild und das Glaetten mittels Gaussfilter klappt wunderbar. Die einzelnen Sobelbilder sehen nicht aus wie Bilder in denen horizontale oder vertikale Kannten betont sind sondern eher wie helle oder dunkle Bilder mit stark veraendertem Gammawert. Im finalen Sobelbild ist quasi nichts zu sehen, d.h. mehr oder weniger einfarbige Flaeche.

Gutelo

Edit: Code ist mit Lazarus erstellt, sollte aber 1:1 auf Delphi laufen

[jfheins]

Hast du mal ein Eingabebild und die zwei Ausgaben, die dein Programm produziert?