// compile as:
// cc -std=c99 -o bmp2laser `pkg-config --cflags --libs MagickWand` bmp2laser.c


// A simple example equivalent to
// convert logo: logo.jpg

//#include <windows.h>
#include <math.h>
#include <string.h>
#include <wand/magick_wand.h>




// minimum width of black features 
//#define MINBLACK 12
int MINBLACK=12;
// minimum width of white features under which they should not be reduced
//#define MINWHITE 20
int MINWHITE=20;
// size of features to erase and ignore
#define SKIPMIN 2

// width of the image-edge areas where there is no copper-growth to be done
// must be bigger than MINBLACK and MINWHITE
#define EDGECELL_X 20
#define EDGECELL_Y 20

#define DEF_ITERATIONS 6



// maximum width and size of images the program can handle
#define MAXWIDTH 20000
#define MAXSIZE 10000*2000

#define MAXSIZE 10000*2000
//unsigned char imgarr[MAXSIZE]; // FIXED SIZE! FIXME
//unsigned char*imgarr1;
//unsigned char*imgarr2;
unsigned char*imgarrI;
unsigned char*imgarrO;
unsigned char*lineconv;
int img_h,img_w,img_wb;

void setimgwb()
{ img_wb=(img_w+7)/8;}

size_t getimgsize()
{ return img_wb*img_h;}

unsigned char getbit(int n)
{ return 1<<n;}
unsigned char getibit(int n)
{ return 0xFF-(1<<n);}

void swapbuffers()
{ unsigned char *c;
  c=imgarrI;
  imgarrI=imgarrO;
  imgarrO=c;
}

void setpixel(int x,int y,int p)
{ size_t i;
  i=img_wb*y+(x>>3);
  if(p)imgarrI[i]|=getbit (x&0x07);
  else imgarrI[i]&=getibit(x&0x07);
}

void setpixelO(int x,int y,int p)
{ size_t i;
  i=img_wb*y+(x>>3);
  if(p)imgarrO[i]|=getbit (x&0x07);
  else imgarrO[i]&=getibit(x&0x07);
}

int getpixel(int x,int y)
{ size_t i,i1;
//  if(y>=img_h)return 0;
//  if(x>=img_w)return 0;
  i=img_wb*y+(x>>3);
  if(imgarrI[i]&getbit(x&0x07))return 1;else return 0;
}

int getpixelO(int x,int y)
{ size_t i,i1;
  i=img_wb*y+(x>>3);
  if(imgarrO[i]&getbit(x&0x07))return 1;else return 0;
}


// set pixel including its neighbors
int setbigpixelO(int x,int y,int c)
{ setpixelO(x,y,c);
  // in "+" shape
  setpixelO(x-1,y,c);
  setpixelO(x+1,y,c);
  setpixelO(x,y-1,c);
  setpixelO(x,y+1,c);
  // corners
//  setpixelO(x-1,y-1,c);
//  setpixelO(x+1,y+1,c);
//  setpixelO(x+1,y-1,c);
//  setpixelO(x-1,y+1,c);
}


// debug function for making text lookalike representation of the data
void dumpimgarr()
{ FILE*f;
  int x,x1,y,t,c;
  printf("Dumping img.txt...\n");
  f=fopen("img.txt","w");if(!f)return;
  for(y=0;y<img_h;y+=8){
    t=0;
    for(x=0;x<img_w;x++){t+=getpixel(x,y);if(x%8==7){c='0'+t;if(t==0)c=' ';else if(c=='8')c='X';fprintf(f,"%c",c);t=0;}}
//    for(x=0;x<img_wb;x++)fprintf(f,"%02x",getpixeloctet(x,y));
  fprintf(f,"\n");
  }
  fclose(f);
}




// refactor this
int laa[MAXWIDTH][2];//color;length
int laacnt;

// ugly
int analyzeline(int y)
{ int t,t1,t2,minl;
  laacnt=0;
  laa[0][0]=getpixel(0,y);laa[0][1]=1;
  for(t=1;t<img_w;t++)
   {//check min length
    minl=0;
    if(getpixel(t,y)!=getpixel(t-1,y)){for(t2=1;t2<SKIPMIN;t2++)if(getpixel(t2+t,y)!=getpixel(t,y)){minl=1;break;}}
    if(minl){for(t2=0;t2<SKIPMIN;t2++){if(getpixel(t2+t,y)!=getpixel(t,y))break;setpixel(t2+t,y,getpixel(t-1,y));}}

    if(getpixel(t,y)==getpixel(t-1,y))laa[laacnt][1]++;
    else{//if(la[laacnt][1]<=SKIPMIN)la[laacnt-1][1]+=la[laacnt][1];laacnt--;
         laacnt++;laa[laacnt][0]=getpixel(t,y);laa[laacnt][1]=1;
        }
   }
  return laacnt;
}

int analyzecolumn(int x)
{ int t,t1,t2,minl;
  laacnt=0;
  laa[0][0]=getpixel(x,0);laa[0][1]=1;
  for(t=1;t<img_h;t++)
   {//check min length
    minl=0;
    if(getpixel(x,t)!=getpixel(x,t-1)){for(t2=1;t2<SKIPMIN;t2++)if(getpixel(x,t2+t)!=getpixel(x,t)){minl=1;break;}}
    if(minl){for(t2=0;t2<SKIPMIN;t2++){if(getpixel(x,t2+t)!=getpixel(x,t))break;setpixel(x,t2+t,getpixel(x,t-1));}}

    if(getpixel(x,t)==getpixel(x,t-1))laa[laacnt][1]++;
    else{//if(la[laacnt][1]<=SKIPMIN)la[laacnt-1][1]+=la[laacnt][1];laacnt--;
         laacnt++;laa[laacnt][0]=getpixel(x,t);laa[laacnt][1]=1;
        }
   }
  return laacnt;
}

void printanalyzeline()
{ int t;
  for(t=0;t<laacnt;t++)
   {printf("%i:%i ",laa[t][0],laa[t][1]);}
  printf("\n");
}

// also ugly
int transformline()
{ int t,diff,tr;
  tr=0;
  for(t=1;t<laacnt-1;t++)
  { if(laa[t][0])continue; // white, skip it
    if(laa[t][1]>=MINBLACK)continue; // wider than limit black, skip it
    diff=MINBLACK-laa[t][1];

    if(laa[t-1][1]<=MINWHITE)if(laa[t+1][1]<=MINWHITE)continue; // nowhere to go
    if(laa[t-1][1]-diff/2>MINWHITE)if(laa[t+1][1]-diff/2>MINWHITE){laa[t-1][1]-=diff/2;laa[t][1]+=diff;diff-=diff/2;laa[t+1][1]-=diff;tr++;continue;} // subtract from both ends
    if(laa[t-1][1]<=MINWHITE){laa[t+1][1]-=diff;laa[t][1]+=diff;tr++;continue;}
                        else {laa[t-1][1]-=diff;laa[t][1]+=diff;tr++;continue;}

 }
 return tr;
}

// ugly way to widen the minimum horizonal black width so the laser can "bite"
// todo: refactor using the cellular automaton
void enforcehorizontalspaces()
{ int tx,ty,t2,t3;
  printf("Enforcing minimal horizontal spacing...                 \r");fflush(stdout);
  // dirty code - refactor!
  for(ty=1;ty<img_h-1;ty++){
      analyzeline(ty);
      transformline();
      t3=0;  
      for(tx=0;tx<laacnt;tx++){
         for(t2=0;t2<laa[tx][1];t2++){if(tx>1)if(laa[tx][1]<2)laa[tx][0]=laa[tx-1][0];//filter 1-pixel dirt
                                      if(laa[tx][0])setpixel(t3,ty,1);
                                        else // setpixel(t3,ty,0);
                                         // avoid introduction of vertical orphans
                                         if((getpixel(t3,ty-1)==1)&&(getpixel(t3,ty+1)==1))setpixel(t3,ty,1);else setpixel(t3,ty,0);
                                      t3++;
                                     }
          }
    }
}

void enforceverticalspaces()
{ int tx,ty,t2,t3;
  printf("Enforcing minimal vertical spacing...                 \r");fflush(stdout);
  // dirty code - refactor!
  for(tx=1;tx<img_w-1;tx++){
      analyzecolumn(tx);
      transformline();
      t3=0;  
      for(ty=0;ty<laacnt;ty++){
         for(t2=0;t2<laa[ty][1];t2++){if(ty>1)if(laa[ty][1]<2)laa[ty][0]=laa[ty-1][0];//filter 1-pixel dirt
                                      if(laa[ty][0])setpixel(tx,t3,1);
                                        else // setpixel(t3,ty,0);
                                         // avoid introduction of vertical orphans
                                         if((getpixel(tx-1,t3)==1)&&(getpixel(tx+1,t3)==1))setpixel(tx,t3,1);else setpixel(tx,t3,0);
                                      t3++;
                                     }
          }
    }
}




// definitions for the polar-coordinates lookup array
// number of steps per circle - 24=15 deg each, 36=10 deg each
#define ANGLECOUNT 24
#define ANGLESTEP 360/ANGLECOUNT
// 
#define RADIUS MINBLACK
#define RADIUSHALF MINBLACK/2
// which of the subarrays is sine and cosine
#define a_SIN 0
#define a_COS 1
int arr_sincos[ANGLECOUNT][2];
int arr_sincoshalf[ANGLECOUNT][2];
// populate the arr_sincos array
void precalculatearr_sincos()
{ int t;
  double d,di;
  for(t=0;t<ANGLECOUNT;t++)
   {di=(double)ANGLESTEP*t*3.1415/180;
    d=sin(di);
    arr_sincos[t][a_SIN]=(int)round(d*RADIUS);
    arr_sincoshalf[t][a_SIN]=(int)round(d*RADIUSHALF);
    d=cos(di);
    arr_sincos[t][a_COS]=(int)round(d*RADIUS);
    arr_sincoshalf[t][a_COS]=(int)round(d*RADIUSHALF);
    // printf("%2i %2i %2i\n",t,arr_sincos[t][a_SIN],arr_sincos[t][a_COS]);
   }
}

// return angle array index, modulo array size
int modang(int ang)
{ if(ang<0)ang+=ANGLECOUNT;
  ang%=ANGLECOUNT;
  return ang;
}

int isfreearound(int c,int x,int y,int n)
{ int t,t1,xx=0;
  int tx,ty;
  int m=n-1;
  int cmp=1;
  int l=1;
  int range=5;
  int start,end,mid;

  int ang[ANGLECOUNT];

  // find the surroundings, calculate the vector at which to look
  for(t=0;t<ANGLECOUNT;t++)
   { ang[t]=getpixel(x+arr_sincoshalf[t][a_COS],y+arr_sincoshalf[t][a_SIN]);
//    printf("%c",'0'+angles[t]);
  }  

  start=-1;end=-1;
  //go through the circle, look for start
  for(t=0;t<ANGLECOUNT+4;t++)
   {  if(ang[modang(t)]==1)continue;
      if(ang[modang(t-1)]==0)continue;
      if(ang[modang(t-1)]==1)
        if((ang[modang(t+1)]==1)||(ang[modang(t+2)]==1)||(ang[modang(t+3)]==1)) continue; // skip shorter sequences <=60 deg
      start=modang(t);
      for(t1=0;t1<ANGLECOUNT;t1++){end=modang(t+t1);if(ang[modang(t+t1+1)]==1)break;}
    }      

//debug: show the detected edges
/*
  for(t=0;t<ANGLECOUNT;t++)
   {if(ang[t]==0){if(ang[modang(t-1)]==1)ang[t]=2;}
    else
    if(ang[t]==1){if(ang[modang(t-1)]==0)ang[modang(t-1)]=3;}
   }
*/

  if(start>end)end+=ANGLECOUNT;
  mid=modang(((end-start)/2)+start);

  if(start<0)return 0;
  if(end<0)return 0;
 
  if(start>=0)ang[start]=2;
  if(end>=0)ang[modang(end)]=3;
  ang[mid]=9;

//debug: display the pixel surrounding data
/* 
  for(t=0;t<ANGLECOUNT;t++)
    printf("%c",'0'+ang[t]);
  printf(" %2i %2i",start,modang(end));
  printf(" %2i  ",mid);
  printf(" %04i %04i",x,y);
  printf("\n");
*/
  
  for(t=start+2;t<end-2;t++)
   { if(getpixel(x+arr_sincos[modang(t)][a_COS],y+arr_sincos[modang(t)][a_SIN])==1)return 0; 
     if(getpixel(x+arr_sincoshalf[modang(t)][a_COS],y+arr_sincoshalf[modang(t)][a_SIN])==1)return 0; 
     if(getpixel(x+arr_sincos[modang(t)][a_COS]/3,y+arr_sincos[modang(t)][a_SIN]/3)==1)return 0; 
   }
  return 1;
}

// for showing how convex the area -range..+range around the pixel is - the higher ratio the less space around the pixel
// total pixels=(2*range+1)^2 - 1=>9, 2=>25, 3=>49
int sumsurrounding(int x,int y,int range)
{ 
  int n=0,tx,ty;
  for(tx=x-range;tx<=x+range;tx++)for(ty=y-range;ty<=y+range;ty++)
   n+=getpixel(tx,ty);
  return n;
}





// decide in a cellular-automaton like where the copper areas should grow
void handlecellpixel(int x,int y,int doconvex,int dogrow)
{ int c;
  int n;
  int tx,ty,range;

  c=getpixel(x,y);

  //just copy the border areas, do not grow there
  if( (x<EDGECELL_X) || (y<EDGECELL_Y) 
    || (x>(img_w-EDGECELL_X)) || (y>(img_h-EDGECELL_Y)) )
      {setpixelO(x,y,c);return;}

/*
  // example for isolating edges
  if(getpixel(x,y-1)==c)if(getpixel(x,y+1)==c)if(getpixel(x-1,y)==c)if(getpixel(x+1,y)==c)
    if(getpixel(x-1,y-1)==c)if(getpixel(x-1,y+1)==c)if(getpixel(x+1,y-1)==c)if(getpixel(x+1,y+1)==c)
      {setpixelO(x,y,0);return;}
  setpixelO(x,y,1);return;
*/

  // filter vertical orphans
  if(getpixel(x,y-1)!=c)if(getpixel(x,y+1)!=c){if(c)setpixelO(x,y,0);else setpixelO(x,y,1);}
  
  // skip white pixels, we are only growing copper and ignoring spaces
  if(c){setpixelO(x,y,c);return;}

  // skip pixels that are not on the edges
  if(getpixel(x,y-1)==c)if(getpixel(x,y+1)==c)if(getpixel(x-1,y)==c)if(getpixel(x+1,y)==c)
    if(getpixel(x-1,y-1)==c)if(getpixel(x-1,y+1)==c)if(getpixel(x+1,y-1)==c)if(getpixel(x+1,y+1)==c)
      {setpixelO(x,y,c);return;}

  // set pixels in convex places
  // total pixels=(2*range+1)^2-1 (assuming the center pixel is empty) - 1=>8, 2=>24, 3=>48
  if(doconvex)if(sumsurrounding(x,y,2)>20){setbigpixelO(x,y,1);return;}

  // grow edge by pixel if enough space forward
  if(dogrow)if(isfreearound(c,x,y,MINBLACK)){setbigpixelO(x,y,1);return;}

  // nothing to do, just copy pixel
  setpixelO(x,y,c);
}



void convertbmp(char*fni,char*fno,int iternum)
{
    MagickWand *mw = NULL;
    int t,t1,t2,t3,tr,i,tx,ty;

    MagickWandGenesis();

    /* Create a wand */
    mw = NewMagickWand();

    /* Read the input image */
  printf("Reading file...                 \n");
    MagickReadImage(mw,fni);

    // adjust resolution to 1000 dpi
    MagickSetImageUnits(mw,PixelsPerInchResolution);
    MagickSetImageResolution(mw,1000,1000);

    img_w=MagickGetImageWidth(mw);
    img_h=MagickGetImageHeight(mw);
    //img_h=h;img_w=w;
    setimgwb();

    imgarrI=(unsigned char*)malloc(img_wb*img_h*sizeof(char)+2);
    imgarrO=(unsigned char*)malloc(img_wb*img_h*sizeof(char)+2);
    //imgarrI=imgarr1;
    //imgarrO=imgarr2;
    lineconv=(unsigned char*)malloc(img_w*sizeof(char)+2);

    precalculatearr_sincos();

  printf("Allocated %li+%li bytes      \n",img_wb*img_h*sizeof(char)+2,img_w*sizeof(char)+2);

    if(img_w>=MAXWIDTH){fprintf(stderr,"ERROR: image wider than %i pixels, aborting.\n",MAXWIDTH);return;}
    if(img_wb*img_h/8>=MAXSIZE){fprintf(stderr,"ERROR: image data bigger than %i bytes, aborting.\n",MAXSIZE);return;}
    //imgarr=(unsigned char*)malloc(getimgsize());

  printf("Quantizing...                 \n");
    MagickQuantizeImage(mw,2,GRAYColorspace,0,MagickFalse,MagickFalse);

//    xlim=img_w;//ylim=540;
//    ylim=img_h;


  printf("Reading...                 \n");
    for(ty=0;ty<img_h;ty++){
      printf("%i/%i %i                      \r",ty,img_h,tr);fflush(stdout);

      MagickExportImagePixels(mw,0,ty,img_w,1,"R",CharPixel,lineconv);
      for(tx=0;tx<img_w;tx++){setpixel(tx,ty,lineconv[tx]);}
    }

  enforcehorizontalspaces();
  enforceverticalspaces();

  for(t=0;t<iternum;t++){
    printf("Iterating cell [%i/%i]...                                 \r",t+1,iternum);fflush(stdout);
    if(t)swapbuffers();
      //if(t<3)t1=1;else t1=0;
      t1=1;
      for(ty=0;ty<img_h;ty++)for(tx=0;tx<img_w;tx++)handlecellpixel(tx,ty,t1,1);
  }

  if(iternum)swapbuffers();
  enforcehorizontalspaces();


  printf("\nWriting image...                                 \n");
  for(ty=0;ty<img_h;ty++){
    printf("%i/%i %i                      \r",ty,img_h);fflush(stdout);
    for(tx=0;tx<img_w;tx++){
      lineconv[tx]=255*getpixel(tx,ty);
      }
    MagickImportImagePixels(mw,0,ty,img_w,1,"I",CharPixel,lineconv);
   }

  printf("almost done                                  \r");fflush(stdout);
  /* write it */
  MagickQuantizeImage(mw,2,GRAYColorspace,0,MagickFalse,MagickFalse);
  printf("Saving file...                                  \n");fflush(stdout);
  MagickWriteImage(mw,fno);
  printf("done                                  \n");

  /* Tidy up */
  if(mw)mw = DestroyMagickWand(mw);
  MagickWandTerminus();
}

void help()
{ printf("Parameters: [-i<iterations>] [-b<MINBLACK>] [-w<MINWHITE>] <input-filename> <output-filename>\n");
  printf("Defaults: -i%i -b%i -w%i\n",DEF_ITERATIONS,MINBLACK,MINWHITE);
}

int main(int argc,char*argv[])
{ int n,t,iter=DEF_ITERATIONS;
  if(argc<3){help();return 0;}
  // input file, output file, number of iterations
  n=1;
  for(t=1;t<argc;t++)
   { 
    if(!strncmp(argv[t],"-h",2)){help();return 0;}
    if(!strncmp(argv[t],"-?",2)){help();return 0;}
    if(!strncmp(argv[t],"--h",3)){help();return 0;}
    if(!strncmp(argv[t],"-i",2)){iter=atoi(argv[t]+2);n++;continue;}
    if(!strncmp(argv[t],"-b",2)){MINBLACK=atoi(argv[t]+2);n++;continue;}
    if(!strncmp(argv[t],"-w",2)){MINWHITE=atoi(argv[t]+2);n++;continue;}
   }
  convertbmp(argv[n],argv[n+1],iter); 
  return 0;
}