Science 版 (精华区)
发信人: qpcwth (独翅鸟), 信区: Science
标 题: 《分形艺术》45
发信站: 哈工大紫丁香 (2001年11月03日18:18:09 星期六), 站内信件
第六章 复平面上的迭代
6.3 朱丽亚集
芒德勃罗集M记录的是整个区域上的c值情况，而朱丽亚集是取一固定的c 值后，观
察复平面上每一点(x，y)在迭代中的表现，并把结果记录下来。朱丽 亚集简称J集。因
此对于同一个迭代函数，既有M集又有J集。而且这 两种集合的确有密切联系 。
一个M集可以对应无数种J集，实际上M集就是所有J集的浓缩 。M集不同部位的形状
，反映了对应于该处的J集的形状，于是用M集 可以对J集进行分类，至少在计算机图形
学的层次上可以这样说。
计算J集时，初始迭代点就不能总取(0，0)了，而是要根据实际位置取实际的 (x，
y)坐标值。仍以迭代z→z^2+c为例说明。先取 定一个c值，比如c=(1.0221,0.2433)，迭
代关系化为
x→x^2-y^2+1.0221，
y→2xy+0.2433.
从屏幕左上角算起，逐行计算，一直算到屏幕右下角。当然，也可以不取整个屏幕
那么大，只选一个200×200的小区域做。
标色的原理与上面讲M集合时完全类似，此略。改变常数c的取值，可以得到 各式各
样的J集。
在源程序中，M集与J集的计算方法十分相似，只需改变两处语句就可以互换 为对方
，具体说明见表6.1。
表6.1从计算机程序上看M集与J集的区别与联系
初始迭代点 迭代关系 说明
M集 x0:=0;y0:=0； x1:=f(x,y)+p0；y1:=g(x,y)+q0； p_0,q_0不断变化
J集 x0:=p0;y0:=q0； x1:=f(x,y)+c_X；y1:=g(x,y)+c_Y；< /TD> c_(X)，c_(Y)固定不
变
下面给出作者编写的通用M类集合与J类集合逐级放大程序，此程序没有考虑 优化问
题，读者可仿照它们编写出更好的程序。
{Mandelbrot Set and Julia Set, Huajie, Fart1995.PAS}
uses graph,crt,vgafont,shelib;
label 10; label 20;
const
MaxIteNum=300;{实际上可以取得小些，如120}
Distance=40;
scalewidthheight=1；
var
bps,i,PixelNumber,gd,gm: integer;
pal : array［0..256*3-1］ of char;
f:file;
p:pointer;
b:char;
IteNumber,np,nq,kcolor,scale:integer;
r,p0,q0,x0,y0,x1,y1,temp,eps,lins1,lins2:real;
pmin,qmin,pmax,qmax,deltap,deltaq:real;
ppx,qqy,s:string;
Procedure PrepareFile(var Filename:string);
BEGIN
writeln('Please input image filename');
write('Filename = ');
readln(s);
IF s='' THEN REPEAT
writeln('You must input a true DOS filename!');
write('Filename = ');
readln(s) UNTIL s<>'';
Writeln('Press any key two times!');
Writeln('( IF you want TO abort press key: ESC )');
assign(f, s);
rewrite(f, 1);
blockwrite(f, win.width, sizeof(win.width));
blockwrite(f, win.height, sizeof(win.height));
blockwrite(f, bps, sizeof(bps));
blockwrite(f, pal, sizeof(pal));
END;
procedure Abort(Msg : string);
begin
Writeln(Msg, ': ', GraphErrorMsg(GraphResult));
Halt(1);
end;
BEGIN
Clrscr;
If RegisterBGIDriver(@VGADriver) <0 then Abort('EGA/VGA');
If RegisterBGIfont(@SansFont) < 0 then Abort('SansSerif'); {3}
If RegisterBGIfont(@TripFont) < 0 then Abort('Triplex'); {1}
writeln('Dr. Liu Huajie');
writeln('Research Center for Science & Society');
writeln('Peking University,100871');
writeln('Beijing,CHINA');
writeln('Email: liuhj@mccux0.mech.pku.eddu.cn');
writeln;
writeln('Now we start...');
writeln('Please press any key to continue');
readln;
Clrscr;TextColor(RED);
writeln('Hint: -3.0 < pmin < +2.0');
write('Please input pmin= ');
readln(pmin);
writeln;
writeln('Hint: -3.0 < qmin < +2.0');
Write('Please input qmin= '); readln(qmin);
qmin:=-qmin;{将屏幕坐标转化成正常坐标}
writeln;
writeln('Hint: 50 < number < 5000');
write('Please input the number of pixels on one direction=');
Readln(PixelNumber);
win.width :=PixelNumber;
win.height :=round(PixelNumber/scalewidthheight);
bps:=8;
scale:=PixelNumber div 400+1;
win.xstart:=0; win.ystart:=0;
pmax:=pmin+5.8;
qmax:=qmin+5.8*(win.height-1)/(win.width-1);
{Set Palette}
FOR i:=0 TO 20 do BEGIN
pal［3*i］:=chr(255-10*i);
pal［3*i+1］:=chr(255-10*i);
pal［3*i+2］:=chr(55+10*i) END;
FOR i:=21 TO 30 do BEGIN
pal［3*i］:=chr(11);
pal［3*i+1］:=chr(255-21*(i-21));
pal［3*i+2］:=chr(11) END;
FOR i:=31 TO 43 do BEGIN
pal［3*i］:=chr(255-11*(i-31));
pal［3*i+1］:=chr(20);
pal［3*i+2］:=chr(20) END;
FOR i:=44 TO 53 do BEGIN
pal［3*i］:=chr(60);
pal［3*i+1］:=chr(255-21*(i-44));
pal［3*i+2］:=chr(60) END;
…
FOR i:=201 TO 255 do BEGIN
pal［i*3］:=chr(255-20*(201-i));
pal［i*3+1］:=chr(0);
pal［i+3+2］:=chr(i) END;
{call procedure to write the initial part of image file}
PrepareFile(s);
gd:=VGA; gm:=VGAHI;
initgraph(gd,gm,'');
SetColor(6);SetTextStyle(1,0,2);
OutTextXY(400,5,'Fractal Image!');
SetColor(4);
OutTextXY(400,50,'Please wait...');
REPEAT
deltap:=(pmax-pmin)/(win.width-1);
deltaq:=(qmax-qmin)/(win.height-1);
FOR nq:=0 TO win.height-1 do
FOR np:=0 TO win.width-1 do BEGIN
p0:=pmin+np*deltap;
q0:=qmin+nq*deltaq;
IteNumber:=0;
x0:=0;y0:=0;
{for Julia set，x0:=p0;y0:=q0;}
10:
x1:=x0*x0-y0*y0+p0;
y1:=2*x0*y0+q0;
{for Julia set，use other values in the place of p0,q0}
{x1:=x0*x0*x0-3*x0*y0*y0-0.431183125;
y1:=3*x0*x0*y0-y0*y0*y0+0.01712875;}
{x1:=sin(x0*x0)-y0*y0/3+0.0136;
y1:=3*x0*y0+0.21; }
{lins1:=(1+eps)*(2*x0*cos(pi/10)-y0*sin(pi/10));
lins2:=(1+eps)*(y0*cos(pi/10)-x0*sin(pi/10));
x1:=x0*x0-y0*y0+lins1;
y1:=2*x0*y0+lins2+0.7;}
inc(IteNumber);
r:=x1*x1+y1*y1;
IF r>Distance THEN BEGIN
FOR i:=1 TO IteNumber do BEGIN
IF (IteNumber>i)and(IteNumber<=i+1) THEN
kcolor:=i END;
putpixel(np div scale,nq div scale,kcolor);
b:=chr(kcolor and 255);
blockwrite(f,b,sizeof(b)) END
ELSE IF (r<=Distance) and (IteNumber
grOK THEN Halt(1);
SetTextStyle(1,0,4);
str(pmin,ppx); str(-qmin,qqy);
OutTextXY(10,380,'px=');
Moveto(80,380);
OutText(ppx);
OutTextXY(10,420,'qy=');
Moveto(80,420);
OutText(qqy);
settextstyle(1,0,2);
OutTextXY(10,200,'Press ESC to EXIT!!');
SetViewPort(win.xstart,win.ystart,
win.width,win.height,ClipOn);
ch:=ReadKey;
if ch=esc then exit;
UNTIL (ch=esc);
Sound(500);Delay(200);nosound;
ClearViewPort;CloseGraph；
END.
上述主程序调用了一个单元Shelib.TPU，此单元的源程序如下：
{SHELIB.PAS,Huajie,1993}
Unit shelib;
interface
uses graph,crt,dos;
type view=Record
xstart,ystart,width,height:integer;
end;
var
gd,gm:integer;
c,d,i,j,xofinit0,xofinit1,yofinit0,yofinit1:integer;
K,Ratio,coef1,coef2,xs,ys,xe,ye,xsub,ysub:real;
xc,yc,step,dx,dy,dot:integer;
color,size:word;
x,y:integer;
s:string;
ch:char;
kofflag:boolean;
win:view;
const
esc=#27;up=#72;down=#80;left=#75;
right=#77;enter=#13;space=#32;
function FuncKey(var ch:char):boolean;
procedure cross(xc,yc:integer);
procedure movecross(var xc,yc:integer;step:integer);
procedure area(var xc,yc,dx,dy:integer);
implementation
function FuncKey(var ch:char):boolean;
begin
ch:=ReadKey;
if ch<>#0 then FuncKey:=False else begin
FuncKey:=True;
Ch:=ReadKey end;
end;
procedure cross(xc,yc:integer);
begin
if xc<0 then xc:=0;
if xc>GetMaxX then xc:=GetMaxX;
if yc<0 then yc:=0;
if yc>GetMaxY then yc:=GetMaxY;
if (GetMaxColor>1) then
Setcolor(GetMaxColor-1);
SetLineStyle(SolidLn,0,NormWidth);
SetWriteMode(Xorput);
Line(xc-5,yc,xc+5,yc);
Line(xc,yc-4,xc,yc-1);
Line(xc,yc+1,xc,yc+4);
end;
procedure movecross(var xc,yc:integer;step:integer);
var
x,y:integer;
begin
x:=xc;y:=yc;
repeat
kofflag:=FuncKey(ch);
if ch=esc then exit;
if (ch=up) or (ch=down) or (ch=left) or (ch=right) then
begin
case ch of
up: y:=y-step;
down: y:=y+step;
left: x:=x-step;
right: x:=x+step
end;
cross(xc,yc);
if x<0 then x:=0;
if x>GetMaxX then x:=GetMaxX;
if y<0 then y:=0;
if y>GetMaxY then y:=GetMaxY;
xc:=x;yc:=y;
cross(xc,yc)
end;
until (ch=enter)
end;
procedure area(var xc,yc,dx,dy:integer);
begin
setwritemode(Xorput);
rectangle(xc,yc,xc+dx,yc+dy);
step:=5;
repeat
kofflag:=FuncKey(ch);
if (ch=up) or (ch=down) or (ch=left) or (ch=right) then
begin
rectangle(xc,yc,xc+dx,yc+dy);
case ch of
up: dy:=dy-step;
down: dy:=dy+step;
left: dx:=dx-step;
right: dx:=dx+step
end;
end;
if (ch=up) or (ch=down) then dx:=round(dy*Ratio);
if (ch=left) or (ch=right) then dy:=round(dx/Ratio);
if ((xc+dx)>win.width) then
dx:=dx-step;
if ((xc+dx)<0) then
dx:=dx+step;
if ((yc+dy)>win.height) then
dy:=dy-step;
if ((yc+dy)<0) then
dy:=dy+step;
rectangle(xc,yc,xc+dx,yc+dy)
until (ch=enter);
rectangle(xc,yc,xc+dx,yc+dy);
end;
begin
end.
使用时先将SHELIB.PAS编译成SHELIB.TPU，然后再编译Fart1995.PAS得到Fart1995
.EXE。你会注意到，Fart1995.PAS中的迭代函数实际上可以取任何你能想象得出的函数
。距离函数可以用其他方式定义。因而上述程序完全适用于下一节要说的广义M集和广义
J 集的计算。
使用上述程序的好处是，M集与J集统一于一个程序，对任意函数都如此。此 外，所
作图形的大小不受显示卡类型限制，因为我们直接往文件中写数据，屏幕只起一个监 视
进程的作用。
采用上述方法得到的文件只是初级的非标准图形文件，一般软件还不能直接读取。
使用我们编写的RAWTOGIF.EXE工具软件，可以将它们转化成标准的GIF格式的文件，用法
是：
RAWTOGIF myfile newfile.GIF
其中myfile是初级文件，newfile.GIF是最终生成的标准GIF格式文件。RAWTOGIF.EXE可
从与本书配套的光盘上找到，如果没买到光盘也可以自己编译一下如下源程序(本程序由
罗军编写)。
//此程序用C语言编写，共包括三部分：GIF.H,GIFSAVE.C和GIFMAIN.C
//以下是GIF.H
#ifndef XGIFDATA
#include
#include
#include
#include
#include
#include
#include
typedef struct tagGLOBAL{
char version［3］;
char subversion［3］;
int screen_width;
int screen_height;
unsigned char global_flag;
unsigned char bg_color;
unsigned char zero;
}GLOBAL;
typedef struct tagLOCAL{
unsigned char comma;
int image_top;
int image_left;
int image_width;
int image_height;
unsigned char local_flag;
unsigned char first_char;
}LOCAL;
extern void pack_gif(FILE * , FILE *);
#endif // XGIFDATA
//以下是GIFSAVE.C
#include "gif.h"
void pack_gif(FILE * , FILE *);
void pack_image(FILE *, FILE *, int, int, int);
void pack_init(int clear);
short lookup_ct(short code, unsigned char thischar);
void flush(void);
void putcode(short code);
FILE *fpout;
int index = 0, pixmask = 0xff;
int first_ch;
int code, oldcode;
int nextcode, nextlim;
int remcnt, reqcnt;
unsigned char rem;
int EOI, CLEAR;
unsigned char DataBuff［4096 * 5］;
int *CTlink = (int*)DataBuff;
unsigned char *CTfirst = DataBuff + 4096 * 2;
int *CTnext = (int *)(DataBuff + 4096 * 3);
unsigned char strbuf［256］;
void pack_gif(FILE * SourceRaw, FILE * DestGif)
{
char palette［256*3］;
int width, height, bps;
GLOBAL global;
LOCAL local;
if (fread(&width, sizeof(width), 1, SourceRaw) != 1){
fclose(SourceRaw);
fclose(DestGif);
fprintf(stderr, "Reading from source file failed 1.\n");
exit(-1);
}
if (fread(&height, sizeof(height), 1, SourceRaw) != 1){
fclose(SourceRaw);
fclose(DestGif);
fprintf(stderr, "Reading from source file failed 2.\n");
exit(-1);
}
if (fread(&bps, sizeof(bps), 1, SourceRaw) != 1){
fclose(SourceRaw);
fclose(DestGif);
fprintf(stderr, "Reading from source file failed 3.\n");
exit(-1);
}
if (fread(&palette, sizeof(palette), 1, SourceRaw) != 1){
fclose(SourceRaw);
fclose(DestGif);
fprintf(stderr, "Reading from source file failed 4.\n");
exit(-1);
}
strncpy(global.version, "GIF", 3);
strncpy(global.subversion, "87a", 3);
global.screen_width = width;
global.screen_height = height;
global.bg_color = 0;
global.zero = 0;
global.global_flag = 0x80 | (bps - 1);
local.comma = ',';
local.image_top = 0;
local.image_left = 0;
local.image_width = width;
local.image_height = height;
local.local_flag = global.global_flag & 0x07;
local.first_char = bps;
if (fwrite(&global, 1L, sizeof(global), DestGif) != (size_t)sizeof(global)
){
fclose(SourceRaw);
fclose(DestGif);
fprintf(stderr, "Save GIF: error writing file.\n");
return;
}
if (fwrite(palette, 1L, (1<
nextlim){
if (reqcnt == 12){
putcode(oldcode);
putcode(CLEAR);
pack_init(CLEAR);
}else{
reqcnt++;
nextlim = nextlim << 1;
if (reqcnt == 12) nextlim--;
}
}
} // for j
} // for i
putcode(oldcode);
putcode(EOI);
flush();
fputc(0, fpout);
fputc(0x3b, fpout);
}
void pack_init(int clear)
{
register int i;
nextcode = clear + 2;
oldcode = -1;
for(i = 0; i < clear; i++){
CTfirst［i］ = i;
for(i = 0; i < 4096; i++){
CTlink［i］ = -2;
CTnext［i］ = -1;
}
reqcnt = first_ch + 1; // first_ch == bps -- BitsPerPixel.
nextlim = 1 << reqcnt;
}
short lookup_ct(short code, unsigned char thischar)
{
if (code == -1) return(unsigned int) thischar;
if (CTlink［code］ == -2){
CTlink［code］ = nextcode;
}else{
code = CTlink［code］;
while((CTnext［code］ != -1) && (CTfirst［code］ != thischar))
code = CTnext［code］;
if (CTfirst［code］ == thischar) return code;
CTnext［code］ = nextcode;
}
CTfirst［nextcode++］ = thischar;
return -1;
}
void putcode(short code)
{
int allcnt, used, codelen;
codelen = reqcnt;
allcnt = remcnt + reqcnt;
while(allcnt >= 8){
if (remcnt > 0){
used = 8 - remcnt;
strbuf［index++］ = rem | (unsigned char)
code = code >> used;
codelen -= used;
remcnt = 0;
}else{
strbuf［index++］ = (unsigned char)(code & pixmask);
codelen -= 8;
code = code >> 8;
}
allcnt -= 8;
if (index == 256){
strbuf［0］ = 0xff;
fwrite(strbuf, index, 1, fpout);
index = 1;
}
}
if (remcnt == 0){
rem = code;
remcnt = codelen;
}else{
rem |= (code << remcnt);
remcnt += codelen;
}
}
void flush()
{
if (remcnt != 0) strbuf［index++］ = rem;
if (index == 1) return;
strbuf［0］ = (unsigned char) (index - 1);
fwrite(strbuf, index, 1, fpout);
index = 1;
}
//以下是GIFMAIN.C
#include "gif.h"
int main(int argc, char *argv［］)
{
char SourceName［256］, DestName［256］;
FILE *SourceRaw, *DestGif;
if (argc < 3){
printf("Usage : %s SOURCE-RAW DEST-GIF\n", argv［0］);
exit(0);
}
SourceRaw = fopen(argv［1］, "rb");
if (SourceRaw == NULL){
printf("Can not open source file %s.\n", argv［1］);
exit(0);
}
DestGif = fopen(argv［2］, "wb");
if (DestGif == NULL){
fclose(SourceRaw);
printf("Can not create destination file %s.\n", argv［2］);
exit(0);
}
pack_gif(SourceRaw, DestGif);
return 0;
}
--
心事浩茫连广宇,于无声处听惊雷
※ 来源:·哈工大紫丁香 bbs.hit.edu.cn·[FROM: 202.118.229.154]
Powered by KBS BBS 2.0 (http://dev.kcn.cn)
页面执行时间：211.411毫秒