Entwickler-Ecke

Algorithmen, Optimierung und Assembler - Buddhabrot-Algorithmus


Jakob_Ullmann - Mo 22.11.10 19:19
Titel: Buddhabrot-Algorithmus
Hi!

Ich wollte jetzt mal den Buddhabrot-Algorithmus zur Darstellung der Mandelbrotmenge ausprobieren. Folgenden Quellcode habe ich dafür geschrieben (Lazarus):


Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
 
unit unit1;

{$mode objfpc}{$H+}

interface

uses
  Classes, SysUtils, FileUtil, LResources, Forms, Controls, Graphics, Dialogs,
  StdCtrls, ExtCtrls;

type
  TKomplex = record
    Re, Im: Double; // Real- und Imaginärteil
  end;

  { TForm1 }

  TForm1 = class(TForm)
    Bevel1: TBevel;
    Button1: TButton;
    Button2: TButton;
    Edit1: TEdit;
    Label1: TLabel;
    Label2: TLabel;
    PaintBox1: TPaintBox;
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
  private
    { private declarations }
  public
    { public declarations }
  end; 

var
  Form1: TForm1;
  cando: Boolean;

implementation

{ TForm1 }

function Komplex(a, b: Single): TKomplex;
begin
  Result.Re := a;
  Result.Im := b;
end;

// Betrag einer komplexen Zahl
function ComplexAbs(x: TKomplex): Double;
begin
  Result := Sqrt(x.Im * x.Im + x.Re * x.Re);
end;

// Zwei komplexe Zahlen addieren:
function ComplexAdd(a, b: TKomplex): TKomplex;
begin
  Result.Re := a.Re + b.Re;
  Result.Im := a.Im + b.Im;
end;

// eine kompl. Zahl quadrieren
function ComplexSqr(x: TKomplex): TKomplex;
begin
  Result.Re := x.Re * x.Re - x.Im * x.Im;
  Result.Im := 2 * x.Re * x.Im;
end;

function RGB(r, g, b: Byte): TColor;
begin
  Result := r * $000001 + g * $000100 + b * $010000;
end;

procedure TForm1.Button1Click(Sender: TObject);
var
  x, y, xx, yy, i: Integer;
  c, z, zz: TKomplex;
  a, w, h: Single;
  pfad: array[1..500] of TKomplex;
  ctr: array[0..399, 0..399] of Single;
  maxctr: Single;
begin
  cando := true;
  for x := 0 to PaintBox1.Width - 1 do
    for y := 0 to PaintBox1.Height - 1 do
      ctr[x, y] := 0;
  maxctr := 0;
  for i := 1 to 500 do
    pfad[i] := Komplex(0, 0);
  w := 1.5;
  h := 1.5;
  // für alle Bildpunkte:
  for x := 0 to PaintBox1.Width - 1 do
    for y := 0 to PaintBox1.Height - 1 do
  begin
    Application.ProcessMessages;
    if not cando then
      Break;
    z := Komplex(0, 0); // z = 0
    c := Komplex((x - PaintBox1.Width  div 2) * w / (PaintBox1.Width  div 2),
                -(y - PaintBox1.Height div 2) * h / (PaintBox1.Height div 2));
    i := 0;
    a := ComplexAbs(z); // sollte 0 sein
    while (a < 100) and (i < 500) do
    begin
      zz := ComplexSqr(z);
      z := ComplexAdd(zz, c); // z' = z² + c
      Inc(i);
      pfad[i] := z;
      Application.ProcessMessages;
      if not cando then
        Break;
      a := ComplexAbs(z);
    end;
    Application.ProcessMessages;
    if not cando then
      Break;
    if a < 99 then
      for i := 1 to 500 do
      begin
        xx := Trunc(PaintBox1.Width  div 2 + pfad[i].Re *
                      (PaintBox1.Width  div 2) / w);
        yy := Trunc(PaintBox1.Height div 2 - pfad[i].Im *
                      (PaintBox1.Height div 2) / h);
        ctr[xx, yy] := ctr[xx, yy] + 1;
        Application.ProcessMessages;
        if not cando then
          Break;
        if (ctr[xx, yy] > maxctr) then
          maxctr := ctr[xx, yy];
      end;
  end;
  if maxctr > 0 then
    maxctr := ln(maxctr + 1);
  for x := 0 to PaintBox1.Width - 1 do
    for y := 0 to PaintBox1.Height - 1 do
  begin
    Application.ProcessMessages;
    if ctr[x, y] >= 1 then
      ctr[x, y] := ln(ctr[x, y] + 1);
    PaintBox1.Canvas.Pixels[x, y] :=
      RGB(0, 0, Trunc(ctr[x, y] / maxctr * 255));
  end;
end;

procedure TForm1.Button2Click(Sender: TObject);
begin
  cando := false;
end;

initialization
  {$I unit1.lrs}

end.


Quellcode sollte ja klar sein, was er macht. Ich habe mir gedacht, wenn die Werte, wie oft ein Punkt getroffen wird, ist es sinnvoll, nicht den Wert selbst, sondern seinen Logarithmus (hier einfach spontan den natürlichen Logarithmus genommen) zur Farbwertbestimmung heranzuziehen.

Hier wird Buddhabrot beschrieben und ein Bild gezeigt:

http://de.wikipedia.org/wiki/Mandelbrot-Menge#Buddhabrot

Klar, dass meins 90° gedreht ist. Aber das sieht ja auch sonst total anders aus.

Mache ich was falsch? Vielleicht hat das ja jemand schon gemacht oder versteht sofort, wie das Bild entstanden ist.

Über die vielen Application.ProcessMessages und die Prüfungen auf cando bitte hinwegsehen. Die habe ich eingebaut, nachdem ich Lazarus zum tausendsten Mal killen musste.


Kha - Mo 22.11.10 20:04
Da ich den Output deines Programms nicht vor Augen habe, kann ich da wenig beurteilen ;) , aber in einem etwas älteren Thread [http://www.delphipraxis.net/54146-buddhabrot.html] findest du unter anderem eine Implementierung von mir - sogar noch in Delphi :) .


bummi - Mo 22.11.10 20:12
Hallo Jakob, ich habe Deinen mal etwas angefummelt, muß aber jetzt weg

Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
 
unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls, StdCtrls;

const
C_W=600;
C_H=600;

type
  TKomplex = record
    Re, Im: Double; // Real- und Imaginärteil
  end;

  TForm1 = class(TForm)
    Button1: TButton;
    Button2: TButton;
    Image1: TImage;
    procedure Button1Click(Sender: TObject);
    procedure Button2Click(Sender: TObject);
  private
    { Private-Deklarationen }
    pfad: array[1..1000] of TKomplex;
    ctr: array[0..C_W, 0..C_H] of Single;
  public
    { Public-Deklarationen }
  end;

var
  Form1: TForm1;
  cando:Boolean;
implementation

{$R *.dfm}


function Komplex(a, b: Single): TKomplex;
begin
  Result.Re := a;
  Result.Im := b;
end;

// Betrag einer komplexen Zahl
function ComplexAbs(x: TKomplex): Double;
begin
  Result := Sqrt(x.Im * x.Im + x.Re * x.Re);
end;

// Zwei komplexe Zahlen addieren:
function ComplexAdd(a, b: TKomplex): TKomplex;
begin
  Result.Re := a.Re + b.Re;
  Result.Im := a.Im + b.Im;
end;

// eine kompl. Zahl quadrieren
function ComplexSqr(x: TKomplex): TKomplex;
begin
  Result.Re := x.Re * x.Re - x.Im * x.Im;
  Result.Im := 2 * x.Re * x.Im;
end;



procedure TForm1.Button1Click(Sender: TObject);
CONST
  PixelCountMax = MaxInt / 3;
TYPE
  pRGBTripleArray = ^TRGBTripleArray;
  TRGBTripleArray = ARRAY[0..$effffff] OF TRGBTriple;
var
  pscanLine : pRGBTripleArray;
  x, y, xx, yy, i: Integer;
  c, z, zz: TKomplex;
  a, w, h: Single;
  maxctr: Single;
  bmp:TBitMap;
  WH,HH:Integer;
begin
  Image1.Width := C_W;
  Image1.Height := C_H;
  WH := C_W div 2;
  HH := C_H div 2;
  bmp:=TBitMap.Create;
  try
  bmp.PixelFormat := pf24bit;
  bmp.Width := C_W;
  bmp.Height := C_H;
  cando := true;
  for x := 0 to C_W - 1 do
    for y := 0 to C_H - 1 do
      ctr[x, y] := 0;
  maxctr := 0;
  for i := 1 to 500 do pfad[i] := Komplex(0, 0);
  w := 1.5;
  h := 1.5;

  // für alle Bildpunkte:
  for x := 0 to C_W - 1 do
    for y := 0 to C_W - 1 do
  begin
    z := Komplex(0, 0); // z = 0
    c := Komplex((x - WH) * w / WH,
                -(y - HH) * h / HH);
    i := 0;
    a := ComplexAbs(z); // sollte 0 sein
    while (a < 100) and (i < 500) do
    begin
      zz := ComplexSqr(z);
      z := ComplexAdd(zz, c); // z' = z² + c
      Inc(i);
      pfad[i] := z;
      a := ComplexAbs(z);
    end;

    if a < 99 then
      for i := 1 to 100 do
      begin
        xx := Trunc(WH + pfad[i].Re * HH / w);
        yy := Trunc(HH - pfad[i].Im * (HH) / h);

//            ASSERT(xx > -1);
//            ASSERT(xx < (C_W - 1));
//            ASSERT(yy > -1);
//            ASSERT(yy < (C_H - 1));


        if (xx > -1) and (xx < (C_W - 1)) and (yy > -1) and  (yy < (C_H - 1)) then
            begin
            ctr[xx, yy] := ctr[xx, yy] + 1;
            if (ctr[xx, yy] > maxctr) then maxctr := ctr[xx, yy];
            end;
      end;
  end;
  if maxctr > 0 then maxctr := ln(maxctr + 1);
  for y := 0 to C_H - 1 do
       begin
       pscanLine := bmp.Scanline[y];
        for x := 0 to C_W -1 do
           begin
             if ctr[x, y] >= 1 then ctr[x, y] := ln(ctr[x, y] + 1);
             pscanLine[x].rgbtBlue   := Trunc(ctr[x, y] / maxctr * 255);
             pscanLine[x].rgbtRed    := 0;
             pscanLine[x].rgbtGreen    := 0;
           end;
       end;
  image1.Picture.bitmap.Assign(bmp);
  finally
  bmp.Free;
  end;

end;


end.


EDIT, ich habe die Arrays in einem 2. Edit mal vom Stack genommen (heißen daher noch nicht Fxxx damit ist das Stackoverflowprobelem ab 800 * 800 erst mal weg)


Jakob_Ullmann - Mo 22.11.10 20:31
Also danke schonmal euch beiden. Ich habe beide Quelltexte versucht, zu verstehen, bin aber noch nicht mal mit Kha's durch (aber ich kann mich an einen Thread erinnern, wo du das TMathImage mal vorgestellt hast. :oops: Habe ihn aber nie wiedergefunden).

Hier die Ausgabe des Programms:

buddhabrot1

Ich habe nochmal den englischen Artikel durchgelesen. Sieht ein wenig wie Antibrot aus (also die "falschen" Punkte verfolgt). Ich schaue mal, ob das der Fehler sein könnte.

EDIT: Scheint der Fehler zu sein. Allerdings schaffe ich das irgendwie nicht, den zu korrigieren. Ich versuchs morgen nochmal.


bummi - Mo 22.11.10 20:40
http://www.delphipraxis.net/51154-tmathbitmap.html


Jakob_Ullmann - Di 23.11.10 16:15
Ich habe es jetzt geschafft, den echten Buddhabrot (und nicht Antibrot oder wie das hieß) zum Laufen zu bringen. Allerdings hat man dann eine Art Kästchenmuster. Deshalb habe ich die Größe geviertelt, sodass jetzt ein Punkt als Mittelwert von vier Punkten berechnet wird. Das Bild ist jetzt bei massiver Rechenzeit total klein und die Linien sind immernoch da. Ich werde das Gefühl nicht los, dass ich hier einen systematischen Fehler mache...

buddha2


Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
 
procedure TForm1.Button1Click(Sender: TObject);
var
  x, y, xx, yy, i, j: Integer;
  c, z, zz: TKomplex;
  a, w, h: Single;
  pfad: array[1..100] of TKomplex;
  ctr: array[0..399, 0..399] of Single;
  mycol: Single;
  maxctr: Single;
begin
  cando := true;
  for x := 0 to PaintBox1.Width - 1 do
    for y := 0 to PaintBox1.Height - 1 do
      ctr[x, y] := 0;
  maxctr := 0;
  for i := 1 to 100 do
    pfad[i] := Komplex(0, 0);
  w := 1.5;
  h := 1.5;
  // für alle Bildpunkte:
  for x := 0 to PaintBox1.Width - 1 do
    for y := 0 to PaintBox1.Height - 1 do
  begin
    Application.ProcessMessages;
    if not cando then
      Break;
    z := Komplex(0, 0); // z = 0
    c := Komplex((x - PaintBox1.Width  div 2) * w / (PaintBox1.Width  div 2),
                -(y - PaintBox1.Height div 2) * h / (PaintBox1.Height div 2));
    i := 0;
    a := ComplexAbs(z); // sollte 0 sein
    while (a < 2) and (i < 100) do
    begin
      zz := ComplexSqr(z);
      z := ComplexAdd(zz, c); // z' = z² + c
      Inc(i);
      pfad[i] := z;
      Application.ProcessMessages;
      if not cando then
        Break;
      a := ComplexAbs(z);
    end;
    Application.ProcessMessages;
    if not cando then
      Break;
    if a > 2 then
    begin
      if i > 0 then
      for j := 1 to i do
      begin
        if (pfad[j].Re = 0) and (pfad[j].Im = 0) then
          Break
        else begin
          xx := Trunc(PaintBox1.Width  div 2 + pfad[j].Re *
                        (PaintBox1.Width  div 2) / w);
          yy := Trunc(PaintBox1.Height div 2 - pfad[j].Im *
                        (PaintBox1.Height div 2) / h);
          if (xx >= 0) and (yy >= 0) and (xx < PaintBox1.Width)
             and (yy < PaintBox1.Height) then
          begin
            ctr[xx, yy] := ctr[xx, yy] + 1;
          //else
            //ShowMessage(IntToStr(xx) + '  ' + IntToStr(yy));
            if (ctr[xx, yy] > maxctr) then
              maxctr := ctr[xx, yy];
          end;
          Application.ProcessMessages;
          if not cando then
            Break;
        end;
      end;
      for i := 1 to 100 do
        pfad[i] := Komplex(0, 0);
    end;
  end;
  if maxctr > 0 then
    maxctr := ln(maxctr + 1);
  for x := 1 to PaintBox1.Width - 1 do
    for y := 1 to PaintBox1.Height - 1 do
      if ctr[x, y] >= 1 then
        ctr[x, y] := ln(ctr[x, y] + 1);
  for x := 1 to PaintBox1.Width div 2 - 2 do
    for y := 1 to PaintBox1.Height div 2 - 2 do
  begin
    Application.ProcessMessages;
    if ctr[x, y] >= 1 then
      ctr[x, y] := ln(ctr[x, y] + 1);
    mycol := 1/4 *
        (ctr[2 * x - 1, 2 * y - 1] + ctr[2 * x, 2 * y - 1] +
         ctr[2 * x - 1, 2 * y]     + ctr[2 * x, 2 * y]);
    PaintBox1.Canvas.Pixels[x, y] :=
      RGB(Trunc(mycol / maxctr * 155), Trunc(mycol / maxctr * 155), Trunc(mycol / maxctr * 255));
  end;
end;

end.


elundril - Di 23.11.10 17:38
Mal meine Vermutung in den Raum geschossen: Du hast ein Application.ProcessMessages; drin, was ziemlich verlangsamen kann; Du arbeitest mit PaintBox.Pixels[]; statt mit Scanline.

Erstes kannst du umgehen indem du ev. mit einem TBitmap und einem TThread arbeitest, glaube ich. Zweites kannst du mit Scanline umgehen, was auch dann deutlich schneller sein sollte.

lg elundril


Jakob_Ullmann - Di 23.11.10 17:43
Also das mit dem Application.ProcessMessages ist klar. Das könnte das Teil ordentlich verlangsamen. Scanline muss aber nicht wirklich sein. Das Zeichnen dauert weniger als eine Sekunde, das Rendern mindestens zwanzig Sekunden. Das ist dann eher was, wenn das Teil mal ordentlich aussieht. Aber das tut es ja noch nicht.

EDIT:
Das mit dem ProcessMessages war ein guter Tipp. So kann ich in derselben Zeit, wie vorher 300 Iterationen gingen, 15000 (fünfzehntausend!) Iterationen machen. Sieht schonmal nicht schlecht aus:

[Vergleich]
user defined image

buddhabrot3

Aber da sind immer noch ein paar "Störpixel". Die werden immer dunkler, je mehr Iterationen ich mache. Aber ich kann mit meinem "kleinen" PC ja schlecht 100.000 Iterationen machen. :shock: Außerdem werden die wieder sichtbar, sobald ich wieder den Log. nat. als Maß für die Farbe nehme. Und das muss ich irgendwann, sonst sieht man gar nichts mehr.


Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
 
procedure TForm1.Button1Click(Sender: TObject);
var
  x, y, xx, yy, i, j: Integer;
  c, z, zz: TKomplex;
  a, w, h: Single;
  pfad: array[1..15000] of TKomplex;
  ctr: array[0..399, 0..399] of Single;
  mycol: Single;
  maxctr: Single;
begin
  cando := true;
  for x := 0 to PaintBox1.Width - 1 do
    for y := 0 to PaintBox1.Height - 1 do
      ctr[x, y] := 0;
  maxctr := 0;
  for i := 1 to 15000 do
    pfad[i] := Komplex(0, 0);
  w := 2;
  h := 2;
  // für alle Bildpunkte:
  for x := 0 to PaintBox1.Width - 1 do
    for y := 0 to PaintBox1.Height - 1 do
  begin
    z := Komplex(0, 0); // z = 0
    c := Komplex((x - PaintBox1.Width  div 2) * w / (PaintBox1.Width  div 2),
                -(y - PaintBox1.Height div 2) * h / (PaintBox1.Height div 2));
    i := 0;
    a := ComplexAbs(z); // sollte 0 sein
    while (a < 3) and (i < 15000) do
    begin
      zz := ComplexSqr(z);
      z := ComplexAdd(zz, c); // z' = z² + c
      Inc(i);
      pfad[i] := z;
      a := ComplexAbs(z);
    end;
    Application.ProcessMessages;
    if not cando then
      Break;
    if a > 3 then
    begin
      if i > 0 then
      for j := 1 to i do
      begin
        if (pfad[j].Re = 0) and (pfad[j].Im = 0) then
          Break
        else begin
          xx := Trunc(PaintBox1.Width  div 2 + pfad[j].Re *
                        (PaintBox1.Width  div 2) / w);
          yy := Trunc(PaintBox1.Height div 2 - pfad[j].Im *
                        (PaintBox1.Height div 2) / h);
          if (xx >= 0) and (yy >= 0) and (xx < PaintBox1.Width)
             and (yy < PaintBox1.Height) then
          begin
            ctr[xx, yy] := ctr[xx, yy] + 1;
          //else
            //ShowMessage(IntToStr(xx) + '  ' + IntToStr(yy));
            if (ctr[xx, yy] > maxctr) then
              maxctr := ctr[xx, yy];
          end;
        end;
      end;
      for i := 1 to 15000 do
        pfad[i] := Komplex(0, 0);
    end;
  end;
  if maxctr > 0 then
    maxctr := maxctr / 8;
  for x := 1 to PaintBox1.Width - 1 do
    for y := 1 to PaintBox1.Height - 1 do
      //if ctr[x, y] <= maxctr / 5 then
      //  ctr[x, y] := 0
      if ctr[x, y] >= maxctr then
        ctr[x, y] := maxctr;
  //maxctr := ln(maxctr + 1);
  for x := 1 to PaintBox1.Width - 1 do
    for y := 1 to PaintBox1.Height - 1 do
  begin
    Application.ProcessMessages;
    //if ctr[x, y] >= 1 then
    //  ctr[x, y] := ln(ctr[x, y] + 1);
    mycol := ctr[x, y];
    PaintBox1.Canvas.Pixels[x, y] :=
      RGB(Trunc(mycol / maxctr * 155), Trunc(mycol / maxctr * 200), Trunc(mycol / maxctr * 255));
  end;
end;


bummi - Di 23.11.10 18:16
Sieht zwar immer noch nicht perfekt aus, aber das Raster ist weg

Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
 
unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls, StdCtrls;
const
C_W=400;
C_H=400;

type
  TKomplex = record
    Re, Im: Double; // Real- und Imaginärteil
  end;

  TForm1 = class(TForm)
    Button1: TButton;
    Image1: TImage;
    procedure Button1Click(Sender: TObject);
  private
    { Private-Deklarationen }
  public
  pfad: array[1..100] of TKomplex;
  ctr: array[0..C_W, 0..C_H] of Single;

    { Public-Deklarationen }
  end;
var
  Form1: TForm1;
  cando:Boolean;
implementation

{$R *.dfm}


function Komplex(a, b: Single): TKomplex;
begin
  Result.Re := a;
  Result.Im := b;
end;

// Betrag einer komplexen Zahl
function ComplexAbs(x: TKomplex): Double;
begin
  Result := Sqrt(x.Im * x.Im + x.Re * x.Re);
end;

// Zwei komplexe Zahlen addieren:
function ComplexAdd(a, b: TKomplex): TKomplex;
begin
  Result.Re := a.Re + b.Re;
  Result.Im := a.Im + b.Im;
end;

// eine kompl. Zahl quadrieren
function ComplexSqr(x: TKomplex): TKomplex;
begin
  Result.Re := x.Re * x.Re - x.Im * x.Im;
  Result.Im := 2 * x.Re * x.Im;
end;



procedure TForm1.Button1Click(Sender: TObject);
CONST
  PixelCountMax = MaxInt / 3;
TYPE
  pRGBTripleArray = ^TRGBTripleArray;
  TRGBTripleArray = ARRAY[0..$effffff] OF TRGBTriple;
var
  pscanLine : pRGBTripleArray;
  x, y, xx, yy, i, j: Integer;
  c, z, zz: TKomplex;
  a, w, h: Single;
  maxctr: Single;
  bmp:TBitMap;
  WH,HH:Integer;
begin
  Image1.Width := C_W;
  Image1.Height := C_H;
  WH := C_W div 2;
  HH := C_H div 2;
  bmp:=TBitMap.Create;
  try
  bmp.PixelFormat := pf24bit;
  bmp.Width := C_W;
  bmp.Height := C_H;
  cando := true;
  for x := 0 to C_W - 1 do
    for y := 0 to C_H - 1 do
      ctr[x, y] := 0;
  maxctr := 0;
  for i := LOW(pfad) to High(pfad) do pfad[i] := Komplex(0, 0);
  w := 1.5;
  h := 1.5;

  // für alle Bildpunkte:
  for x := 0 to C_W - 1 do
    for y := 0 to C_W - 1 do
  begin
    z := Komplex(0, 0); // z = 0
    c := Komplex((x - WH) * w / WH,
                -(y - HH) * h / HH);
    i := 0;
    a := ComplexAbs(z); // sollte 0 sein
    while (a < 2) and (i < High(pfad)) do
    begin
      zz := ComplexSqr(z);
      z := ComplexAdd(zz, c); // z' = z² + c
      Inc(i);
      pfad[i] := z;
      a := ComplexAbs(z);
    end;


    if (a > 2)   then
    begin
      for j := 1 to i do
      begin
        if not ( (pfad[j].Re = 0) and (pfad[j].Im = 0)) then
          begin
            xx := Round(WH + pfad[j].Re * WH / w);
            yy := Round(HH - pfad[j].Im * (HH) / h);

            if (xx > -1) and (xx < (C_W - 1)) and (yy > -1) and  (yy < (C_H - 1)) then
                begin
                ctr[xx, yy] := ctr[xx, yy] + 1;
                if (ctr[xx, yy] > maxctr) then maxctr := ctr[xx, yy];
                end;
          end;
      end;
    end;
  end;
  if maxctr > 0 then maxctr := ln(maxctr + 1);
  for y := 0 to C_H - 1 do
       begin
       pscanLine := bmp.Scanline[y];
        for x := 0 to C_W -1 do
           begin
             if ctr[x, y] >= 1 then ctr[x, y] := ln(ctr[x, y] + 1);
            // pscanLine[x].rgbtBlue   := 255 - Trunc(ctr[x, y] / maxctr * 255);
             pscanLine[x].rgbtBlue   :=  Trunc(ctr[x, y] / maxctr * 255);
             pscanLine[x].rgbtRed    := 0;
             pscanLine[x].rgbtGreen    := 0;
           end;
       end;
  image1.Picture.bitmap.Assign(bmp);
  finally
  bmp.Free;
  end;

end;

end.


Jakob_Ullmann - Di 23.11.10 18:27
Aber eigentlich machst du doch nichts anders als das TBitmap32 und das Scanline... OK, du hast es z. B. auf High() und Low() umgeschrieben und es dadurch eleganter gelöst, aber im Grunde ist doch die eigentliche Funktionsweise dieselbe...


bummi - Di 23.11.10 18:34
Ist ein normales Bitmap und ich habe diese Stelle umgebaut

Delphi-Quelltext 
1:
2:
 
            xx := Round(WH + pfad[j].Re * WH / w);
            yy := Round(HH - pfad[j].Im * (HH) / h);


Jakob_Ullmann - Di 23.11.10 19:20
Ich hätte gar nicht gedacht, dass Round und Trunc so einen großen Unterschied macht. Aber die Störpixel sind leider trotzdem noch da. Eigenartigerweise scheint es jetzt ein Raster zu geben, dass nicht kartesisch, sondern polar (also vom Mittelpunkt in alle Richtungen ausgehend) ausgerichtet ist. Aber das kann auch nur Einbildung sein (denn ich sehe dafür keine logische Erklärung).

Es wird ja abgebrochen, sobald der Bailout-Radius 2 überschritten ist. Daher lasse ich jetzt pfad[j] nur bis pfad[i - 1] laufen. Dann hat man außerhalb von |z| < 2 keine Störpixel mehr. Aber irgendwie müssen doch auch die vorhandenen weggehen. :(


Delphi-Laie - Di 23.11.10 20:15
user profile iconJakob_Ullmann hat folgendes geschrieben Zum zitierten Posting springen:
Ich habe mir gedacht, wenn die Werte, wie oft ein Punkt getroffen wird, ist es sinnvoll, nicht den Wert selbst, sondern seinen Logarithmus (hier einfach spontan den natürlichen Logarithmus genommen) zur Farbwertbestimmung heranzuziehen.


Dieser Satz ist leider wirr. Es scheint darin etwas verlorengegangen zu sein. Was meintest Du?


bummi - Mi 24.11.10 16:26
so sieht es bei mir eigentlich gut aus

Delphi-Quelltext 
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:
54:
55:
56:
57:
58:
59:
60:
61:
62:
63:
64:
65:
66:
67:
68:
69:
70:
71:
72:
73:
74:
75:
76:
77:
78:
79:
80:
81:
82:
83:
84:
85:
86:
87:
88:
89:
90:
91:
92:
93:
94:
95:
96:
97:
98:
99:
100:
101:
102:
103:
104:
105:
106:
107:
108:
109:
110:
111:
112:
113:
114:
115:
116:
117:
118:
119:
120:
121:
122:
123:
124:
125:
126:
127:
128:
129:
130:
131:
132:
133:
134:
135:
136:
137:
138:
139:
140:
141:
142:
143:
144:
145:
146:
147:
148:
149:
150:
151:
152:
153:
154:
155:
156:
157:
158:
159:
160:
161:
162:
163:
164:
165:
166:
167:
168:
169:
170:
171:
172:
173:
174:
175:
176:
177:
178:
179:
180:
181:
182:
183:
184:
185:
186:
187:
188:
189:
190:
191:
192:
193:
194:
195:
196:
197:
198:
199:
200:
201:
202:
203:
204:
205:
206:
207:
208:
209:
210:
211:
212:
213:
214:
215:
216:
217:
218:
219:
220:
221:
222:
223:
224:
225:
226:
227:
228:
229:
230:
231:
232:
233:
234:
235:
236:
237:
238:
239:
240:
241:
 
unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls, StdCtrls,math;
const
C_W=500;
C_H=500;

type
  TKomplex = record
    Re, Im: Double; // Real- und Imaginärteil
  end;

  TForm1 = class(TForm)
    Button1: TButton;
    Image1: TImage;
    procedure Button1Click(Sender: TObject);
    procedure Image1Click(Sender: TObject);
  private
    { Private-Deklarationen }
  public
  pfad: array[1..20000] of TKomplex;
  ctr: array[0..C_W, 0..C_H] of Single;

    { Public-Deklarationen }
  end;
var
  Form1: TForm1;
  cando:Boolean;
implementation

{$R *.dfm}


function Komplex(a, b: Single): TKomplex;
begin
  Result.Re := a;
  Result.Im := b;
end;

// Betrag einer komplexen Zahl
function ComplexAbs(x: TKomplex): Double;
begin
  Result := Sqrt(x.Im * x.Im + x.Re * x.Re);
end;

// Zwei komplexe Zahlen addieren:
function ComplexAdd(a, b: TKomplex): TKomplex;
begin
  Result.Re := a.Re + b.Re;
  Result.Im := a.Im + b.Im;
end;

// eine kompl. Zahl quadrieren
function ComplexSqr(x: TKomplex): TKomplex;
begin
  Result.Re := x.Re * x.Re - x.Im * x.Im;
  Result.Im := 2 * x.Re * x.Im;
end;



procedure TForm1.Button1Click(Sender: TObject);
CONST
  PixelCountMax = MaxInt / 3;
TYPE
  pRGBTripleArray = ^TRGBTripleArray;
  TRGBTripleArray = ARRAY[0..$effffff] OF TRGBTriple;
var
  pscanLine : pRGBTripleArray;
  x, y, xx, yy, i, j: Integer;
  c, z, zz: TKomplex;
  a, w, h: Single;
  maxctr: Single;
  bmp:TBitMap;
  WH,HH:Integer;
begin
  Image1.Width := C_W;
  Image1.Height := C_H;
  WH := C_W div 2;
  HH := C_H div 2;
  bmp:=TBitMap.Create;
  try
  bmp.PixelFormat := pf24bit;
  bmp.Width := C_W;
  bmp.Height := C_H;

  for x := 0 to C_W - 1 do
    for y := 0 to C_H - 1 do
      ctr[x, y] := 0;
  maxctr := 0;
  for i := LOW(pfad) to High(pfad) do pfad[i] := Komplex(0, 0);
  w := 1.5;
  h := 1.5;

  // für alle Bildpunkte:
  for x := 0 to C_W - 1 do
    for y := 0 to C_W - 1 do
  begin
    z := Komplex(0, 0); // z = 0
    c := Komplex((x - WH) * w / WH,
                -(y - HH) * h / HH);
    i := 0;
    a := ComplexAbs(z); // sollte 0 sein
    while (a <= 2) and (i < High(pfad)) do
    begin
      zz := ComplexSqr(z);
      z := ComplexAdd(zz, c); // z' = z² + c
      Inc(i);
      pfad[i] := z;
      a := ComplexAbs(z);
    end;


    if (a > 2)   then
    begin
      for j := 1 to i do
      begin
        if not ( (pfad[j].Re = 0) and (pfad[j].Im = 0)) then
          begin
            xx := Round(WH + pfad[j].Re * WH / w);
            yy := Round(HH - pfad[j].Im * (HH) / h);

            if (xx > -1) and (xx < (C_W - 1)) and (yy > -1) and  (yy < (C_H - 1)) then
                begin
                ctr[xx, yy] := ctr[xx, yy] + 1;
                if (ctr[xx, yy] > maxctr) then maxctr := ctr[xx, yy];
                end;
          end;
      end;
    end;
  end;
  if maxctr > 0 then maxctr := LnXP1(maxctr );
  for y := 0 to C_H - 1 do
       begin
       pscanLine := bmp.Scanline[y];
        for x := 0 to C_W -1 do
           begin

             if ctr[x, y] >= 1 then ctr[x, y] := LnXP1(ctr[x, y] );
             pscanLine[x].rgbtBlue   :=  Trunc(255 * ctr[x, y] / maxctr );
             pscanLine[x].rgbtRed    := 0;
             pscanLine[x].rgbtGreen    := 0;
           end;
       end;
  image1.Picture.bitmap.Assign(bmp);
  finally
  bmp.Free;
  end;

end;





Procedure InvertBitMap24(bmp:TBitMap);
CONST
  PixelCountMax = MaxInt / 3;
TYPE
  pRGBTripleArray = ^TRGBTripleArray;
  TRGBTripleArray = ARRAY[0..$effffff] OF TRGBTriple;
var
  pscanLine : pRGBTripleArray;
  x,y:Integer;
begin
    for Y := 0 to bmp.Height -1 do
     begin
         pscanLine := bmp.Scanline[y];
         for x := 0 to bmp.Width -1 do
           begin
           pscanLine[x].rgbtBlue   := pscanLine[x].rgbtBlue XOR 255;
           pscanLine[x].rgbtGreen   := pscanLine[x].rgbtGreen XOR 255;
           pscanLine[x].rgbtRed   := pscanLine[x].rgbtRed XOR 255;
           end;
     end;
end;

Procedure InvertBitMap32(bmp:TBitMap);
CONST
  PixelCountMax = MaxInt / 3;
TYPE
  pRGBQuadArray = ^TRGBQuadArray;
  TRGBQuadArray = ARRAY[0..$effffff] OF TRGBQuad;
var
  pscanLine : pRGBQuadArray;
  x,y:Integer;
begin
    for Y := 0 to bmp.Height -1 do
     begin
         pscanLine := bmp.Scanline[y];
         for x := 0 to bmp.Width -1 do
           begin
           pscanLine[x].rgbBlue   := pscanLine[x].rgbBlue XOR 255;
           pscanLine[x].rgbGreen   := pscanLine[x].rgbGreen XOR 255;
           pscanLine[x].rgbRed   := pscanLine[x].rgbRed XOR 255;
           end;
     end;
end;

Procedure InvertBitMap(bmp:TBitMap);
begin
  if bmp.PixelFormat=pf32Bit then  InvertBitMap32(bmp)
  else if bmp.PixelFormat=pf24Bit then  InvertBitMap24(bmp);
end;



procedure MirrorBitmap(Bmp, MBmp: TBitmap;Horizonal:Boolean=true);
var
  x1, x2, y1, y2: integer;
begin
  MBmp.Width := Bmp.Width;
  MBmp.Height := Bmp.Height;
  if Horizonal then
    begin
        x1 := MBmp.Width - 1;
        x2 := - 1;
        y1 := 0;
        y2 := MBmp.Height;
    end
  else
    begin
        x1 := 0;
        x2 := MBmp.Width;
        y1 := MBmp.Height - 1;
        y2 := -1;
    end;
  MBmp.Canvas.CopyRect(Rect(x1, y1, x2, y2), Bmp.Canvas, Rect(0, 0, MBmp.Width, MBmp.Height));
end;
procedure TForm1.Image1Click(Sender: TObject);
begin
   // MirrorBitmap(Image1.Picture.Bitmap,Image1.Picture.Bitmap)
   InvertBitMap(Image1.Picture.Bitmap);
   Image1.Invalidate;
end;

end.


Jakob_Ullmann - Mi 24.11.10 18:03
user profile iconJakob_Ullmann hat folgendes geschrieben Zum zitierten Posting springen:
Ich habe mir gedacht, wenn die Werte, wie oft ein Punkt getroffen wird, sehr stark variieren ist es sinnvoll, nicht den Wert selbst, sondern seinen Logarithmus (hier einfach spontan den natürlichen Logarithmus genommen) zur Farbwertbestimmung heranzuziehen.


Das meinte ich eigentlich. Es gibt halt ein paar Punkte, die sehr oft getroffen werden, aber die meisten werden wesentlich seltener getroffen. Das hätte einen riesengroßen Kontrast zur Folge. Einige Punkte würden sehr hell leuchten, während andere untergingen. Wenn man nicht x nimmt, sondern ln(x+1), so hat man eine ausgewogene Spannweite, aber leider werden auch die störenden Pixel am Rand wesentlich heller.

EDIT: Ich habe mal Nachforschungen angestellt. Für gute Bilder werden schon teilweise Renderzeiten von bis zu 1h in Kauf genommen. Dann kann man natürlich problemlos ein paar Treffer abziehen, ohne das Gesamtergebnis zu beeinträchtigen. Ich bin jetzt bis auf 35000 Iterationen hochgegangen. Die Störpixel sind noch da, aber die Punkte werden mit mehr Iterationen immer dichter und irgendwann wird fast jeder signifikante Punkt mal ein paar Treffer abbekommen haben und man kann abziehen. Des weiteren enthalten die Wiki-Bilder natürlich auch mehr Details. Einfacher wird es sein, ich versuche mich mit Nebulabrot.


Entwickler-Ecke.de based on phpBB
Copyright 2002 - 2011 by Tino Teuber, Copyright 2011 - 2025 by Christian Stelzmann Alle Rechte vorbehalten.
Alle Beiträge stammen von dritten Personen und dürfen geltendes Recht nicht verletzen.
Entwickler-Ecke und die zugehörigen Webseiten distanzieren sich ausdrücklich von Fremdinhalten jeglicher Art!