''' MorphMapper.py

This file defines a class, called MorphMapper,
whose instances are callable objects that compute
a kind of affine mappings of the plane.  They are
components of the warping method used for morphing
by Beier and Neely.

'''

import math

class MorphMapper:
  # static constants used in the weighting formula:
  a = 0.01
  b = 1
  p = 0.5

  def __init__(self, Pd, Qd, Ps, Qs):
    # unpack the points
    (self.xpd, self.ypd) = Pd
    (self.xqd, self.yqd) = Qd
    (self.xps, self.yps) = Ps
    (self.xqs, self.yqs) = Qs
    # make vectors from the segments:
    (self.destDX, self.destDY) = (self.xqd - self.xpd, self.yqd - self.ypd)
    (self.srcDX,  self.srcDY)  = (self.xqs - self.xps, self.yqs - self.yps)
    # find parameters for rotating dest. vec. to align with Y axis:
    self.destHypot = math.sqrt(self.destDX**2 + self.destDY**2)
    if self.destHypot < 1.0: self.destHypot = 1.0 # to prevent zeroDivision errors.
    self.costheta = self.destDY / self.destHypot
    self.sintheta = self.destDX / self.destHypot
    perpX = self.srcDY
    perpY = - self.srcDX
    self.srcHypot = math.sqrt(self.srcDX**2 + self.srcDY**2)
    if self.srcHypot < 1.0: self.srcHypot = 1.0 # to prevent zeroDevision errors
    self.unitPerpX = perpX / self.srcHypot
    self.unitPerpY = perpY / self.srcHypot
    #print "(destHypot, srcHypot)+"+str((self.destHypot, srcHypot))

  def __call__(self, x,y):
    # translate (x,y) by -Pd
    xt = x - self.xpd
    yt = y - self.ypd
    # rotate
    xr = self.costheta*xt - self.sintheta*yt
    yr = self.sintheta*xt + self.costheta*yt
    u = yr / self.destHypot
    v = xr
    #print "(u,v)="+str((u,v)),
    xs = self.xps + (u*self.srcDX) + (v*self.unitPerpX)
    ys = self.yps + (u*self.srcDY) + (v*self.unitPerpY)
    # compute distance from (x,y) to the dest. line segment.
    dist = abs(v) # default case
    if   u<0: dist = math.sqrt((x - self.xpd)**2 + (y - self.ypd)**2)
    elif u>1: dist = math.sqrt((x - self.xqd)**2 + (y - self.yqd)**2)
    # compute weight according to Beier & Neely method
    weight = (self.destHypot**MorphMapper.p / (MorphMapper.a + dist))**MorphMapper.b
    return (xs, ys, dist, weight)

def test():
  Pd = (0,0); Qd=(1,0)
  Ps = (5,5); Qs=(6,5)
  X = (13,18)

  map1 = MorphMapper(Pd, Qd, Ps, Qs)

  for x in range(3):
    for y in range(3):
      (xs,ys,dist,weight) = map1(x,y)
      print str((x,y))+" maps to "+str((xs,ys))+"; distance="+str(dist)+\
        "; weight="+str(weight)