extern void cbar(int sx, int sy, int wx, int wy);
extern void screen_range(double xs, double xe, double ys, double ye);
extern void draw_img_line(double *p1, double *p2);
extern void draw_img_line3(double *p1, double *p2);
extern void draw_img_triangle(double tr[3][2], int col, bool line);
extern void draw_img_rect(double sx, double ex, double h, int col);
extern void init_img();
extern void write_img(char *fn);
extern void draw_img_ball(double x0, double y0, double z0, double r, int hue, double lig, double amb);
extern void draw_img_triangle_tp(double tr[3][3], int col, double lig, double tra, double damp_l);
#define IMG_COL_WHITE (-1)
#define IMG_COL_GRAY (-2)
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "vis.h"
#define XSIZE 800
#define YSIZE 600
static unsigned char buf[3*XSIZE*YSIZE];
static double zbuf[XSIZE*YSIZE];
static double xmin,xmax,ymin,ymax;
static void putcol(unsigned char *buf, int col, double lig)
{
double d=lig;
if(col == IMG_COL_GRAY){
buf[0]=(unsigned char)(d*128);
buf[1]=(unsigned char)(d*128);
buf[2]=(unsigned char)(d*128);
return;
}
else
if(col == IMG_COL_WHITE){
buf[0]=(unsigned char)(d*255);
buf[1]=(unsigned char)(d*255);
buf[2]=(unsigned char)(d*255);
return;
}
if(col<64)
{
buf[0]=0;
buf[1]=(unsigned char)(col*4*d);
buf[2]=(unsigned char)(255*d);
return;
}
col-=64;
if(col<64){
buf[0]=0;
buf[1]=(unsigned char)(255*d);
buf[2]=(unsigned char)((63-col)*4*d);
return;
}
col-=64;
if(col<64)
{
buf[0]=(unsigned char)(col*4*d);
buf[1]=(unsigned char)(255*d);
buf[2]=0;
return;
}
col-=64;
if(col<64)
{
buf[0]=(unsigned char)(255*d);
buf[1]=(unsigned char)((63-col)*4*d);
buf[2]=0;
return;
}
col-=64;
if(col<64)
{
buf[0]=(unsigned char)(255*d);
buf[1]=0;
buf[2]=(unsigned char)(col*4*d);
return;
}
col-=64;
if(col<64)
{
buf[0]=(unsigned char)((63-col)*4*d);
buf[1]=0;
buf[2]=(unsigned char)(255*d);
return;
}
}
static void puthl(unsigned char *buf, double lig)
{
int h = (int)(255*lig);
if (buf[0]<h) buf[0]=h;
if (buf[1]<h) buf[1]=h;
if (buf[2]<h) buf[2]=h;
}
extern void cbar(int sx, int sy, int wx, int wy)
{
int x,y;
for(x=0;x<wx;x++)
{
int col=(int)(255.0*x/wx);
for(y=0;y<wy;y++)
putcol(&buf[3*(x+y*XSIZE)], col, 1.0);
}
}
void draw_line(int x1, int y1, double z1, int x2, int y2, double z2, double d)
{
int i,x,y;
int sx=0;
int sy=0;
int dx = abs(x1-x2);
int dy = abs(y1-y2);
int xx1=0,xx2=0,yy1=0,yy2=0;
if (x2<x1) {sx = -1; xx1=x2;xx2=x1;}
else if(x2>x1) {sx = 1; xx1=x1;xx2=x2;}
if (y2<y1) { sy=-1;yy1=y2;yy2=y1;;}
else if(y2>y1) { sy=1; yy1=y1;yy2=y2;}
if(dx==0 && dy==0)return;
if (dx>dy)
{
for (i=0;i<=dx; i++){
x = xx1 + i;
if (sx == 1)
y = y1 + (dy*sy*i)/dx;
else
y = y2 - (dy*sy*i)/dx;
if(x<0 || y<0 || x>=XSIZE || y>=YSIZE)continue;
double c= (x-x2)*1.0/(x1-x2);
double z0 = c*z1+(1-c)*z2;
if ( zbuf[(x+y*XSIZE)] < z0) continue;
putcol(&buf[3*(x+y*XSIZE)], IMG_COL_WHITE, d);
zbuf[(x+y*XSIZE)]=z0;
}
}
else
{
for (i=0;i<=dy; i++){
y = yy1 + i;
if (sy == 1)
x = x1 + (dx*sx*i)/dy;
else
x = x2 - (dx*sx*i)/dy;
if(x<0 || y<0 || x>=XSIZE || y>=YSIZE)continue;
double c= (y-y2)*1.0/(y1-y2);
double z0 = c*z1+(1-c)*z2;
if ( zbuf[(x+y*XSIZE)] < z0) continue;
putcol(&buf[3*(x+y*XSIZE)], IMG_COL_WHITE, d);
}
}
}
void draw_line(int x1, int y1, int x2, int y2, double d)
{
draw_line(x1, y1, 0.0, x2, y2, 0.0, d);
}
void screen_range(double xs, double xe, double ys, double ye)
{
double xlen = xe-xs;
double ylen = ye-ys;
double xm = (xs+xe)/2;
double ym = (ys+ye)/2;
if (xlen/XSIZE < ylen/YSIZE)
{
xlen = ylen/YSIZE*XSIZE;
}
else
{
ylen = xlen/XSIZE*YSIZE;
}
xmin= xm-xlen/2;
xmax= xm+xlen/2;
ymin= ym-ylen/2;
ymax= ym+ylen/2;
}
void screen_cord(double dx, double dy, double *sx, double *sy)
{
*sx=XSIZE*( (dx-xmin)/(xmax-xmin));
*sy=YSIZE*( (dy-ymin)/(ymax-ymin));
}
void real_cord(int sx, int sy, double *rx, double *ry)
{
*rx = xmin + (xmax-xmin)/XSIZE * sx;
*ry = ymin + (ymax-ymin)/YSIZE * sy;
}
static void clr_buf()
{
int i,j;
for(i=0;i<XSIZE*YSIZE*3;i++)buf[i]=64;
for(i=0;i<XSIZE*YSIZE;i++)zbuf[i]=9999999.9;
}
void draw_img_line(double *p1, double *p2)
{
double sx1,sy1;
double sx2,sy2;
screen_cord(p1[0],p1[1], &sx1, &sy1);
screen_cord(p2[0],p2[1], &sx2, &sy2);
draw_line((int)sx1,(int)sy1, (int)sx2,(int)sy2, 1.0);
}
void draw_img_line3(double *p1, double *p2)
{
double sx1,sy1;
double sx2,sy2;
screen_cord(p1[0],p1[1], &sx1, &sy1);
screen_cord(p2[0],p2[1], &sx2, &sy2);
draw_line((int)sx1,(int)sy1, p1[2], (int)sx2,(int)sy2, p2[2], 1.0);
}
void draw_img_triangle(double tr[3][2], int col, bool line)
{
int i,x,y;
double sx[3],sy[3];
double cx[3],cy[3],c0[3];
double dd;
int xmax,xmin,ymax,ymin;
xmax=ymax=0;
xmin=XSIZE-1;
ymin=YSIZE-1;
for(i=0;i<3;i++)
{
screen_cord(tr[i][0], tr[i][1], &sx[i], &sy[i]);
sx[i]=1.0*(int)sx[i];
sy[i]=1.0*(int)sy[i];
int xx=(int)sx[i];
int yy=(int)sy[i];
if(xmax<xx)xmax=xx;
if(ymax<yy)ymax=yy;
if(xmin>xx)xmin=xx;
if(ymin>yy)ymin=yy;
}
for(i=0;i<3;i++)
{
int i2=(i+1)%3;
int i3=(i+2)%3;
double cx0= sx[i2]-sx[i];
double cy0= sy[i2]-sy[i];
double c00 = -sx[i]*cy0 + sy[i]*cx0;
double d= sx[i3]*cy0-sy[i3]*cx0 +c00;
if(d>0){ cx[i]=cy0; cy[i]=-cx0; c0[i]=c00;}
else { cx[i]=-cy0; cy[i]=cx0; c0[i]=-c00;}
}
if(xmin<0)xmin=0;
if(ymin<0)ymin=0;
if(xmax>XSIZE-1)xmax=XSIZE-1;
if(ymax>YSIZE-1)ymax=YSIZE-1;
for(x=xmin;x<=xmax;x++){
for(y=ymin;y<=ymax;y++){
int flag=1;
for(i=0;i<3;i++)if (x*cx[i]+y*cy[i]+c0[i]<0.0) flag=0;
if(flag==0) continue;
{
putcol(&buf[(x+y*XSIZE)*3], col, 1.0);
}
}}
if(!line)return;
draw_line((int)sx[1],(int)sy[1], (int)sx[2],(int)sy[2], 1.0);
draw_line((int)sx[0],(int)sy[0], (int)sx[2],(int)sy[2], 1.0);
draw_line((int)sx[0],(int)sy[0], (int)sx[1],(int)sy[1], 1.0);
}
void draw_img_triangle2(double tr[4][2], int col, bool line)
{
double tr2[3][2];
int i;
for(i=0;i<3;i++) { tr2[i][0] = tr[i][0]; tr2[i][1]=tr[i][1];}
draw_img_triangle(tr2,col, line);
for(i=0;i<3;i++) { tr2[i][0] = tr[i+1][0]; tr2[i][1]=tr[i+1][1];}
draw_img_triangle(tr2, col, line);
}
void draw_img_rect(double sx, double ex, double h, int col)
{
int i,x,y;
int xmax,xmin,ymax,ymin;
double xx,yy;
screen_cord(sx,0, &xx, &yy);
xmin = (int)xx;
screen_cord(ex,0, &xx, &yy);
xmax = (int)xx;
if(xmin<0)xmin=0;
if(xmax>=XSIZE)xmax=XSIZE-1;
ymin = (int)((YSIZE-1)*(1-h));
ymax = YSIZE-1;
for(x=xmin;x<=xmax;x++){
for(y=ymin;y<=ymax;y++){
putcol(&buf[(x+y*XSIZE)*3], col, 1.0);
}}
}
void init_img()
{
char buf[99];
int i;
double theta;
clr_buf();
}
void write_img(char *fn)
{
FILE *fp;
fp=fopen(fn,"wb");
fprintf(fp,"P6\n%d %d\n255\n",XSIZE,YSIZE);
fwrite(buf,1,3*XSIZE*YSIZE,fp);
fclose(fp);
}
void draw_img_triangle_tp(double tr[3][3], int col, double lig, double tra, double damp_l)
{
int i,j,x,y;
double tmp[3];
double sx[3],sy[3];
double cx[3],cy[3],c0[3];
double sz[3];
double dd;
unsigned char rgb[3];
putcol(rgb, col, 0.2+lig*0.8);
double a,b,c,d,det;
int xmax,xmin,ymax,ymin;
xmax=ymax=0;
xmin=XSIZE-1;
ymin=YSIZE-1;
dd=lig;
for(i=0;i<3;i++)
{
int xx,yy;
screen_cord(tr[i][0], tr[i][1], &sx[i], &sy[i]);
sx[i]=1.0*(int)sx[i];
sy[i]=1.0*(int)sy[i];
sz[i] = tr[i][2];
xx=(int)sx[i];
yy=(int)sy[i];
if(xmax<xx)xmax=xx;
if(ymax<yy)ymax=yy;
if(xmin>xx)xmin=xx;
if(ymin>yy)ymin=yy;
}
a = sx[1]-sx[0];
b = sx[2]-sx[0];
c = sy[1]-sy[0];
d = sy[2]-sy[0];
det = a*d - b*c;
for(i=0;i<3;i++)
{
int i2=(i+1)%3;
int i3=(i+2)%3;
double cx0= sx[i2]-sx[i];
double cy0= sy[i2]-sy[i];
double c00 = -sx[i]*cy0 + sy[i]*cx0;
double d= sx[i3]*cy0-sy[i3]*cx0 +c00;
if(d>0){ cx[i]=cy0; cy[i]=-cx0; c0[i]=c00;}
else { cx[i]=-cy0; cy[i]=cx0; c0[i]=-c00;}
}
if(xmin<0)xmin=0;
if(ymin<0)ymin=0;
if(xmax>XSIZE-1)xmax=XSIZE-1;
if(ymax>YSIZE-1)ymax=YSIZE-1;
for(x=xmin;x<=xmax;x++){
for(y=ymin;y<=ymax;y++){
int flag=0;
for(i=0;i<3;i++)if (x*cx[i]+y*cy[i]+c0[i]<0.0) flag=1;
if(flag==1)continue;
#if 1
double x0 = x-sx[0];
double y0 = y-sy[0];
double s = ( d*x0 -b*y0)/det;
double t = (-c*x0 +a*y0)/det;
double zz = sz[0] + s*(sz[1]-sz[0]) + t*(sz[2]-sz[0]);
#endif
double znow=zbuf[x+y*XSIZE];
if(zz<znow)
{
double dz = znow - zz;
if (dz>20) dz=20;
unsigned char *rgb0 =&buf[(x+y*XSIZE)*3];
int k;
for(k=0;k<3;k++)
{
double dat = rgb[k] + (1-tra)*rgb0[k] * exp(-dz/damp_l);
rgb0[k] = (unsigned char)dat;
}
double hl = 1.0-lig*lig;
hl = exp(-hl*10);
puthl(rgb0, hl);
zbuf[x+y*XSIZE]=zz;
}
}}
draw_line( (int)sx[0], (int)sy[0], tr[0][2]-1e-2, (int)sx[1], (int)sy[1], tr[1][2]-1e-2, lig);
draw_line( (int)sx[1], (int)sy[1], tr[1][2]-1e-2, (int)sx[2], (int)sy[2], tr[2][2]-1e-2, lig);
draw_line( (int)sx[0], (int)sy[0], tr[0][2]-1e-2, (int)sx[2], (int)sy[2], tr[2][2]-1e-2, lig);
}
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "btree.h"
#include "vis.h"
#include "mymath.h"
#define RMAX 2.4
#define DR (0.12)
#define MAXV 1000000
int nv;
class Vtx *vtx[MAXV];
#if 0
#endif
#define AddT(v1,v2,v3) { Tri *t = new Tri(v1,v2,v3);btt->Insert(t->vidx, t);}
class Vtx
{
public:
int idx;
int maxf;
double Pos[3];
double nv0;
int nn;
class Vtx *nei[6];
double nvec[3];
int nnvec;
void AddNei(Vtx *v2)
{
if(nn==6)
{
printf("too much nn\n");exit(1);
}
nei[nn++]=v2;
}
void NormPos()
{
double x=Pos[0];
double y=Pos[1];
double z=Pos[2];
double l=sqrt(x*x+y*y+z*z);
Pos[0]/=l;
Pos[1]/=l;
Pos[2]/=l;
}
Vtx(double x, double y, double z, bool norm)
{
Pos[0]=x;
Pos[1]=y;
Pos[2]=z;
if (norm)NormPos();
nn=0;
};
Vtx(double *pos0, bool norm)
{
Pos[0]=pos0[0];
Pos[1]=pos0[1];
Pos[2]=pos0[2];
if (norm)NormPos();
nn=0;
};
Vtx(Vtx *v1, Vtx *v2, bool norm)
{
int k;
for(k=0;k<3;k++)
Pos[k]= (v1->Pos[k] + v2->Pos[k])*0.5;
if (norm)NormPos();
nn=0;
};
void dump()
{
printf("%f %f %f %d\n", Pos[0],Pos[1],Pos[2],idx);
}
};
void AddV(class Vtx *v)
{
if (nv==MAXV)
{
printf("vmax\n");exit(1);
}
v->idx = nv;
vtx[nv++]=v;
}
class Edge
{
public:
int idx;
Vtx *v1;
Vtx *v2;
Vtx *nV;
int vidx[2];
double len;
void CalLen()
{
len = sqrt(V_DIST3(v1->Pos, v2->Pos));
}
Edge(Vtx *v10, Vtx *v20)
{
v1=v10;
v2=v20;
nV=NULL;
CalLen();
vidx[0]=v10->idx;
vidx[1]=v20->idx;
};
};
Btree *bte;
Btree *btt;
class Tri
{
public:
int idx;
int Lv;
double lmax;
double nvec[3];
class Vtx *pVtx[3];
class Tri *Child[4];
class Vtx *nVtx[3];
int vidx[3];
class Edge *pEdge[3];
Tri(Vtx *v1, Vtx *v2, Vtx *v3)
{
pVtx[0]=v1;
pVtx[1]=v2;
pVtx[2]=v3;
int k;
for(k=0;k<3;k++)vidx[k]=pVtx[k]->idx;
int vv[2];
for(k=0;k<3;k++)
{
int k1=k;
int k2=(k+1)%3;
vv[0]= vidx[k1];
vv[1]= vidx[k2];
Edge *e = (Edge *)bte->Find(vv);
if(e==NULL)
{
e = new Edge( pVtx[k1], pVtx[k2]);
bte->Insert(e->vidx, e);
}
pEdge[k]=e;
}
double nv1[3],nv2[3];
for(k=0;k<3;k++)
{
nv1[k]= v2->Pos[k] - v1->Pos[k];
nv2[k]= v3->Pos[k] - v1->Pos[k];
V_OPROD3(nv1, nv2, nvec);
V_NORMALIZE3(nvec);
}
};
~Tri()
{
#if 0
#endif
};
void CalcNv()
{
int k;
double nv1[3],nv2[3];
for(k=0;k<3;k++)
{
nv1[k]= pVtx[1]->Pos[k] - pVtx[0]->Pos[k];
nv2[k]= pVtx[2]->Pos[k] - pVtx[0]->Pos[k];
V_OPROD3(nv1, nv2, nvec);
}
V_NORMALIZE3(nvec);
double p = V_IPROD3(nvec, pVtx[0]->Pos);
if(p<0.0)
{
for(k=0;k<3;k++) nvec[k]=-nvec[k];
}
}
bool Div(bool norm)
{
Vtx *v1 = pVtx[0];
Vtx *v2 = pVtx[1];
Vtx *v3 = pVtx[2];
int k;
int ndiv=0;
int ks=-1, kl=-1;
for(k=0;k<3;k++)
{
Edge *e = pEdge[k];
if(e->nV != NULL)
{
kl=k;
ndiv++;
}
else
ks=k;
}
if (ndiv==0) return false;
if (ndiv==3)
{
for(k=0;k<3;k++)
{
Edge *e = pEdge[k];
if (e->nV == NULL)
{
exit(1);
e->nV = new Vtx( e->v1, e->v2, norm);
AddV(e->nV);
}
}
Vtx *v12 = pEdge[0]->nV;
Vtx *v23 = pEdge[1]->nV;
Vtx *v31 = pEdge[2]->nV;
AddT(v1,v12, v31);
AddT(v2,v12, v23);
AddT(v3,v23, v31);
AddT(v12,v23,v31);
return true;
}
else if (ndiv==1)
{
Edge *e = pEdge[kl];
if (e->nV == NULL)
{
exit(1);
e->nV = new Vtx( e->v1, e->v2, norm);
AddV(e->nV);
}
int k1 = kl;
int k2 = (kl+1)%3;
int k3 = (kl+2)%3;
Vtx *v1 = pVtx[k1];
Vtx *v2 = pVtx[k2];
Vtx *v3 = pVtx[k3];
Vtx *vn = e->nV;
AddT(v3, v1, vn);
AddT(v3, v2, vn);
return true;
}
else if (ndiv==2)
{
int k1 = ks;
int k2 = (ks+1)%3;
int k3 = (ks+2)%3;
Edge *e1 = pEdge[k2];
Edge *e2 = pEdge[k3];
Vtx *v1 = pVtx[k1];
Vtx *v2 = pVtx[k2];
Vtx *v3 = pVtx[k3];
Vtx *vn1 = e1->nV;
Vtx *vn2 = e2->nV;
AddT(v1, v2, vn1);
AddT(vn1, vn2, v1);
AddT(v3, vn1, vn2);
return true;
}
return false;
}
};
void RndDir(double v[3])
{
double z = rand()*2.0/RAND_MAX-1;
double th = rand()*2*M_PI/RAND_MAX;
double r=sqrt(1-z*z);
v[2]=z;
v[0]=r*cos(th);
v[1]=r*sin(th);
}
void Grow(double dw, double dl)
{
double nvec[3];
RndDir(nvec);
double zmax=0.0;
int i;
int imin;
for(i=0;i<nv;i++)
{
Vtx *v = vtx[i];
double zz = nvec[0]*v->Pos[0] + nvec[1]*v->Pos[1] +nvec[2]*v->Pos[2] ;
v->nv0=zz;
if(zz>zmax)
{
zmax=zz;
imin=i;
}
}
double *p0 = vtx[imin]->Pos;
double r2 = V_NORM3(p0);
if(r2>RMAX*RMAX) return;
for(i=0;i<nv;i++)
{
Vtx *v = vtx[i];
int k;
#if 0
#else
double r2 = V_DIST3(v->Pos, p0)/(dl*dl);
if (r2< 4.0)
{
double rr = V_NORM3(v->Pos);
if (rr>RMAX*RMAX)continue;
double dw2 = dw*exp(-r2);
int k;
double nx=v->nvec[0];
double ny=v->nvec[1];
double nz=v->nvec[2];
nx=nx*nx;
ny=ny*ny;
nz=nz*nz;
double r2 = nx+ny+nz;
double cub = 0.4+(nx*ny*ny +ny*nx*nx + ny*nz*nz + nz*ny*ny + nz*nx*nx + nx*nz*nz)/(r2*r2*r2)*0.4;
dw2 *= cub;
for(k=0;k<3;k++)
v->Pos[k] += dw2* v->nvec[k];
}
#endif
}
}
void Cut(double r0)
{
int i;
for(i=0;i<nv;i++)
{
Vtx *v = vtx[i];
double r2 = V_NORM3(v->Pos);
if (r2>RMAX*RMAX)
{
double rr = RMAX/sqrt(r2);
int k;
for(k=0;k<3;k++) v->Pos[k] *=rr;
}
}
}
void check_max()
{
int i;
for(i=0;i<nv;i++)
{
vtx[i]->maxf=1;
}
BTFOR(bte)
{
Edge *e = (Edge *)(bte->LoopVar());
double r1 = V_NORM3(e->v1->Pos);
double r2 = V_NORM3(e->v2->Pos);
if(r1>r2) e->v2->maxf = 0;
else e->v1->maxf=0;
}
for(i=0;i<nv;i++)
{
if (vtx[i]->maxf ==1) printf("%f %d\n", sqrt(V_NORM3(vtx[i]->Pos)), i);
}
}
int tstart, tend;
void CalcNv()
{
int i;
for(i=0;i<nv;i++)
{
Vtx *v = vtx[i];
v->nnvec=0;
int k;
for(k=0;k<3;k++) v->nvec[k]=0.0;
}
BTFOR(btt)
{
Tri *t = (Tri *)(btt->LoopVar());
t->CalcNv();
int j;
for(j=0;j<3;j++)
{
Vtx *v = t->pVtx[j];
v->nnvec++;
int k;
for(k=0;k<3;k++) v->nvec[k] += t->nvec[k];
}
}
for(i=0;i<nv;i++)
{
Vtx *v = vtx[i];
int k;
for(k=0;k<3;k++) v->nvec[k]/= v->nnvec;
}
}
void InitMesh()
{
bte = new Btree(2, true);
btt = new Btree(3, true);
Vtx *vx1 = new Vtx(1.0,0.0,0.0, true);
Vtx *vx2 = new Vtx(-1.0,0.0,0.0, true);
Vtx *vy1 = new Vtx(0.0, 1.0,0.0, true);
Vtx *vy2 = new Vtx(0.0, -1.0,0.0, true);
Vtx *vz1 = new Vtx(0.0, 0.0, 1.0, true);
Vtx *vz2 = new Vtx(0.0, 0.0, -1.0, true);
AddV(vx1);
AddV(vx2);
AddV(vy1);
AddV(vy2);
AddV(vz1);
AddV(vz2);
AddT(vx1, vy1, vz1);
AddT(vx1, vy1, vz2);
AddT(vx1, vy2, vz1);
AddT(vx1, vy2, vz2);
AddT(vx2, vy1, vz1);
AddT(vx2, vy1, vz2);
AddT(vx2, vy2, vz1);
AddT(vx2, vy2, vz2);
}
int RefineMesh1(double dx, bool norm)
{
int inc=0;
BtList tlist(btt);
int i, k;
printf("#begin edge len check\n");
BTFOR(bte)
{
Edge *e = (Edge *)(bte->LoopVar());
e->CalLen();
e->nV=NULL;
if (e->len > dx)
{
e->nV = new Vtx( e->v1, e->v2, norm);
AddV(e->nV);
}
}
printf("#begin tetra div\n");
for(i=0;i<tlist.NData;i++)
{
Tri *t = (Tri *)(tlist.Ptr[i]);
if (t->Div(dx))
{
btt->Delete(t->vidx);
delete t;
inc++;
};
}
BTFOR(bte)
{
Edge *e = (Edge *)(bte->LoopVar());
if (e->nV)
{
bte->Delete(e->vidx);
delete e;
}
}
return inc;
}
void RefineMesh(double dx, bool norm)
{
while(1)
{
int inc=RefineMesh1(dx, norm);
if(inc==0) break;
}
CalcNv();
}
void w_img(int st)
{
double tr[3][3];
init_img();
screen_range( -2.,2., -2., 2.);
double th = st*2 *3.14159*2 / 2800.0;
double cc=cos(th);
double ss=sin(th);
BTFOR(btt)
{
Tri *t = (Tri *)(btt->LoopVar());
int k,kk;
for(k=0;k<3;k++)
{
double x0 = t->pVtx[k]->Pos[0];
double y0 = t->pVtx[k]->Pos[1];
double z0 = t->pVtx[k]->Pos[2];
tr[k][1] = y0;
tr[k][0] = cc*x0 + ss*z0;
tr[k][2] = cc*z0 - ss*x0;
}
double lig = fabs(cc*t->nvec[2] -ss*t->nvec[0]);
draw_img_triangle_tp( tr, 64*2, lig, 1.0, 100.0);
}
char fn[32];
sprintf(fn,"ppm/%04d.ppm", st);
write_img(fn);
}
int main(int argc, char **argv)
{
nv=0;
srandom(12345);
double lw = atof(argv[1]);
int i, j;
InitMesh();
CalcNv();
RefineMesh( DR, true);
printf("# Begin growth\n");
double cdep = lw*lw;
double mul=atof(argv[2]);
int maxst = 3000;
#if 1
for(i=0;i<maxst ;i++)
{
w_img(i);
Grow(1e-2, lw);
if (i%10==0)RefineMesh( DR, false);
}
#endif
printf("#Step %d\n",maxst);
check_max();
}