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 :D)

[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?

[Gutelo]

Hallo,

Ich habe die Bilder an dieses Post gehaengt

Gutelo

[Medium]

Zwar sehen deine "gesobelten" Bilder in der Tat etwas komisch aus (eher wie kontrastverminderte und teilinvertierte Versionen des Originals), aber ich sehe ein großes Problem bei deiner Vorgehensweise: Du nutzt als Zwischenspeicher für deine Schritte immer ein TBitmap. Du brauchst aber auch die negativen Werte der Sobel-Komponenten um nachher die korrekten Magnituden und Richtungen zu berechnen, da kommst du nicht drum rum!

Ich hatte für mein Programm damals eine eigenen FloatBitmap Klasse gemacht, und habe diese für sämtliche Operationen genutzt - ausser Laden und Speichern, dafür konnte diese Klasse ein .CreateFromBitmap() und ein .ToBitmap() durchführen. Ich habe die Bilder darin, also zumindest die, die zur Darstellung genutzt wurden, auf die Range 0..1 skaliert, und alle Farbkomponenten als Single hinterlegt. (Hatte auch den Hintergrund, dass ich das nachher so umgebaut habe, dass Gauss und Sobel in einem Shader von der Graka gemacht wurden, wo man das eh so braucht.)
Für die Farben hatte ich eine eigene FloatColor Klasse, die Methoden wie .ToYCC() und .ToRGB() oder .ToColor() gleich mit gebracht hat. Das FloatBitmap hat davon einfach ein 2D Array verwaltet, und ein paar Methoden zum Umwandeln eben.

Nutzt man sowas als "Mittler" zwischen den Schritten, sind Negative Farben und denormalisierte "Bilder" kein Problem mehr (streng genommen sind es ja schon fast keine Bilder mehr dann). Das würde ich als aller erstes machen, und dann mal schauen wie es danach aussieht.

Sowohl die Normaierung auf 0..1, als auch die deutlich erhöhte Genauigkeit durch Floats sind gerade bei den Zwischenschritten zum einen eine Vereinfachung, und das Zweite je nach Anwendungsfall schon fast eine Notwendigkeit finde ich. Wenn die Geschwindigkeit nicht hin haut, könnte man auch noch mit auf Integer-Range skalierten Farben arbeiten, muss dann aber aufpassen mit Ausreissern in den Zwischenschritten die wie beim Sobel dann auch mal über "1" (bzw. MaxInt) liegen können. Und es ist eine ziemliche Rumrechnerei. Wenn wirklich sehr schnelle Verarbeitung sein soll, würde ich letztlich wieder zum Shader greifen. Da das aber eine ganze Ecke Infrastruktur und potenziell Anpassungen im gesamten Programm braucht wenn man es nicht von Anfang an darauf ausgelegt hat, will ich das nicht uneingeschränkt als Heiland preisen.

Aber deine negativen Werte brauchst du. Definitiv. Wenn du die Zwischenschritte anzeigen willst, würde ich dafür das FloatBitmap von (-1..1) auf (0..255) mappen, und Werte unter/über -1/1 clampen. Dann bleibt genügend "Kontrast" um vernünftig was sehen zu können. (Diese Bilder taugen dann aber nicht mehr wirklich zum weiter damit arbeiten, nur zur Kontrollanzeige.)

[Gutelo]

Hallo Medium

Stimmt, ich sollte erstmal die Minuswerte beibehalten.

Was meinst du mit "würde ich dafür das FloatBitmap von (-1..1) auf (0..255) mappen, und Werte unter/über -1/1 clampen". Wenn ich nicht normiere dann varrieren die Werte von minimal -4*256 bis maximal 4*256 in den beiden "Teilbildern". Im naechsten Schritt - sqrt(Gx*Gx + Gy*Gy) - sollten aehnliche Werte rauskommen. Am Ende muesste ich doch diese Bereiche auf (0..255) mappen oder nicht?

Also ich rechne jetzt mal ohne Normierungsfaktor fuer die Sobel Matrizen die Werte fuer beide Sobel operatoren aus und anschliessend das finale Sobel Ergebnis. Danach bestimme ich die maximalen negativen und positiven Werte. Dann skaliere ich, so dass alle Werte zwischen 0 und 255 liegen. Richtig?

Gutelo

[Medium]

Zitat:

Zitat von Gutelo (Beitrag 1266757)

Wenn ich nicht normiere dann varrieren die Werte von minimal -4*256 bis maximal 4*256 in den beiden "Teilbildern".

Ja, deswegen meinte ich ja, dass die so gemappten Bilder nur zum Angucken taugen. Werte jenseits von -1 und 1 sind recht seltene Extrema in Fotos, und wenn man etwa -4 bis 4 auf 0 bis 255 mapped wird das Bild überwiegend wie 50% Grau aussehen. Geht natürlich auch, fand ich zum zwischen drin mal beschauen aber nicht so prima.

Zitat:

Im naechsten Schritt - sqrt(Gx*Gx + Gy*Gy) - sollten aehnliche Werte rauskommen. Am Ende muesste ich doch diese Bereiche auf (0..255) mappen oder nicht?

Nur, wenn du das Gradienten-"Bild" anzeigen willst. Letztlich willst du doch eigentlich nur deine Kanten haben, oder? Der Sobel-Schritt wird also am Ende nie angezeigt. Warum sollte man ihn dann auf Bitmap-Werte umrechnen? Einfach so lassen wie es aus der Filterung kommt, und damit die Non-Maxima-Unterdrückung machen. Was davon übrig bleibt dann als S/W Bild nehmen (Kanten weiss, Rest schwarz) oder wie auch immer weiter damit arbeiten. Aber löse dich davon, dass die Zwischenschritte "Bilder" sind. Man kann sie zur Anzeige als solche interpretieren und daraufhin anpassen, aber an sich sind das nur noch große Matrizen mit "irgendwelchen" Werten. Deswegen plädiere ich auch dafür, dass man das TBitmap als Container für diese nicht nimmt.

Zitat:

Also ich rechne jetzt mal ohne Normierungsfaktor fuer die Sobel Matrizen die Werte fuer beide Sobel operatoren aus und anschliessend das finale Sobel Ergebnis. Danach bestimme ich die maximalen negativen und positiven Werte. Dann skaliere ich, so dass alle Werte zwischen 0 und 255 liegen. Richtig?

Jain, s.o. :)

[Gutelo]

Soderle, hab es jetzt mal umgeschrieben. Funktioniert soweit so gut -> Siehe angehaengte Bilder.
Als naechstes nehme ich NMS, Hysteresis und Hough in Angriff.


Falls jemand etwas Aehnliches machen moechte, hier mein Code bis hierhin. An einigen Stellen sicher noch etwas suboptimal und verbesserungsfaehig (z.B. dass das gesammte Array extra Durchlaufen wird um die minimalen und maximalen Werte zu bestimmen.)

Delphi-Quellcode:

unit Canny;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, Graphics,  LCLIntf, Dialogs, Math;

type // type to store a filter
     TFilter = Array of Array of Integer;
     // types for fast access to TBitmap image pixels
     TRGB = Array[1..3] of Byte;
     PixArray = Array[1..MaxInt div SizeOf(TRGB)-1] of TRGB;
     PPixArray = ^PixArray;
     PixArray2D = Array of PPixArray;
     // type to store the pixels of grayscale TBitmaps in an Array
     // and to perform filter operations on the array values
     TPixelArray = Array of Array of integer;


// Export
Procedure ConvertToGrayScale(ImgIn : TBitmap; var ImgOut : TBitmap);
Procedure CalcSobel(ImgIn : TBitmap; var ImgOut : TBitmap);

implementation

// ############################################################################
// Convert Color Image to Grayscale Image
// ############################################################################

Procedure ConvertToGrayScale(ImgIn : TBitmap; var ImgOut : TBitmap);
//
var ImgIn2D : PixArray2D;     // Access to ImgIn Pixels
    ImgOut2D : PixArray2D;    // Access to ImgOut Pixels
    rows,cols : Integer;      // rows and cols of an image
    GrayShade : Byte;         // Gray value
//
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;

// ############################################################################
// Make Gauss Filter Matrix
// ############################################################################

Procedure MakeGaussFilter(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;

// ############################################################################
// Make Sobel-X Filter Matrix
// ############################################################################

Procedure MakeSobelXFilter(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;

// ############################################################################
// Make Sobel-Y Filter Matrix
// ############################################################################

Procedure MakeSobelYFilter(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;

// ############################################################################
// Convert "Image -> Pixel Array" and "Pixel Array -> Image"
// ############################################################################

Procedure ImgToPixArray(ImgIn : TBitmap; var PixArr : TPixelArray);
//
var H,W : Integer;
    ImgIn2D : PixArray2D;
    rows, cols : Integer;
//
begin
  // Store height and width of input image
  H := ImgIn.Height;
  W := ImgIn.Width;
  // Fast access to input image:
  SetLength(ImgIn2D,H);
  For rows := 0 to H-1 do ImgIn2D[rows] := ImgIn.ScanLine[rows];
  // Set dimensions of pixel array
  SetLength(PixArr,H);
  For rows := 0 to H-1 do SetLength(PixArr[rows],W);
  // Copy pixel values of input image to pixel array
  For rows := 0 to H-1 do
  For cols := 0 to W-1 do
  begin
    PixArr[rows,cols] := ImgIn2D[rows]^[cols][1];
  end;
end;

Procedure PixArrayToImg(PixArr : TPixelArray; var ImgOut : TBitmap);
//
var H,W : Integer;
    ImgOut2D : PixArray2D;
    rows, cols : Integer;
    MinVal, MaxVal : integer;
    Pixel : Integer;
//
begin
   // Store height and width of input image
   H := Length(PixArr);
   W := Length(PixArr[0]);
   // Initialize output image
   ImgOut := TBitmap.Create;
   ImgOut.Width := W;
   ImgOut.Height := H;
   // Fast access to output Image:
   SetLength(ImgOut2D,H);
   For rows := 0 to H-1 do ImgOut2D[rows] := ImgOut.ScanLine[rows];
   // Determine minimum and maximum values in pixel array
   MinVal := PixArr[0,0];
   MaxVal := PixArr[0,0];
   For rows := 0 to H-1 do
   For cols := 0 to W-1 do
   begin
      if PixArr[rows,cols] < MinVal then MinVal := PixArr[rows,cols];
      if PixArr[rows,cols] > MaxVal then MaxVal := PixArr[rows,cols];
   end;
   // Rescale value range to 0..255 and copy values to output image
   For rows := 0 to H-1 do
   For cols := 0 to W-1 do
   begin
     // only scale if necessary
     if (MinVal < 0) or (MaxVal > 255)
     then
         Pixel := Round(255*((PixArr[rows,cols] - MinVal)/ (MaxVal-MinVal)))
     else
         Pixel := PixArr[rows,cols];
     // Copy values to output image
     ImgOut2D[rows]^[cols][1] := Pixel;
     ImgOut2D[rows]^[cols][2] := Pixel;
     ImgOut2D[rows]^[cols][3] := Pixel;
   end;
end;

// ############################################################################
// Apply Filter on Pixel Array
// ############################################################################

Procedure ApplyFilterOnPixArray(PA_In : TPixelArray; var PA_Out : TPixelArray; Filter : TFilter; DoWeight : Boolean);
//
var H,W : Integer;        // Height and width of ImgIn
    rows,cols : Integer;  // rows and cols in the image
    weight : Integer;     // Weight for Filter Matrix
    FilterBT : Integer;   // Border Thickness of Filter around central element
    x, y : Integer;       // Variables for Convolution
    PixelSum : Integer;   // Convolution Value
    Pixel : Integer;      // Final Pixel Value
//
begin
   // Border thickness of filter around central matrix element:
   // Examples: 3x3 Matrix -> BT=1 , 5x5 Matrix -> BT=3
   FilterBT := Round(0.5 * (Length(Filter)-1));
   // Store height and width of input pixel array (PA_In):
   H := Length(PA_In);
   W := Length(PA_In[0]);
   // Set dimensions of output pixel array (PA_Out), "removes Filter Border":
   SetLength(PA_Out, H-(2*FilterBT));
   For rows := 0 to High(PA_Out) do SetLength(PA_Out[rows],W-(2*FilterBT));
   // Determine filter matrix weight:
   weight := 0;
   For y := 0 to High(Filter[0]) do
   For x := 0 to High(Filter) do
   begin
      weight := weight + Abs(Filter[x,y]);
   end;
   // Calculate filtered Image, i.e. perform convolution of Filter and PA_In -> PA_Out
   For rows := FilterBT to High(PA_In)-FilterBT do
   For cols := FilterBT to High(PA_In[0])-FilterBT do
   begin
      // Perform convolution
      PixelSum := 0;
      For y := -FilterBT to FilterBT do
      For x := -FilterBT to FilterBT do
      begin
        PixelSum := PixelSum + Filter[x+FilterBT,y+FilterBT]* PA_In[rows+y,cols+x];
      end;
      // Apply weight when required
      if DoWeight then Pixel := Round(PixelSum/weight)
                  else Pixel := PixelSum;
      // write output pixel array
      PA_Out[rows-FilterBT,cols-FilterBT] := Pixel;
   end;
end;

// ############################################################################
// Make final Sobel pixel array based on Sobel-X and Sobel-Y
// ############################################################################

Procedure MakeSobel(PA_Sx, PA_Sy : TPixelArray; var PA_Out : TPixelArray);
//
var rows, cols : Integer;
    Sx, Sy : Integer;
//
begin
   // Set dimensions of output pixel array
   SetLength(PA_Out,Length(PA_Sx));
   For rows := 0 to High(PA_Out) do SetLength(PA_Out[rows], Length(PA_Sx[0]));
   // calculate output pixel array:
   For rows := 0 to High(PA_Sx) do
   For cols := 0 to High(PA_Sx[0]) do
   begin
      Sx := PA_Sx[rows,cols];
      Sy := PA_Sy[rows,cols];
      // write output pixel array
      PA_Out[rows,cols] := Round(sqrt(Sx*Sx + Sy*Sy));
   end;
end;

// ############################################################################
// Procedure to perform complete Sobel calculation
// ############################################################################

Procedure CalcSobel(ImgIn : TBitmap; var ImgOut : TBitmap);
var ImgGray : TBitmap;
    PA_Gray, PA_Gauss, PA_SobelX, PA_SobelY, PA_Sobel : TPixelArray;
    GaussFilter, SobelX, SobelY : TFilter;
begin
  ImgGray := TBitmap.Create;
  ImgGray.Width := ImgIn.Width;
  ImgGray.Height := ImgIn.Height;
  // Convert color input image (ImgIn) to grayscale image (ImgGray)
  // and convert ImgGray to pixel array format (PA_Gray)
  ConvertToGrayScale(ImgIn,ImgGray);
  ImgToPixArray(ImgGray, PA_Gray);
  // Apply Gauss smoothening filter to PA_Gray -> PA_Gauss
  MakeGaussFilter(GaussFilter);
  ApplyFilterOnPixArray(PA_Gray, PA_Gauss, GaussFilter, true);
  // Apply SobelX filter to PA_Gauss -> PA_SobelX
  MakeSobelXFilter(SobelX);
  ApplyFilterOnPixArray(PA_Gauss, PA_SobelX, SobelX, false);
  // Apply SobelY filter to PA_Gauss -> PA_SobelY
  MakeSobelYFilter(SobelY);
  ApplyFilterOnPixArray(PA_Gauss, PA_SobelY, SobelY, false);
  // Combine PA_SobelX and PA_SobelY -> PA_Sobel
  MakeSobel(PA_SobelX, PA_SobelY, PA_Sobel);
  // Convert PA_Sobel to output image (ImgOut)
  PixArrayToImg(PA_Sobel, ImgOut);
end;

end.

Verwendung aus Hauptprogramm:

Delphi-Quellcode:

procedure TForm1.Button1Click(Sender: TObject);
var ImgIn, ImgOut : TBitmap;
begin
  ImgIn := TBitmap.Create;
  ImgOut := TBitmap.Create;
  ImgIn.LoadFromFile('Test.bmp');
  CalcSobel(ImgIn, ImgOut);
  Image3.Picture.Bitmap := ImgOut;
end;

Gutelo

[Medium]

Das sieht doch schon sehr richtig aus! Nett.

Darf man fragen, was du mit den Kanten nachher vor hast? So von Nachtschwärmer zu Nachtschwärmer :drunken:

[Gutelo]

:) Jepp

Ich schreib ein Programm welches folgendes machen soll:

1.) Mittels Kamera ein Dokument auf einem Tisch aufnehmen (Webcam oder Handy)
2.) Perspektivische Verzerrung rausrechnen, so dass quasi eine Aufnahme "von oben" erhalten wird
3.) Einlesen von Zahlen mittels OCR

1. und 2. sind soweit fertig.

Den Canny+Hough brauche ich um die Eckpunkte des Dokuments im verzerrten Bild automatisch zu finden.

Fuer die OCR werde ich wohl ein einfaches Pixel-matching verwenden: Einmal perfekt mit Flachbettscanner einscannen und fuer jede Zahl von 0 bis 9 eine Referenz anfertigen. Dann die gescannten Zahlen ueber die Referenzen durchschieben und den maximalen Match bestimmen. Das sollte eigentlich ganz gut funktionieren wenn man Hoehe/Breite (Pixelaufloesung) einigermassen waehlt. Vielleicht kommen auch Neuronale Netze in Frage, obwohl mich die Treffsicherheit der simplen OCR Implementationen mittels NN die ich im Internet gefunden habe nicht wirklich ueberzeugt haben.

Eigentlich kam mir bei der ganzen Canny Geschichte folgendes in den Sinn bezueglich der OCR: Man koennte doch ueber jede gescannte Zahl den Canny laufen lassen und anschliessend die Gradientenrichtungen (Statistik ueber Gradienten in 0,45,90,.. Grad Richtung) der Kannten benutzen um die Zahl zu bestimmen. Was haelst du von dieser Idee? Koennte man auch mit dem Pixel-matching kombinieren.

Gutelo

Warum das in der urpruenglichen Version nicht funktioniert hat ist mir aber immer noch raetselhaft, da ich dort im Grunde ja das gleiche mache bis zu den Sx und Sy. Ich glaube aber dass es eventuell an der Normierung liegt. Wenn ich bei der aktuellen Version bereits die Pixel Arrays der Sx und Sy normiere, dann ist der Effekt nicht ganz so stark ausgepraegt als wenn ich erst am Ende nach dem Sobel normiere. Ein anderes Problem der alten Version ist der nicht behandelte Rand. Dieser verfaelscht die Skala. In der neuen Version werden die Raender abgeschnitten nach der Filteranwendung.

[Medium]

Das klingt verdammt spannend! Und beeindruckend oben drauf.
Bei OCR bin ich nicht so sehr im Thema, aber die Idee das anhand der "Winkellisten" zu vergleichen finde ich zumindest schon mal höchst kreativ - und gar nicht mal so abwegig. Da könnte es dann spannend werden einen Vergleichsmechanismus zu finden, der 2 Listen potenziell unterschiedlicher Länge ein geeignetes Ähnlichkeitsmaß zuordnet. Das wirklich interessante bei dem Ansatz ist vor allem, dass man im Prinzip ohne ROI auskommen kann. Einfach alle Gradientenlisten durchgehen und gut. Fies wird es da allerdings wohl bei Zahlen, die in mehere Listen zerfallen können. Also die beiden, die gekreuzte Linien haben können: 4 und 8. Da müsste man sich was überlegen.

Zum Sobel bzw. dessen Nicht-Funktion vorher: Da ich nicht so ganz durch deine Normierungsmethode steige gerade, weiss ich auch nicht genau was da schief lief. Grundsätzlich bin ich aber wegen solchen potenziellen Problemen eher dazu geneigt, die Container den Daten anzupassen als umgekehrt :cyclops:

[Gutelo]

Oh ich dachte eigentlich nur an ein Haufigkeitsdiagram, also Häufigkeit über Richtung. Aber deine "Liste" ist noch viel besser. Man nimmt einfach die äusserste Kante die zusammenhängend die Zahl umschliesst und läuft die Kannte Pixel fûr Pixel ab und macht ein Array mit den Gradientrichtungen. Die Sprûnge in den Richtungen und die Länge der Liste sind charakteristisch für die Zahl. Probleme mit 4 und 8 gibt es dann auch nicht. Gute Idee.

Wofür steht ROI?

Edit: ah. Region of interest. Das Problem auch die Zahlen zu separieren liegt nicht an der Kreuzung sondern dass 4 und 8 mehrere geschlossene Kanten haben. Aber die Zahlen zu separieren ist nicht so wild, das hab ich schon fertig.

[Medium]

Ha, in dieser Form hatte ich an die Listen jetzt auch nicht gedacht :D Aber das klingt gut, im Grunde eine Hüllkurvenauswertung. Könnte nur bei 6, 8 und 9 je nach Font knifflig werden, aber ausprobierenswert. Vor allem wenn du die Separierung schon hast - das ist der wichtigste Schritt überhaupt. Alles danach lässt sich dann ja recht einfach ausprobieren und verfeinern.

[Gutelo]

Mit 8 und 6/9 sollte es kein Problrm geben. 6 und 9 haben etwa gleiche Aenderungen. Hier muss man dann einen festen Startpunkt wählen, z.B. immer beim nördlichsten Punkt anfangen.

Ich muss erst noch die NMS, Hysteresis und Hough fertig machen dann kümmer ich mich um die OCR idee. Halte dich auf dem Laufenden

[Gutelo]

Hallo,

Non-Maximum Suppression (NMS) und Hysteresis arbeiten jetzt auch. Ich haenge ein kleines Demo-Projekt fuer Lazarus an dieses Post. Laeuft alles schnell genug, kann aber auch noch reichlich optimiert werden. Soll nur als ein Grundgeruest dienen fuer Leute die Aehnliches vor haben. Ich habe den Code reichlich kommentiert. Das Binary ist leider elendig gross geworden (28mb). Lazarus scheint noch mehr aufzublaehen als Delphi. Ferner nicht wundern wenn nach dem Verlassen des kompilierten Programms der Debugger meckert. Ein leidiges Lazarus Problem bei Verwendung von OpenDialog. Am besten einmal kompilieren und dann das Program ausserhalb der IDE aufrufen.

WICHTIG: >>> geht NUR mit *.bmp bitmaps, keine jpg, keine gif, ... <<<

Naechste Schritte:

1) Ich muss alle kurzen Kanten entfernen und nur die Langen ueberleben lassen, d.h. Kannten reduzieren.
2) Ich brauche Hough um Geraden und Ecken zu identifizieren.

Gutelo

[Gutelo]

Mist gerade bei Hough gemerkt, dass etwas nicht 100% stimmt.

Erstmal muss in Canny.pas der Typ PixArray mit dem Array bei 0 beginnen und nicht bei 1. Bis zum Sobel ist alles okay.

Die NMS Berechnung und/oder Berechnung der Kantenrichtungen sind fehlerhaft. Bei exakt horizontalen Kanten verschwinden die Kanten fast komplett aus dem Bild bei der NMS Berechnung. Woran kann das liegen? An definition von Atan2 ?

Gleiches Problem bei exakt vertikalen Kanten.

Gutelo

[Medium]

Ohne den Code zu deiner NMS (inkl. der Richtungsbestimmungen, falls die ausgelagert sein sollte) zu sehen schwer zu sagen. Ich hatte da eigentlich kein Problem. Du beachtest aber, dass du senkrecht zu den tatsächlichen Gradientenrichtungen vorgehen musst? Das würde mir so ad hoc als schnell gemachter Fehler einfallen.

[Gutelo]

Hallo Medium,

der Code ist in dem Zip-file des vorletzten Posts

Gutelo

[Medium]

Oh, dachte da wäre die NMS noch nicht bei gewesen. Scusi! Werde ich heute und vermutlich auch morgen aber nicht mehr zu kommen leider :?

[Gutelo]

Kein Problem, ich poste es hier auch noch mal:

Delphi-Quellcode:

// ############################################################################
// Build pixel array of gradient directions
// ############################################################################

Procedure MakeEdgeDirections(PA_Sx, PA_Sy : TPixelArray; var PA_Out : TPixelArray);
//
var rows, cols : Integer;
    Th : Double;         // Angle Theta
//
begin
   // Set dimensions of output pixel array
   SetLength(PA_Out, Length(PA_Sx));
   For rows := 0 to High(PA_Out) do SetLength(PA_Out[rows], Length(PA_Sx[0]));
   // calculate output pixel array:
   For rows := 0 to High(PA_Sx) do
   For cols := 0 to High(PA_Sx[0]) do
   begin
      Th := (180/Pi) * ArcTan2(PA_Sy[rows,cols],PA_Sx[rows,cols]);
      // write output pixel array
      PA_Out[rows,cols] := Round(Th);
   end;
end;

Function FMod(x,y : Double) : Double;
begin
   If y = 0 then
   begin
     Showmessage('Error in function FMod: Cannot divide through 0');
     FMod := 0;
   end
   else
   begin
     FMod := x-Int(x/y)*y;
   end;
end;

// ############################################################################
// Perform Non-maximum Suppression (quantisized way)
// ############################################################################

// Symmetry => 180 degree cover the four main directions: 0,45,90,and 135 deg
// Corresponding equivalent directions: 180,225,270,and 315 deg

// Divide the 180 degree region into four quadrant sections
//  (0) to (Pi/8) and (Pi-(Pi/8)) to (Pi)  <= quadrant around 0 (or symmetry equivalent
//                                                               180 deg,respectively)
//  (1*Pi/4) +- Pi/8    <= quadrant around  45 deg
//  (2*Pi/4) +- Pi/8    <= quadrant around  90 deg
//  (3*Pi/4) +- Pi/8    <= quadrant around 135 deg

// Introduce 'dir' value: dir := (FMod(PA_ED[rows,cols]+180, 180) /180) * 8
// Used to classify the member quadrant for a particular gradient direction
//
// How 'dir' works:
// Example: Input angle at point (x,y) : PA_ED[x,y] = -135
//          Shift to get positive angle: PA_ED[x,y] + 180 = 45
//          Value on 'dir' scale: (FMod(45,180)/180)*8= (45/180)*8 = 2
//
// Integer 'dir' values correspond to the following angles:
//         dir = 0  <=> Th =   0           = 0      { center of 1st quadrant
//         dir = 1  <=> Th =   Pi/8        = 1*Pi/8 <= border from 1st to 2nd quadrant
//         dir = 2  <=> Th =   Pi/4        = 2*Pi/8 { center of 2nd quadrant
//         dir = 3  <=> Th =   Pi/4 + Pi/8 = 3*Pi/8 <= border from 2nd to 3rd quadrant
//         dir = 4  <=> Th =   Pi/2        = 4*Pi/8 { center of 3rd quadrant
//         dir = 5  <=> Th =   Pi/2 + Pi/8 = 5*Pi/8 <= border from 3rd to 4th quadrant
//         dir = 6  <=> Th = 3*Pi/4        = 6*Pi/8 { center of 4th quadrant
//         dir = 7  <=> Th = 3*Pi/4 + Pi/8 = 7*Pi/8 <= border from 4th to 1st quadrant
//         dir = 8  <=> Th =   Pi          = 8*Pi/8 { center of 1st quadrant

Procedure MakeNMS(PA_SO, PA_ED : TPixelArray; var PA_Out : TPixelArray);
//
var dir : Double;
    rows, cols : Integer;
//
begin
   // Set dimensions of output pixel array
   SetLength(PA_Out, Length(PA_SO));
   For rows := 0 to High(PA_Out) do SetLength(PA_Out[rows], Length(PA_SO[0]));
   // Calculate output pixel array:
   For rows := 1 to High(PA_SO)-1 do
   For cols := 1 to High(PA_SO[0])-1 do
   begin
      // Calculate 'dir' value to determine the right quadrant
      dir := (FMod(PA_ED[rows,cols]+180, 180) /180) * 8;
      //
      // Check membership in quadrant, keep maxima,or suppress minima
      //
      // 0 degree quadrant (EE <-> WW)
      If  ((dir <= 1) or (dir > 7))
          and (PA_SO[rows,cols] > PA_SO[rows,cols+1])  // EE
          and (PA_SO[rows,cols] > PA_SO[rows,cols-1])  // WW
      then PA_Out[rows,cols] := PA_SO[rows,cols]
      //
      else
      // 45 degree quadrant (NW <-> SE)
      If  ((dir > 1) or (dir <= 3))
          and (PA_SO[rows,cols] > PA_SO[rows-1,cols-1]) // NW
          and (PA_SO[rows,cols] > PA_SO[rows+1,cols+1]) // SE
      then PA_Out[rows,cols] := PA_SO[rows,cols]
      //
      else
      // 90 degree quadrant (NN <-> SS)
      If  ((dir > 3) or (dir <= 5))
          and (PA_SO[rows,cols] > PA_SO[rows-1,cols]) // NN
          and (PA_SO[rows,cols] > PA_SO[rows+1,cols]) // SS
      then PA_Out[rows,cols] := PA_SO[rows,cols]
      //
      else
      // 135 degree quadrant (NE <-> SW)
      If  ((dir > 5) or (dir <= 7))
          and (PA_SO[rows,cols] > PA_SO[rows-1,cols+1]) // NE
          and (PA_SO[rows,cols] > PA_SO[rows+1,cols-1]) // SW
      then PA_Out[rows,cols] := PA_SO[rows,cols]
      //
      else
      // If no maximum -> Suppress Value
          PA_Out[rows,cols] := 0;
   end;
end;

Also ein Copy+Paste Fehler ist schonmal, dass es

If ((dir > 1) and (dir <= 3))
If ((dir > 3) and (dir <= 5))
((dir > 5) and (dir <= 7))

heissen muss. Da es alles verschlimmert stimmt wohl auch irgendwas mit den Richtungen nicht.

[Gutelo]

Medium, bist du dir denn sicher, dass du dieses Problem nicht auch im code hast?

Bei grossen Bildern faellt es nicht besonders auf, vorallem wenn die meisten Kanten verdreht sind gegenueber der horizontalen oder vertikalen Richtung. Probier mal das Testbild im Anhand (Cross.bmp) mit deinem Code aus.

Ich habe im Internet bei einigen Implementationen einen Kommentar gefunden ala "Problem mit vertikalen und horizontalen Kanten korrigiert"

Zitat:

"Du beachtest aber, dass du senkrecht zu den tatsächlichen Gradientenrichtungen vorgehen musst? "

Das verstehe ich nicht. Die Gradientenrichtung steht schon senkrecht auf der Kante und das ist genau die Richtung in die zusaetzliche Pixel entfernt werden muessen.

Gutelo

[Medium]

Das wirklich tragische ist gerade, dass ich mein Projekt nicht mehr finde. Ist ja auch nur meine Bachelorarbeit, die in >1 Jahr als Leidenschaftsprojekt erstellt wurde :pale:. Ich glaube fast, dass mir der liebe Ukash Virus Anfang des Jahres die letzte Sicherung davon geschrottet hat. Ich hoffe meine Eltern haben die CD noch, ich glaube denen hatte ich die gegeben. Wie ärgerlich :(

Edit: Das mit dem senkrecht vergiss übrigens gleich wieder, da hatte ich akuten Hirnknoten. Irgendwas war damit beim Tracing der Kanten in meiner Arbeit, das habe ich wohl vermischt gestern.

[Gutelo]

Oh das ist nicht gut. Hast du in der Uni keine Kopie hinterlassen?

Ich vermute dass der Atan2 probleme macht. Das Richtungsbild sieht etwas komisch aus um die horizontale Kante.

[Gutelo]

Am Atan2 liegt es nicht. Das Problem ist dass man die Aenderungen bei der NMS am Original durchführen muss. Blöder Fehler. Kann man sich auch einfach überlegen: Angenommen es liegen bei einer horizontalen linie zwei maxima senkrecht übereinander. Beide Pixel würden unterdrückt da der Wert nicht grösser ist als die nord und süd pixel. Ändert man jedoch am Original so ist der eine bereits unterdrückt und der zweite wird behalten.

Hab auch noch einen anderen Fehler gefunden der dazu führte dass SobelX und SobelY vertauscht waren. Überarbeite alles und poste den code nochmal

Gutelo

Edit: Die Problematik ist doch etwas komplizierter. Es reicht auch nicht aus die NMS auf dem Originalbild auszufuehren, da dann andere Pixel stehen bleiben die eigentlich unterdrueckt werden muessen. Viele Implementationen im Internet aendern nicht am Original und haben das Problem dass horizontale oder vertikale Linien bei der NMS komplett unterdrueckt werden.

[Medium]

Ui, der fällt aber auch schon in die Kategorie "Muss man erstmal drauf kommen". Dann hat sich das Problem bei mir vermutlich nicht gezeigt, weil es bei Floats deutlich unwahrscheinlicher (wenn auch nicht unmöglich) ist genau gleiche Werte nebeneinander zu bekommen. (Ich habe nämlich auch in einen frischen Puffer gemalt und nicht im Original.) Dazu kommt noch, dass im Endprodukt zwei Mal gesobelt wird, weil ich die "Kanten der Kanten" haben wollte, mit einem kleinen Zwischenschritt (Potenzierung des 1. Gradienten), wodurch recht schmale Strukturen entstanden sind.

Freut mich, dass es jetzt läuft! Werde übers Wochenende sicherlich auch dazu kommen, mir das genauer anzuschauen!

[Gutelo]

Aaaaargs. Naetuerlich liegt es daran dass ich vorher runde. Ich Depp. Wahrscheinlich kann dieser Fall ohne Runden nicht auftreten. Falls doch muesste man zweimal durchlaufen:

1) Lese aus Original und Schreibe in Kopie. Loesche alle Pixel die einen kleineren Wert haben als mindestens einer von beiden Nachbarpixeln in Gradientenrichtung. -> Elimination aller Pixel die garantiert kein Maximum sind

2) Nehme Kopie und loesche dort alle Pixel die nicht groesser sind als ihre beiden Nachbarpixel in Gradientenrichtung. -> Eine breite Kante wird von einer Seite her ausgeduennt bis nur noch die letzte Pixelzeile uebrig bleibt.

Hoechstwahrscheinlich ist das aber nicht noetig wenn man vorher alles bei floats belaesst.

Danke Medium

Edit: Oh, anscheinend tritt das Problem auch gelegentlich auf wenn man mit floats arbeitet. Medium kannst du das mal checken? Einige horizontale Kanten verschwinden vollstaendig, andere ueberleben. Ich werde die oben beschriebene Methode implementieren die zweimal das Sobelbild durchgeht um diesen Effekt zu korrigieren bzw zu vermeiden.

[Medium]

Ich muss gestehen, dass ich das bei mir noch nicht beobachtet habe, kann es aber mangels Code auch nicht prüfen. Meine FH hat leider immer nur zu meinen Arbeitszeiten geöffnet und ist etwas weiter weg :?. Ich muss noch mal alle USB Sticks fein säuberlich durchgehen.
Da ich aber nach einem Sobel bei mir nicht aufhöre, sondern noch etliche Dinge daran anschließen bis zur NMS, müsste ich auch noch ein wenig umbauen. Mir fällt allerdings ein: Ich glaube, dass ich am Ende tatsächlich die Werte im Original auch auf Null gesetzt habe bei der NMS. Oder es war erst beim Tracing. Da bin ich gerade unsicher, ist aber ja nun auch 4-5 Jahre her =). Das war aber, soweit ich mich erinnere, nicht wegen Problemen mit Senk- und Waagerechten. Was genau weiss ich aber leider auch nicht mehr.
Ich schiebe zudem auch noch einen Schritt nach, der eine 8-Nachbarschaft mit exakt 3 Punkten pro 3x3 Fenster herstellt. Dabei könnten eventuell verbleibende Doppelkanten, die blieben wenn man nicht mit ">" sondern ">=" arbeitete, auch noch erodiert werden. Ich hab da so lange so viel verschiedenes ausprobiert... ich muss es wirklich sehen um sinnvoll antworten zu können fürchte ich :?

[Gutelo]

Die ungenauigkeit von float werten bewirkt dass Teile der horizontalen und vertikalen Kanten die NMSüberleben. Aber Teile verschwinden.

[Gutelo]

Das zweistufige Verfahren ist auch murks. Dann sind zwar die horizontalen und vertikalen Werte okay aber er verballhornt die anderen Richtungen. Hingegen funktioniert es gut wenn man für eine Gradientenrichtung >= verwendet und für die andere >. Ich lasse es erstmal so. Hough funktioniert auch prima :) Allerdings ist Hough naturgegeben sau lahm. Werde mir mal die optimierten Hough Algoritmen anschauen.