Regionen erkennen aus einer Matrix

[Peter666]

Hi,

ich hab hier eine Aufgabe über der ich schon seit einiger Zeit grübele:

Als Vorgabe besitze ich eine eindimensionale Matrix (bestehend aus Breite*Höhe):

Map: array[0..256 * 256 - 1] of Boolean;

Die Aufgabe ist nun sämtliche gesetzte Punkte zu zählen und Regionen zuzuordnen. Mein Grundgerüst schaut wie folgt aus:

Code:

procedure FindRegions(const Width, Height: Integer; const Map: array of Boolean; var Regions: array of Integer);
var tmp: integer;
  i, j: integer;
  x, y: integer;
begin
  // Regionen löschen
  Fillchar(Regions, sizeof(regions), 0);

  for y := 0 to Height - 1 do
    for x := 0 to Width - 1 do
    begin
     
    end;

  // Bubblesort von Groß, zu Klein
  for i := 0 to high(regions) - 2 do
    for j := 0 to high(regions) - 2 do
      if (i <> j) and (regions[i] > regions[j]) then
      begin
        tmp := regions[i];
        regions[i] := regions[j];
        regions[j] := tmp;
      end;
end;

Als Ausgabe sollte eine Integer Liste mit der Pixelanzahl der jeweiligen Regionen sein. Meine Idee ist anstelle des Map: array[0..256 * 256 - 1] of Boolean;
ein Map: array[0..256 * 256 - 1] of Integer; zu machen, dem bei der Initialisierung eine fortlaufende Nummer zu geben - der Einfachheit der Offset und die Punkte drumherum zu zählen und ggf. zusammenzuführen, also Punkt links, Punkt Oben, Punkt rechts <> 0, dann setze alle Punkte auf den Wert des ersten Punktes.
Das ganze hat aber Macken, weiß jemand von euch eventuell Rat?

Peter

[Aphton]

Was ist eine "Region"?
Edit: Ist es eine Menge von Punkten, die sich benachbarn?

[Peter666]

Genau alle Punkte die benachbart sind, sollten in eine Region zusammengefasst werden. Ich brauche allerdings nur eine Liste mit wie viele Punkte pro Region zu finden sind.

Peter

[jfheins]

Das grundsätzliche Vorgehen das ich wählen würde sieht so aus:

Wenn du nicht nur die Anzahl der Regionen nund deren Felder haben willst, sondern auch die Zuordnung Punkt => Region dann musst du eine Integer-Array hernehmen. Ich definiere mir jetzt mal folgende Werte:
-1 : Zu zählender Punkt
-2 : Zu ignorierender Punkt
>=0: Punkt gehört zu der Region mit dem Index

Zuerst: Schreib eine rekursive Funktion, die dir die umliegenden Felder abklappert und zu einer Region verbindet. Klingt kompliziert, ist es aber nicht.
Dann: Alle Felder durchgehen und für jedes Feld die Funktion aufrufen.

Zu Punkt 1:

Delphi-Quellcode:

function MakeRegions(const Map: Array of Integer; index, regionindex, width: Integer): Integer;
Result := 0;
if (index >= 0) and (index < length(Map)) then
if Map[index] == -1 then
begin
  Map[index] = regionindex;
  Result := 1;

  Result += MakeRegions(Map, index + 1, regionindex, width);
  Result += MakeRegions(Map, index - 1, regionindex, width);
  Result += MakeRegions(Map, index + width, regionindex, width);
  Result += MakeRegions(Map, index - width, regionindex, width);
end;

Sowas in der Richtung. Ist jetzt ein wenig gemischter Code aber der Gedanke soll klar werden.

Aufruf dann für jeden Punkt:

Delphi-Quellcode:

  for o := 0 to length(Map) - 1 do
    begin
     MakeRegions(Map, i, i, width);
    end;

Falls der Rückgabewert > 0 ist, wurde eine neue Region gefunden, mit der angegebenen Zahl von Feldern.
Das Array sollte nur mit -1 und -2 initialisiert sein. Danach kannst du anhand des Werts im Array die Region ablesen, zu der der Punkt gehört.

[Aphton]

Wie sieht die Region aus? Muss sie rechteckig, quadratisch und noch wichtiger - müssen alle gleich sein (breite, höhe)?
Falls nicht, dann würde ich dir ein Quadtree vorschlagen, bei der du dann versuchen solltest, die Regionen durch das umgebende Quad zu bestimmen.

MfG

[Peter666]

Die ist nicht rechteckig.

@jfheins: So ähnlich hab ich mir das auch schon gedacht - ich probier das mal. Hmm das geht glaub ich nicht wirklich, denn bei nem Array aus 100 Feldern dauert das ja schon knapp ne Minute.

Peter

[BoolString]

Möglicherweise meinst du eine Implementation des Hoshen-Kopelman Algorithmus? Das ist ein recht komplexes Forschungsfeld; zumindest sobald es darum geht mehr als ein paar einzelne Pixel auf einer binären Maske zu erkennen.

Ich finde hier allerdings meine Delphi-Implementierung nicht mehr. Aber es ist nicht endlos kompliziert (Englisch vorausgesetzt). Man kann das wohl in ein-zwei Nachmittagen schaffen...

Hintergrundinformationen zum Hoshen-Kopelman Algorithmus.
Ein Vortrag zur Problematik der 'Musterfindung'.
Der Original-Artikel von Hoshen über den verbesserten Algorithmus von 1997.


Beste Grüße

Jan

[jfheins]

Zitat:

Zitat von Peter666 (Beitrag 1073438)

Die ist nicht rechteckig.

@jfheins: So ähnlich hab ich mir das auch schon gedacht - ich probier das mal. Hmm das geht glaub ich nicht wirklich, denn bei nem Array aus 100 Feldern dauert das ja schon knapp ne Minute.

Peter

Da sollte noch was rauszuholen sein. Wichtig ist, dass das Array nicht ständig umherkopiert wird, also entweder nen Pointer übergeben oder das Array global machen.
Außerdem kannst du den letzten rekursiven Aufruf weglassen, wenn du immer von niedrigen zu hohen indexen iterierst.

[Peter666]

Genial, danke der Algorithmus scheint das wohl lösen zu können :)

[Medium]

Irgendwie erinnert mich das Problem an FloodFilling, nur dass man sich wegschreiben muss, welche Idizes man denn "gefüllt" hat.

[Peter666]

@Medium: Ja mich auch :)

Der HK Algorithmus sollte wie folgt aussehen:

Code:

var labels: array of integer = nil;
  n_labels: integer = 0; (* length of the labels array *)

(*  uf_find returns the canonical label for the equivalence class containing x *)

function uf_find(x: integer): integer;
var y, z: integer;
begin
  y := x;
  while (labels[y] <> y) do
    y := labels[y];

  while (labels[x] <> x) do
  begin
    z := labels[x];
    labels[x] := y;
    x := z;
  end;
  result := y;
end;

(*  uf_union joins two equivalence classes and returns the canonical label of the resulting class. *)

function uf_union(x, y: integer): Integer;
begin
  result := uf_find(y);
  labels[uf_find(x)] := result;
end;


(*  uf_make_set creates a new equivalence class and returns its label *)

function uf_make_set: integer;
begin
  inc(labels[0]);
  labels[labels[0]] := labels[0];
  result := labels[0];
end;

(*  uf_intitialize sets up the data structures needed by the union-find implementation. *)

procedure uf_initialize(max_labels: integer);
begin
  n_labels := max_labels;
  SetLength(labels, n_labels);
  labels[0] := 0;
end;

(*  uf_done frees the memory used by the union-find data structures *)

procedure uf_done;
begin
  n_labels := 0;
  setlength(labels, 0);
  labels := nil;
end;

(* End Union-Find implementation *)

function max(const a, b: integer): integer;
begin
  if a > b then result := a else result := b;
end;

function min(const a, b: integer): integer;
begin
  if a > b then result := b else result := a;
end;

type
  TMatrix = array of array of integer;

(* Label the clusters in "matrix".  Return the total number of clusters found. *)

function hoshen_kopelman(matrix: TMatrix; m, n: integer): integer;
var i, j: integer;
  up, left: integer;
  new_labels: array of integer;
  x: integer;
begin
  uf_initialize(m * n div 2);

  (* scan the matrix *)

  for i := 0 to m - 1 do
    for j := 0 to n - 1 do
      if (matrix[i][j] <> 0) then // if occupied ...
      begin
        if i = 0 then up := 0 else up := matrix[i - 1][j]; //  look up
        if j = 0 then left := 0 else left := matrix[i][j - 1]; //  look left

        case up + left of //switch (!!up + !!left) {
          0: matrix[i][j] := uf_make_set(); // a new cluster
          1: matrix[i][j] := max(up, left); // part of an existing cluster
                                            // whichever is nonzero is labelled
          2: matrix[i][j] := uf_union(up, left); // this site binds two clusters
        end;
      end;

  (* apply the relabeling to the matrix *)

  (* This is a little bit sneaky.. we create a mapping from the canonical labels
     determined by union/find into a new set of canonical labels, which are
     guaranteed to be sequential. *)

  Setlength(new_labels, n_labels); // allocate array, initialized to zero

  for i := 0 to m - 1 do
    for j := 0 to n - 1 do
      if (matrix[i][j] <> 0) then
      begin
        x := uf_find(matrix[i][j]);
        if (new_labels[x] = 0) then
        begin
          inc(new_labels[0]);
          new_labels[x] := new_labels[0];
        end;
        matrix[i][j] := new_labels[x];
      end;

  result := new_labels[0];

  setlength(new_labels, 0);
  uf_done();
end;

Was mir nen bisschen unbekannt vorkommt ist dieses Konstrukt:
switch (!!up + !!left)

[Peter666]

Ich hab das Beispiel (http://www.ocf.berkeley.edu/~fricke/...nkopelman.html) - danke BoolString dafür - komplett nachgebaut.
So ganz funktioniert das aber nicht, da die case Anweisung wohl nicht richtig übersetzt wurde.
Wie muss ich mir "switch (!!Operator1 + !!Operator2)" vorstellen? "case (byte(Operator1>0) + byte(Operator2>0)) of" ist ja wohl doch was anderes, oder?

Peter
PS: Danke nochmal für die hilfreichen Antworten :)
PPS: Nun funktioniert es doch - der Fehler lag in uf_union. Manchmal hilft auch lesen...