#include <stdio.h>
#include <stdlib.h>
#include <strings.h>

#define NDIMS 3

typedef struct xyi XYI;
typedef struct xyzi XYZI;
typedef struct facept FACEPT;
typedef struct face FACE;
typedef struct tri TRI;
typedef struct loc LOC;
typedef struct cell CELL;

struct xyi {
  int x;
  int y;
  } ;

struct xyzi {
  int x;
  int y;
  int z;
  } ;

struct loc {
  LOC *link;
  XYZI p;
  int n;
  } ;

struct tri {
  TRI *link;
  LOC *a;
  LOC *b;
  LOC *c;
  int n;
  } ;

struct cell {
  char filled;
  char flag1;
  char flag2;
  } ;

struct facept {
  FACEPT *flink;
  FACEPT *blink;
  LOC *l;
  } ;

struct face {
  int npts;
  FACEPT *pts;
  FACEPT **pttail;
  XYZI normal;
  } ;

#define XHP ((XYZI){1,0,0})
#define YHP ((XYZI){0,1,0})
#define ZHP ((XYZI){0,0,1})
#define XHN ((XYZI){-1,0,0})
#define YHN ((XYZI){0,-1,0})
#define ZHN ((XYZI){0,0,-1})

static int pcsz[NDIMS];
static int pcm[NDIMS];
static int pcnv;
static CELL *pcv;
#define pc(x,y,z) pcv[((x)*pcm[2])+((y)*pcm[1])+((z)*pcm[0])]
static LOC *locs;
static int nlocs;
static TRI *tris;
static int ntris;
static FACE face;

static XYI xyi_add(XYI, XYI) __attribute__((__const__));
static XYI xyi_add(XYI a, XYI b)
{
 return((XYI){a.x+b.x,a.y+b.y});
}

static XYI xyi_sub(XYI, XYI) __attribute__((__const__));
static XYI xyi_sub(XYI a, XYI b)
{
 return((XYI){a.x-b.x,a.y-b.y});
}

static XYI xyi_rotcw(XYI) __attribute__((__const__));
static XYI xyi_rotcw(XYI a)
{
 return((XYI){a.y,-a.x});
}

static XYI xyi_rotccw(XYI) __attribute__((__const__));
static XYI xyi_rotccw(XYI a)
{
 return((XYI){-a.y,a.x});
}

static int xyi_zerop(XYI) __attribute__((__const__));
static int xyi_zerop(XYI a)
{
 return((a.x == 0) && (a.y == 0));
}

static int xyzi_zerop(XYZI) __attribute__((__const__));
static int xyzi_zerop(XYZI v)
{
 return((v.x == 0) && (v.y == 0) && (v.z == 0));
}

static int xyzi_dot(XYZI a, XYZI b)
{
 return((a.x * b.x) + (a.y * b.y) + (a.z * b.z));
}

static XYZI xyzi_sub(XYZI, XYZI) __attribute__((__const__));
static XYZI xyzi_sub(XYZI a, XYZI b)
{
 return( (XYZI) { .x = a.x - b.x,
		  .y = a.y - b.y,
		  .z = a.z - b.z } );
}

static XYZI xyzi_cross(XYZI, XYZI) __attribute__((__const__));
static XYZI xyzi_cross(XYZI a, XYZI b)
{
 return( (XYZI) { .x = (a.y * b.z) - (a.z * b.y),
		  .y = (a.z * b.x) - (a.x * b.z),
		  .z = (a.x * b.y) - (a.y * b.x) } );
}

static void baddescr(const char *how)
{
 fprintf(stderr,"Bad piece (%s)\n",how);
 exit(1);
}

static void read_pc(void)
{
 int ch;
 int dims[NDIMS];
 int lens[NDIMS];
 int muls[NDIMS];
 int x;
 int y;
 int z;
 char *bits;
 int bithave;
 int bitlen;
 int i;
 int depth;
 int m;
 CELL c;

 bits = 0;
 bitlen = 0;
 bithave = 0;
 for (i=NDIMS-1;i>=-0;i--) dims[i] = -1;
 depth = 0;
 lens[0] = 0;
 while <"readloop"> (1)
  { ch = getchar();
    switch (ch)
     { default:
	  baddescr("bad character in bitmap");
	  break;
       case EOF:
	  break <"readloop">;
	  break;
       case ' ': case '\t': case '\n':
	  break;
       case '(':
	  depth ++;
	  if (depth >= NDIMS) baddescr("too many (s in bitmap");
	  lens[depth] = 0;
	  break;
       case ')':
	  if (dims[depth] < 0)
	   { dims[depth] = lens[depth];
	   }
	  else if (lens[depth] != dims[depth])
	   { baddescr("length mismatch in bitmap");
	   }
	  depth --;
	  if (depth < 0) baddescr("unmatched ) in bitmap");
	  lens[depth] ++;
	  break;
       case '0': case '1':
	  if (depth != NDIMS-1) baddescr("bit not deep enough");
	  lens[depth] ++;
	  if (bitlen >= bithave) bits = realloc(bits,bithave=bitlen+64);
	  bits[bitlen++] = ('1'==ch);
	  break;
     }
  }
 if (depth != 0) baddescr("misbalanced ()s in bitmap");
 dims[0] = lens[0];
 m = 1;
 for (i=NDIMS-1;i>=0;i--)
  { muls[i] = m;
    m *= dims[i];
  }
 m = 1;
 for (i=NDIMS-1;i>=0;i--)
  { pcsz[i] = dims[i] + 2;
    pcm[i] = m;
    m *= pcsz[i];
  }
 pcnv = m;
 pcv = malloc(m*sizeof(CELL));
 for (x=pcsz[2]-1;x>=0;x--)
  { for (y=pcsz[1]-1;y>=0;y--)
     { for (z=pcsz[0]-1;z>=0;z--)
	{ if ( (x == 0) || (x == pcsz[2]-1) ||
	       (y == 0) || (y == pcsz[1]-1) ||
	       (z == 0) || (z == pcsz[0]-1) ||
	       !bits[((x-1)*muls[2])+((y-1)*muls[1])+((z-1)*muls[0])] )
	   { c.filled = 0;
	   }
	  else
	   { c.filled = 1;
	   }
	  pc(x,y,z) = c;
	}
     }
  }
 free(bits);
}

static LOC *loc_cache(XYZI p)
{
 LOC *l;

 for (l=locs;l;l=l->link)
  { if ((l->p.x == p.x) && (l->p.y == p.y) && (l->p.z == p.z)) return(l);
  }
 l = malloc(sizeof(LOC));
 l->p = p;
 l->link = locs;
 locs = l;
 nlocs ++;
 return(l);
}

static void save_tri(LOC *a, LOC *b, LOC *c)
{
 TRI *t;

 printf("@ triangle (%d,%d,%d) - (%d,%d,%d) - (%d,%d,%d)\n",
	a->p.x, a->p.y, a->p.z,
	b->p.x, b->p.y, b->p.z,
	c->p.x, c->p.y, c->p.z );
 t = malloc(sizeof(TRI));
 t->a = a;
 t->b = b;
 t->c = c;
 t->link = tris;
 tris = t;
 ntris ++;
}

static int point_in_triangle(XYZI p, XYZI a, XYZI b, XYZI c)
{
 XYZI ab;
 XYZI bc;
 XYZI ca;

 a = xyzi_sub(a,p);
 b = xyzi_sub(b,p);
 c = xyzi_sub(c,p);
 ab = xyzi_cross(a,b);
 bc = xyzi_cross(b,c);
 ca = xyzi_cross(c,a);
 if (xyzi_zerop(ab) && (xyzi_dot(a,b) <= 0)) return(1);
 if (xyzi_zerop(bc) && (xyzi_dot(b,c) <= 0)) return(1);
 if (xyzi_zerop(ca) && (xyzi_dot(c,a) <= 0)) return(1);
 return( (xyzi_dot(ab,bc) > 0) &&
	 (xyzi_dot(bc,ca) > 0) );
}

static void face_start(XYZI normal)
{
 face.normal = normal;
 face.npts = 0;
 face.pttail = &face.pts;
}

static void face_pt(XYZI p)
{
 FACEPT *fp;

 printf("@ face_pt(%d,%d,%d)\n",p.x,p.y,p.z);
 fp = malloc(sizeof(FACEPT));
 fp->l = loc_cache(p);
 *face.pttail = fp;
 face.pttail = &fp->flink;
 face.npts ++;
}

static void face_finish(void)
{
 FACEPT *p;
 FACEPT *q;
 FACEPT *r;
 FACEPT *s;
 int i;
 XYZI p_q;
 XYZI r_q;
 XYZI c;

 *face.pttail = 0;
 q = 0;
 for (p=face.pts;p;p=p->flink)
  { p->blink = q;
    q = p;
  }
 q->flink = face.pts;
 face.pts->blink = q;
 printf("@ face %d:\n",face.npts);
 for (p=face.pts,i=face.npts;i>0;p=p->flink,i--)
  { printf("@ %d %d %d\n",p->l->p.x,p->l->p.y,p->l->p.z);
  }
 if (face.npts < 3)
  { fprintf(stderr,"face with point count %d\n",face.npts);
    abort();
  }
 while <"lopping"> (1)
  { p = face.pts;
    q = p->flink;
    r = q->flink;
    if (face.npts < 4) break;
    printf("@ p=(%d,%d,%d) q=(%d,%d,%d) r=(%d,%d,%d)\n",
	p->l->p.x, p->l->p.y, p->l->p.z,
	q->l->p.x, q->l->p.y, q->l->p.z,
	r->l->p.x, r->l->p.y, r->l->p.z );
    p_q = xyzi_sub(p->l->p,q->l->p);
    r_q = xyzi_sub(r->l->p,q->l->p);
    c = xyzi_cross(p_q,r_q);
    printf("@ p_q=(%d,%d,%d) r_q=(%d,%d,%d) c=(%d,%d,%d)\n",
	p_q.x, p_q.y, p_q.z,
	r_q.x, r_q.y, r_q.z,
	c.x, c.y, c.z );
    if (xyzi_dot(c,face.normal) < 0)
     { printf("@ concavity\n");
       face.pts = q;
       continue;
     }
    if (xyzi_zerop(c))
     { printf("@ snipping\n");
     }
    else
     { for (s=r->flink,i=face.npts-3;i>0;s=s->flink,i--)
	{ if (point_in_triangle(s->l->p,p->l->p,q->l->p,r->l->p))
	   { printf("@ contains (%d,%d,%d)\n",s->l->p.x,s->l->p.y,s->l->p.z);
	     face.pts = q;
	     continue <"lopping">;
	   }
	}
       printf("@ lopping\n");
       save_tri(p->l,q->l,r->l);
     }
    p->flink = r;
    r->blink = p;
    free(q);
    face.npts --;
  }
 save_tri(p->l,q->l,r->l);
 free(p);
 free(q);
 free(r);
}

static void setup_mesh(void)
{
 int x;
 int y;
 int z;

 void find_face(XYZI nrm, int (*filled)(XYI), void (*flag)(XYI), XYZI (*ptmap)(XYI))
  { XYI p;
    XYI c;
    XYI f;
    XYI r;
    XYI cf;
    XYI cl;
    p = (XYI){0,0};
    c = (XYI){0,0};
    f = (XYI){1,0};
    r = (XYI){0,-1};
    printf("@ find_face, cell at (%d,%d,%d), normal (%d,%d,%d)\n",
	x,y,z,nrm.x,nrm.y,nrm.z);
    face_start(nrm);
    do
     { printf("@ c=(%d,%d) p=(%d,%d) f=(%d,%d) r=(%d,%d)\n",
		c.x,c.y, p.x,p.y, f.x,f.y, r.x,r.y );
       (*flag)(c);
       cf = xyi_add(c,f);
       cl = xyi_sub(cf,r);
       p = xyi_add(p,f);
       if (! (*filled)(cf))
	{ printf("@ right, pt (%d,%d)\n",p.x,p.y);
	  face_pt((*ptmap)(p));
	  f = r;
	  r = xyi_rotcw(r);
	}
       else if ((*filled)(cl))
	{ printf("@ left, pt (%d,%d)\n",p.x,p.y);
	  face_pt((*ptmap)(p));
	  r = f;
	  f = xyi_rotccw(f);
	  (*flag)(cf);
	  c = cl;
	}
       else
	{ printf("@ straight\n");
	  c = cf;
	}
     } while (! xyi_zerop(p));
    face_finish();
  }

 printf("@size: (%d,%d,%d)\n",pcsz[2],pcsz[1],pcsz[0]);
 for (z=0;z<pcsz[0];z++)
  { printf("@\n");
    for (y=0;y<pcsz[1];y++)
     { printf("@ ");
       for (x=0;x<pcsz[2];x++)
	{ printf(" %d",pc(x,y,z).filled);
	}
       printf("\n");
     }
  }
 locs = 0;
 nlocs = 0;
 tris = 0;
 ntris = 0;
 for (x=pcnv-1;x>=0;x--)
  { pcv[x].flag1 = 0;
    pcv[x].flag2 = 0;
  }
 for (x=pcsz[2]-1;x>0;x--)
  { for (y=pcsz[1]-1;y>0;y--)
     { for (z=pcsz[0]-1;z>0;z--)
	{ if (pc(x,y,z).filled && !pc(x+1,y,z).filled && !pc(x,y,z).flag1)
	   { find_face( XHP,
			({ int filled_xp(XYI p)
			    { return( pc(x,y+p.y,z-p.x).filled &&
				      !pc(x+1,y+p.y,z-p.x).filled );
			    }
			   &filled_xp;
			 }),
			({ void flag_xp(XYI p)
			    { pc(x,y+p.y,z-p.x).flag1 = 1;
			    }
			   &flag_xp;
			 }),
			({ XYZI ptmap_xp(XYI p)
			    { return((XYZI){x+1,y+1+p.y,z+1-p.x});
			    }
			   &ptmap_xp;
			 }) );
	   }
	  if (pc(x,y,z).filled && !pc(x-1,y,z).filled && !pc(x,y,z).flag2)
	   { find_face( XHN,
			({ int filled_xn(XYI p)
			    { return( pc(x,y-p.x,z+p.y).filled &&
				      !pc(x-1,y-p.x,z+p.y).filled );
			    }
			   &filled_xn;
			 }),
			({ void flag_xn(XYI p)
			    { pc(x,y-p.x,z+p.y).flag2 = 1;
			    }
			   &flag_xn;
			 }),
			({ XYZI ptmap_xn(XYI p)
			    { return((XYZI){x,y+1-p.x,z+1+p.y});
			    }
			   &ptmap_xn;
			 }) );
	   }
	}
     }
  }
 for (x=pcnv-1;x>=0;x--)
  { pcv[x].flag1 = 0;
    pcv[x].flag2 = 0;
  }
 for (y=pcsz[1]-1;y>0;y--)
  { for (z=pcsz[0]-1;z>0;z--)
     { for (x=pcsz[2]-1;x>0;x--)
	{ if (pc(x,y,z).filled && !pc(x,y+1,z).filled && !pc(x,y,z).flag1)
	   { find_face( YHP,
			({ int filled_yp(XYI p)
			    { return( pc(x-p.x,y,z+p.y).filled &&
				      !pc(x-p.x,y+1,z+p.y).filled );
			    }
			   &filled_yp;
			 }),
			({ void flag_yp(XYI p)
			    { pc(x-p.x,y,z+p.y).flag1 = 1;
			    }
			   &flag_yp;
			 }),
			({ XYZI ptmap_yp(XYI p)
			    { return((XYZI){x+1-p.x,y+1,z+1+p.y});
			    }
			   &ptmap_yp;
			 }) );
	   }
	  if (pc(x,y,z).filled && !pc(x,y-1,z).filled && !pc(x,y,z).flag2)
	   { find_face( YHN,
			({ int filled_yn(XYI p)
			    { return( pc(x+p.y,y,z-p.x).filled &&
				      !pc(x+p.y,y-1,z-p.x).filled );
			    }
			   &filled_yn;
			 }),
			({ void flag_yn(XYI p)
			    { pc(x+p.y,y,z-p.x).flag2 = 1;
			    }
			   &flag_yn;
			 }),
			({ XYZI ptmap_yn(XYI p)
			    { return((XYZI){x+1+p.y,y,z+1-p.x});
			    }
			   &ptmap_yn;
			 }) );
	   }
	}
     }
  }
 for (x=pcnv-1;x>=0;x--)
  { pcv[x].flag1 = 0;
    pcv[x].flag2 = 0;
  }
 for (z=pcsz[0]-1;z>0;z--)
  { for (x=pcsz[2]-1;x>0;x--)
     { for (y=pcsz[1]-1;y>0;y--)
	{ if (pc(x,y,z).filled && !pc(x,y,z+1).filled && !pc(x,y,z).flag1)
	   { find_face( ZHP,
			({ int filled_yp(XYI p)
			    { return( pc(x+p.y,y-p.x,z).filled &&
				      !pc(x+p.y,y-p.x,z+1).filled );
			    }
			   &filled_yp;
			 }),
			({ void flag_yp(XYI p)
			    { pc(x+p.y,y-p.x,z).flag1 = 1;
			    }
			   &flag_yp;
			 }),
			({ XYZI ptmap_yp(XYI p)
			    { return((XYZI){x+1+p.y,y+1-p.x,z+1});
			    }
			   &ptmap_yp;
			 }) );
	   }
	  if (pc(x,y,z).filled && !pc(x,y,z-1).filled && !pc(x,y,z).flag2)
	   { find_face( ZHN,
			({ int filled_yn(XYI p)
			    { return( pc(x-p.x,y+p.y,z).filled &&
				      !pc(x-p.x,y+p.y,z-1).filled );
			    }
			   &filled_yn;
			 }),
			({ void flag_yn(XYI p)
			    { pc(x-p.x,y+p.y,z).flag2 = 1;
			    }
			   &flag_yn;
			 }),
			({ XYZI ptmap_yn(XYI p)
			    { return((XYZI){x+1-p.x,y+1+p.y,z});
			    }
			   &ptmap_yn;
			 }) );
	   }
	}
     }
  }
}

static void dump_output(void)
{
 LOC *l;
 TRI *t;
 int i;
 double cx;
 double cy;
 double cz;

 cx = 0;
 cy = 0;
 cz = 0;
 for (i=0,l=locs;l;l=l->link,i++)
  { l->n = i;
    cx += l->p.x;
    cy += l->p.y;
    cz += l->p.z;
  }
 cx /= nlocs;
 cy /= nlocs;
 cz /= nlocs;
 if (i != nlocs) abort();
 for (i=0,t=tris;t;t=t->link,i++) t->n = i;
 if (i != ntris) abort();
 printf("@\nModelV1\n");
 printf("1 %d %d %d\n",nlocs,nlocs,ntris);
 printf("White 0 0 0 .7 .7 .7 10 .5\n");
 for (l=locs;l;l=l->link) printf("%g %g %g\n",l->p.x-cx,l->p.y-cy,l->p.z-cz);
 for (l=locs;l;l=l->link) printf("%d c\n",l->n);
 for (t=tris;t;t=t->link) printf("0 %d %d %d\n",t->a->n,t->b->n,t->c->n);
}

int main(void);
int main(void)
{
 read_pc();
 setup_mesh();
 dump_output();
 return(0);
}