Random Midpoint Displacement Algorithm

Algorithms to generate fractal clouds are described in detail in the book The Science of Fractal Images by Peitgen and Saupe.
Fractint's algorithm is from Bret Mulvey. I reconverted the C source of Fractint back to a Pascal program.

The program, as it is here, draws a plasma tile 100x100. It can be modified to draw a simple fullscreen plasma, if the constants xmax, ymax are set to 639 resp. 479, and if the six (blue) lines before "readln" are deleted.

program plasma;
uses Graph,Dos;

const grad=700;
   xmax=50; ymax=50;
var
   GM,GD,x,y,level: integer;
   pal: array[0..255] of array[0..2] of Byte;
   reg: registers;

function adjust(xa,ya,x,y,xb,yb:integer): integer;
var rand:integer;
begin
   rand:=random(grad) shr level;
   if (random(2)=0) then rand:=-rand;
   rand:=((getpixel(xa,ya)+getpixel(xb,yb)+1) shr 1)+rand;
   if (rand<1) then rand:=1;
   if (rand>250) then rand:=250;
   putpixel(x,y,rand);
   adjust:=rand
end;

procedure divide(x1,y1,x2,y2:integer);
var
   x,y,i,v:integer;
begin
   if ((x2-x1<2) and (y2-y1<2)) then exit;
   inc(level);
   x:=(x1+x2) shr 1;
   y:=(y1+y2) shr 1;
   v:=getpixel(x,y1);
   if(v=0) then v:=adjust(x1,y1,x,y1,x2,y1);
   i:=v;
   v:=getpixel(x2,y);
   if(v=0) then v:=adjust(x2,y1,x2,y,x2,y2);
   i:=i+v;
   v:=getpixel(x,y2);
   if(v=0) then v:=adjust(x1,y2,x,y2,x2,y2);
   i:=i+v;
   v:=getpixel(x1,y);
   if(v=0) then v:=adjust(x1,y1,x1,y,x1,y2);
   i:=i+v;
   if (getpixel(x,y)=0) then putpixel(x,y,(i+2) shr 2);
   divide(x1,y1,x,y);
   divide(x,y1,x2,y);
   divide(x,y,x2,y2);
   divide(x1,y,x,y2);
   dec(level)
end;

begin
   randomize;
   GD:= InstallUserDriver('svga256',nil);GM:= 2;InitGraph(GD,GM,'bgi');
   for y:=0 to 255 do for x:=0 to 2 do pal[y][x]:=y div 4;
   reg.ax := $1012;
   reg.bx := 0;
   reg.cx := 256;
   reg.es := Seg(Pal);
   reg.dx := Ofs(Pal);
   intr($10,reg);
   putpixel(0,0,1+random(255));
   putpixel(xmax,0,1+random(255));
   putpixel(0,ymax,1+random(255));
   putpixel(xmax,ymax,1+random(255));
   level:=0;divide(0,0,xmax,ymax);
   for y:=0 to ymax do putpixel(2*xmax,y,getpixel(0,y));
   level:=0;divide(xmax,0,2*xmax,ymax);
   for x:=0 to 2*xmax do putpixel(x,2*ymax,getpixel(x,0));
   level:=0;divide(xmax,ymax,2*xmax,2*ymax);
   for y:=ymax to 2*ymax do putpixel(0,y,getpixel(2*xmax,y));
   level:=0;divide(0,ymax,xmax,2*ymax);
   readln;
   CloseGraph
end.

Graphic Driver: Plasma Fractals should have at least 256 colors; to get an adequate video mode, I use a non-standard Borland driver, svga256.bgi ((c) Jordan Hargraphix); I put that file in both the directory (c:\tp) containing my pascal programs, and the subdirectory c:\tp\bgi.

Variables:
xmax,ymax are the dimensions of the initial rectangle [A,B,D,C];
grad is the "graininess": the random part in the midpoint calculation. This random part is divided by the size of the currently calculated sub-rectangle, which is given by level.

Modules:
the function Adjust draws the Mean Point between (xa,ya) and (xb,yb);
the procedure Divide (calls itself recursively) draws a given rectangle;
the main program first generates a grey color palette, then calls routine Divide to draw the rectangle [A,B,D,C]; and finally (blue lines) draws 3 other rectangles arround [A,B,D,C] in order to get a periodic tile.

Optimisation: the program is optimised for speed, not for clearliness ;-)