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

#define LINEPIX 2048
#define PX 1023
#define PY 1023

/* Maximum recursive subdivision depth. 20 is ample for very large local stretches. */
#define MAX_DEPTH 20

void printerror( int status);

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);
}

static int deposit_quad(double x1,double y1,double x2,double y2,
                        double x3,double y3,double x4,double y4,
                        double weight,double value,double zo[LINEPIX][LINEPIX],
                        int depth,int src_i,int src_j)
{
  int i1,i2,i3,i4,j1,j2,j3,j4;
  int ii1,ii2,ii3,ii4,jj1,jj2,jj3,jj4;
  int ic,jc;
  int imin,imax,jmin,jmax;
  double ff1,ff2,ff3,ff4,ds;
  double m12x,m12y,m23x,m23y,m34x,m34y,m41x,m41y,cx,cy;

  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--;

  imin=i1;
  if(i2<imin) imin=i2;
  if(i3<imin) imin=i3;
  if(i4<imin) imin=i4;
  imax=i1;
  if(i2>imax) imax=i2;
  if(i3>imax) imax=i3;
  if(i4>imax) imax=i4;
  jmin=j1;
  if(j2<jmin) jmin=j2;
  if(j3<jmin) jmin=j3;
  if(j4<jmin) jmin=j4;
  jmax=j1;
  if(j2>jmax) jmax=j2;
  if(j3>jmax) jmax=j3;
  if(j4>jmax) jmax=j4;

  if(imax-imin>1 || jmax-jmin>1) {
    if(depth>=MAX_DEPTH) {
      fprintf(stderr,"Adaptive subdivision exceeded depth %d at source (%d,%d):\n",MAX_DEPTH,src_i,src_j);
      fprintf(stderr,"  (%g,%g) (%g,%g) (%g,%g) (%g,%g)\n",x1,y1,x2,y2,x3,y3,x4,y4);
      return 0;
    }
    m12x=(x1+x2)/2.; m12y=(y1+y2)/2.;
    m23x=(x2+x3)/2.; m23y=(y2+y3)/2.;
    m34x=(x3+x4)/2.; m34y=(y3+y4)/2.;
    m41x=(x4+x1)/2.; m41y=(y4+y1)/2.;
    cx=(x1+x2+x3+x4)/4.; cy=(y1+y2+y3+y4)/4.;

    if(!deposit_quad(x1,y1,m12x,m12y,cx,cy,m41x,m41y,weight/4.,value,zo,depth+1,src_i,src_j)) return 0;
    if(!deposit_quad(m12x,m12y,x2,y2,m23x,m23y,cx,cy,weight/4.,value,zo,depth+1,src_i,src_j)) return 0;
    if(!deposit_quad(cx,cy,m23x,m23y,x3,y3,m34x,m34y,weight/4.,value,zo,depth+1,src_i,src_j)) return 0;
    if(!deposit_quad(m41x,m41y,cx,cy,m34x,m34y,x4,y4,weight/4.,value,zo,depth+1,src_i,src_j)) return 0;
    return 1;
  }

	      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;
      if(fabs(ds)<1.e-15 || ds<0.) {
        fprintf(stderr,"Invalid mapped area at source (%d,%d): ds=%g\n",src_i,src_j,ds);
        fprintf(stderr,"  (%g,%g) (%g,%g) (%g,%g) (%g,%g)\n",x1,y1,x2,y2,x3,y3,x4,y4);
        return 0;
      }
	      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 at source (%d,%d): (%d,%d)(%d,%d)(%d,%d)(%d,%d)(%g,%g)(%g,%g)(%g,%g)(%g,%g)\n",src_j,src_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 at source (%d,%d): (%d,%d)(%d,%d)(%d,%d)(%d,%d)(%g,%g)(%g,%g)(%g,%g)(%g,%g) %g %g %g %g\n",src_j,src_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... source (%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",src_j,src_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);getchar();
	      }

  if(ii1<0) ii1=0; else if(ii1>LINEPIX-1) ii1=LINEPIX-1;
  if(ii2<0) ii2=0; else if(ii2>LINEPIX-1) ii2=LINEPIX-1;
  if(ii3<0) ii3=0; else if(ii3>LINEPIX-1) ii3=LINEPIX-1;
  if(ii4<0) ii4=0; else if(ii4>LINEPIX-1) ii4=LINEPIX-1;
  if(jj1<0) jj1=0; else if(jj1>LINEPIX-1) jj1=LINEPIX-1;
  if(jj2<0) jj2=0; else if(jj2>LINEPIX-1) jj2=LINEPIX-1;
  if(jj3<0) jj3=0; else if(jj3>LINEPIX-1) jj3=LINEPIX-1;
  if(jj4<0) jj4=0; else if(jj4>LINEPIX-1) jj4=LINEPIX-1;
  if(isnan(zo[ii1][jj1])) zo[ii1][jj1]=0.;
  if(isnan(zo[ii2][jj2])) zo[ii2][jj2]=0.;
  if(isnan(zo[ii3][jj3])) zo[ii3][jj3]=0.;
  if(isnan(zo[ii4][jj4])) zo[ii4][jj4]=0.;
  zo[ii1][jj1]+=value*weight*ff1/ds;
  zo[ii2][jj2]+=value*weight*ff2/ds;
  zo[ii3][jj3]+=value*weight*ff3/ds;
  zo[ii4][jj4]+=value*weight*ff4/ds;
  return 1;
}

int main(int argc,char *argv[])
{
  FILE *fp;
  int i,j,ni[3],nj[3],iq,jq,nzi,nzo,nan;
  float fac;
  static float zx[LINEPIX][LINEPIX],zy[LINEPIX][LINEPIX],zi[LINEPIX][LINEPIX];
  static double zo[LINEPIX][LINEPIX];
  double gx[3][3],gy[3][3],x1,x2,x3,x4,y1,y2,y3,y4;
  double szi,szo;

  if(argc!=6){
    printf("Usage: %s warp_x warp_y factor input output \n",argv[0]);
    exit(1);
  }

  fprintf(stderr,"Read %s ... ",argv[1]);
  if((fp = fopen(argv[1],"rb")) == NULL){ printf("Cannot open %s\n",argv[1]); exit(0); }
  fread(zx,4,LINEPIX*LINEPIX,fp); fclose(fp); fprintf(stderr,"OK\n");
  fprintf(stderr,"Read %s ... ",argv[2]);
  if((fp = fopen(argv[2],"rb")) == NULL){ printf("Cannot open %s\n",argv[2]); exit(0); }
  fread(zy,4,LINEPIX*LINEPIX,fp); fclose(fp); fprintf(stderr,"OK\n");
  sscanf(argv[3],"%f",&fac);
  fprintf(stderr,"factor=%g\n",fac);
  fprintf(stderr,"Read %s ... ",argv[4]);
  if((fp = fopen(argv[4],"rb")) == NULL){ printf("Cannot open %s\n",argv[4]); exit(0); }
  fread(zi,4,LINEPIX*LINEPIX,fp); fclose(fp); fprintf(stderr,"OK\n");

  for(i=0;i<LINEPIX;i++) for(j=0;j<LINEPIX;j++) zo[i][j]=0./0.;

  for(i=0;i<LINEPIX;i++){
    ni[1]=i;
    if(i==0) ni[0]=i; else ni[0]=i-1;
    if(i==LINEPIX-1) ni[2]=i; else ni[2]=i+1;
    for(j=0;j<LINEPIX;j++){
      if(isnan(zi[i][j])||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==LINEPIX-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==LINEPIX-1&&iq==2) gy[iq][jq]+=1.;
          if(j==0&&jq==0) gx[iq][jq]-=1.;
          else if(j==LINEPIX-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();
          }

          x1=gx[iq][jq+1]+0.5;       y1=gy[iq+1][jq+1]+0.5;
          x2=gx[iq][jq]+0.5;         y2=gy[iq+1][jq]+0.5;
          x3=gx[iq][jq]+0.5;         y3=gy[iq][jq]+0.5;
          x4=gx[iq][jq+1]+0.5;       y4=gy[iq][jq]+0.5;

          if(!deposit_quad(x1,y1,x2,y2,x3,y3,x4,y4,0.25,zi[i][j],zo,0,i,j)){
            fprintf(stderr,"Adaptive remapping failed at source pixel (%d,%d).\n",i,j);
            exit(1);
          }
        }
      }
    }
  }

  szi=0.;nzi=0; szo=0.;nzo=0;
  for(i=0;i<LINEPIX;i++){
    for(j=0;j<LINEPIX;j++){
      if(isnan(zi[i][j])==0&&isnan(zx[i][j])==0&&isnan(zy[i][j])==0){ szi+=zi[i][j]; nzi++; }
      if(isnan(zo[i][j])==0){
        szo+=zo[i][j]; nzo++;
        if(isnan(szo)){fprintf(stderr,"zo[%d][%d]=%g szo=%g nzo=%d\n",i,j,zo[i][j],szo,nzo);getchar();}
      }
      zi[i][j]=zo[i][j];
    }
  }
  nzi=LINEPIX*LINEPIX-nzi; nzo=LINEPIX*LINEPIX-nzo;
  fprintf(stderr,"Sum : %lf => %lf\n",szi,szo);
  fprintf(stderr,"NAN : %d => %d\n",nzi,nzo);

  fprintf(stderr,"Write %s ... ",argv[5]);
  fp=fopen(argv[5],"wb");
  fwrite(zi,4,LINEPIX*LINEPIX,fp);
  fclose(fp);
  fprintf(stderr,"OK\n");
  return 0;
}