AEP.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from vector_utils import sph2cart, cart2sph"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def expmap_aep_smith(base_vectors, tangent_vectors):\n",
      "    # If we've been passed a single vector to map, then add the extra axis\n",
      "    # Number of sample first\n",
      "    if len(tangent_vectors.shape) == 3:\n",
      "        vector_count = tangent_vectors.shape[0]\n",
      "    else:\n",
      "        tangent_vectors = tangent_vectors[None, ...]\n",
      "        vector_count = 1\n",
      "\n",
      "    N = tangent_vectors.shape[1]\n",
      "    ns = np.zeros([vector_count, N, tangent_vectors.shape[2] + 1])\n",
      "    eps = np.spacing(1)\n",
      "    \n",
      "    kset = tangent_vectors\n",
      "    \n",
      "    long0lat1 = cart2sph(-base_vectors[:, 1], base_vectors[:, 0], base_vectors[:, 2])\n",
      "    long0lat1[:, 0] += np.pi\n",
      "    long0lat1 = long0lat1[None, ...]\n",
      " \n",
      "    rho = np.sqrt(kset[..., 0] ** 2 + kset[..., 1] ** 2)\n",
      "    Z = np.exp(np.arctan2(kset[..., 1], kset[..., 0]) * 1j)\n",
      "    V1 = rho * np.real(Z)\n",
      "    V2 = rho * np.imag(Z)\n",
      "    \n",
      "    # To prevent divide by 0 warnings when rho is 0\n",
      "    ir = rho == 0\n",
      "    rho[ir] = eps\n",
      "    c = rho\n",
      "    \n",
      "    c[ir] = eps\n",
      "    \n",
      "    Y = np.real(np.arcsin(np.cos(c) * np.sin(long0lat1[..., 1]) + np.cos(long0lat1[..., 1]) * np.sin(c) * V2 / rho))\n",
      "    X = np.real(long0lat1[..., 0] + np.arctan2(V1 * np.sin(c), np.cos(long0lat1[..., 1]) * np.cos(c) * rho - np.sin(long0lat1[..., 1]) * V2 * np.sin(c)))\n",
      "    \n",
      "    ns = sph2cart(X - np.pi, Y, np.ones_like(Y))\n",
      "    # Swap x and y axes\n",
      "    ns[..., [0, 1]] = ns[..., [1, 0]]\n",
      "    ns[..., 1] = -ns[..., 1]\n",
      "        \n",
      "    return ns\n",
      " \n",
      "def logmap_aep_smith(base_vectors, sd_vectors):\n",
      "    # If we've been passed a single vector to map, then add the extra axis\n",
      "    # Number of sample first\n",
      "    if len(sd_vectors.shape) == 3:\n",
      "        vector_count = sd_vectors.shape[0]\n",
      "    else:\n",
      "        sd_vectors = sd_vectors[None, ...]\n",
      "        vector_count = 1\n",
      "    \n",
      "    N = sd_vectors.shape[1]    \n",
      "    vs = np.zeros([vector_count, N, sd_vectors.shape[2] - 1])\n",
      "    \n",
      "    kset = sd_vectors\n",
      "    \n",
      "    longlat = cart2sph(-kset[..., 1], kset[..., 0], kset[..., 2])\n",
      "    longlat[..., 0] += np.pi\n",
      "    \n",
      "    long0lat1 = cart2sph(-base_vectors[:, 1], base_vectors[:, 0], base_vectors[:, 2]) \n",
      "    long0lat1[:, 0] += np.pi\n",
      "    long0lat1 = long0lat1[None, ...]\n",
      "    \n",
      "    c = np.arccos(np.sin(long0lat1[..., 1]) * np.sin(longlat[..., 1]) + np.cos(long0lat1[..., 1]) * np.cos(longlat[..., 1]) * np.cos(longlat[..., 0] - long0lat1[..., 0]))\n",
      "    \n",
      "    k = 1.0 / np.sin(c)\n",
      "    k = np.nan_to_num(k)\n",
      "    k = c * k\n",
      "    \n",
      "    vs[..., 0] = k * (np.cos(longlat[..., 1]) * np.sin(longlat[..., 0] - long0lat1[..., 0]))\n",
      "    vs[..., 1] = k * (np.cos(long0lat1[..., 1]) * np.sin(longlat[..., 1]) - np.sin(long0lat1[..., 1]) * np.cos(longlat[..., 1]) * np.cos(longlat[..., 0] - long0lat1[..., 0]))\n",
      "                                                                                                           \n",
      "    return vs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def logmap_aep_snape(base_vectors, sd_vectors):\n",
      "    # If we've been passed a single vector to map, then add the extra axis\n",
      "    # Number of sample first\n",
      "    if len(sd_vectors.shape) < 3:\n",
      "        sd_vectors = sd_vectors[None, ...]\n",
      "\n",
      "    thetaav = elevation(base_vectors[:, 2])\n",
      "    phiav = azimuth(base_vectors[:, 0], base_vectors[:, 1])\n",
      "    phiav, thetaav = clip_angles(phiav, thetaav)\n",
      "\n",
      "    # find any zero normals as they present a real problem\n",
      "    # the column indicies are the same in both sets of data\n",
      "    zero_indices = np.sum(np.abs(sd_vectors), axis=-1) == 0.0\n",
      "\n",
      "    thetak = elevation(sd_vectors[..., 2])\n",
      "    phik = azimuth(sd_vectors[..., 0], sd_vectors[..., 1])\n",
      "    phik, thetak = clip_angles(phik, thetak)\n",
      "\n",
      "    # cos(c) = sin(thetaav) * sin(thetak) +\n",
      "    #          cos(thetaav) * cos(thetak) * cos[phik - phiav]\n",
      "    cosc = (np.sin(thetaav) * np.sin(thetak) +\n",
      "            np.cos(thetaav) * np.cos(thetak) * np.cos(phik - phiav))\n",
      "    c = np.arccos(cosc)\n",
      "    # kprime = c / sin(c)\n",
      "    kprime = c / np.sin(c)\n",
      "\n",
      "    # xs = kprime * cos(thetak) * sin[phik - phiav]\n",
      "    xs = kprime * np.cos(thetak) * np.sin(phik - phiav)\n",
      "    # ys = kprime * (cos(thetaav) * sin(phik) -\n",
      "    #      sin(thetaav) * cos(thetak) * cos[phik - phiav]\n",
      "    ys = kprime * (np.cos(thetaav) * np.sin(thetak) -\n",
      "                   np.sin(thetaav) * np.cos(thetak) * np.cos(phik - phiav))\n",
      "    vs = np.dstack([xs[..., None], ys[..., None]])\n",
      "\n",
      "    # reset the zero normals back to 0\n",
      "    vs[zero_indices, :] = 0.0\n",
      "\n",
      "    return vs\n",
      "\n",
      "\n",
      "def expmap_aep_snape(base_vectors, tangent_vectors):\n",
      "    # If we've been passed a single vector to map, then add the extra axis\n",
      "    # Number of sample first\n",
      "    if len(tangent_vectors.shape) < 3:\n",
      "        tangent_vectors = tangent_vectors[None, ...]\n",
      "\n",
      "    thetaav = elevation(base_vectors[:, 2])\n",
      "    phiav = azimuth(base_vectors[:, 0], base_vectors[:, 1])\n",
      "    phiav, thetaav = clip_angles(phiav, thetaav)\n",
      "\n",
      "    # find any zero normals as they present a real problem\n",
      "    # the column indicies are the same in both sets of data\n",
      "    zero_indices = np.sum(np.abs(tangent_vectors), axis=-1) == 0.0\n",
      "\n",
      "    xs = tangent_vectors[..., 0]\n",
      "    ys = tangent_vectors[..., 1]\n",
      "\n",
      "    c = np.sqrt(xs ** 2 + ys ** 2)\n",
      "    recipc = 1.0 / c\n",
      "\n",
      "    # thetas = asin[cos(c) * sin(thetaav) - (1/c) * yk * sin(c) * cos(thetav)]\n",
      "    s = np.cos(c) * np.sin(thetaav) + recipc * ys * np.sin(c) * np.cos(thetaav)\n",
      "\n",
      "    # Elevation\n",
      "    el = np.arcsin(s)\n",
      "    # Azimuth\n",
      "    # phis = phiav + atan(psi)\n",
      "    phis = psi_new(c, thetaav, xs, ys)\n",
      "    azi = phiav + np.arctan2(phis[..., 0], phis[..., 1])\n",
      "    azi, el = clip_angles(azi, el)\n",
      "\n",
      "    # convert angles to coordinates\n",
      "    xs = np.cos(azi) * np.sin(el)\n",
      "    ys = np.sin(azi) * np.sin(el)\n",
      "    zs = np.cos(el)\n",
      "    ns = np.dstack([xs[..., None], ys[..., None], zs[..., None]])\n",
      "\n",
      "    # reset zero projections back to zero\n",
      "    ns[zero_indices, :] = 0.0\n",
      "\n",
      "    return ns\n",
      "\n",
      "\n",
      "def clip_angles(phi, theta):\n",
      "    # round thetas back in to the range [-pi/2, pi/2]\n",
      "    pi_over_2 = np.pi / 2.0\n",
      "    theta[theta > pi_over_2] -= np.pi\n",
      "    theta[theta < -pi_over_2] += np.pi\n",
      "\n",
      "    phi[phi > np.pi] -= 2.0 * np.pi\n",
      "    phi[phi < -np.pi] += 2.0 * np.pi\n",
      "\n",
      "    return phi, theta\n",
      "\n",
      "\n",
      "# psi = thetaav != (pi / 2) ->\n",
      "#           xk * sin(c) / c * cos(thetaav) * cos(c) - yk * sin(thetav) * sin(c)\n",
      "#       thetaav == (pi/2)   -> -(xk/yk)\n",
      "#       thetaav == -(pi/2)  -> xk/yk\n",
      "def psi_new(c, thetaav, xs, ys):\n",
      "    eps = np.spacing(1)\n",
      "    out = np.zeros([xs.shape[0], thetaav.size, 2])\n",
      "\n",
      "    neg_pi_ind = np.abs(thetaav - (np.pi / 2.0)) < eps\n",
      "    if np.any(neg_pi_ind):\n",
      "        out[:, neg_pi_ind, :] = [-xs[:, neg_pi_ind], -ys[:, neg_pi_ind]]\n",
      "\n",
      "    pos_pi_ind = np.abs(thetaav + (np.pi / 2.0)) < eps\n",
      "    if np.any(pos_pi_ind):\n",
      "        out[:, pos_pi_ind, :] = [xs[:, pos_pi_ind], ys[:, pos_pi_ind]]\n",
      "\n",
      "    ind = ~np.logical_and(neg_pi_ind, pos_pi_ind)\n",
      "    indexed_c = c[:, ind]\n",
      "    indexed_thetaav = thetaav[ind]\n",
      "    out[:, ind, 0] = xs[:, ind] * np.sin(indexed_c)\n",
      "    out[:, ind, 1] = (indexed_c * np.cos(indexed_thetaav) * np.cos(indexed_c) -\n",
      "                      ys[:, ind] * np.sin(indexed_thetaav) * np.sin(indexed_c))\n",
      "\n",
      "    return out\n",
      "\n",
      "\n",
      "# theta = (pi / 2) - asin(nz)\n",
      "def elevation(zs):\n",
      "    return (np.pi / 2.0) - np.arcsin(zs)\n",
      "\n",
      "\n",
      "# phi = atan(ny,nx)\n",
      "def azimuth(xs, ys):\n",
      "     return np.arctan2(ys, xs)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from vector_utils import normalise_vector\n",
      "\n",
      "N = normalise_vector(np.random.uniform(low=-1.0, high=1.0, size=(100, 1000, 3)))\n",
      "bp = normalise_vector(np.random.uniform(low=-1.0, high=1.0, size=(1000, 3)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "snape = expmap_aep_snape(bp, logmap_aep_snape(bp, N))\n",
      "smith = expmap_aep_smith(bp, logmap_aep_smith(bp, N))\n",
      "print np.allclose(snape, N)\n",
      "print np.allclose(smith, N)\n",
      "np.allclose(snape, smith)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit expmap_aep_snape(bp, logmap_aep_snape(bp, N))\n",
      "%timeit expmap_aep_smith(bp, logmap_aep_smith(bp, N))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}