Konpeito

// vis.h
extern void cbar(int sx, int sy, int wx, int wy);
//extern void draw_line(int x1, int y1, int x2, int y2, double d);
//extern void draw_line(int x1, int y1, double z1, int x2, int y2, double z2, double d);
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);

// transparent tri
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)

//vis.cc

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "vis.h"

//image size
#define XSIZE 800
#define YSIZE 600
//image pixel buffer
static unsigned char buf[3*XSIZE*YSIZE];
static double zbuf[XSIZE*YSIZE];
static double xmin,xmax,ymin,ymax;

// col=hue (0 to 64*6) , lig=brightness (0.0 to 1.0)
static void putcol(unsigned char *buf, int col, double lig)
{
    double d=lig;
    //double d=1.0;
    //if(d<0.0) d=0.0;
    //d=(0.3+d)/1.31;

    if(col == IMG_COL_GRAY){
	//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){
	    //White
	    buf[0]=(unsigned char)(d*255);
	    buf[1]=(unsigned char)(d*255);
	    buf[2]=(unsigned char)(d*255);
	    return;
	}

    // B to C
    if(col<64)
    {
	buf[0]=0;
	buf[1]=(unsigned char)(col*4*d);
	buf[2]=(unsigned char)(255*d);
	return;
    }

    col-=64;
    //C to G
    if(col<64){
	buf[0]=0;
	buf[1]=(unsigned char)(255*d);
	buf[2]=(unsigned char)((63-col)*4*d);
	return;
    }

    col-=64;
    //G to Y
    if(col<64)
    {
	buf[0]=(unsigned char)(col*4*d);
	buf[1]=(unsigned char)(255*d);
	buf[2]=0;
	return;
    }

    col-=64;
    // Y to R
    if(col<64)
    {
	buf[0]=(unsigned char)(255*d);
	buf[1]=(unsigned char)((63-col)*4*d);
	buf[2]=0;
	return;
    }
    col-=64;
    // R to M
    if(col<64)
    {
	buf[0]=(unsigned char)(255*d);
	buf[1]=0;
	buf[2]=(unsigned char)(col*4*d);
	return;
    }
    col-=64;
    // M to B
    if(col<64)
    {
	buf[0]=(unsigned char)((63-col)*4*d);
	buf[1]=0;
	buf[2]=(unsigned char)(255*d);
	return;
    }

}

// hi lite
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); 
    }
}


// point 1,2 screen x,y and zbuffer d=brightness
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;

    //printf("%d %d\n%d %d\n\n",x1,y1,x2,y2);

    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); 
		//putcol(&buf[3*(x+y*XSIZE)], -2, 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;
}

//画面上の座標  vr:座標 vx, vy ベクトル
void screen_cord(double dx, double dy, double *sx, double *sy)
{
    *sx=XSIZE*( (dx-xmin)/(xmax-xmin));
    *sy=YSIZE*( (dy-ymin)/(ymax-ymin));
}

// inverse conversion
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;
    //gray BG
    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();

}

//flush uffer
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);
}


// lig: surface brightness
// tra: transparency
// damp_l: damping length
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);

  // for Matrix operation
  double a,b,c,d,det;

  int xmax,xmin,ymax,ymin;
  xmax=ymax=0;
  xmin=XSIZE-1;
  ymin=YSIZE-1;

  dd=lig;

  // calc screen draw range
  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;
  }

  // Matrix for color interpolation
  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;

  //calc 3 equations for triangle
  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
	  // calc Z on the plane
	  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];
// tra: transparency
// damp_l: damping length

	  if(zz<znow) 
	  {
	      double dz = znow - zz;
	      if (dz>20) dz=20;
              // calc color
              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; //(dx*dx+dy*dy)/r/r;
	    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);

}

//konpei.cc

#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; //, nt;

//class Tri *tri[MAXT];
class Vtx *vtx[MAXV];

#if 0
void AddT(class Tri *t)
{
    tri[nt++]=t;
    if(nt==MAXT) nt=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);
	    }

	}; // constructor

	~Tri()
	{
#if 0
	    int k;
	    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)
		{
		    bte->Delete(vv);
		    delete e;
		}
	    }
#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;
	    }//k
	    if (ndiv==0) return false;

	    // triple div
	    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;

		//printf("#Div %x %x %x  %x %x %x\n", v1,v2,v3,v12,v23,v31);

		AddT(v1,v12, v31);
		AddT(v2,v12, v23);
		AddT(v3,v23, v31);
		AddT(v12,v23,v31);
		return true;
	    }//ndiv=3
	    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;

		//printf("v1 v2 v3 vn %x %x %x %x\n",v1,v2,v3,vn);

		AddT(v3, v1, vn);
		AddT(v3, v2, vn);
		return true;
	    }//ndiv!=3
	    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;
	}//div
};
//Tri

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
	double dz = (zmax - (v->nv0));
	if (dz < dw*2)
	{
	    for(k=0;k<3;k++)
		v->Pos[k] += dw*nvec[k];
	}
#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];
    // write img
    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; //(int)(nv / cdep/20*mul);

#if 1
    for(i=0;i<maxst ;i++)
    {
        w_img(i);
	Grow(1e-2, lw);

	//if (i%20==0)Cut(RMAX);
	if (i%10==0)RefineMesh( DR, false);
    }
#endif



	printf("#Step %d\n",maxst);

    check_max();
}