morph.cpp

Go to the documentation of this file.

#include <qstring.h>
#include <math.h>

inline double min(double a, double b)
{
  if (a<b) return a;
  return b;
}

extern double dist(double, double, double, double);

int sign(double b)
{
  if (b>0) return 1;
  if (b<0) return -1;
  return 0;
}

void pos(double x1, double y1, double x2, double y2, double pct, double& x, double &y)
{
  double dx = x2-x1;
  double dy = y2-y1;
  x = x1 + dx*pct;
  y = y1 + dy*pct;
}

class NetPoint{
  public:
  double x;
  double y;
  double shouldbe_x;
  double shouldbe_y;
  double a;
  double b;
  double nx;
  double ny;
  bool free;
  NetPoint * up;
  NetPoint * down;
  NetPoint * left;
  NetPoint * right;
  NetPoint()
  {
    x = y = shouldbe_x = shouldbe_y = 0;
    free = true;
    up = down = left = right = NULL;
  }
  void setNotFree()
  {
    shouldbe_x=x;
    shouldbe_y=y;
    free = false;
  }
  void mark_as_tx()
  {
    a = x;
    b = y;
    x = shouldbe_x;
    y = shouldbe_y;
  }
  void move(bool goto_shouldbe)
  {
    double bx = x;
    double by = y;
    if (free)
      {
        int i = 0;
        x = 0;
        y = 0;
        if (left)
          {
            x+=left->x;
            y+=left->y;
            i++;
          }
        if (right)
          {
            x+=right->x;
            y+=right->y;
            i++;
          }
        if (up)
          {
            x+=up->x;
            y+=up->y;
            i++;
          }
        if (down)
          {
            x+=down->x;
            y+=down->y;
            i++;
          }
        x/=i;
        y/=i;
      }
    else if (goto_shouldbe && (shouldbe_x != x || shouldbe_y != y))
      {
        double dx=0;
        double dy=0;
        dx = sign(shouldbe_x-x);
        dy = sign(shouldbe_y-y);
        x+=dx;
        y+=dy;
      }
    nx = x;
    ny =y;
    x = bx;
    y = by;
  }
  
  void fixx()
  {
    if (!free) return;
    // where are the current endpoints
    NetPoint * L;
    L = left;
    while(L)
      {
        if (!L->free) break;
        L=L->left;
      }

    NetPoint * R;
    R = right;
    while(R)
      {
        if (!R->free) break;
        R=R->right;
      }
    assert(L&&R);
    double pct = (a-L->a)/(R->a-L->a);
    printf("%g\n",pct);
    x = pct*(R->x-L->x)+L->x;
  }

  void fixy()
  {
    if (!free) return;
    // where are the current endpoints
    NetPoint * U;
    U = up;
    while(U)
      {
        if (!U->free) break;
        U=U->up;
      }
    
    NetPoint * D;
    D = down;
    while(D)
      {
        if (!D->free) break;
        D=D->down;
      }
    assert(U&&D);
    double pct = (b-U->b)/(D->b-U->b);
    y = pct*(D->y-U->y)+U->y;
  }
  
  void realize()
  {
    x = nx;
    y = ny;
  }
  
  void set_bound(int mx, int my)
  {
    if (x>mx-1) x = mx-1;
    if (y>my-1) y = my-1;
    if (a>mx-1) a = mx-1;
    if (b>my-1) b = my-1;
    if (x<0) x = 0;
    if (y<0) y = 0;
    if (a<0) a = 0;
    if (b<0) b = 0;
  }
};

extern QString type_name(NetPoint);

#include "value.h"
#include "runtime.h"

/*
Value * Runtime::morph(Value * image, Value * source, Value * target)
{
  // input checks
  if (image==NULL)  return NULL;
  if (source==NULL) return NULL;
  if (target==NULL) return NULL;
  if (typeid(*image)!=typeid(ImageValue)) return NULL;
  if (typeid(*source)!=typeid(ValueArray)) return NULL;
  if (typeid(*target)!=typeid(ValueArray)) return NULL;
  ImageValue * I = (ImageValue*)image;
  ValueArray * S = (ValueArray*)source;
  ValueArray * T = (ValueArray*)target;
  if (S->sx != T->sx) return NULL;
  if (S->sy != T->sy) return NULL;
  if (S->sy != 1 && S->sx != 1) return NULL;
  int marks = S->sx;
  if (marks==1) marks=S->sy;
  int xs = I->width();
  int ys = I->height();
  int iteraties=dist(xs,ys,0,0);
  int ch  = I->channelcount();
  double target_x[marks];
  double target_y[marks];
  double source_x[marks];
  double source_y[marks];
  if (S->sx==1)
    for(int i = 0 ; i < marks ; i++)
      {
        Value * sv = S->get(i,0);
        Value * tv = T->get(i,0);
        if (sv==NULL) return NULL;
        if (tv==NULL) return NULL;
        if (typeid(*sv)!=typeid(PointValue)) return NULL;
        if (typeid(*tv)!=typeid(PointValue)) return NULL;
        PointValue * s = (PointValue*)sv;
        PointValue * t = (PointValue*)tv;
        target_x[i]=(t->x + 1)*xs/2.0;
        target_y[i]=(t->y + 1)*ys/2.0;
        source_x[i]=(s->x + 1)*xs/2.0;
        source_y[i]=(s->y + 1)*ys/2.0;
      }
  else
    for(int i = 0 ; i < marks ; i++)
      {
        Value * sv = S->get(0,i);
        Value * tv = T->get(0,i);
        if (sv==NULL) return NULL;
        if (tv==NULL) return NULL;
        if (typeid(*sv)!=typeid(PointValue)) return NULL;
        if (typeid(*tv)!=typeid(PointValue)) return NULL;
        PointValue * s = (PointValue*)sv;
        PointValue * t = (PointValue*)tv;
        target_x[i]=(t->x + 1)*xs/2.0;
        target_y[i]=(t->y + 1)*ys/2.0;
        source_x[i]=(s->x + 1)*xs/2.0;
        source_y[i]=(s->y + 1)*ys/2.0;
      }
  // we will work with a superimposed net 
  // each going through the given points.
  Array<2,NetPoint> net(marks+2,marks+2);
  Array<1,NetPoint> toprow=net(Select<1>(0));
  int nsx = 1;
  double equality_distance = 9;
  double px = toprow[0].x = 0;
  while(px!=xs-1)
    {
      double md = xs-1-px;
      for(int n = 0 ; n < marks ; n++)
        {
          double d = target_x[n]-px;
          if (d>equality_distance && d < md) md = d;
        }
      px+=md;
      toprow[nsx++].x=px;
    }
  // do the same for the y coordinates
  Array<1,NetPoint> leftcol=net(Select<1>(1));
  int nsy = 1;
  double py = leftcol[0].y = 0;
  while(py!=ys-1)
    {
      // find the next closest
      double md = ys-1-py;
      for(int n = 0 ; n < marks ; n++)
        {
          double d = target_y[n]-py;
          if (d>equality_distance && d < md) md = d;
        }
      py+=md;
      leftcol[nsy++].y=py;
    }
  printf("%d:%d\n",nsx,nsy);
  nsx-=2;
  nsy-=2;
  Array<2,NetPoint> N = net(From<2>(1,1),Size<2>(nsx-2,nsy-2));
  Array<2,NetPoint>::positions it(N);
  while(it.more())
    {
      // take the x and y value of the resptive
      NetPoint & np = it;
      NetPoint & left = N[Position<2>(0,it[1])];
      NetPoint & top  = N[Position<2>(it[0],0)];
      if (it[0]>0)
        {
          np.y=left.y;
          np.left = &N[Position<2>(it[0]-1,it[1])];
        }
      if (it[0]<nsx-1)
        {
          np.right = &N[Position<2>(it[0]+1,it[1])];
        }
      if (it[1]>0)
        {
          np.x=top.x;
          np.up = &N[Position<2>(it[0],it[1]-1)];
        }
      if (it[1]<nsy-1)
        {
          np.down = &N[Position<2>(it[0],it[1]+1)];
        }
      if (it[0]==0 || it[1]==0 ||
          it[0]==nsx - 1 || it[1]==nsy-1)
        np.setNotFree();
      // is there a matching point ?
      for(int m = 0 ; m < marks ; m++)
        {
          if (fabs(target_x[m]-np.x)<=equality_distance &&
              fabs(target_y[m]-np.y)<=equality_distance)
            {
              assert(np.free);
              np.x=target_x[m];
              np.y=target_y[m];
              np.shouldbe_x=source_x[m];
              np.shouldbe_y=source_y[m];
              np.free = false;
            }
        }
      ++it;
    }
  // now we oscilate the net until we find a fixpoint
  int con = iteraties;
  do
    {
      Array<2,NetPoint>::values v(N);
      while(v.more())
        {
          NetPoint &np=v;
          np.move(false);
          ++v;
        }
      Array<2,NetPoint>::values w(N);
      while(w.more())
        {
          NetPoint &np=w;
          np.realize();
          ++w;
        }
    }
  while(con--);
  
  for(Array<2,NetPoint>::values v(N); v.more(); ++v)
    v.current()->mark_as_tx();
  printf("Finished marking target\n");
  {
    Array<2,NetPoint>::values v(N);
    while(v.more())
      {
        NetPoint &np=v;
        np.fixx();
        ++v;
      }
  }
  {
    Array<2,NetPoint>::values v(N);
    while(v.more())
      {
        NetPoint &np=v;
        np.fixy();
        ++v;
      }
  }
  // now we oscilate the net until we find a fixpoint
    con=iteraties;
  do
    {
      Array<2,NetPoint>::values v(N);
      while(v.more())
        {
          NetPoint &np=v;
          np.move(true);
          ++v;
        }
      Array<2,NetPoint>::values w(N);
      while(w.more())
        {
          NetPoint &np=w;
          np.realize();
          ++w;
        }
    }
  while(--con);
  
  printf("Drawing triangles\n");

  
  // print the result net
  for(Array<2,NetPoint>::values v(N); v.more(); ++v)
    v.current()->set_bound(xs,ys);

  
  ImageValue * A = new ImageValue(xs,ys,3);
  for(Array<2,NetPoint>::values v(N); v.more(); ++v)
    {
      NetPoint &np=v;
      A->set(np.a,np.b,0,255);
      A->set(np.x, np.y, 1,255);
      if (!np.free)
        {
          A->set(np.a,np.b,2,255);
          A->set(np.x,np.y,2,255);
        }
    }
  //  return A;
  
        
  
// create the output image

  ImageValue * O = new ImageValue(xs,ys,ch);
  Array<2,NetPoint>::values boxes(N);
  while(boxes.more())
    {
      NetPoint & tl = boxes;
      ++boxes;
      if (!tl.right) continue;
      if (!tl.down) continue;
      if (!tl.right->down) continue;
      NetPoint & br = *tl.right->down;
      NetPoint & bl = *tl.down;
      NetPoint & tr = *tl.right;
      
      // a,b is the target image
      // x,y is the source image
      // h,v are the virtual positions
      double dh = 1.0/max(dist(tl.a,tl.b,tr.a,tr.b),
                          dist(bl.a,bl.b,br.a,br.b));
      double dv = 1.0/max(dist(tl.a,tl.b,bl.a,bl.b),
                          dist(tr.a,tr.b,br.a,br.b));
      for(double h = 0 ; h < 1.0 ; h+=dh)
        {       
          double Ta, Tb;
          double Ba, Bb;
          pos(tl.a,tl.b,tr.a,tr.b,h,Ta,Tb);
          pos(bl.a,bl.b,br.a,br.b,h,Ba,Bb);

          double Tx, Ty;
          double Bx, By;
          pos(tl.x,tl.y,tr.x,tr.y,h,Tx,Ty);
          pos(bl.x,bl.y,br.x,br.y,h,Bx,By);
          for(double v = 0 ; v < 1.0 ; v+=dv)
            {
              double a,b,x,y;
              pos(Ta,Tb,Ba,Bb,v,a,b);
              pos(Tx,Ty,Bx,By,v,x,y);

              for(int c = 0 ; c < ch ; c++)
                O->set(a,b,c,I->get(x,y,c));
            }
        }
    }
  return O;

}
*/

QString type_name(NetPoint)
{
  return "NetPoint-type";
}

/*
  // we now go through all triangles and point the target pixels
  Array<2,NetPoint>::positions it(net);
  while(it.more())
    {
      // take the x and y value of the resptive
      NetPoint & np = it;
      if (
      NetPoint & left = net[Position<2>(0,it[1])];
      NetPoint & top  = net[Position<2>(it[0],0)];
      if (it[0]>0)
        {
          np.y=left.y;
          np.left = &net[Position<2>(it[0]-1,it[1])];
        }
      if (it[0]<nsx)
        {
          np.right = &net[Position<2>(it[0]+1,it[1])];
        }
      if (it[1]>0)
        {
          np.x=top.x;
          np.up = &net[Position<2>(it[0],it[1]-1)];
        }
      if (it[1]<nsy)
        {
          np.down = &net[Position<2>(it[0],it[1]+1)];
        }
      // is there a mathcing point ?
      for(int m = 0 ; m < marks ; m++)
        {
          if (target_x[m]==np.x &&
              target_y[m]==np.y)
            {
              np.shouldbe_x=source_x[m];
              np.shouldbe_y=source_y[m];
              np.free = false;
            }
        }
      ++it;
    }
  
  // calculate the x addresses
  ImageValue * A = new ImageValue(xs,ys,3);
  for(int a = 0 ; a < xs ; a++)
    {
      for(int m = 0 , m < marks ; m++)
        {
          double tx = target_x[m];
          double sx;
          if (a < tx)
            {
              double pxt = a/tx;
              sx = 
            }
        }
    }
  // we first sort the points according to y
  
  int idx[marks];
  qsort(idx,marks,sizeof(int),compare_element);

  ImageValue * B = new ImageValue(xs,ys,3);
  // calculate the target addresses
  // in every step we set the know positions and then 
  // the field is calculated to fill in the remaining data
  // until a stable point is reached
  for(int a = 0 ; a < xs ; a++)
    for(int b = 0 ; b < ys ; b++)
      {
        A->set(a,b,0,a);
        A->set(a,b,1,b);
        B->set(a,b,0,a);
        B->set(a,b,1,b);
      }
  for(int i = 0 ; i < 100 ; i ++)
    {
      printf("%d\n",i);
      for(int m = 0 ; m < marks ; m++)
        {
          A->set(target_x[m],target_y[m],0,source_x[m]);
          A->set(target_x[m],target_y[m],1,source_y[m]);
        }
      for(int a = 1 ; a < xs-1 ; a++)
        for(int b = 1 ; b < ys-1 ; b++)
          {
            double sum_x=0;
            double sum_y=0;
            for(int x = - 1 ; x <= 1 ; x++)
              for(int y = -1 ; y <= 1 ; y++)
                {
                  // we subtract x and y because
                  // if we take the pixel left then we must assume 
                  // that pixel ab will be one position higher
                  sum_x += A->get(a+x,b+y,0)-x-y;
                  sum_y += A->get(a+x,b+y,1)-x-y;
                }
            sum_x/=9;
            sum_y/=9;
            B->set(a,b,0,sum_x);
            B->set(a,b,1,sum_y);
          }
      ImageValue * t = A;
      A = B;
      B = t;
    }

  ImageValue * O = new ImageValue(xs,ys,ch);
  for(int a = 0 ; a < xs ; a++)
    for(int b = 0 ; b < ys ; b++)
      {
        int x = A->get(a,b,0);
        int y = A->get(a,b,1);
        //      printf("Fetching %d:%d from %d:%d\n",a,b,x,y);
        if (x<0) x = 0;
        if (y<0) y = 0;
        if (x>=xs) x = xs-1;
        if (y>=ys) y = ys-1;
        // copy pixel
        for(int c = 0 ; c < ch ; c++)
          O->set(a,b,c,I->get(x,y,c));
      }
  return O;
}
*/

Value * Runtime::morph(vector<Value*> v)
{
  if (v.size()!=3) return NULL;
  return morph(v[0],v[1],v[2]);
}

int side(double a, double b, double x1, double y1, double x2, double y2)
{
  return sign( (a-x1)*(y2-y1) - (b-y1)*(x2-x1) );
}


/*Value * Runtime::morph(Value * image, Value * source, Value * target)
{
  // input checks
  if (image==NULL)  return NULL;
  if (source==NULL) return NULL;
  if (target==NULL) return NULL;
  if (typeid(*image)!=typeid(ImageValue)) return NULL;
  if (typeid(*source)!=typeid(ValueArray)) return NULL;
  if (typeid(*target)!=typeid(ValueArray)) return NULL;
  ImageValue * I = (ImageValue*)image;
  ValueArray * S = (ValueArray*)source;
  ValueArray * T = (ValueArray*)target;
  if (S->sx != T->sx) return NULL;
  if (S->sy != T->sy) return NULL;
  if (S->sy != 1) return NULL;
  int marks = S->sx;
  int xs = I->width();
  int ys = I->height();
  int ch  = I->channelcount();
  double target_x[marks];
  double target_y[marks];
  double source_x[marks];
  double source_y[marks];
  for(int i = 0 ; i < marks ; i++)
    {
      Value * sv = S->get(i,0);
      Value * tv = T->get(i,0);
      if (sv==NULL) return NULL;
      if (tv==NULL) return NULL;
      if (typeid(*sv)!=typeid(PointValue)) return NULL;
      if (typeid(*tv)!=typeid(PointValue)) return NULL;
      PointValue * s = (PointValue*)sv;
      PointValue * t = (PointValue*)tv;
      I->ocsToDcs(t,target_x[i],target_y[i]);
      i->ocsToDcs(s,source_x[i],source_y[i]);
    }
  // create the output image
  ImageValue * O = new ImageValue(xs,ys,ch);
  for(int a = 0 ; a < xs ; a++)
    for(int b = 0 ; b < ys ; b++)
      {
        // first closest point
        double md = dist(0,0,xs,ys);
        int p1 = -1;
        for(int m = 0 ; m < marks ; m++)
          {
            double d =  dist(a,b,target_x[m],target_y[m]);
            if (d<md)
              {
                md = d;
                p1 = m;
              }
          }
        // second and third point
        md = dist(0,0,xs,ys);
        md*=md;
        int p2 = -1;
        int p3 = -1;
        for(int m = 0 ; m < marks ; m++)
          {
            if (m==p1) continue;
            for(int n = 0 ; n < marks ; n++)
              {
                if (n == p1) continue;
                if (n == m) continue;
                
                // does the triangle contain (a,b) ?
                int s1 = side(a,b,target_x[p1],target_y[p1],target_x[m],target_y[m]);
                int s2 = side(a,b,target_x[m],target_y[m],target_x[n],target_y[n]);
                int s3 = side(a,b,target_x[n],target_y[n],target_x[p1],target_y[p1]);
                bool ok = (s1>=0 && s2>=0 && s3>=0) || (s1<=0 && s2<=0 && s3<=0);
                if (!ok) continue;
                // weigh the triangle. the smaller the d the better
                double d =  dist(target_x[n],target_y[n],target_x[m],target_y[m])
                  * dist(target_x[n],target_y[n],target_x[p1],target_y[p1]);
                if (d<md)
                  {
                    md = d;
                    p2 = m;
                    p3 = n;
                  }
              }
          }
        if (a>30 && a < xs-30 && b > 30 && b < ys-30)
          {
            assert(p3!=-1);
          }
        if (p3==-1) continue;
        // we find two points on two lines of the target triangle
        double d12 = dist(target_x[p1],target_y[p1],target_x[p2],target_y[p2]);
        double d23 = dist(target_x[p3],target_y[p3],target_x[p2],target_y[p2]);
        double d13 = dist(target_x[p1],target_y[p1],target_x[p3],target_y[p3]);
        int q1,q2,q3;
        if (d12>d23)
          {
            if (d13 > d12) 
              {
                // 1 - 3 is maximum side
                q1 = p1;
                q2 = p3;
                q3 = p2;
              }
            else
              {
                // 1 - 2 is maximum side
                q1 = p1;
                q2 = p2;
                q3 = p3;
              }
          }
        else
          {
            if (d13 > d23)
              {
                // 1 - 3 is maximum side
                q1 = p1;
                q2 = p3;
                q3 = p2;
              }
            else
              {
                // 2 - 3 is maximum side
                q1 = p2;
                q2 = p3;
                q3 = p1;
              }
          }
        // here we determine then which point axis provides highest
        // accuracy
        double ddx = fabs(target_x[q1]-target_x[q2]);
        double ddy = fabs(target_y[q1]-target_y[q2]);
        bool swap = false;
        if (ddx>ddy)
          {
            // we take the x axis
            double sq2 = sign(a - target_x[q2]);
            double sq3 = sign(a - target_x[q3]);
            swap = (sq3>=0 && sq2<=0) || (sq3<=0 && sq2>=0);
          }
        else
          {
            // we must take the y axis for best accuacy
            double sq2 = sign(b - target_y[q2]);
            double sq3 = sign(b - target_y[q3]);
            swap = (sq3>=0 && sq2<=0) || (sq3<=0 && sq2>=0);
          }
        // we now swap them
        if (swap)
          {
            int tmp = q1;
            q1 = q2;
            q2 = tmp;
          }
        // determine x or y positions on segments 1-3 & 1-2
        int x, y;
        {
          double d13;
          double d12;
          if (ddx>ddy)
            {
              d13 = (a - target_x[q1]) / (target_x[q3]-target_x[q1]);
              d12 = (a - target_x[q1]) / (target_x[q2]-target_x[q1]);
            }
          else
            {
              d13 = (b - target_y[q1]) / (target_y[q3]-target_y[q1]);
              d12 = (b - target_y[q1]) / (target_y[q2]-target_y[q1]);
            }
          //      assert(d13>=0);
          //      assert(d13<=1.1);
          //      assert(d12>=0);
          //      assert(d12<=1.1);
          // find the line segment position starting from d12   
          double t12x;
          double t12y;
          double t13x;
          double t13y;
          pos(target_x[q1],target_y[q1],target_x[q2],target_y[q2],d12,t12x,t12y);
          pos(target_x[q1],target_y[q1],target_x[q3],target_y[q3],d13,t13x,t13y);
          double ds = dist(a,b,t12x,t12y)/dist(t13x,t13y,t12x,t12y);
          if (ds>1) continue;
          assert(ds>=0);
          assert(ds<=1);
          double s12x;
          double s12y; 
          double s13x;
          double s13y;
          pos(source_x[q1],source_y[q1],source_x[q2],source_y[q2],d12,s12x,s12y);
          pos(source_x[q1],source_y[q1],source_x[q3],source_y[q3],d13,s13x,s13y);
          double fx,fy;
          pos(s12x,s12y,s13x,s13y,ds,fx,fy);
          x = fx;
          y = fy;
        }
        // within bounds ?
        //      printf("Fetching %d:%d from %d:%d\n",a,b,x,y);
        if (x<0) x = 0;
        if (y<0) y = 0;
        if (x>=xs) x = xs-1;
        if (y>=ys) y = ys-1;
        // copy pixel
        for(int c = 0 ; c < ch ; c++)
          O->set(a,b,c,I->get(x,y,c));
      }
  return O;
}

*/
// triangles ?
// if we make a grid of equal sized corners then it should be possible
// to bend various lines, somewhat like a spring. Those lines can then be used to 
// 

int * getRing(double *x, double *y, int & count, bool* enabled, int points)
{
  int * result = new int[points];
  count = 0;
  double largestx = -1;
  int idx = -1;
  for(int i = 0 ; i < points; i++)
    if (enabled[i] && x[i]>largestx)
      largestx=x[idx=i];
  if (idx==-1) 
    {
      printf("No points left\n");
      return NULL;
    }
  result[count++]=idx;
  printf("pointidx %d\n",idx);
  
  while(true)
    {
      bool added_something = false;
      for(int i = 0 ; i < points; i++)
        {
          if (!enabled[i]) continue;
          bool already_in_result = false;
          for(int k = 0 ; k < count ; k++)
            if (result[k]==i) already_in_result = true;
          if (already_in_result) continue;
          int j;
          for(j = 0 ; j < points; j++)
            {
              if (!enabled[j]) continue;
              if (j==result[count-1]) continue;
              if (j==i) continue;
              if (side(x[j],
                       y[j],
                       x[result[count-1]],
                       y[result[count-1]],
                       x[i],
                       y[i])>0) break;
            }
          if (j==points)
            {
              printf("pointidx %d\n",i);
              added_something = true;
              result[count++]=i;
              assert(count < 1000);
              break;
            }
        }
      if (!added_something) return result;
    }
}

class Triangle
{
public:
  int A;
  int B;
  int C;
};

void newpoint(double ax, double ay, double bx, double by, double cx, double cy, 
              double ix, double iy, double &tx, double &ty)
{
  ax-=bx;
  ay-=by;
  cx-=bx;
  cy-=by;
  tx = ax + cx;
  ty = ay + cy;
  double d = dist(0,0,tx,ty);
  tx /= d;
  ty /= d;
  tx *= 40;
  ty *= 40;
  tx = ix - tx;
  ty = iy - ty;
}

/*

Value * Runtime::morph(Value * image, Value * source, Value * target)
{
  // input checks
  if (image==NULL)  return NULL;
  if (source==NULL) return NULL;
  if (target==NULL) return NULL;
  if (typeid(*image)!=typeid(ImageValue)) return NULL;
  if (typeid(*source)!=typeid(ValueArray)) return NULL;
  if (typeid(*target)!=typeid(ValueArray)) return NULL;
  ImageValue * I = (ImageValue*)image;
  ValueArray * S = (ValueArray*)source;
  ValueArray * T = (ValueArray*)target;
  if (S->sx != T->sx) return NULL;
  if (S->sy != T->sy) return NULL;
  if (S->sy != 1 && S->sx != 1) return NULL;
  int marks = S->sx;
  if (marks==1) marks=S->sy;
  int xs = I->width();
  int ys = I->height();
  int ch  = I->channelcount();
  // the extra 1 element is useful to create the last ring
  double *target_x = new double[marks+1];
  double *target_y = new double[marks+1];
  double *source_x = new double[marks+1];
  double *source_y = new double[marks+1];
  if (S->sx==1)
    for(int i = 0 ; i < marks ; i++)
      {
        Value * sv = S->get(i,0);
        Value * tv = T->get(i,0);
        if (sv==NULL) return NULL;
        if (tv==NULL) return NULL;
        if (typeid(*sv)!=typeid(PointValue)) return NULL;
        if (typeid(*tv)!=typeid(PointValue)) return NULL;
        PointValue * s = (PointValue*)sv;
        PointValue * t = (PointValue*)tv;
        target_x[i]=(t->x + 1)*xs/2.0;
        target_y[i]=(t->y + 1)*ys/2.0;
        source_x[i]=(s->x + 1)*xs/2.0;
        source_y[i]=(s->y + 1)*ys/2.0;
      }
  else
    for(int i = 0 ; i < marks ; i++)
      {
        Value * sv = S->get(0,i);
        Value * tv = T->get(0,i);
        if (sv==NULL) return NULL;
        if (tv==NULL) return NULL;
        if (typeid(*sv)!=typeid(PointValue)) return NULL;
        if (typeid(*tv)!=typeid(PointValue)) return NULL;
        PointValue * s = (PointValue*)sv;
        PointValue * t = (PointValue*)tv;
        target_x[i]=(t->x + 1)*xs/2.0;
        target_y[i]=(t->y + 1)*ys/2.0;
        source_x[i]=(s->x + 1)*xs/2.0;
        source_y[i]=(s->y + 1)*ys/2.0;
      }

  int triangle_count = 0;
  Triangle triangles[1000];
  bool enabled[marks];
  for(int i = 0 ; i < marks ; i++)
    enabled[i]=true;
  int outer_count; 
  int * outer_index = getRing(target_x,target_y,outer_count,enabled,marks);
  // fix the marks array with outer_count extra fields
  if (outer_count > 2)
    {
      double *ntarget_x=new double[marks+outer_count+1];
      double *ntarget_y=new double[marks+outer_count+1];
      double *nsource_x=new double[marks+outer_count+1];
      double *nsource_y=new double[marks+outer_count+1];
      for(int i = 0 ; i < marks ; i ++)
        {
          ntarget_x[i]=target_x[i];
          ntarget_y[i]=target_y[i];
          nsource_x[i]=source_x[i];
          nsource_y[i]=source_y[i];
        }
      target_x=ntarget_x;
      target_y=ntarget_y;
      source_x=nsource_x;
      source_y=nsource_y;
      for(int c = 0 ; c < outer_count ; c++)
        {
          int p = c - 1;
          if (p<0) p=outer_count-1;
          int n = (c + 1)%outer_count;
          newpoint(target_x[outer_index[p]],target_y[outer_index[p]],
                   target_x[outer_index[c]],target_y[outer_index[c]],
                   target_x[outer_index[n]],target_y[outer_index[n]],
                   target_x[outer_index[c]],target_y[outer_index[c]],
                   target_x[c+marks],       target_y[c+marks]);
          newpoint(target_x[outer_index[p]],target_y[outer_index[p]],
                   target_x[outer_index[c]],target_y[outer_index[c]],
                   target_x[outer_index[n]],target_y[outer_index[n]],
                   source_x[outer_index[c]],source_y[outer_index[c]],
                   source_x[c+marks],       source_y[c+marks]);
        }
      marks+=outer_count;
      // restart with extra ring...
      for(int i = 0 ; i < marks ; i++) enabled[i]=true;
      outer_index = getRing(target_x,target_y,outer_count,enabled,marks);
    }

  // here we add a boundary to the outer layer;
  bool last_ring_processed = false;
  while (!last_ring_processed)
    {
      // we disable all these positions
      for(int i = 0 ; i < outer_count ; i++)
        enabled[outer_index[i]]=false;
      // and find the inner ring
      int inner_count;
      int * inner_index = getRing(target_x,target_y,inner_count,enabled,marks);
      // we then create the triangles
      if (!inner_index)
        {
          // find the mean
          double a,x;
          double b,y;
          a=b=x=y=0;
          for(int i = 0 ; i < outer_count; i++)
            {
              a+=target_x[outer_index[i]];
              b+=target_y[outer_index[i]];
              x+=source_x[outer_index[i]];
              y+=source_y[outer_index[i]];
            }
          a/=outer_count;
          b/=outer_count;
          x/=outer_count;
          y/=outer_count;
          target_x[marks]=a;
          target_y[marks]=b;
          source_x[marks]=x;
          source_y[marks]=y;
          last_ring_processed = true;
          inner_count = 1;
          inner_index = new int[1];
          inner_index[0] = marks;
        }
      int oi = 0;
      int ii = 0;
      while(true)
        {
          // the current segment is oi -- ii
          // the next one is either from o to i2
          // or from i to o2
          int oi2 = (oi+1)%outer_count;
          int ii2 = (ii+1)%inner_count;
          double do1i2=dist(target_x[outer_index[oi]],target_y[outer_index[oi]],
                            target_x[inner_index[ii2]],target_y[inner_index[ii2]]);
          double do2i1=dist(target_x[outer_index[oi2]],target_y[outer_index[oi2]],
                            target_x[inner_index[ii]],target_y[inner_index[ii]]);
          int tA;
          int tB;
          int tC;
          if (do1i2<do2i1 && ii!=ii2)
            {
              tA = outer_index[oi];
              tB = inner_index[ii];
              tC = inner_index[ii2];
              oi = oi;
              ii = ii2;
            }
          else
            {
              tA = outer_index[oi];
              tB = inner_index[ii];
              tC = outer_index[oi2];
              oi = oi2;
              ii = ii;
            }
          // do we have the triangle ?
          bool have_triangle = false;
          for(int i = 0 ; i < triangle_count ; i++)
            {
              if (triangles[i].A==tA
                  && triangles[i].B==tB
                  && triangles[i].C==tC)
                have_triangle = true;
            }
          // we add the triangle
          if (have_triangle)
            break;
          triangles[triangle_count].A=tA;
          triangles[triangle_count].B=tB;
          triangles[triangle_count].C=tC;
          triangle_count++;
          assert(triangle_count<1000);
          //      if (last_ring_processed)
            printf("%d - %d\n",oi,ii);
        }
      printf("triangle_count %d\n",triangle_count);
      outer_index = inner_index;
      outer_count = inner_count;
    }
  
  // create the output image
  ImageValue * O = new ImageValue(xs,ys,ch);
  for(int triangle = 0 ; triangle< triangle_count ; triangle++)
    {
      Triangle t = triangles[triangle];
      double Aa = target_x[t.A];
      double Ab = target_y[t.A];
      double Ax = source_x[t.A];
      double Ay = source_y[t.A];

      double Ba = target_x[t.B];
      double Bb = target_y[t.B];
      double Bx = source_x[t.B];
      double By = source_y[t.B];

      double Ca = target_x[t.C];
      double Cb = target_y[t.C];
      double Cx = source_x[t.C];
      double Cy = source_y[t.C];

      // a,b is the target image
      // x,y is the source image
      // h,v are the virtual positions
      double dh = 1.0/max(dist(Aa,Ab,Ba,Bb),
                          dist(Ca,Cb,Ca,Cb));
      double dv = 1.0/max(dist(Aa,Ab,Ca,Cb),
                          dist(Ba,Bb,Ca,Cb));
      for(double h = 0 ; h < 1.0 ; h+=dh)
        {       
          double Ta, Tb;
          pos(Aa,Ab,Ba,Bb,h,Ta,Tb);

          double Tx, Ty;
          pos(Ax,Ay,Bx,By,h,Tx,Ty);
          for(double v = 0 ; v < 1.0 ; v+=dv)
            {
              double a,b,x,y;
              pos(Ta,Tb,Ca,Cb,v,a,b);
              pos(Tx,Ty,Cx,Cy,v,x,y);

              if (a<0) continue;
              if (b<0) continue;
              if (x<0) continue;
              if (y<0) continue;
              if (a>=xs) continue;
              if (b>=ys) continue;
              if (x>=xs) continue;
              if (y>=ys) continue;

              for(int c = 0 ; c < ch ; c++)
                O->set(a,b,c,I->get(x,y,c));
            }
        }
    }
  return O;
}

*/

bool Q1(double a)
{
  while (a<0) a+=2*M_PI;
  while (a>=2*M_PI) a-=2*M_PI;
  bool r = a<=M_PI/2;
  return r;
}

bool Q3(double a)
{
  while (a<0) a+=2*M_PI;
  while (a>=2*M_PI) a-=2*M_PI;
  bool r = a>=3*M_PI/2;
  return r;
}

double calc_step(double x1, double y1, double x2, double y2)
{
  double d = dist(x1,y1,x2,y2);
  if (d==0) return 1;
  return 0.5/d;
}

Value * Runtime::morph(Value * image, Value * source, Value * target)
{
  // input checks
  if (image==NULL)  return NULL;
  if (source==NULL) return NULL;
  if (target==NULL) return NULL;
  if (typeid(*image)!=typeid(ImageValue)) return NULL;
  if (typeid(*source)!=typeid(ValueArray)) return NULL;
  if (typeid(*target)!=typeid(ValueArray)) return NULL;
  ImageValue * I = (ImageValue*)image;
  ValueArray * S = (ValueArray*)source;
  ValueArray * T = (ValueArray*)target;
  if (S->sx != T->sx) return NULL;
  if (S->sy != T->sy) return NULL;
  if (S->sy != 1 && S->sx != 1) return NULL;
  int marks = S->sx;
  if (marks==1) marks=S->sy;
  int xs = I->width_direct();
  int ys = I->height_direct();
  int ch  = I->channelcount();
  // the extra 1 element is useful to create the last ring
  double *target_x = new double[marks+1];
  double *target_y = new double[marks+1];
  double *source_x = new double[marks+1];
  double *source_y = new double[marks+1];
  if (S->sx==1)
    for(int i = 0 ; i < marks ; i++)
      {
        Value * sv = S->get(i,0);
        Value * tv = T->get(i,0);
        if (sv==NULL) return NULL;
        if (tv==NULL) return NULL;
        if (typeid(*sv)!=typeid(PointValue)) return NULL;
        if (typeid(*tv)!=typeid(PointValue)) return NULL;

        PointValue * s = (PointValue*)sv;
        PointValue * t = (PointValue*)tv;
        I->ocsToDcs(*t,target_x[i],target_y[i]);
        I->ocsToDcs(*s,source_x[i],source_y[i]);
      }
  else
    for(int i = 0 ; i < marks ; i++)
      {
        Value * sv = S->get(0,i);
        Value * tv = T->get(0,i);
        if (sv==NULL) return NULL;
        if (tv==NULL) return NULL;
        if (typeid(*sv)!=typeid(PointValue)) return NULL;
        if (typeid(*tv)!=typeid(PointValue)) return NULL;
        PointValue * s = (PointValue*)sv;
        PointValue * t = (PointValue*)tv;
        I->ocsToDcs(*t,target_x[i],target_y[i]);
        I->ocsToDcs(*s,source_x[i],source_y[i]);
      }

  bool enabled[marks];
  for(int i = 0 ; i < marks ; i++)
    enabled[i]=true;
  int outer_count; 
  int * outer_index = getRing(target_x,target_y,outer_count,enabled,marks);
  ImageValue * O = new ImageValue(xs,ys,ch);
  while (true)
    {
      // find the middle of the outer circle...
      double ca,cx;
      double cb,cy;
      ca=cb=cx=cy=0;
      for(int i = 0 ; i < outer_count; i++)
        {
          ca+=target_x[outer_index[i]];
          cb+=target_y[outer_index[i]];
          cx+=source_x[outer_index[i]];
          cy+=source_y[outer_index[i]];
        }
      ca/=outer_count;
      cb/=outer_count;
      cx/=outer_count;
      cy/=outer_count;
      // we disable all the outer positions
      for(int i = 0 ; i < outer_count ; i++)
        enabled[outer_index[i]]=false;
      // and find the inner ring
      int inner_count;
      int * inner_index = getRing(target_x,target_y,inner_count,enabled,marks);
      // go through 360 degree and fill the gap between the outer and inner circle

      // we don't go through all possible angles, instead we iterate over every outer segment
      // and try to figure out a matching inner segment
      for(int s = 0 ; s < outer_count ; s++)
        {
          int n = (s+1)%outer_count;
          double ab = target_x[outer_index[s]];
          double bb = target_y[outer_index[s]];
          double ae = target_x[outer_index[n]];
          double be = target_y[outer_index[n]];
          double xb = source_x[outer_index[s]];
          double yb = source_y[outer_index[s]];
          double xe = source_x[outer_index[n]];
          double ye = source_y[outer_index[n]];
          double pct_out_delta = calc_step(ab,bb,ae,be);
          for(double pct_out = 0 ; pct_out<=1.0 ; pct_out+=pct_out_delta)
            {
              double a_out,b_out;
              double x_out,y_out;
              double angle = atan2(b_out-cb,a_out-ca);
              pos(ab,bb,ae,be,pct_out,a_out,b_out);
              pos(xb,yb,xe,ye,pct_out,x_out,y_out);
              
              // we need to find the correct inner segment
              double a_in=ca,b_in=cb;
              double x_in=cx, y_in=cy;
              for(int s = 0 ; s < inner_count ; s++)
                {
                  int n = (s+1)%inner_count;
                  double ab = target_x[inner_index[s]];
                  double bb = target_y[inner_index[s]];
                  double ae = target_x[inner_index[n]];
                  double be = target_y[inner_index[n]];
                  double angle_b = atan2(bb-cb,ab-ca);
                  double angle_e = atan2(be-cb,ae-ca);
                  angle_b-=angle;
                  angle_e-=angle;
                  if ( (Q1(angle_b) && Q3(angle_e)) ||
                       (Q3(angle_b) && Q1(angle_e)))
                    {
                      double db = dist(ab,bb,ca,cb);
                      double de = dist(ae,be,ca,cb);
                      // we subtract the angle
                      double x1 = cos(angle_b)*db;
                      double x2 = cos(angle_e)*de;
                      double y1 = sin(angle_b)*db;
                      double y2 = sin(angle_e)*de;
                      double din = (-y1*(x2-x1)/(y2-y1))+x1;
                      a_in = cos(angle)*din + ca;
                      b_in = sin(angle)*din + cb;
                      double pct_in = dist(a_in,b_in,ab,bb)/dist(ab,bb,ae,be);
                      double xb = source_x[inner_index[s]];
                      double yb = source_y[inner_index[s]];
                      double xe = source_x[inner_index[n]];
                      double ye = source_y[inner_index[n]];
                      pos(xb,yb,xe,ye,pct_in,x_in,y_in);
                      break;
                    }
                }
              // iterate through the positions
              double pct_delta = calc_step(a_out,b_out,a_in,b_in);
              for(double pct = 0 ; pct<=1.0 ; pct+=pct_delta)
                {
                  double a,b;
                  double x,y;
                  pos(a_out,b_out,a_in,b_in,pct,a,b);
                  pos(x_out,y_out,x_in,y_in,pct,x,y);
                  
                  if (a<0) continue;
                  if (b<0) continue;
                  if (x<0) continue;
                  if (y<0) continue;
                  if (a>=xs) continue;
                  if (b>=ys) continue;
                  if (x>=xs) continue;
                  if (y>=ys) continue;
                  
                  for(int c = 0 ; c < ch ; c++)
                    O->set_direct(a,b,c,I->get_direct(x,y,c));
                }
            }
        }
      if (inner_count==0 || inner_index==NULL) break;
      for(int s = 0 ; s < inner_count ; s++)
        {
          int n = (s+1)%inner_count;
          double ab = target_x[inner_index[s]];
          double bb = target_y[inner_index[s]];
          double ae = target_x[inner_index[n]];
          double be = target_y[inner_index[n]];
          double xb = source_x[inner_index[s]];
          double yb = source_y[inner_index[s]];
          double xe = source_x[inner_index[n]];
          double ye = source_y[inner_index[n]];
          double pct_out_delta = calc_step(ab,bb,ae,be);
          for(double pct_out = 0 ; pct_out<=1.0 ; pct_out+=pct_out_delta)
            {
              double a_out,b_out;
              double x_out,y_out;
              double angle = atan2(b_out-cb,a_out-ca);
              pos(ab,bb,ae,be,pct_out,a_out,b_out);
              pos(xb,yb,xe,ye,pct_out,x_out,y_out);
              
              // we need to find the correct outer segment
              double a_in,b_in;
              double x_in, y_in;
              for(int s = 0 ; s < outer_count ; s++)
                {
                  int n = (s+1)%outer_count;
                  double ab = target_x[outer_index[s]];
                  double bb = target_y[outer_index[s]];
                  double ae = target_x[outer_index[n]];
                  double be = target_y[outer_index[n]];
                  double angle_b = atan2(bb-cb,ab-ca);
                  double angle_e = atan2(be-cb,ae-ca);
                  angle_b-=angle;
                  angle_e-=angle;
                  if ( (Q1(angle_b) && Q3(angle_e)) ||
                       (Q3(angle_b) && Q1(angle_e)))
                    {
                      double db = dist(ab,bb,ca,cb);
                      double de = dist(ae,be,ca,cb);
                      // we subtract the angle
                      double x1 = cos(angle_b)*db;
                      double x2 = cos(angle_e)*de;
                      double y1 = sin(angle_b)*db;
                      double y2 = sin(angle_e)*de;
                      double din = (-y1*(x2-x1)/(y2-y1))+x1;
                      a_in = cos(angle)*din + ca;
                      b_in = sin(angle)*din + cb;
                      double pct_in = dist(a_in,b_in,ab,bb)/dist(ab,bb,ae,be);
                      double xb = source_x[outer_index[s]];
                      double yb = source_y[outer_index[s]];
                      double xe = source_x[outer_index[n]];
                      double ye = source_y[outer_index[n]];
                      pos(xb,yb,xe,ye,pct_in,x_in,y_in);
                      break;
                    }
                }
              // iterate through the positions
              double pct_delta = calc_step(a_out,b_out,a_in,b_in);
              for(double pct = 0 ; pct<=1.0 ; pct+=pct_delta)
                {
                  double a,b;
                  double x,y;
                  pos(a_out,b_out,a_in,b_in,pct,a,b);
                  pos(x_out,y_out,x_in,y_in,pct,x,y);
                  
                  if (a<0) continue;
                  if (b<0) continue;
                  if (x<0) continue;
                  if (y<0) continue;
                  if (a>=xs) continue;
                  if (b>=ys) continue;
                  if (x>=xs) continue;
                  if (y>=ys) continue;
                  
                  for(int c = 0 ; c < ch ; c++)
                    O->set_direct(a,b,c,I->get_direct(x,y,c));
                }
            }
        }
      outer_index = inner_index;
      outer_count = inner_count;
    }
return O;
}

---