Anteil eines Werts in einem Umkreis bestimmen

[Frida]

Moin,

ich möchte für jeden Punkt einer Karte den Anteil eines Wertes innerhalb eines Radius um den Punkt bestimmen. Um das ganze zu veranschaulichen: Ich habe eine Landkarte, in der eigetragen ist, ob ein Pixel bebaut ist oder eben nicht. Nun möchte ich für jedes Pixel berechnen, wie hoch der Anteil bebauter Fläche innerhalb eines Kreises um diesen Punkt herum ist.
Das funktioniert auch gut, dauert auch gar nicht so lange, aber da ich mit recht großen Arrays arbeite und das ganze für viele verschiedene Radien mache, summiert sich die Rechenzeit doch zu einer beträchtlichen Länge. Hat jemand eine Idee, wie man das ganze schneller machen kann?
Hier kommt der Code, mit dem ich es bisher berechen:

Delphi-Quellcode:

procedure CalcLandContext(radius : integer; const InData : TIntArray; var OutData : TDoubleArray);
var
  dr2                        : integer;
  i, j, ii, jj, nx, ny, n, a : integer;
  ncol, nrow                 : integer;
begin
  // InData ist ein Array mit 0 für unbebaut und 1 für bebaut, -9999 steht für NoDataValue
  ncol := Length(InData);
  nrow := Length(InData[0]);
  SetLength(OutData, ncol, nrow);
  for j:=0 to High(OutData[0]) do
    begin
    for i:=0 to High(OutData) do
      begin
      if InData[i,j] = -9999 then
        OutData[i,j]:=-9999
      else
        begin
        n:=0;
        a:=0;
        for  nx:=-radius-1 to radius+1 do
          begin
          for ny:=-radius-1 to radius+1 do
            begin
            ii := i+nx;
            jj := j+ny;
            if (ii>=0) and (ii<ncol) and (jj>=0) and (jj<nrow)  //nicht über Ränder
              then begin
              dr2 := nx*nx+ny*ny;
              if dr2 <= sqr(radius)
                then begin
                if InData[ii,jj] = 1 then
                  inc(n);
                inc(a);
                end;
              end;
            end;
          end;
          OutData[i,j] := 100*n/a;
        end;
      end;
    end;
end;

Ist es vielleicht besser, das Array erst in eine Bitmap umzuwandel und dann Farbwerte zu vergleichen? Oder gibt's eine Funktion ähnlich FlodFill, die auch mitzählt, wieviele Pixel sie einfärbt?

Wäre super, wenn irgendjemand ein paar Tipps für mich hätte.

[ibp]

zum einen...

Delphi-Quellcode:

...
  rquadr : integer;
  ..
  rquadr=sqr(radius);  // dann brauchs nicht jdesmal neu berechnet zu werden daher ausserhalb der schleife berechnen
  ..
              if nx*nx+ny*ny <= rquadr then

zum anderen, wenn du das für verschiedene radien machen willst, dann mache dies gleichzeitig für den jeweiligen mittelpunkt! d.h. du solltest dein OutData[i,j] um eine ebene der radien erweitern!


ps.: herzlich willkommen bei dp!

[Frida]

Vielen Dank, das hilft schon ein wenig.

[Sidorion]

Soll keine Besserwisserei sein, aber die Hypothenuse ist die Wurzel der Summe der Kathethenquadrate.
Also musst Du Sqrt(nx*nx+ny*ny) mit radius vergleichen.
Hierbei kann man allerdings FETT performance sparen, wenn man die Wurzeln weglässt, da gilt a<b <=> a^2<b^2.
@idb: Den Radius nochmal zu wurzeln is natürlich grund falsch, da somit a^2 mit sqrt(b) verglichen wid.(a ist die Hypothenuse, b der Radius)

Die richtge Lösung st also, den Radius zu quadrieren, da man sich dann in jedem Schleifendurchlauf die Wurzel spart.

[Khabarakh]

@Sidoron:
Hier hat niemand die Hypothenuse anders definiert, hier hat niemand einen Radius gewurzelt und all deine Verbesserungen sind im Code schon längst enthalten.

@Frida: Die Überprüfung, ob man sich außerhalb des Randes befindet, kannst du wegoptimieren, indem du die Schleifen erst am Rand starten lässt.

[Sidorion]

Zitat von ibp:

Delphi-Quellcode:

...
  rquadr : integer;
  ..
  rquadr=sqr(radius);  // dann brauchs nicht jdesmal neu berechnet zu werden daher ausserhalb der schleife berechnen
  ..
              if nx*nx+ny*ny <= rquadr then

q.e.d.

[Khabarakh]

c² wird außerhalb der Schleife gespeichert und dann mit a² + b² verglichen. Genauso, wie du es danach vorgeschlagen hast. So what?

PS: Könnte es sein, dass du ein "t" hindichtest, wo keines ist ?

[Frida]

Ich glaube, das mit dem Herrn Pythagoras haben wir geklärt; sollte so stimmen.

@Khabarakh
Wie meinst Du das?

Zitat:

@Frida: Die Überprüfung, ob man sich außerhalb des Randes befindet, kannst du wegoptimieren, indem du die Schleifen erst am Rand starten lässt.

Die äußeren Schleifen müssen ja schon ganz außen anfangen, also z.B. in einer Ecke bei InData[0,0], die dann nur von einem Viertelkreis umschlagen wird.

[ibp]

was ich nicht verstehe, ist warum du ...

Delphi-Quellcode:

..
        for  nx:=-radius-1 to radius+1 do 
          begin 
          for ny:=-radius-1 to radius+1 do
..

hier zum radius noch |1| dazuzählst?
d.h. du hast den durchmesser um 2 vergrößert!

folgendes reicht auch und du hast damit 2*2*(1+radius) schleifendurchläufe weniger!

Delphi-Quellcode:

..
        for  nx:=-radius to radius do 
          begin 
          for ny:=-radius to radius do
..

[marabu]

Hallo Frida,

du besuchst jeden Punkt p(x,y) deiner Karte mehrmals - häufiger mit zunehmendem Radius. Eine mögliche Optimierung stelle ich mir etwa so vor: Du betrachtest das deinem Kreis k(p,r) umschriebene Quadrat q(k). Dann iterierst du zeilenweise über alle Punkte der Karte so, dass gerade und ungerade Zeilen gegenläufig besucht werden und berechnest den benötigten Summenwert s. Die Fläche a(k) kannst du erstmal wie bisher mitzählen.

Durch die Wegvorschrift ändert sich die betrachtete Punktmenge in q(k), es kommen maximal etwa 4r neue Punkte hinzu und maximal die gleiche Anzahl an Punkten fällt weg. Nur für diese Punkte musst du prüfen, ob sie zur Kreisfläche k(p,r) gehören. Das errechnete Delta addierst du auf den zwischengespeicherten Wert s und bildest deinen Funktionswert f(x,y) = s/a.

Freundliche Grüße vom marabu