#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define LINEPIX 4096
#define NLINE 100
#define LHEAD 69
#define HW 10
#define NX 1
#define NY 1
#define WT 160
#define BK 96

  static float zx[LINEPIX][LINEPIX],zy[LINEPIX][LINEPIX];

double tr1(double a,double b,double c,double d,double e,double f,int md)
{
  double s;
  if(md==1) s=(a/(a-c)*(d-b)+a/(a-e)*(b-f))*a/2;
  else if(md==2) s=(b/(b-d)*(a-c)+b/(b-f)*(e-a))*b/2;
  return (fabs(s));
}

double tr2(double a,double b,double c,double d)
{
  a=fabs(a);b=fabs(b);c=fabs(c);d=fabs(d);
  return((a*d-b*c)*(a*d-b*c)/(a+c)/(b+d)/2);
}

double rc1(double a,double b,double c,double d,double e,double f)
{
  return(fabs(a*d-b*c)/fabs(a-c)*fabs(a)/2+fabs(b*e-a*f)/fabs(b-f)*fabs(b)/2);
}

double rc2(double a,double b,double c,double d,double e,double f,double g,double h,int md)
{
  double s;
  if(md==1) s=(a+g)*(b-h)/2+a*a*(d-b)/(a-c)/2+g*g*(h-f)/(g-e)/2;
  else if(md==2) s=(b+d)*(a-c)/2+b*b*(g-a)/(b-h)/2+d*d*(c-e)/(d-f)/2;
  else if(md==3) s=(c+e)*(f-d)/2+c*c*(b-d)/(a-c)/2+e*e*(f-h)/(g-e)/2;
  else if(md==4) s=(f+h)*(e-g)/2+h*h*(a-g)/(b-h)/2+f*f*(e-c)/(d-f)/2;
  if(s<0.){fprintf(stderr,"%g %d\n(%g,%g)\n(%g,%g)\n(%g,%g)\n(%g,%g)\n%g %g %g\n",s,md,a,b,c,d,e,f,g,h,(c+e)*(f-d)/2,c*c*(b-d)/(a-c)/2,e*e*(f-h)/(g-e)/2);getchar();}
  return (s);
}

void conv(int jmax, int imax, unsigned char image[LINEPIX][LINEPIX])
{
  int i,j,ni[3],nj[3],i1,i2,i3,i4,j1,j2,j3,j4,iq,jq,is,js,ic,jc,ii1,ii2,ii3,ii4,jj1,jj2,jj3,jj4,nan;
  float fac=1.,fx1,fx2,fy1,fy2;
  static float z[LINEPIX][LINEPIX];
  double gx[3][3],gy[3][3],x1,x2,x3,x4,y1,y2,y3,y4,ff1,ff2,ff3,ff4,ds;

  for(i=0;i<LINEPIX;i++){
    for(j=0;j<LINEPIX;j++) z[i][j]=0./0.;
  }
  for(i=0;i<imax;i++){
    ni[1]=i;
    if(i==0) ni[0]=i; else ni[0]=i-1;
    if(i==imax-1) ni[2]=i; else ni[2]=i+1;
    for(j=0;j<jmax;j++){
      image[i][j]=255-image[i][j];
      if(isnan(zx[i][j])||isnan(zy[i][j])) continue;
      nj[1]=j;
      if(j==0) nj[0]=j; else nj[0]=j-1;
      if(j==jmax-1) nj[2]=j; else nj[2]=j+1;
      nan=0;
      for(iq=0;iq<3;iq++){
	for(jq=0;jq<3;jq++){
	  if(isnan(zy[ni[iq]][nj[jq]])) {gy[iq][jq]=zy[i][j]*fac+ni[iq];nan++;} else gy[iq][jq]=zy[ni[iq]][nj[jq]]*fac+ni[iq];
	  if(isnan(zx[ni[iq]][nj[jq]])) {gx[iq][jq]=zx[i][j]*fac+nj[jq];nan++;} else gx[iq][jq]=zx[ni[iq]][nj[jq]]*fac+nj[jq];
	  if(i==0&&iq==0)  gy[iq][jq]-=1.;
	  else if(i==imax-1&&iq==2)  gy[iq][jq]+=1.;
	  if(j==0&&jq==0)  gx[iq][jq]-=1.;
	  else if(j==jmax-1&&jq==2)  gx[iq][jq]+=1.;
	}
      }
      if(nan) continue;
      for(iq=0;iq<3;iq++){
	for(jq=0;jq<3;jq++){
	  gx[iq][jq]=(gx[iq][jq]+gx[1][1])/2;
	  gy[iq][jq]=(gy[iq][jq]+gy[1][1])/2;
	}
      }
      for(iq=0;iq<2;iq++){
	for(jq=0;jq<2;jq++){
	  if((gx[iq][jq]>gx[iq+1][jq+1]&&gy[iq][jq]>gy[iq+1][jq+1])||(gx[iq+1][jq]>gx[iq][jq+1]&&gy[iq][jq+1]>gy[iq+1][jq])){
	    fprintf(stderr,"Crushed grid: %d %d %d %d %d %d\n",nj[0],nj[1],nj[2],ni[0],ni[1],ni[2]);
	    fprintf(stderr,"(%g,%g)(%g,%g)(%g,%g)(%g,%g)\n",gx[iq][jq],gy[iq][jq],gx[iq][jq+1],gy[iq][jq+1],gx[iq+1][jq],gy[iq+1][jq],gx[iq+1][jq+1],gy[iq+1][jq+1]);getchar();
	  }
	  for(is=0;is<NY;is++){
	    fy1=(float)is/NY;
	    fy2=(float)(is+1)/NY;
	    for(js=0;js<NX;js++){
	      fx1=(float)js/NX;
	      fx2=(float)(js+1)/NX;
	      y1=gy[iq][jq]*(1-fy2)*(1-fx2)+gy[iq][jq+1]*(1-fy2)*fx2+gy[iq+1][jq]*fy2*(1-fx2)+gy[iq+1][jq+1]*fy2*fx2+0.5;
	      x1=gx[iq][jq]*(1-fy2)*(1-fx2)+gx[iq][jq+1]*(1-fy2)*fx2+gx[iq+1][jq]*fy2*(1-fx2)+gx[iq+1][jq+1]*fy2*fx2+0.5;
	      y2=gy[iq][jq]*(1-fy2)*(1-fx1)+gy[iq][jq+1]*(1-fy2)*fx1+gy[iq+1][jq]*fy2*(1-fx1)+gy[iq+1][jq+1]*fy2*fx1+0.5;
	      x2=gx[iq][jq]*(1-fy2)*(1-fx1)+gx[iq][jq+1]*(1-fy2)*fx1+gx[iq+1][jq]*fy2*(1-fx1)+gx[iq+1][jq+1]*fy2*fx1+0.5;
	      y3=gy[iq][jq]*(1-fy1)*(1-fx1)+gy[iq][jq+1]*(1-fy1)*fx1+gy[iq+1][jq]*fy1*(1-fx1)+gy[iq+1][jq+1]*fy1*fx1+0.5;
	      x3=gx[iq][jq]*(1-fy1)*(1-fx1)+gx[iq][jq+1]*(1-fy1)*fx1+gx[iq+1][jq]*fy1*(1-fx1)+gx[iq+1][jq+1]*fy1*fx1+0.5;
	      y4=gy[iq][jq]*(1-fy1)*(1-fx2)+gy[iq][jq+1]*(1-fy1)*fx2+gy[iq+1][jq]*fy1*(1-fx2)+gy[iq+1][jq+1]*fy1*fx2+0.5;
	      x4=gx[iq][jq]*(1-fy1)*(1-fx2)+gx[iq][jq+1]*(1-fy1)*fx2+gx[iq+1][jq]*fy1*(1-fx2)+gx[iq+1][jq+1]*fy1*fx2+0.5;
	      i1=y1;i2=y2;i3=y3;i4=y4;
	      if(y1<0) i1--;
	      if(y2<0) i2--;
	      if(y3<0) i3--;
	      if(y4<0) i4--;
	      j1=x1;j2=x2;j3=x3;j4=x4;
	      if(x1<0) j1--;
	      if(x2<0) j2--;
	      if(x3<0) j3--;
	      if(x4<0) j4--;
	      if(i1-i3>1||i2-i4>1){
		fprintf(stderr,"NY is too small: (%d,%d)(%d,%d)(%d,%d)(%d,%d)\n",j1,i1,j2,i2,j3,i3,j4,i4);
		fprintf(stderr,"%d %d:(%g,%g)(%g,%g)(%g,%g)(%g,%g)\n",jq,iq,x1,y1,x2,y2,x3,y3,x4,y4);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq],ni[iq],zx[ni[iq]][nj[jq]],zy[ni[iq]][nj[jq]]);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq+1],ni[iq],zx[ni[iq]][nj[jq+1]],zy[ni[iq]][nj[jq+1]]);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq+1],ni[iq+1],zx[ni[iq+1]][nj[jq+1]],zy[ni[iq+1]][nj[jq+1]]);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq],ni[iq+1],zx[ni[iq+1]][nj[jq]],zy[ni[iq+1]][nj[jq]]);
		exit(0);
	      }
	      if(j1-j3>1||j4-j2>1){
		fprintf(stderr,"NX is too small: (%d,%d)(%d,%d)(%d,%d)(%d,%d)\n",j1,i1,j2,i2,j3,i3,j4,i4);
		fprintf(stderr,"%d %d:(%g,%g)(%g,%g)(%g,%g)(%g,%g)\n",jq,iq,x1,y1,x2,y2,x3,y3,x4,y4);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq],ni[iq],zx[ni[iq]][nj[jq]],zy[ni[iq]][nj[jq]]);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq+1],ni[iq],zx[ni[iq]][nj[jq+1]],zy[ni[iq]][nj[jq+1]]);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq+1],ni[iq+1],zx[ni[iq+1]][nj[jq+1]],zy[ni[iq+1]][nj[jq+1]]);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq],ni[iq+1],zx[ni[iq+1]][nj[jq]],zy[ni[iq+1]][nj[jq]]);
		exit(0);
	      }
	      if(x1+x2+x3+x4<-2.) jc=(x1+x2+x3+x4)/4-0.5;else jc=(x1+x2+x3+x4)/4+0.5;
	      if(y1+y2+y3+y4<-2.) ic=(y1+y2+y3+y4)/4-0.5;else ic=(y1+y2+y3+y4)/4+0.5;
	      x1-=jc;x2-=jc;x3-=jc;x4-=jc;y1-=ic;y2-=ic;y3-=ic;y4-=ic;
	      ds=(x1*y2-x2*y1+x2*y3-x3*y2+x3*y4-x4*y3+x4*y1-x1*y4)/2;
	      ii1=i1;ii2=i2;ii3=i3;ii4=i4;jj1=j1;jj2=j2;jj3=j3;jj4=j4;ff1=0.;ff2=0.;ff3=0.;ff4=0.;
	      //1234
	      if(i1==i2&&i2==i3&&i3==i4&&j1==j2&&j2==j3&&j3==j4){ff1=ds;}
	      //123-4
	      else if(i1==i2&&i2==i3&&i3==i4&&j1==j2&&j2==j3){ff4=tr1(x4,y4,x1,y1,x3,y3,1);ff1=ds-ff4;}
	      else if(i1==i2&&i2==i3&&j1==j2&&j2==j3&&j3==j4){ff4=tr1(x4,y4,x3,y3,x1,y1,2);ff1=ds-ff4;}
	      else if(i1==i2&&i2==i3&&j1==j2&&j2==j3){
		if(x3*y4-x4*y3<0.){jj1++;ff1=tr2(x4,y4,x1,y1)-tr2(x3,y3,x4,y4);ff4=tr1(x4,y4,x3,y3,x1,y1,2);ff2=ds-ff1-ff4;}
		else if(x4*y1-x1*y4<0.){ii3--;ff3=tr2(x3,y3,x4,y4)-tr2(x4,y4,x1,y1);ff4=tr1(x4,y4,x1,y1,x3,y3,1);ff2=ds-ff3-ff4;}
		else {jj1++;ii3--;ff3=tr2(x3,y3,x4,y4);ff1=tr2(x4,y4,x1,y1);ff4=rc1(x4,y4,x3,y3,x1,y1);ff2=ds-ff1-ff3-ff4;}
	      }
	      //124-3
	      else if(i1==i2&&i2==i3&&i3==i4&&j1==j2&&j2==j4){ff3=tr1(x3,y3,x2,y2,x4,y4,1);ff1=ds-ff3;}
	      else if(i1==i2&&i2==i4&&j1==j2&&j2==j3&&j3==j4){ff3=tr1(x3,y3,x4,y4,x2,y2,2);ff1=ds-ff3;}
	      else if(i1==i2&&i2==i4&&j1==j2&&j2==j4){
		if(x2*y3-x3*y2<0.){ii4--;ff4=tr2(x3,y3,x4,y4)-tr2(x2,y2,x3,y3);ff3=tr1(x3,y3,x2,y2,x4,y4,1);ff1=ds-ff3-ff4;}
		else if(x3*y4-x4*y3<0.){jj2--;ff2=tr2(x2,y2,x3,y3)-tr2(x3,y3,x4,y4);ff3=tr1(x3,y3,x4,y4,x2,y2,2);ff1=ds-ff2-ff3;}
		else {jj2--;ii4--;ff2=tr2(x2,y2,x3,y3);ff4=tr2(x3,y3,x4,y4);ff3=rc1(x3,y3,x4,y4,x2,y2);ff1=ds-ff2-ff3-ff4;}
	      }
	      //134-2
	      else if(i1==i2&&i2==i3&&i3==i4&&j1==j3&&j3==j4){ff2=tr1(x2,y2,x3,y3,x1,y1,1);ff1=ds-ff2;}
	      else if(i1==i3&&i3==i4&&j1==j2&&j2==j3&&j3==j4){ff2=tr1(x2,y2,x1,y1,x3,y3,2);ff1=ds-ff2;}
	      else if(i1==i3&&i3==i4&&j1==j3&&j3==j4){
		if(x1*y2-x2*y1<0.){jj3--;ff3=tr2(x2,y2,x3,y3)-tr2(x1,y1,x2,y2);ff2=tr1(x2,y2,x1,y1,x3,y3,2);ff4=ds-ff2-ff3;}
		else if(x2*y3-x3*y2<0.){ii1++;ff1=tr2(x1,y1,x2,y2)-tr2(x2,y2,x3,y3);ff2=tr1(x2,y2,x3,y3,x1,y1,1);ff4=ds-ff1-ff2;}
		else {ii1++;jj3--;ff1=tr2(x1,y1,x2,y2);ff3=tr2(x2,y2,x3,y3);ff2=rc1(x2,y2,x1,y1,x3,y3);ff4=ds-ff1-ff2-ff3;}
	      }
	      //234-1
	      else if(i1==i2&&i2==i3&&i3==i4&&j2==j3&&j3==j4){ff1=tr1(x1,y1,x4,y4,x2,y2,1);ff2=ds-ff1;}
	      else if(i2==i3&&i3==i4&&j1==j2&&j2==j3&&j3==j4){ff1=tr1(x1,y1,x2,y2,x4,y4,2);ff2=ds-ff1;}
	      else if(i2==i3&&i3==i4&&j2==j3&&j3==j4){
		if(x4*y1-x1*y4<0.){ii2++;ff2=tr2(x1,y1,x2,y2)-tr2(x4,y4,x1,y1);ff1=tr1(x1,y1,x4,y4,x2,y2,1);ff3=ds-ff1-ff2;}
		else if(x1*y2-x2*y1<0.){jj4++;ff4=tr2(x4,y4,x1,y1)-tr2(x1,y1,x2,y2);ff1=tr1(x1,y1,x2,y2,x4,y4,2);ff3=ds-ff1-ff4;}
		else {ii2++;jj4++;ff4=tr2(x4,y4,x1,y1);ff2=tr2(x1,y1,x2,y2);ff1=rc1(x1,y1,x2,y2,x4,y4);ff3=ds-ff1-ff2-ff4;}
	      }
	      //14-23
	      else if(i1==i2&&i2==i3&&i3==i4&&j1==j4&&j2==j3){ff1=rc2(x1,y1,x2,y2,x3,y3,x4,y4,1);ff3=rc2(x1,y1,x2,y2,x3,y3,x4,y4,3);}
	      else if(i1==i4&&i2==i3&&i2<i4&&j1==j2&&j2==j3&&j3==j4){ff1=rc2(x4,y4,x1,y1,x2,y2,x3,y3,2);ff3=rc2(x4,y4,x1,y1,x2,y2,x3,y3,4);}
	      else if(i1==i4&&i2==i3&&i2<i4&&j1==j4&&j2==j3){
		if(x1*y2-x2*y1<0.){ii4--;ff3=rc2(x1,y1,x2,y2,x3,y3,x4,y4,3);ff1=rc2(x4,y4,x1,y1,x2,y2,x3,y3,2);ff4=tr2(x3,y3,x4,y4)-tr2(x1,y1,x2,y2);}
		else if(x3*y4-x4*y3<0.){ii2++;ff1=rc2(x1,y1,x2,y2,x3,y3,x4,y4,1);ff3=rc2(x4,y4,x1,y1,x2,y2,x3,y3,4);ff2=tr2(x1,y1,x2,y2)-tr2(x3,y3,x4,y4);}
		else {ii2++;ii4--;ff2=tr2(x1,y1,x2,y2);ff4=tr2(x3,y3,x4,y4);ff1=rc2(x1,y1,x2,y2,x3,y3,x4,y4,1)-ff4;ff3=rc2(x1,y1,x2,y2,x3,y3,x4,y4,3)-ff2;}
	      }
	      else if(i1==i4&&i2==i3&&i1<i3&&j1==j2&&j2==j3&&j3==j4){ff2=rc2(x2,y2,x3,y3,x4,y4,x1,y1,2);ff4=rc2(x2,y2,x3,y3,x4,y4,x1,y1,4);}
	      else if(i1==i4&&i2==i3&&i1<i3&&j1==j4&&j2==j3){
		if(x1*y2-x2*y1<0.){ii3--;ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,1);ff2=rc2(x2,y2,x3,y3,x4,y4,x1,y1,2);ff3=tr2(x3,y3,x4,y4)-tr2(x1,y1,x2,y2);}
		else if(x3*y4-x4*y3<0.){ii1++;ff2=rc2(x1,y1,x2,y2,x3,y3,x4,y4,3);ff4=rc2(x2,y2,x3,y3,x4,y4,x1,y1,4);ff1=tr2(x1,y1,x2,y2)-tr2(x3,y3,x4,y4);}
		else {ii1++;ii3--;ff1=tr2(x1,y1,x2,y2);ff3=tr2(x3,y3,x4,y4);ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,1)-ff1;ff2=rc2(x1,y1,x2,y2,x3,y3,x4,y4,3)-ff3;}
	      }
	      //12-34
	      else if(i1==i2&&i3==i4&&j1==j2&&j2==j3&&j3==j4){ff2=rc2(x1,y1,x2,y2,x3,y3,x4,y4,2);ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4);}
	      else if(i1==i2&&i2==i3&&i3==i4&&j1==j2&&j3==j4&&j4<j2){ff1=rc2(x2,y2,x3,y3,x4,y4,x1,y1,1);ff3=rc2(x2,y2,x3,y3,x4,y4,x1,y1,3);}
	      else if(i1==i2&&i3==i4&&j1==j2&&j3==j4&&j4<j2){
		if(x2*y3-x3*y2<0.){jj4++;ff1=rc2(x1,y1,x2,y2,x3,y3,x4,y4,2);ff3=rc2(x2,y2,x3,y3,x4,y4,x1,y1,3);ff4=tr2(x4,y4,x1,y1)-tr2(x2,y2,x3,y3);}
		else if(x4*y1-x1*y4<0.){jj2--;ff3=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4);ff1=rc2(x2,y2,x3,y3,x4,y4,x1,y1,1);ff2=tr2(x2,y2,x3,y3)-tr2(x4,y4,x1,y1);}
		else {jj2--;jj4++;ff2=tr2(x2,y2,x3,y3);ff4=tr2(x4,y4,x1,y1);ff1=rc2(x1,y1,x2,y2,x3,y3,x4,y4,2)-ff2;ff3=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4)-ff4;}
	      }
	      else if(i1==i2&&i2==i3&&i3==i4&&j1==j2&&j3==j4&&j1<j3){ff4=rc2(x4,y4,x1,y1,x2,y2,x3,y3,1);ff2=rc2(x4,y4,x1,y1,x2,y2,x3,y3,3);}
	      else if(i1==i2&&i3==i4&&j1==j2&&j3==j4&&j1<j3){
		if(x2*y3-x3*y2<0.){jj1++;ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4);ff2=rc2(x4,y4,x1,y1,x2,y2,x3,y3,3);ff1=tr2(x4,y4,x1,y1)-tr2(x2,y2,x3,y3);}
		else if(x4*y1-x1*y4<0.){jj3--;ff2=rc2(x1,y1,x2,y2,x3,y3,x4,y4,2);ff4=rc2(x4,y4,x1,y1,x2,y2,x3,y3,1);ff3=tr2(x2,y2,x3,y3)-tr2(x4,y4,x1,y1);}
		else {jj1++;jj3--;ff3=tr2(x2,y2,x3,y3);ff1=tr2(x4,y4,x1,y1);ff2=rc2(x1,y1,x2,y2,x3,y3,x4,y4,2)-ff1;ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4)-ff3;}
	      }
	      //12-3-4
	      else if(i1==i2&&i3==i4&&j1==j2&&j2==j4){
		if(x2*y3-x3*y2>0.) {jj2--;ff2=tr2(x2,y2,x3,y3);ff3=rc1(x3,y3,x4,y4,x2,y2);ff4=rc1(x4,y4,x3,y3,x1,y1);ff1=rc2(x1,y1,x2,y2,x3,y3,x4,y4,2)-ff2;}
		else {ff1=rc2(x1,y1,x2,y2,x3,y3,x4,y4,2);ff3=tr1(x3,y3,x2,y2,x4,y4,1);ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4)-ff3;}
	      }
	      else if(i1==i2&&i2==i3&&j1==j2&&j2==j4){
		if(x3*y4-x4*y3>0.) {jj2--;ii3--;ff3=tr2(x3,y3,x4,y4);ff2=rc1(x3,y3,x2,y2,x4,y4);ff4=rc1(x4,y4,x3,y3,x1,y1);ff1=ds-ff2-ff3-ff4;}
		else {ff4=tr1(x4,y4,x3,y3,x1,y1,2);ff3=tr1(x3,y3,x2,y2,x4,y4,1);ff1=ds-ff3-ff4;}
	      }
	      else if(i1==i2&&i2==i3&&j1==j2&&j3==j4&&j4<j2){
		if(x4*y1-x1*y4>0.) {jj2--;ii3--;jj4++;ff4=tr2(x4,y4,x1,y1);ff2=rc1(x3,y3,x2,y2,x4,y4);ff3=rc1(x4,y4,x1,y1,x3,y3);ff1=rc2(x2,y2,x3,y3,x4,y4,x1,y1,1)-ff4;}
		else {ff4=tr1(x4,y4,x3,y3,x1,y1,2);ff3=rc1(x3,y3,x2,y2,x4,y4)-tr2(x4,y4,x1,y1);ff1=rc2(x2,y2,x3,y3,x4,y4,x1,y1,1);}
	      }
	      else if(i1==i2&&i3==i4&&j1==j2&&j2==j3){
		if(x4*y1-x1*y4>0.) {jj1++;ff1=tr2(x4,y4,x1,y1);ff4=rc1(x4,y4,x3,y3,x1,y1);ff3=rc1(x3,y3,x4,y4,x2,y2);ff2=rc2(x1,y1,x2,y2,x3,y3,x4,y4,2)-ff1;}
		else {ff2=rc2(x1,y1,x2,y2,x3,y3,x4,y4,2);ff4=tr1(x4,y4,x1,y1,x3,y3,1);ff3=rc1(x3,y3,x4,y4,x2,y2)-tr2(x4,y4,x1,y1);}
	      }
	      else if(i1==i2&&i2==i4&&j1==j2&&j2==j3){
		if(x3*y4-x4*y3>0.) {jj1++;ii4--;ff4=tr2(x3,y3,x4,y4);ff1=rc1(x4,y4,x1,y1,x3,y3);ff3=rc1(x3,y3,x4,y4,x2,y2);ff2=ds-ff1-ff3-ff4;}
		else {ff3=tr1(x3,y3,x4,y4,x2,y2,2);ff4=tr1(x4,y4,x1,y1,x3,y3,1);ff2=ds-ff3-ff4;}
	      }
	      else if(i1==i2&&i2==i4&&j1==j2&&j3==j4&&j3>j1){
		if(x2*y3-x3*y2>0.) {jj1++;ii4--;jj3--;ff3=tr2(x2,y2,x3,y3);ff1=rc1(x4,y4,x1,y1,x3,y3);ff4=rc1(x3,y3,x2,y2,x4,y4);ff2=rc2(x4,y4,x1,y1,x2,y2,x3,y3,3)-ff3;}
		else {ff3=tr1(x3,y3,x4,y4,x2,y2,2);ff4=rc1(x4,y4,x1,y1,x3,y3)-tr2(x2,y2,x3,y3);ff2=rc2(x4,y4,x1,y1,x2,y2,x3,y3,3);}
	      }
	      //23-1-4
	      else if(i1==i2&&i2==i3&&j1==j4&&j2==j3){
		if(x3*y4-x4*y3>0.) {ii3--;ff3=tr2(x3,y3,x4,y4);ff4=rc1(x4,y4,x3,y3,x1,y1);ff1=rc1(x1,y1,x2,y2,x4,y4);ff2=rc2(x1,y1,x2,y2,x3,y3,x4,y4,3)-ff3;}
		else {ff2=rc2(x1,y1,x2,y2,x3,y3,x4,y4,3);ff4=tr1(x4,y4,x3,y3,x1,y1,2);ff1=rc1(x1,y1,x2,y2,x4,y4)-tr2(x3,y3,x4,y4);}
	      }
	      else if(i1==i2&&i2==i3&&j2==j3&&j3==j4){
		if(x4*y1-x1*y4>0.) {jj4++;ii3--;ff4=tr2(x4,y4,x1,y1);ff3=rc1(x4,y4,x1,y1,x3,y3);ff1=rc1(x1,y1,x2,y2,x4,y4);ff2=ds-ff1-ff3-ff4;}
		else {ff1=tr1(x1,y1,x4,y4,x2,y2,1);ff4=tr1(x4,y4,x3,y3,x1,y1,2);ff2=ds-ff1-ff4;}
	      }
	      else if(i1==i4&&i2==i3&&j2==j3&&j3==j4&&i1<i3){
		if(x1*y2-x2*y1>0.) {jj4++;ii3--;ii1++;ff1=tr2(x1,y1,x2,y2);ff3=rc1(x4,y4,x1,y1,x3,y3);ff4=rc1(x1,y1,x4,y4,x2,y2);ff2=rc2(x2,y2,x3,y3,x4,y4,x1,y1,2)-ff1;}
		else {ff1=tr1(x1,y1,x4,y4,x2,y2,1);ff4=rc1(x4,y4,x1,y1,x3,y3)-tr2(x1,y1,x2,y2);ff2=rc2(x2,y2,x3,y3,x4,y4,x1,y1,2);}
	      }
	      else if(i2==i3&&i3==i4&&j1==j4&&j2==j3){
		if(x1*y2-x2*y1>0.) {ii2++;ff2=tr2(x1,y1,x2,y2);ff1=rc1(x1,y1,x2,y2,x4,y4);ff4=rc1(x4,y4,x3,y3,x1,y1);ff3=rc2(x1,y1,x2,y2,x3,y3,x4,y4,3)-ff2;}
		else {ff3=rc2(x1,y1,x2,y2,x3,y3,x4,y4,3);ff1=tr1(x1,y1,x2,y2,x4,y4,2);ff4=rc1(x4,y4,x3,y3,x1,y1)-tr2(x1,y1,x2,y2);}
	      }
	      else if(i2==i3&&i3==i4&&j1==j2&&j2==j3){
		if(x4*y1-x1*y4>0.) {jj1++;ii2++;ff1=tr2(x4,y4,x1,y1);ff2=rc1(x1,y1,x4,y4,x2,y2);ff4=rc1(x4,y4,x3,y3,x1,y1);ff3=ds-ff1-ff2-ff4;}
		else {ff4=tr1(x4,y4,x1,y1,x3,y3,1);ff1=tr1(x1,y1,x2,y2,x4,y4,2);ff3=ds-ff1-ff4;}
	      }
	      else if(i1==i4&&i2==i3&&j1==j2&&j2==j3&&i4>i2){
		if(x3*y4-x4*y3>0.) {jj1++;ii2++;ii4--;ff4=tr2(x3,y3,x4,y4);ff2=rc1(x1,y1,x4,y4,x2,y2);ff1=rc1(x4,y4,x1,y1,x3,y3);ff3=rc2(x4,y4,x1,y1,x2,y2,x3,y3,4)-ff4;}
		else {ff4=tr1(x4,y4,x1,y1,x3,y3,1);ff1=rc1(x1,y1,x4,y4,x2,y2)-tr2(x3,y3,x4,y4);ff3=rc2(x4,y4,x1,y1,x2,y2,x3,y3,4);}
	      }
	      //34-1-2
	      else if(i1==i2&&i3==i4&&j2==j3&&j3==j4){
		if(x4*y1-x1*y4>0.) {jj4++;ff4=tr2(x4,y4,x1,y1);ff1=rc1(x1,y1,x2,y2,x4,y4);ff2=rc1(x2,y2,x1,y1,x3,y3);ff3=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4)-ff4;}
		else {ff3=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4);ff1=tr1(x1,y1,x4,y4,x2,y2,1);ff2=rc1(x2,y2,x1,y1,x3,y3)-tr2(x4,y4,x1,y1);}
	      }
	      else if(i1==i3&&i3==i4&&j2==j3&&j3==j4){
		if(x1*y2-x2*y1>0.) {jj4++;ii1++;ff1=tr2(x1,y1,x2,y2);ff4=rc1(x1,y1,x4,y4,x2,y2);ff2=rc1(x2,y2,x1,y1,x3,y3);ff3=ds-ff1-ff2-ff4;}
		else {ff2=tr1(x2,y2,x1,y1,x3,y3,2);ff1=tr1(x1,y1,x4,y4,x2,y2,1);ff3=ds-ff1-ff2;}
	      }
	      else if(i1==i3&&i3==i4&&j1==j2&&j3==j4&&j2>j4){
		if(x2*y3-x3*y2>0.) {jj4++;ii1++;jj2--;ff2=tr2(x2,y2,x3,y3);ff4=rc1(x1,y1,x4,y4,x2,y2);ff1=rc1(x2,y2,x3,y3,x1,y1);ff3=rc2(x2,y2,x3,y3,x4,y4,x1,y1,3)-ff2;}
		else {ff2=tr1(x2,y2,x1,y1,x3,y3,2);ff1=rc1(x1,y1,x4,y4,x2,y2)-tr2(x2,y2,x3,y3);ff3=rc2(x2,y2,x3,y3,x4,y4,x1,y1,3);}
	      }
	      else if(i1==i2&&i3==i4&&j1==j3&&j3==j4){
		if(x2*y3-x3*y2>0.) {jj3--;ff3=tr2(x2,y2,x3,y3);ff2=rc1(x2,y2,x1,y1,x3,y3);ff1=rc1(x1,y1,x2,y2,x4,y4);ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4)-ff3;}
		else {ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,4);ff2=tr1(x2,y2,x3,y3,x1,y1,1);ff1=rc1(x1,y1,x2,y2,x4,y4)-tr2(x2,y2,x3,y3);}
	      }
	      else if(i2==i3&&i3==i4&&j1==j3&&j3==j4){
		if(x1*y2-x2*y1>0.) {jj3--;ii2++;ff2=tr2(x1,y1,x2,y2);ff3=rc1(x2,y2,x3,y3,x1,y1);ff1=rc1(x1,y1,x2,y2,x4,y4);ff4=ds-ff1-ff2-ff3;}
		else {ff1=tr1(x1,y1,x2,y2,x4,y4,2);ff2=tr1(x2,y2,x3,y3,x1,y1,1);ff4=ds-ff1-ff2;}
	      }
	      else if(i2==i3&&i3==i4&&j1==j2&&j3==j4&&j1<j3){
		if(x4*y1-x1*y4>0.) {jj3--;ii2++;jj1++;ff1=tr2(x4,y4,x1,y1);ff3=rc1(x2,y2,x3,y3,x1,y1);ff2=rc1(x1,y1,x4,y4,x2,y2);ff4=rc2(x4,y4,x1,y1,x2,y2,x3,y3,1)-ff1;}
		else {ff1=tr1(x1,y1,x2,y2,x4,y4,2);ff2=rc1(x2,y2,x3,y3,x1,y1)-tr2(x4,y4,x1,y1);ff4=rc2(x4,y4,x1,y1,x2,y2,x3,y3,1);}
	      }
	      //14-2-3
	      else if(i1==i3&&i3==i4&&j1==j4&&j2==j3){
		if(x1*y2-x2*y1>0.) {ii1++;ff1=tr2(x1,y1,x2,y2);ff2=rc1(x2,y2,x1,y1,x3,y3);ff3=rc1(x3,y3,x4,y4,x2,y2);ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,1)-ff1;}
		else {ff4=rc2(x1,y1,x2,y2,x3,y3,x4,y4,1);ff2=tr1(x2,y2,x1,y1,x3,y3,2);ff3=rc1(x3,y3,x4,y4,x2,y2)-tr2(x1,y1,x2,y2);}
	      }
	      else if(i1==i3&&i3==i4&&j1==j2&&j2==j4){
		if(x2*y3-x3*y2>0.) {jj2--;ii1++;ff2=tr2(x2,y2,x3,y3);ff1=rc1(x2,y2,x3,y3,x1,y1);ff3=rc1(x3,y3,x4,y4,x2,y2);ff4=ds-ff1-ff2-ff3;}
		else {ff3=tr1(x3,y3,x2,y2,x4,y4,1);ff2=tr1(x2,y2,x1,y1,x3,y3,2);ff4=ds-ff2-ff3;}
	      }
	      else if(i1==i4&&i2==i3&&j1==j2&&j2==j4&&i3>i1){
		if(x3*y4-x4*y3>0.) {jj2--;ii1++;i3--;ff3=tr2(x3,y3,x4,y4);ff1=rc1(x2,y2,x3,y3,x1,y1);ff2=rc1(x3,y3,x2,y2,x4,y4);ff4=rc2(x2,y2,x3,y3,x4,y4,x1,y1,4)-ff3;}
		else {ff3=tr1(x3,y3,x2,y2,x4,y4,1);ff2=rc1(x2,y2,x3,y3,x1,y1)-tr2(x3,y3,x4,y4);ff4=rc2(x2,y2,x3,y3,x4,y4,x1,y1,4);}
	      }
	      else if(i1==i2&&i2==i4&&j1==j4&&j2==j3){
		if(x3*y4-x4*y3>0.) {ii4--;ff4=tr2(x3,y3,x4,y4);ff3=rc1(x3,y3,x4,y4,x2,y2);ff2=rc1(x2,y2,x1,y1,x3,y3);ff1=rc2(x1,y1,x2,y2,x3,y3,x4,y4,1)-ff4;}
		else {ff1=rc2(x1,y1,x2,y2,x3,y3,x4,y4,1);ff3=tr1(x3,y3,x4,y4,x2,y2,2);ff2=rc1(x2,y2,x1,y1,x3,y3)-tr2(x3,y3,x4,y4);}
	      }
	      else if(i1==i2&&i2==i4&&j1==j3&&j3==j4){
		if(x2*y3-x3*y2>0.) {jj3--;ii4--;ff3=tr2(x2,y2,x3,y3);ff4=rc1(x3,y3,x2,y2,x4,y4);ff2=rc1(x2,y2,x1,y1,x3,y3);ff1=ds-ff2-ff3-ff4;}
		else {ff2=tr1(x2,y2,x3,y3,x1,y1,1);ff3=tr1(x3,y3,x4,y4,x2,y2,2);ff1=ds-ff2-ff3;}
	      }
	      else if(i1==i4&&i2==i3&&j1==j3&&j3==j4&&i2<i4){
		if(x1*y2-x2*y1>0.) {jj3--;ii4--;ii2++;ff2=tr2(x1,y1,x2,y2);ff4=rc1(x3,y3,x2,y2,x4,y4);ff3=rc1(x2,y2,x3,y3,x1,y1);ff1=rc2(x4,y4,x1,y1,x2,y2,x3,y3,2)-ff2;}
		else {ff2=tr1(x2,y2,x3,y3,x1,y1,1);ff3=rc1(x3,y3,x2,y2,x4,y4)-tr2(x1,y1,x2,y2);ff1=rc2(x4,y4,x1,y1,x2,y2,x3,y3,2);}
	      }
	      //13-2-4
	      else if(i1==i2&& i2==i3&&j1==j3&&j3==j4){ff2=tr1(x2,y2,x3,y3,x1,y1,1);ff4=tr1(x4,y4,x3,y3,x1,y1,2);ff1=ds-ff2-ff4;}
	      else if(i1==i3&& i3==i4&&j1==j2&&j2==j3){ff2=tr1(x2,y2,x1,y1,x3,y3,2);ff4=tr1(x4,y4,x1,y1,x3,y3,1);ff3=ds-ff2-ff4;}
	      //24-1-3
	      else if(i1==i2&& i2==i4&&j2==j3&&j3==j4){ff1=tr1(x1,y1,x4,y4,x2,y2,1);ff3=tr1(x3,y3,x4,y4,x2,y2,2);ff2=ds-ff1-ff3;}
	      else if(i2==i3&& i3==i4&&j1==j2&&j2==j4){ff1=tr1(x1,y1,x2,y2,x4,y4,2);ff3=tr1(x3,y3,x2,y2,x4,y4,1);ff4=ds-ff1-ff3;}
	      //1-2-3-4
	      else if(i1==i2&&i3==i4&&j1==j4&&j2==j3){
		ff1=rc1(x1,y1,x2,y2,x4,y4);ff2=rc1(x2,y2,x1,y1,x3,y3);ff3=rc1(x3,y3,x4,y4,x2,y2);ff4=rc1(x4,y4,x3,y3,x1,y1);
	      }
	      else{
		fprintf(stderr,"Undefined configuration: %d %d: (%d,%d)(%d,%d)(%d,%d)(%d,%d)(%g,%g)(%g,%g)(%g,%g)(%g,%g)\n",j,i,j1,i1,j2,i2,j3,i3,j4,i4,x1,y1,x2,y2,x3,y3,x4,y4);
		exit(0);
	      }
	      if(isinf(ff1)||isinf(ff2)||isinf(ff3)||isinf(ff4)){
		fprintf(stderr,"Inf factor: %d %d: (%d,%d)(%d,%d)(%d,%d)(%d,%d)(%g,%g)(%g,%g)(%g,%g)(%g,%g) %g %g %g %g\n",j,i,j1,i1,j2,i2,j3,i3,j4,i4,x1,y1,x2,y2,x3,y3,x4,y4,ff1,ff2,ff3,ff4);getchar();
	      }
	      if(fabs(ff1+ff2+ff3+ff4-ds)>1e-6||ff1<0.||ff2<0.||ff3<0.||ff4<0.){
		fprintf(stderr,"Something Wrong...: %d %d: (%d,%d)(%d,%d)(%d,%d)(%d,%d)(%g,%g)(%g,%g)(%g,%g)(%g,%g) %g %g %g %g %g %g\n",j,i,j1,i1,j2,i2,j3,i3,j4,i4,x1,y1,x2,y2,x3,y3,x4,y4,ff1,ff2,ff3,ff4,ff1+ff2+ff3+ff4,ds);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq],ni[iq],zx[ni[iq]][nj[jq]],zy[ni[iq]][nj[jq]]);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq+1],ni[iq],zx[ni[iq]][nj[jq+1]],zy[ni[iq]][nj[jq+1]]);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq+1],ni[iq+1],zx[ni[iq+1]][nj[jq+1]],zy[ni[iq+1]][nj[jq+1]]);
		fprintf(stderr,"%d %d:(%g,%g)\n",nj[jq],ni[iq+1],zx[ni[iq+1]][nj[jq]],zy[ni[iq+1]][nj[jq]]);
		getchar();
	      }
	      if(ii1<0) ii1=0; else if(ii1>imax-1) ii1=imax-1;
	      if(ii2<0) ii2=0; else if(ii2>imax-1) ii2=imax-1;
	      if(ii3<0) ii3=0; else if(ii3>imax-1) ii3=imax-1;
	      if(ii4<0) ii4=0; else if(ii4>imax-1) ii4=imax-1;
	      if(jj1<0) jj1=0; else if(jj1>jmax-1) jj1=jmax-1;
	      if(jj2<0) jj2=0; else if(jj2>jmax-1) jj2=jmax-1;
	      if(jj3<0) jj3=0; else if(jj3>jmax-1) jj3=jmax-1;
	      if(jj4<0) jj4=0; else if(jj4>jmax-1) jj4=jmax-1;
	      if(isnan(z[ii1][jj1])) z[ii1][jj1]=0.;
	      if(isnan(z[ii2][jj2])) z[ii2][jj2]=0.;
	      if(isnan(z[ii3][jj3])) z[ii3][jj3]=0.;
	      if(isnan(z[ii4][jj4])) z[ii4][jj4]=0.;
	      z[ii1][jj1]+=(float)image[i][j]/NX/NY/4*ff1/ds;
	      z[ii2][jj2]+=(float)image[i][j]/NX/NY/4*ff2/ds;
	      z[ii3][jj3]+=(float)image[i][j]/NX/NY/4*ff3/ds;
	      z[ii4][jj4]+=(float)image[i][j]/NX/NY/4*ff4/ds;
	    }
	  }
	}
      }
    }
  }
  for(i=0;i<imax;i++){
    for(j=0;j<jmax;j++){
      if(isnan(z[i][j])) image[i][j]=255;
      else if(z[i][j]>255.) image[i][j]=0;
      else image[i][j]=(char)(255.5-z[i][j]);
    }
  }
}

int main(int argc,char *argv[])
{
  FILE *fp;
  unsigned short header[LHEAD];
  static unsigned char image[LINEPIX][LINEPIX][4],buf[4][LINEPIX][LINEPIX];
  static float tn[LINEPIX][LINEPIX];
  float y,zyt[LINEPIX],zy0[LINEPIX],zy0t;
  double a,b,sx,sy,sx2,sxy,sy2,zys,zxs,sg,zys0;
  int i,j,im,jm,cm,k,km,ii,jj,n,nm,nn,i0,j0,zxt[LINEPIX],p,pp,pm,p1,p2,j1,j2,i1,i2,ex[2][NLINE],cy[NLINE][LINEPIX],nw,lh,nzero;
  unsigned char zero[3];

  if(argc != 3 && argc !=4){
    fprintf(stderr,"Usage : %s input output [input2]\n",argv[0]);
    exit(0);
  }
  for(i=0;i<LINEPIX;i++){
    for(j=0;j<LINEPIX;j++){
      buf[0][i][j]=255;
      buf[1][i][j]=255;
      buf[2][i][j]=255;
      zx[i][j]=0./0.;
      zy[i][j]=0./0.;
      tn[i][j]=0./0.;
    }
  }
  for(p=0;p<NLINE;p++){
    for(j=0;j<LINEPIX;j++) cy[p][j]=0;
  }
  fprintf(stderr,"Read %s ...",argv[1]);
  if((fp = fopen(argv[1],"rb")) == NULL){
    fprintf(stderr,"Cannot open %s\n",argv[1]);
    exit(0);
  }
  fread(header,2,6,fp);
  lh=header[5]/2;
  if(lh>LHEAD){
    fprintf(stderr,"Check header length (%d)\n",header[5]);
    exit(0);
  }
  fseek(fp, 0L, SEEK_SET);
  fread(header,2,lh,fp);
  jm=header[9];im=header[11];cm=header[14]/8;
  fprintf(stderr," Image size (%d x %d) ",jm,im);
  nzero=3-(jm*cm-1)%4;
  if(jm>LINEPIX || im>LINEPIX){
    fprintf(stderr," must be smaller than %d x %d.\n",LINEPIX,LINEPIX);
    exit(0);
  }
  for(i=0;i<im;i++){
    for(j=0;j<jm;j++) fread(image[i][j],1,cm,fp);
    fread(zero,1,nzero,fp);
  }
  fclose(fp);
  fprintf(stderr,"OK\n");
  for(i=0;i<im;i++){
    for(j=0;j<jm;j++) buf[3][i][j]=image[i][j][1];
  }
  for(i=HW;i<im-HW;i++){
    for(j=0;j<jm;j++){
      //if(image[i][j][3]!=255) fprintf(stderr,"%d %d\n",j,i);
      if(image[i-HW][j][1]>WT && image[i+HW][j][1]>WT){
	for(k=-HW+1;k<HW;k++){
	  if(k==0) continue;
	  if(image[i+k][j][1]<=image[i][j][1]){
	    //if(abs(im-1-i-283)<2 && abs(j-183)<2) fprintf(stderr,"%d %d: %d %d %d %d\n",j,im-1-i,image[i-1][j][1],image[i][j][1],image[i+1][j][1],k);
	    if(abs(k)==1 && image[i+k][j][1]==image[i][j][1]) continue;
	    else break;
	  }
	}
	//if(abs(im-1-i-283)<2 && abs(j-183)<2) fprintf(stderr,"%d %d: %d %d %d %d\n",j,im-1-i,image[i-1][j][1],image[i][j][1],image[i+1][j][1],k);
	if(k==HW&&image[i][j][1]<BK){
	  for(k=-HW+1;k<HW;k++) buf[0][i+k][j]=image[i+k][j][1];
	  if(buf[0][i-1][j]+buf[0][i+1][j]==2*buf[0][i][j]) zy[i][j]=0.+i;
	  else zy[i][j]=(double)(buf[0][i-1][j]-buf[0][i+1][j])/(buf[0][i-1][j]+buf[0][i+1][j]-2*buf[0][i][j])/2.+i;
	  //if(abs(im-1-i-283)<2 && abs(j-183)<2) fprintf(stderr,"%d %d: %d %d %d %g %d\n",j,im-1-i,buf[0][i-1][j],buf[0][i][j],buf[0][i+1][j],zy[i][j],k);
	}
      }
//if(j==113&&im-1-i==1483) fprintf(stderr,"%d: %d %g\n",j,im-1-i,zy[i][j]);
    }
  }
  for(km=1;km<=HW/3;km++){
    for(i=0;i<im;i++){
      for(j=km;j<jm-km;j++){
	if(buf[0][i][j-km]>WT && buf[0][i][j+km]>WT){
	  for(k=-km;k<=km;k++){buf[0][i][j+k]=255;zy[i][j+k]=0./0.;}
	}
      }
    }
  }
  j1=jm;j2=0;
  for(i=HW;i<im-HW;i++){
    for(j=10*HW;j<jm-10*HW;j++){
//if(j==1417&&im-1-i==1446) fprintf(stderr,"%d: %d %g\n",j,im-1-i,zy[i][j]);
      if(isnan(zy[i][j])==0){
	if(j1>j) j1=j;
	if(j2<j) j2=j;
	a=0.;b=i;sg=1000.;
	for(k=1;k<=4;k++){
	  sx2=0.;
	  sxy=0.;
	  sx=0.;
	  sy2=0.;
	  sy=0.;
	  n=0;
	  for(jj=-10*HW*k/4;jj<=10*HW*k/4;jj++){
	    for(ii=-HW+a*jj;ii<=HW+a*jj;ii++){
	      if(isnan(zy[i+ii][j+jj])==0){
		if(fabs(a*jj+b-zy[i+ii][j+jj])<2*sg){
	          sx2+=jj*jj;
	          sxy+=jj*(double)zy[i+ii][j+jj];
	          sx+=jj;
	          sy2+=(double)zy[i+ii][j+jj]*zy[i+ii][j+jj];
	          sy+=zy[i+ii][j+jj];
	          n++;
//if(j==1417&&im-1-i==1446) fprintf(stderr,"%d: %d %d %d %lf %lf %lf %lf %lf %lf\n",k,jj,ii,n,zy[i+ii][j+jj],sx2,sxy,sx,sy2,sy);
		}
	      }
	    }
	  }
	  sx2/=n;
	  sxy/=n;
	  sx/=n;
	  sy2/=n;
	  sy/=n;
	  a=(sxy-sx*sy)/(sx2-sx*sx);
	  b=(sx2*sy-sxy*sx)/(sx2-sx*sx);
	  sg=sqrt(a*a*sx2+b*b+sy2+2*a*b*sx-2*a*sxy-2*b*sy);
	  if(sg<0.2) sg=0.2;
//if(j==1417&&im-1-i==1446) fprintf(stderr,"%d %d %d: %d %g %g %g\n",k,j,im-1-i,n,a,b,a*a*sx2+b*b+sy2+2*a*b*sx-2*a*sxy-2*b*sy); 
	  if(a>2.) a=2.;
	  else if(a<-2.) a=-2.;
	}
	tn[i][j]=a;
	//printf("%d %d: %d %g %g\n",j,im-1-i,n,a,zy[i][j]); 
      }
    }
  }
  j1-=HW;j2+=HW;
  //fprintf(stderr,"j1=%d j2=%d",j1,j2);getchar();

  ii=HW;p=0;
  while(1){
    for(i=ii;i<im-HW;i++){
      for(j=(j1+j2)/2-HW*10;j<=(j1+j2)/2+HW*10;j++){
	if(isnan(zy[i][j])==0&&buf[1][i][j]==255) break;
      }
      if(isnan(zy[i][j])==0&&buf[1][i][j]==255) break;
    }
    //fprintf(stderr,"%d %d %g %d\n",j,im-1-i,zy[i][j],buf[1][i][j]);
    if(i>=im-HW-1) break;
    ii=i+2*HW;
    y=zy[i][j];
    a=tn[i][j];
    y-=a;
    i=y+0.5;
    j--;
    //fprintf(stderr,"%d %d %g %g\n",j,im-1-i,y,a);
    while(j>=j1&&i>0&&i<im-1){
      for(k=-HW+1;k<HW;k++){
	if(isnan(tn[i+k][j])==0){y=zy[i+k][j];a=tn[i+k][j];break;}
      }
      y-=a;
      i=y+0.5;
      j--;
    }
    i0=i;j0=j;
    //fprintf(stderr,"%d %d %g %g\n",j,im-1-i,y,a);
    y+=a;
    i=y+0.5;
    j++;
    while(j<=j2&&i>0&&i<im-1){
      for(k=-HW+1;k<HW;k++){
	if(isnan(tn[i+k][j])==0){y=zy[i+k][j];a=tn[i+k][j];break;}
      }
//if(p==15) fprintf(stderr,"%d: %d %d %g %g\n",p,j,im-1-i,y,a);
      y+=a;
      i=y+0.5;
      j++;
    }
    //fprintf(stderr,"%d %d %g %g\n",j,im-1-i,y,a);
    y-=a;
    i=y+0.5;
    j--;
    n=0;
    while(j>=j1&&i>0&&i<im-1){
      for(k=-HW+1;k<HW;k++){
	if(isnan(tn[i+k][j])==0){
	  y=zy[i+k][j];
	  a=tn[i+k][j];
	  zxt[n]=j;
	  zyt[n]=y;
	  //if(p==15) fprintf(stderr,"%d %d %g %g\n",n,zxt[n],zyt[n],a);
	  n++;
	  break;
	}
      }
      y-=a;
      i=y+0.5;
      j--;
    }
    nm=n;
    //fprintf(stderr,"%d %d %g %g %d\n",j,im-1-i,y,a,n);getchar();
    if(i==i0&&j==j0&&nm>(j2-j1)/4){
      zys=0.;
      ex[1][p]=zxt[0];cy[p][ex[1][p]]=(int)(zyt[0]+0.5);
      ex[0][p]=zxt[nm-1];cy[p][ex[0][p]]=(int)(zyt[nm-1]+0.5);
      nw=1;
      for(pp=0;pp<p;pp++){
	if(abs(ex[0][p]-ex[0][pp])<5&&abs(cy[p][ex[0][p]]-cy[pp][ex[0][pp]])<5){
	  if(abs(ex[1][p]-ex[1][pp])<5&&abs(cy[p][ex[1][p]]-cy[pp][ex[1][pp]])<5) nw=0;
	  else fprintf(stderr,"Something wrong... %d:(%d,%d)-(%d,%d) %d:(%d,%d)-(%d,%d)\n",p,ex[0][p],im-1-cy[p][ex[0][p]],ex[1][p],im-1-cy[p][ex[1][p]],pp,ex[0][pp],im-1-cy[pp][ex[0][pp]],ex[1][pp],im-1-cy[pp][ex[1][pp]]);
	}
      }
      if(nw){
	for(n=0;n<nm;n++){
	  if(n==0) zys+=zyt[n]*(j2-zxt[n]);
	  else zys+=(zyt[n]+zyt[n-1])/2.*(zxt[n-1]-zxt[n]);
	  if(n==nm-1) zys+=zyt[n]*(zxt[n]-j1);
	  //if(p==15){fprintf(stderr,"%d %d %d %d %g %g\n",n,j1,zxt[n],j2,zyt[n],zys);}
	}
	zys/=j2-j1;
	//fprintf(stderr,"%d: %d %d",p,j1,j2);getchar();
	if(p==0) zys0=zys;
	n=10*HW;
	for(jj=j1;jj<=j2;jj++){
	  if(jj<=(j1+j2)/2){
	    j=(j1+j2)/2-jj+j1;
	    while(zxt[n]>j){
	      if(n==nm-1-10*HW) break;
	      n++;
	    }
	  }
	  else{
	    j=jj;
	    while(zxt[n]<j){
	      if(n==10*HW) break;
	      n--;
	    }
	  }
	  //fprintf(stderr,"%d %d %d",j,n,zxt[n]);getchar();
	  if(p!=0 && zys<zys0+7*HW && (j<zxt[nm-1]||j>zxt[0])) b+=zy0t-zy0[j];
	  else{
	    sx2=0.;
	    sxy=0.;
	    sx=0.;
	    sy=0.;
	    for(nn=n-10*HW;nn<=n+10*HW;nn++){
	      sx2+=(double)(zxt[nn]-j)*(zxt[nn]-j);
	      sxy+=(double)(zxt[nn]-j)*zyt[nn];
	      sx+=(zxt[nn]-j);
	      sy+=zyt[nn];
	      //if(j==142){fprintf(stderr,"%d %d %g",nn,zxt[nn],zyt[nn]);getchar();}
	    }
	    sx2/=20*HW+1;
	    sxy/=20*HW+1;
	    sx/=20*HW+1;
	    sy/=20*HW+1;
	    b=(sx2*sy-sxy*sx)/(sx2-sx*sx);
          }
	  i=b+0.5;
	  zx[i][j]=zys-b;
	  zy0t=zy0[j];
          zy0[j]=zx[i][j];
	  //fprintf(stderr,"j=%d\n",j);
	  //if(j==j1){fprintf(stderr,"%d: %d %d %g %g",p,j,im-i,im-b,zx[i][j]);getchar();}
	  for(k=-HW;k<=HW;k++) buf[1][i+k][j]=buf[0][i+k][j];
	  buf[2][i][j]=0;
	  cy[p][j]=i;
	}
	zys0=zys;
	//fprintf(stderr,"p=%d (%d,%d)-(%d,%d), %g\n",p,ex[0][p],im-1-cy[p][ex[0][p]],ex[1][p],im-1-cy[p][ex[1][p]],zys);
	if(p>0){
	  if(abs(ex[0][p]-ex[0][p-1])<5&&abs(cy[p][ex[0][p]]-cy[p-1][ex[0][p-1]])<5&&abs(ex[1][p]-ex[1][p-1])<5&&abs(cy[p][ex[1][p]]-cy[p-1][ex[1][p-1]])<5) p--;
	}
	p++;
	if(p==NLINE){
	  fprintf(stderr,"Number of detected lines reaches NLINE\n");
	  exit(0);
        }
      }
    }
  }
  pm=p;
  //for(i=0;i<im;i++){
    //for(j=0;j<jm;j++) printf("%lf\n",zx[i][j]);
  //}
  for(j=j1;j<=j2;j++){
    i1=0;
    for(i2=0;i2<im;i2++){
      if(isnan(zx[i2][j])==0){
	if(i1==0){
	  for(i=i1;i<i2;i++) zy[i][j]=zx[i2][j];
	}
	else{
	  for(i=i1;i<i2;i++){
	    zy[i][j]=zx[i1][j]*(i2-i)/(i2-i1)+zx[i2][j]*(i-i1)/(i2-i1);
	    //fprintf(stderr,"%d: %d %d %d %g %g %g",j,i1,i,i2,zx[i1][j],zy[i][j],zx[i2][j]);getchar();
	  }
	}
	i1=i2;
      }
      if(i2==im-1){
	for(i=i1;i<=i2;i++) zy[i][j]=zx[i1][j];
      }
    }
  }
  for(j=0;j<j1;j++){
    for(i=0;i<im;i++) zy[i][j]=zy[i][j1];
  }
  for(j=j2+1;j<jm;j++){
    for(i=0;i<im;i++) zy[i][j]=zy[i][j2];
  }
  for(i=0;i<im;i++){
    for(j=0;j<jm;j++){
      for(k=0;k<3;k++) image[i][j][k]=buf[k][i][j];
      zx[i][j]=0./0.;
      if(isnan(zy[i][j])){
	fprintf(stderr,"zy[%d][%d]=NaN, No Lines ?\n",i,j);
	exit(0);
      }
    }
  }
  for(p=0;p<pm;p++){
    //if(p==28) fprintf(stderr,"%d %d => ",ex[0][p],im-cy[p][ex[0][p]]);
    for(j=ex[0][p];j>HW;j--){
      i=cy[p][j];
      nw=1;
      for(ii=i;ii<=i+2*HW;ii++){
	if(buf[3][ii][j]>BK) nw=0;
        //if(p==28) fprintf(stderr,"\n %d %d %d %d",j,ii,buf[3][ii][j],nw);
      }
      if(nw==0){
	nw=1;
	for(ii=i;ii>=i-2*HW;ii--){
	  if(buf[3][ii][j]>BK) nw=0;
          //if(p==28) fprintf(stderr,"\n %d %d %d %d",j,ii,buf[3][ii][j],nw);
	}
      }
      //if(p==28) fprintf(stderr,"%d-%d ",j,nw);
      if(nw){
	ex[0][p]=j;
	if(buf[3][i][j-HW]>WT&&buf[3][i-1][j-HW]>WT&&buf[3][i+1][j-HW]>WT&&
           buf[3][i][j-HW-1]>WT&&buf[3][i-1][j-HW-1]>WT&&buf[3][i+1][j-HW-1]>WT) break;
      }
    }
    //fprintf(stderr,"%d: %d %d => ",p,ex[0][p],im-cy[p][ex[0][p]]);
    j=ex[0][p];
    while(1){
      if(buf[3][i-HW][j+HW]>WT&&buf[3][i+HW][j+HW]>WT) break;
      j++;
    }
    //fprintf(stderr,"%d %d",ex[0][p],im-cy[p][ex[0][p]]);getchar();
    ex[0][p]=j;
    for(ii=-5;ii<=5;ii++){
      image[cy[p][ex[0][p]]+ii][ex[0][p]][0]=255;
      image[cy[p][ex[0][p]]+ii][ex[0][p]][1]=0;
      image[cy[p][ex[0][p]]+ii][ex[0][p]][2]=0;
    }
    for(j=ex[1][p];j<jm-HW;j++){
      i=cy[p][j];
      nw=1;
      for(ii=i;ii<=i+2*HW;ii++){
	if(buf[3][ii][j]>BK) nw=0;
      }
      if(nw==0){
	nw=1;
	for(ii=i;ii>=i-2*HW;ii--){
	  if(buf[3][ii][j]>BK) nw=0;
	}
      }
      if(nw){
	ex[1][p]=j;
	if(buf[3][i][j+HW]>WT&&buf[3][i-1][j+HW]>WT&&buf[3][i+1][j+HW]>WT&&
	   buf[3][i][j+HW+1]>WT&&buf[3][i-1][j+HW+1]>WT&&buf[3][i+1][j+HW+1]>WT) break;
      }
    }
    j=ex[1][p];
    while(1){
      if(buf[3][i-HW][j-HW]>WT&&buf[3][i+HW][j-HW]>WT) break;
      j--;
    }
    ex[1][p]=j;
    for(ii=-5;ii<=5;ii++){
      image[cy[p][ex[1][p]]+ii][ex[1][p]][0]=255;
      image[cy[p][ex[1][p]]+ii][ex[1][p]][1]=0;
      image[cy[p][ex[1][p]]+ii][ex[1][p]][2]=0;
    }
  }
  /*
  fprintf(stderr,"Write stavealign.bmp ...");
  fopen("stavealign.bmp","wb");
  fwrite(header,2,lh,fp);
  for(i=0;i<im;i++){
    for(j=0;j<jm;j++) fwrite(image[i][j],1,cm,fp);
    fwrite(zero,1,nzero,fp);
  }
  fclose(fp);
  fprintf(stderr,"OK\n");
  */
  if(pm<4){
    fprintf(stderr,"Number of the detected line is too small (%d)\n",pm);
    exit(0);
  }
  for(k=0;k<2;k++){//端補正(左：k=0、右：k=1)
    for(p=0;p<pm;p++){
      a=0.;b=0.;sg=1e10;
      for(n=3;n>=0;n--){
	sx2=0.;
	sxy=0.;
	sy2=0.;
	sx=0.;
	sy=0.;
	nm=0;
	p1=p;
	while(p1>0){
	  if(cy[p1][ex[k][p1]]-cy[p1-1][ex[k][p1-1]]>HW*8) break;
	  p1--;
	}
	p2=p;
	while(p2<pm-1){
	  if(cy[p2+1][ex[k][p2+1]]-cy[p2][ex[k][p2]]>HW*8) break;
	  p2++;
	}
	//fprintf(stderr,"n:p1/p/p2=%d:%d/%d/%d\n",n,p1,p,p2);
	for(pp=p1;pp<=p2;pp++){
	  if(fabs(a*(cy[pp][ex[k][pp]]-cy[p][ex[k][p]])+b-ex[k][pp])<1.5*sg){
	    sx2+=(cy[pp][ex[k][pp]]-cy[p][ex[k][p]])*(cy[pp][ex[k][pp]]-cy[p][ex[k][p]]);
	    sxy+=(cy[pp][ex[k][pp]]-cy[p][ex[k][p]])*ex[k][pp];
	    sy2+=ex[k][pp]*ex[k][pp];
	    sx+=(cy[pp][ex[k][pp]]-cy[p][ex[k][p]]);
	    sy+=ex[k][pp];
	    nm++;
	  }
	  if(n==0&&p==pp){
	    //fprintf(stderr,"Replaced: (%d,%d) fitted_x=%g ",ex[k][p],im-1-cy[p][ex[k][p]],a*(cy[pp][ex[k][pp]]-cy[p][ex[k][p]])+b);
	    ex[k][p]=(int)(a*(cy[pp][ex[k][pp]]-cy[p][ex[k][p]])+b+0.5);
	    //fprintf(stderr,"=> (%d,%d)\n",ex[k][p],im-1-cy[p][ex[k][p]]);
	    for(ii=-3;ii<=3;ii++){
	      image[cy[p][ex[k][p]]+ii][ex[k][p]][0]=255;
	      image[cy[p][ex[k][p]]+ii][ex[k][p]][1]=0;
	      image[cy[p][ex[k][p]]+ii][ex[k][p]][2]=255;
	    }
	  }
	}
	if(nm<=2){
	  fprintf(stderr,"a=%g b=%g sg=%g nm=%d\n",a,b,sg,nm);
	  n=-1;
	  break;
	}
	sx2/=nm;
	sxy/=nm;
	sy2/=nm;
	sx/=nm;
	sy/=nm;
	a=(sxy-sx*sy)/(sx2-sx*sx);
	b=(sx2*sy-sxy*sx)/(sx2-sx*sx);
	sg=sqrt(a*a*sx2+b*b+sy2+2*a*b*sx-2*a*sxy-2*b*sy);
	if(sg<1 || isnan(sg)) sg=1;
	//fprintf(stderr,"a=%g b=%g sg=%g nm=%d\n",a,b,sg,nm);
	//fprintf(stderr,"%d %d %d %g %g",k,cy[p][ex[k][p]],pp,sy,sg);getchar();
      }
    }
  }
  fprintf(stderr,"Write stavealign.bmp ...");
  fopen("stavealign.bmp","wb");
  fwrite(header,2,lh,fp);
  for(i=0;i<im;i++){
    for(j=0;j<jm;j++) fwrite(image[i][j],1,cm,fp);
    fwrite(zero,1,nzero,fp);
  }
  fclose(fp);
  fprintf(stderr,"OK\n");
/*
  if(cy[0][ex[0][0]]<cy[0][ex[1][0]]) i1=cy[0][ex[0][0]];
  else i1=cy[0][ex[1][0]];
  if(cy[pm-1][ex[0][pm-1]]>cy[pm-1][ex[1][pm-1]]) i2=cy[pm-1][ex[0][pm-1]];
  else i2=cy[pm-1][ex[1][pm-1]];
  for(k=0;k<2;k++){
    zxs=0.;
    for(p=0;p<pm;p++) zxs+=ex[k][p];
    zxs/=pm;
    p=HW;
    for(i=i1;i<=i2;i++){
      while(cy[p][ex[k][p]]<i){
        if(p==pm-1-HW) break;
	p++;
      }
      sx2=0.;
      sxy=0.;
      sy2=0.;
      sx=0.;
      sy=0.;
      for(pp=p-HW;pp<=p+HW;pp++){
	sx2+=(cy[pp][ex[k][pp]]-i)*(cy[pp][ex[k][pp]]-i);
	sxy+=(cy[pp][ex[k][pp]]-i)*ex[k][pp];
	sy2+=ex[k][pp]*ex[k][pp];
	sx+=(cy[pp][ex[k][pp]]-i);
	sy+=ex[k][pp];
      }
      sx2/=2*HW+1;
      sxy/=2*HW+1;
      sy2/=2*HW+1;
      sx/=2*HW+1;
      sy/=2*HW+1;
      b=(sx2*sy-sxy*sx)/(sx2-sx*sx);
      zx[i][(int)(b+0.5)]=zxs-b;
      if(i==i1){
        for(ii=i1-1;ii>=0;ii--) zx[ii][(int)(b+0.5)]=zx[i][(int)(b+0.5)];
      }
      if(i==i2){
        for(ii=i2+1;ii<im;ii++) zx[ii][(int)(b+0.5)]=zx[i][(int)(b+0.5)];
      }
    }
  }
  for(i=0;i<im;i++){
    j1=0;
    for(j2=0;j2<jm;j2++){
      if(isnan(zx[i][j2])==0){
	if(j1==0){
	  for(j=j1;j<j2;j++) zx[i][j]=zx[i][j2];
	}
	else{
	  for(j=j1+1;j<j2;j++) zx[i][j]=zx[i][j1]*(j2-j)/(j2-j1)+zx[i][j2]*(j-j1)/(j2-j1);
	}
	j1=j2;
      }
      if(j2==jm-1){
	for(j=j1+1;j<=j2;j++) zx[i][j]=zx[i][j1];
      }
    }
  }
*/
  for(p=0;p<pm;p++){
    for(k=0;k<2;k++){//端補正(左：k=0、右：k=1)
      p1=p;
      while(p1>0){
	if(cy[p1][ex[k][p1]]-cy[p1-1][ex[k][p1-1]]>HW*8) break;
	p1--;
      }
      p2=p;
      while(p2<pm-1){
	if(cy[p2+1][ex[k][p2+1]]-cy[p2][ex[k][p2]]>HW*8) break;
	p2++;
      }
      p1=0;p2=pm-1;
      zxs=0.;
      for(pp=p1;pp<=p2;pp++) zxs+=ex[k][pp];
      zxs/=p2-p1+1;
      zx[cy[p][ex[k][p]]][ex[k][p]]=zxs-ex[k][p];
    }
    //fprintf(stderr,"%d: %d %d %g - %d %d %g\n",p,ex[0][p],cy[p][ex[0][p]],zx[cy[p][ex[0][p]]][ex[0][p]],ex[1][p],cy[p][ex[1][p]],zx[cy[p][ex[1][p]]][ex[1][p]]);
    for(j=0;j<jm;j++){
      if(j<ex[0][p]){cy[p][j]=cy[p][ex[0][p]];zx[cy[p][j]][j]=zx[cy[p][ex[0][p]]][ex[0][p]];}
      else if(j>ex[1][p]){cy[p][j]=cy[p][ex[1][p]];zx[cy[p][j]][j]=zx[cy[p][ex[1][p]]][ex[1][p]];}
      else zx[cy[p][j]][j]=(zx[cy[p][ex[0][p]]][ex[0][p]]*(ex[1][p]-j)+zx[cy[p][ex[1][p]]][ex[1][p]]*(j-ex[0][p]))/(ex[1][p]-ex[0][p]);
      //if(p==0){fprintf(stderr,"%d,%d; %g",j,ex[0][p],zx[cy[p][ex[0][p]]][j]);getchar();}
    }
  }
  //fprintf(stderr,"%d,%d; %g %g\n",ex[0][0],cy[0][ex[0][0]],zx[cy[0][ex[0][0]]][0],zx[cy[0][ex[0][0]]][ex[0][0]]);
  //fprintf(stderr,"%d,%d; %g %g\n",ex[0][pm-1],cy[pm-1][ex[0][pm-1]],zx[cy[pm-1][ex[0][pm-1]]][0],zx[cy[pm-1][ex[0][pm-1]]][ex[0][pm-1]]);
  for(j=0;j<jm;j++){
    for(p=0;p<=pm;p++){
      if(p==0){i1=0;i2=cy[p][j];}
      else if(p==pm){i1=cy[p-1][j];i2=im-1;}
      else {i1=cy[p-1][j];i2=cy[p][j];}
      for(i=i1;i<=i2;i++){
	if(p==0) zx[i][j]=zx[i2][j];
	else if(p==pm) zx[i][j]=zx[i1][j];
	else zx[i][j]=(zx[i1][j]*(i2-i)+zx[i2][j]*(i-i1))/(i2-i1);
      }
    }
  }
  //fprintf(stderr,"%d,%d; %g %g\n",ex[0][0],cy[0][ex[0][0]],zx[cy[0][ex[0][0]]][0],zx[cy[0][ex[0][0]]][ex[0][0]]);
  //fprintf(stderr,"%d,%d; %g %g\n",ex[0][pm-1],cy[pm-1][ex[0][pm-1]],zx[cy[pm-1][ex[0][pm-1]]][0],zx[cy[pm-1][ex[0][pm-1]]][ex[0][pm-1]]);
  if(argc==4){
    fprintf(stderr,"Read %s ...",argv[3]);
    if((fp = fopen(argv[3],"rb")) == NULL){
      fprintf(stderr,"Cannot open %s\n",argv[3]);
      exit(0);
    }
    fread(header,2,6,fp);
    lh=header[5]/2;
    if(lh>LHEAD){
      fprintf(stderr,"Check header length (%d)\n",header[5]);
      exit(0);
    }
    fseek(fp, 0L, SEEK_SET);
    fread(header,2,lh,fp);
    jm=header[9];im=header[11];cm=header[14]/8;
    fprintf(stderr," Image size (%d x %d) ",jm,im);
    nzero=3-(jm*cm-1)%4;
    if(jm>LINEPIX || im>LINEPIX){
      fprintf(stderr," must be smaller than %d x %d.\n",LINEPIX,LINEPIX);
      exit(0);
    }
    for(i=0;i<im;i++){
      for(j=0;j<jm;j++) fread(image[i][j],1,cm,fp);
      fread(zero,1,nzero,fp);
    }
    fclose(fp);
    fprintf(stderr,"OK\n");
    for(i=0;i<im;i++){
      for(j=0;j<jm;j++) buf[3][i][j]=image[i][j][1];
    }
  }
  for(i=0;i<LINEPIX;i++){
    for(j=0;j<LINEPIX;j++){
      zx[i][j]=0.;//横変形無しの場合
    }
  }
  for(i=0;i<im;i++){
    //for(j=0;j<jm;j++) printf("%lf\n",zx[i][j]);
    //for(j=0;j<jm;j++) printf("%lf\n",zy[i][j]);
  }
  conv(jm,im,buf[3]);
  for(i=0;i<im;i++){
    for(j=0;j<jm;j++){
      for(k=0;k<3;k++) image[i][j][k]=buf[3][i][j];
    }
  }
  fprintf(stderr,"Write %s ...",argv[2]);
  fopen(argv[2],"wb");
  fwrite(header,2,lh,fp);
  for(i=0;i<im;i++){
    for(j=0;j<jm;j++) fwrite(image[i][j],1,cm,fp);
    fwrite(zero,1,nzero,fp);
  }
  fclose(fp);
  fprintf(stderr,"OK\n");
  return 0;
}