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

FILE *fp;
int width,height;
long pbyte,lbyte;
int *histx,*histy;
int aleft,aright,atop,abottom;
int fx,fy;
float *vv;

main(int argc,char *argv[]){
  int c,tx,ty,*hx,i,j;
  float v;
  if(argc!=2){
    fprintf(stderr,"Usage: %s pbmfile\n",argv[0]);
    exit(2);
  }
  fp=fopen(argv[1],"r");
  if(fp==NULL){
    fprintf(stderr,"Cannot open %s\n",argv[1]);
    exit(1);
  }
  c=fgetc(fp);
  if(c!='P'||fgetc(fp)!='4'){
    fprintf(stderr,"%s not raw pbmfile\n",argv[1]);
    exit(1);
  }
  for(c=fgetc(fp);c!=EOF&&c<'0'||c>'9';c=fgetc(fp))
    ;
  width=c-'0';
  for(c=fgetc(fp);c>='0'&&c<='9';c=fgetc(fp))
    width=width*10+c-'0';
  while(c!=EOF&&c<'0'||c>'9')
    c=fgetc(fp);
  height=c-'0';
  for(c=fgetc(fp);c>='0'&&c<='9';c=fgetc(fp))
    height=height*10+c-'0';
  pbyte=ftell(fp);
  lbyte=(width+7)>>3;
  histx=(int*)calloc(width,sizeof(int));
  histy=(int*)calloc(height,sizeof(int));
  makehist(0,0,width,height);
  tx=threshold(histx,0,width);
  ty=threshold(histy,0,height);
  for(aleft=0;histx[aleft]<tx;aleft++)
    ;
  for(aright=width-1;histx[aright]<tx;aright--)
    ;
  aright++;
  for(atop=0;histy[atop]<ty;atop++)
    ;
  for(abottom=height-1;histy[abottom]<ty;abottom--)
    ;
  abottom++;
  fprintf(stderr,"bbox:%d,%d,%d,%d\n",aleft,atop,aright,abottom);
  fx=basicfreq(histx,aleft,aright);
  fy=basicfreq(histy,atop,abottom);
  fprintf(stderr,"basicfreq x:%d y:%d\n",fx,fy);
  hx=(int*)calloc(width,sizeof(int));
  for(i=0;i<aleft;i++)
    hx[i]=abottom-histx[aleft];
  for(i=aleft;i<aright;i++)
    hx[i]=abottom-histx[i];
  for(i=aright;i<width;i++)
    hx[i]=abottom-histx[aright];
  vv=(float*)calloc((fx>fy)?fx:fy,sizeof(float));
  for(i=0;i<fx;i++)
    vv[i]=0;
  for(i=aleft;i<aright;i++)
    vv[i%fx]+=hx[i];
  v=0;
  for(i=0;i<fx;i++){
    if(v<vv[i]){
      v=vv[i];
      j=i;
    }
  }
  i=(aright/fx)*fx+j;
  if(i<aright)
    i+=fx>>1;
  else
    i-=fx>>1;
  while(i>aleft+(fx>>1)){
    j=findgap(hx,i-fx*2,i);
    checkline(j,atop,i,abottom);
    i=j;
  }
  exit(0);
}

makehist(l,t,r,b){
  int x,y,c;
  for(x=0;x<width;x++)
    histx[x]=0;
  for(y=0;y<height;y++)
    histy[y]=0;
  for(y=t;y<b;y++){
    fseek(fp,y*lbyte+(l>>3)+pbyte,SEEK_SET);
    if(l%8>0)
      c=fgetc(fp);
    for(x=l;x<r;x++){
      if(x%8==0)
        c=fgetc(fp);
      if(c&(128>>(x%8))){
        histx[x]++;
        histy[y]++;
      }
    }
  }
}

threshold(int d[],int k,int s){
  int t,l,i;
  float *h,*n,*m,f,v,vmax;
  l=d[k];
  for(i=k+1;i<s;i++){
    if(l<d[i])
      l=d[i];
  }
  h=(float*)calloc(l+1,sizeof(float));
  n=(float*)calloc(l+1,sizeof(float));
  m=(float*)calloc(l+1,sizeof(float));
  for(i=0;i<=l;i++)
    h[i]=0;
  for(i=k;i<s;i++)
    h[d[i]]++;
  n[0]=h[0];
  m[0]=0;
  for(i=0;i<l;i++){
    n[i+1]=n[i]+h[i+1];
    m[i+1]=m[i]+h[i+1]*(i+1);
  }
  vmax=0;
  t=0;
  for(i=0;i<=l;i++){
    f=m[i]/n[i]-(m[l]-m[i])/(n[l]-n[i]);
    v=n[i]*(n[l]-n[i])*f*f;
    if(vmax<v){
      vmax=v;
      t=i;
    }
  }
  free(m);
  free(n);
  free(h);
  return(t);
}

basicfreq(int d[],int k,int s){
  int i,j;
  float m,v,vv;
  m=0;
  for(i=k;i<s;i++)
    m+=d[i];
  m/=s-k;
  v=0;
  for(i=k;i<s-1;i++)
    v+=(d[i]-m)*(d[i+1]-m);
  v/=s-1-k;
  vv=v;
  for(j=2;vv>=v;j++){
    vv=v;
    v=0;
    for(i=k;i<s-j;i++)
      v+=(d[i]-m)*(d[i+j]-m);
    v/=s-j-k;
  }
  for(j++;vv<=v;j++){
    vv=v;
    v=0;
    for(i=k;i<s-j;i++)
      v+=(d[i]-m)*(d[i+j]-m);
    v/=s-j-k;
  }
  return(j-2);
}

findgap(int d[],int k,int s){
  int t,l,i;
  float *n,*m,f,v,vmax;
  l=s-k;
  n=(float*)calloc(l+1,sizeof(float));
  m=(float*)calloc(l+1,sizeof(float));
  n[0]=d[k];
  m[0]=0;
  for(i=0;i<l;i++){
    n[i+1]=n[i]+d[k+i+1];
    m[i+1]=m[i]+d[k+i+1]*(i+1);
  }
  vmax=0;
  t=0;
  for(i=0;i<=l;i++){
    f=m[i]/n[i]-(m[l]-m[i])/(n[l]-n[i]);
    v=n[i]*(n[l]-n[i])*f*f;
    if(vmax<v){
      vmax=v;
      t=i;
    }
  }
  free(m);
  free(n);
  return(k+t);
}

checkline(l,t,r,b){
  int *hy,i,j;
  float v;
  makehist(l,t,r,b);
  hy=(int*)calloc(height,sizeof(int));
  for(i=0;i<t;i++)
    hy[i]=r-l-histy[t];
  for(i=t;i<b;i++)
    hy[i]=r-l-histy[i];
  for(i=b;i<width;i++)
    hy[i]=r-l-histy[b-1];
  for(i=0;i<fy;i++)
    vv[i]=0;
  for(i=t;i<b;i++)
    vv[i%fy]+=hy[i];
  v=0;
  for(i=0;i<fy;i++){
    if(v<vv[i]){
      v=vv[i];
      j=i;
    }
  }
  i=(t/fy)*fy+j;
  if(i<t)
    i+=fy>>1;
  else
    i-=fy>>1;
  while(i<b-(fy>>1)){
    j=findgap(hy,i,i+fy*2);
    printf("%d,%d,%d,%d,NaN\n",l+1,i+1,r-l-2,j-i-2);
    i=j;
  }
}