Ferret Neocortical Folding.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: pylab import has clobbered these variables: ['normal', 'pylab']\n",
      "`%matplotlib` prevents importing * from pylab and numpy\n"
     ]
    }
   ],
   "source": [
    "%pylab inline\n",
    "\n",
    "import numpy\n",
    "import pylab\n",
    "\n",
    "# Open mesh in ply format\n",
    "def open_ply(filepath):\n",
    "    f=open(filepath,'r');\n",
    "    i=0;\n",
    "    ip=0;\n",
    "    it=0;\n",
    "    np=0;\n",
    "    nt=0;\n",
    "    for str in f:\n",
    "        i+=1;\n",
    "        arr=str.split(\" \");\n",
    "        if(arr[0]==\"element\" and arr[1]==\"vertex\"):\n",
    "            np=int(arr[2]);\n",
    "            p=numpy.zeros((np,3));\n",
    "        elif(arr[0]==\"element\" and arr[1]==\"face\"):\n",
    "            nt=int(arr[2]);\n",
    "            t=numpy.zeros((nt,3));\n",
    "        elif(i>=11 and i<11+np):\n",
    "            p[ip,0]=float(arr[0]);\n",
    "            p[ip,1]=float(arr[1]);\n",
    "            p[ip,2]=float(arr[2]);\n",
    "            ip+=1;\n",
    "        elif(i>=11+np and i<11+np+nt):\n",
    "            t[it,0]=int(arr[1]);\n",
    "            t[it,1]=int(arr[2]);\n",
    "            t[it,2]=int(arr[3]);\n",
    "            it+=1;\n",
    "    mesh={};\n",
    "    mesh[\"np\"]=np;\n",
    "    mesh[\"nt\"]=nt;\n",
    "    mesh[\"p\"]=p;\n",
    "    mesh[\"t\"]=t;\n",
    "    return mesh;\n",
    "\n",
    "# Save mesh in ply format\n",
    "def save_ply(mesh,filepath):\n",
    "    p=mesh[\"p\"];\n",
    "    t=mesh[\"t\"];\n",
    "    np=mesh[\"np\"];\n",
    "    nt=mesh[\"nt\"];\n",
    "    f=open(filepath,'w');\n",
    "    f.write(\"%i 1\\n\"%np);\n",
    "    f.write(\"ply\\n\");\n",
    "    f.write(\"format ascii 1.0\\n\");\n",
    "    f.write(\"comment ipython ferret, R. Toro 2015\\n\");\n",
    "    f.write(\"element vertex %i\\n\"%np);\n",
    "    f.write(\"property float x\\n\");\n",
    "    f.write(\"property float y\\n\");\n",
    "    f.write(\"property float z\\n\");\n",
    "    f.write(\"element face %i\\n\"%nt);\n",
    "    f.write(\"property list uchar int vertex_indices\\n\");\n",
    "    f.write(\"end_header\\n\");\n",
    "    for i in range(0,np):\n",
    "        f.write(\"%f %f %f\\n\"%(p[i,0],p[i,1],p[i,2]));\n",
    "    for i in range(0,nt):\n",
    "        f.write(\"3 %i %i %i\\n\"%(t[i,0],t[i,1],t[i,2]));\n",
    "    f.close();\n",
    "\n",
    "# Save a scalar over vertices to a text file\n",
    "def saveScalar(s,filename):\n",
    "    np=len(s);\n",
    "    f=open(filename,\"w\");\n",
    "    f.write(\"%i 1\\n\"%np);\n",
    "    for i in range(0,np):\n",
    "        f.write(\"%f\\n\"%s[i]);\n",
    "    f.close();\n",
    "\n",
    "# Area of triangle p0, p1, p2\n",
    "def A(p0,p1,p2):\n",
    "    #a,b,c    // side lengths\n",
    "    #s        // semiperimeter\n",
    "    a=numpy.linalg.norm(p0-p1);\n",
    "    b=numpy.linalg.norm(p1-p2);\n",
    "    c=numpy.linalg.norm(p2-p0);\n",
    "    s=(a+b+c)/2.0;\n",
    "    return numpy.sqrt(s*(s-a)*(s-b)*(s-c))\n",
    "\n",
    "# Normal of triangle a,b,c\n",
    "def normal(a,b,c):\n",
    "    v=numpy.cross(b-a,c-a)\n",
    "    n=v/numpy.linalg.norm(v)\n",
    "    return n;\n",
    "\n",
    "# Compute the normals for each vertex in mesh\n",
    "def normals(mesh):\n",
    "    np=mesh[\"np\"];\n",
    "    nt=mesh[\"nt\"];\n",
    "    p=mesh[\"p\"];\n",
    "    t=mesh[\"t\"];\n",
    "    n=numpy.zeros((np,3));\n",
    "    for i in range(0,nt):\n",
    "        nn=normal(p[t[i,0]],p[t[i,1]],p[t[i,2]]);\n",
    "        n[t[i,0]]=numpy.add(n[t[i,0]],nn);\n",
    "        n[t[i,1]]=numpy.add(n[t[i,1]],nn);\n",
    "        n[t[i,2]]=numpy.add(n[t[i,2]],nn);\n",
    "    for i in range(0,np):\n",
    "        n[i]=n[i]/numpy.linalg.norm(n[i]);\n",
    "    return n;\n",
    "\n",
    "# Compute number of neighbours for each vertex in mesh\n",
    "def neighbours(mesh):\n",
    "    np=mesh[\"np\"];\n",
    "    nt=mesh[\"nt\"];\n",
    "    t=mesh[\"t\"];\n",
    "    neig=numpy.zeros(np);\n",
    "    for i in range(0,nt):\n",
    "        neig[t[i,0]]+=1;\n",
    "        neig[t[i,1]]+=1;\n",
    "        neig[t[i,2]]+=1;\n",
    "    return neig;\n",
    "\n",
    "# Compute log Jacobian of the deformation of each triangle between two meshes\n",
    "# with the same topology but different geometry. Distribute the value on\n",
    "# the vertices\n",
    "def logJacobian(m1,m2):\n",
    "    np=m1[\"np\"];\n",
    "    nt=m1[\"nt\"];\n",
    "    t=m1[\"t\"];\n",
    "    p1=m1[\"p\"];\n",
    "    p2=m2[\"p\"];\n",
    "    Jmap=numpy.zeros(np);\n",
    "    neig=neighbours(m1);\n",
    "    for i in range(0,nt):\n",
    "        # compute deformation at each triangle\n",
    "        A1=A(p1[t[i,0]],p1[t[i,1]],p1[t[i,2]]);\n",
    "        A2=A(p2[t[i,0]],p2[t[i,1]],p2[t[i,2]]);\n",
    "        if(A1<1e-5):\n",
    "            continue;\n",
    "        J=numpy.log(A2/A1);\n",
    "        # distribute the deformation among the vertices\n",
    "        Jmap[t[i,0]]+=J;\n",
    "        Jmap[t[i,1]]+=J;\n",
    "        Jmap[t[i,2]]+=J;\n",
    "    for i in range(0,np):\n",
    "        if(neig[i]>0):\n",
    "            Jmap[i]/=neig[i];\n",
    "    return Jmap;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 250,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Adapted from intersect_RayTriangle()\n",
    "# Copyright 2001, softSurfer (www.softsurfer.com)\n",
    "# This code may be freely used and modified for any purpose\n",
    "# providing that this copyright notice is included with it.\n",
    "# SoftSurfer makes no warranty for this code, and cannot be held\n",
    "# liable for any real or imagined damage resulting from its use.\n",
    "# Users of this code must verify correctness for their application.\n",
    "#    Input:  vector \"x\", triangle index \"it\"\n",
    "#    Output: *c0, *c1 = the triangle-based coordinates of the intersection (when it exists)\n",
    "#    Return: -1 = triangle is degenerate (a segment or point)\n",
    "#             0 = disjoint (no intersect)\n",
    "#             1 = intersect in unique point I1\n",
    "#             2 = are in the same plane\n",
    "# code from:http://geometryalgorithms.com/Archive/algorithm_0105/algorithm_0105.htm#intersect_RayTriangle()\n",
    "def intersect_VectorTriangle(x, i, c, m):\n",
    "    EPSILON=1e-6;\n",
    "    np=m[\"np\"];\n",
    "    nt=m[\"nt\"];\n",
    "    t=m[\"t\"];\n",
    "    p=m[\"p\"];\n",
    "    T=t[i];\n",
    "    xx=[3];\n",
    "    \n",
    "    u=numpy.subtract(p[T[1]],p[T[0]]);\n",
    "    v=numpy.subtract(p[T[2]],p[T[0]]);\n",
    "    n=numpy.cross(u,v);\n",
    "    if(numpy.sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2])<1e-10):\n",
    "        return -1;\n",
    "\n",
    "    direc=x;\n",
    "    w0=-p[T[0]];\n",
    "    a = numpy.dot(n,w0);\n",
    "    b = numpy.dot(n,direc);\n",
    "    if (b>-EPSILON and b<EPSILON):  # ray is parallel to triangle plane\n",
    "        if (a == 0.0):               # ray lies in triangle plane\n",
    "            return 2;\n",
    "        else:\n",
    "            return 0;              # ray disjoint from plane\n",
    "\n",
    "    # get intersect point of ray with triangle plane\n",
    "    r = -a/b;\n",
    "    # for a segment, also test if (r > 1.0) => no intersect\n",
    "    if (r < 0.0):                   # ray goes away from triangle\n",
    "        return 0;                   # => no intersect\n",
    "\n",
    "    xx=direc*r;    # intersect point of ray and plane\n",
    "\n",
    "    # is I inside T?\n",
    "    uu=numpy.dot(u,u);\n",
    "    uv=numpy.dot(u,v);\n",
    "    vv=numpy.dot(v,v);\n",
    "    w=numpy.subtract(xx,p[T[0]]);\n",
    "    wu=numpy.dot(w,u);\n",
    "    wv=numpy.dot(w,w);\n",
    "    D = uv * uv - uu * vv;\n",
    "\n",
    "    # get and test parametric coords\n",
    "    ss = (uv * wv - vv * wu) / D;\n",
    "    if(ss>-EPSILON and ss<EPSILON):\n",
    "        ss=0;\n",
    "    if((1-ss)>-EPSILON and (1-ss)<EPSILON):\n",
    "        ss=1;\n",
    "    \n",
    "    tt = (uv * wu - uu * wv) / D;\n",
    "    if(tt>-EPSILON and tt<EPSILON):\n",
    "        tt=0;\n",
    "    if((1-tt)>-EPSILON and (1-tt)<EPSILON):\n",
    "        tt=1;\n",
    "    \n",
    "    c[0]=ss;\n",
    "    c[1]=tt;\n",
    "    \n",
    "    if (ss < 0.0 or tt < 0.0 or (ss + tt) > 1.0):  # I is outside T\n",
    "        return 0;\n",
    "\n",
    "    return 1;                      # I is in T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 271,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def barycentric(p,mesh):\n",
    "    # compute barycentric coordinates relative to the topology mesh of the vertices in euclidean space\n",
    "    # (normally, the vertices in euclidean are in the geometry of mesh)\n",
    "    \n",
    "    np_xyz=len(xyz);\n",
    "    p=mesh[\"p\"];\n",
    "    t=mesh[\"t\"];\n",
    "    np=mesh[\"np\"];\n",
    "    nt=mesh[\"nt\"];\n",
    "    \n",
    "    # Interpolate coordinates on reference mesh\n",
    "    b=numpy.zeros((np_xyz,3));   # the new points in barycentric coordinates are stored in b\n",
    "    case_deg=0;\n",
    "    case_disj=0;\n",
    "    case_parl=0;\n",
    "    non_mapped=0;\n",
    "    for i in range(0,np_xyz):\n",
    "        # look for a triangle in the spherical original mesh containing\n",
    "        # point i of the target mesh\n",
    "        \n",
    "        # display progress\n",
    "        if(int(i*100/np_xyz)>int((i-1)*100/np_xyz)):\n",
    "            print(\"%i%% \" % (i*100/np_xyz));\n",
    "        \n",
    "        for j in range(0,nt):\n",
    "            c=[0,0];\n",
    "            result=intersect_VectorTriangle(xyz[i],j,c,mesh);\n",
    "            if(result==1):\n",
    "                b[i]=[j,c[0],c[1]];\n",
    "                #b[i]=(1-c0-c1)*p[t[j,0]];\n",
    "                #b[i]=numpy.add(b[i],c[0]*p[t[j,1]]);\n",
    "                #b[i]=numpy.add(b[i],c[1]*p[t[j,2]]);\n",
    "                break;\n",
    "            if(result==-1):\n",
    "                case_deg+=1;\n",
    "            if(result==0):\n",
    "                case_disj+=1;\n",
    "            if(result==2):\n",
    "                case_parl+=1;\n",
    "        \n",
    "        # if j==nt no triangle was found in the spherical original mesh that contained\n",
    "        # point i. As an approximation, pick the closest point in the spherical original\n",
    "        # mesh.\n",
    "        if(j==nt):\n",
    "            mindist=numpy.linalg.norm(sub3D(xyz[i],p[0]));\n",
    "            for k in range(0,np):\n",
    "                if(numpy.linalg.norm(numpy.subtract(xyz[i],p[k]))<mindist):\n",
    "                    mindist=numpy.linalg.norm(numpy.subtract(xyz[i],p[k]));\n",
    "                    imindist=k;\n",
    "            c=[0,0];\n",
    "            for k in range(0,nt):\n",
    "                if(t[k,0]==imindist):\n",
    "                    break;\n",
    "                if(t[k,1]==imindist):\n",
    "                    c=[1,0];\n",
    "                    break;\n",
    "                if(t[k,2]==imindist):\n",
    "                    c=[0,1];\n",
    "                    break;\n",
    "            b[i]=[imindist,c[0],c[1]];\n",
    "            non_mapped+=1;\n",
    "    if(non_mapped):\n",
    "        print \"WARNING: %i vertices could not be resampled and were mapped to the closest vertex\\n\" % non_mapped;\n",
    "    return b;\n",
    "\n",
    "def euclidean(barycentric,mesh):\n",
    "    # compute vertices in euclidean coordinates from the barycentric coordinates in barycentric and their\n",
    "    # reference mesh\n",
    "    return;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0% \n",
      "1% \n",
      "2% \n",
      "3% \n",
      "4% \n",
      "5% \n",
      "6% \n",
      "7% \n",
      "8% \n",
      "9% \n",
      "10% \n",
      "11% \n",
      "12% \n",
      "13% \n",
      "14% \n",
      "15% \n",
      "16% \n",
      "17% \n",
      "18% \n",
      "19% \n",
      "20% \n",
      "21% \n",
      "22% \n",
      "23% \n",
      "24% \n",
      "25% \n",
      "26% \n",
      "27% \n",
      "28% \n",
      "29% \n",
      "30% \n",
      "31% \n",
      "32% \n",
      "33% \n",
      "34% \n",
      "35% \n",
      "36% \n",
      "37% \n",
      "38% \n",
      "39% \n",
      "40% \n",
      "41% \n",
      "42% \n",
      "43% \n",
      "44% \n",
      "45% \n",
      "46% \n",
      "47% \n",
      "48% \n",
      "49% \n",
      "50% \n",
      "51% \n",
      "52% \n",
      "53% \n",
      "54% \n",
      "55% \n",
      "56% \n",
      "57% \n",
      "58% \n",
      "59% \n",
      "60% \n",
      "61% \n",
      "62% \n",
      "63% \n",
      "64% \n",
      "65% \n",
      "66% \n",
      "67% \n",
      "68% \n",
      "69% \n",
      "70% \n",
      "71% \n",
      "72% \n",
      "73% \n",
      "74% \n",
      "75% \n",
      "76% \n",
      "77% \n",
      "78% \n",
      "79% \n",
      "80% \n",
      "81% \n",
      "82% \n",
      "83% \n",
      "84% \n",
      "85% \n",
      "86% \n",
      "87% \n",
      "88% \n",
      "89% \n",
      "90% \n",
      "91% \n",
      "92% \n",
      "93% \n",
      "94% \n",
      "95% \n",
      "96% \n",
      "97% \n",
      "98% \n",
      "99% \n"
     ]
    }
   ],
   "source": [
    "m1=open_ply(\"/Users/roberto/Desktop/ferret/ferret-sbn/P8-F10/_morph/P8-F10_to_P8-F10.sphere.ply\");\n",
    "m2=open_ply(\"/Users/roberto/Desktop/ferret/ferret-sbn/P8-F10/_morph/P8-F10_to_P8-F11.sphere.ply\");\n",
    "xyz=m1[\"p\"];\n",
    "b=barycentric(xyz,m2);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorithm seems to be working but it's extremely slow. I'll use C code instead.\n",
    "In meshgeometry I could use the function 'resample'. Resample takes 3 inputs:\n",
    "1. a native mesh m\n",
    "2. a spherical version of the native mesh m, sm\n",
    "3. a spherical version of a target (reference) mesh rm.\n",
    "\n",
    "Surface landmarks over the two spherical meshes are supposed to be in correspondance. The algorithm then resamples the native mesh m to have the topology of the spherical mesh rm: It finds for every vertex in rm the corresponding vertex in sm. In general, the vertex from rm will not map directly to a vertex in sm, but will be inside one of its triangles. The algorithm computes the barycentric coordinates of rm's vertices relative to the triangle in sm that contains it, and reproduces the vertex in the native mesh m by applying the barycentric coordinates to the corresponding native m triangle.\n",
    "\n",
    "In our case, rm should contain the vertices of the source mesh mapped into a sphere, and m and sm should be the target sphere. For example, rm could be Adult-F03_to_Adult-F03.sphere.py, m = P32-F16_to_P32-F16.ply, and sm=P32-F16_to_P32-F16.sphere.py.\n",
    "\n",
    "The code for a complete transformation would then be:\n",
    "\n",
    "mg -i P32-F16.ply -resample P32-F16.sphere.ply Adult-F03.sphere.ply -o P32-F16_as_Adult-F03.ply"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Decompose normal and tangential deformation\n",
    "\n",
    "# load meshes\n",
    "\n",
    "#f1=\"P0-F02\";\n",
    "#f1=\"P2-F25\";\n",
    "#f1=\"P4-F06\";\n",
    "f1=\"P8-F10\"\n",
    "\n",
    "dir=\"/Users/roberto/Desktop/ferret/ferret-sbn/\";\n",
    "m1=open_ply(dir+f1+\"/_morph/\"+f1+\"_to_\"+f1+\".ply\");\n",
    "nt=m1[\"nt\"];\n",
    "np=m1[\"np\"];\n",
    "t=m1[\"t\"];\n",
    "p1=m1[\"p\"];\n",
    "\n",
    "# average deformation field, average jacobian\n",
    "avrD=numpy.zeros((np,3));\n",
    "avrJ=numpy.zeros(np);\n",
    "\n",
    "# compute neighbours\n",
    "neig=neighbours(m1);\n",
    "\n",
    "#for f2 in (\"P0-F23\",\"P2-F25\",\"P2-F26\",\"P2-F27\"):\n",
    "#for f2 in (\"P2-F26\",\"P2-F27\",\"P4-F06\",\"P4-F07\",\"P4-F08\",\"P4-F09\"):\n",
    "#for f2 in (\"P4-F07\",\"P4-F08\",\"P4-F09\",\"P8-F10\",\"P8-F11\",\"P8-F13\",\"P8-F14\"):\n",
    "for f2 in (\"P8-F11\",\"P8-F13\",\"P8-F14\",\"P16-F12\",\"P16-F15\",\"P16-F21\"):\n",
    "    m2=open_ply(dir+f2+\"/_morph/\"+f2+\"_to_\"+f1+\".ply\");\n",
    "    p2=m2[\"p\"];\n",
    "    \n",
    "    # compute deformation vector field\n",
    "    d=numpy.zeros((np,3));\n",
    "    for i in range(0,np):\n",
    "        d[i]=numpy.subtract(p2[i,],p1[i,]);\n",
    "\n",
    "    # decompose deformation into normal and tangential\n",
    "    n=normals(m1); # compute normal direction\n",
    "    nd=numpy.zeros((np,3));\n",
    "    td=numpy.zeros((np,3));\n",
    "    for i in range(0,np):\n",
    "        nd[i]=n[i]*numpy.dot(d[i],n[i]); # normal deformation\n",
    "        td[i]=numpy.subtract(d[i],nd[i]);   # tangential deformation\n",
    "\n",
    "    # smooth the tangential part of the deformation\n",
    "    td0=numpy.zeros((np,3));\n",
    "    g=0.99;\n",
    "    niter=5;\n",
    "    for j in range(0,niter):\n",
    "        for i in range(0,nt):\n",
    "            td0[t[i,0]]=numpy.add(td0[t[i,0]],numpy.add(td[t[i,1]],td[t[i,2]]));\n",
    "            td0[t[i,1]]=numpy.add(td0[t[i,1]],numpy.add(td[t[i,2]],td[t[i,0]]));\n",
    "            td0[t[i,2]]=numpy.add(td0[t[i,2]],numpy.add(td[t[i,0]],td[t[i,1]]));\n",
    "        for i in range(0,np):\n",
    "            td[i]=numpy.add(g*td[i],(1-g)*td0[i]/(2.0*neig[i]));\n",
    "\n",
    "    # recompose the deformation with a weighted tangential part, add to average deformation field\n",
    "    w=0.5;\n",
    "    for i in range(0,np):\n",
    "        d[i]=numpy.add(nd[i],td[i]*w);\n",
    "    avrD=numpy.add(avrD,d);\n",
    "    \n",
    "    # replace mesh m2 with m1 + new deformation\n",
    "    for i in range(0,np):\n",
    "        p2[i]=numpy.add(p1[i],d[i]);\n",
    "\n",
    "    # compute jacobian of area deformation, add to average jacobian\n",
    "    J=logJacobian(m1,m2);\n",
    "    avrJ=numpy.add(avrJ,J);\n",
    "    \n",
    "avrJ=(1/7.0)*avrJ;\n",
    "avrD=(1/7.0)*avrD;\n",
    "\n",
    "saveScalar(avrJ,\"/Users/roberto/Desktop/avrJ.txt\");\n",
    "for i in range(0,np):\n",
    "    p2[i]=numpy.add(p1[i],avrD[i]);\n",
    "save_ply(m2,\"/Users/roberto/Desktop/average_P8-F10base.ply\");\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dir=\"/Users/roberto/Desktop/ferret/ferret-sbn/\";\n",
    "\n",
    "t=(0,0, 1,1,1, 2,2,2,2, 3,3,3,3, 4,4,4, 5,5,5,5, 6,6,6);\n",
    "n=((\"P0-F02\",\"P0-F23\"),\n",
    "   (\"P2-F25\",\"P2-F26\",\"P2-F27\"),\n",
    "   (\"P4-F06\",\"P4-F07\",\"P4-F08\",\"P4-F09\"),\n",
    "   (\"P8-F10\",\"P8-F11\",\"P8-F13\",\"P8-F14\"),\n",
    "   (\"P16-F12\",\"P16-F15\",\"P16-F21\"),\n",
    "   (\"P32-F16\",\"P32-F17\",\"P32-F18\",\"P32-F19\"),\n",
    "   (\"Adult-F01\",\"Adult-F03\",\"Adult-F04\"));\n",
    "m=[];\n",
    "for i in range(0,len(n)):\n",
    "    m.append([]);\n",
    "    for j in range(0,len(n[i])):\n",
    "        filename=dir+n[i][j]+\"/_morph/\"+n[i][j]+\"_to_\"+n[i][0]+\".ply\";\n",
    "        m[i].append(open_ply(filename));\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each brain was transformed into each other based on a series of common landmarks.\n",
    "In each transformation a_to_b the geometry is that of a and the topology that of b. Then, a_to_a is brain a with its geometry and its native topology. The geometric transformation between a and b can be achieved either over the topology of a or the topology of b. In general, it may be better to use the topology with the denser sampling of the geometry.\n",
    "\n",
    "In the following problem, I would like to map all vertices of Adult-F03 into the geometry of P02-F02. For this, I will first map the points into the geometry of one P32 brain, then into that of a P16 brain, P8, P4, P2 and finally P02-F02. The choice of the intermediate brains is arbitrary. I will simply choose the first from the list, but later I would like to chose them all, to produce different perspectives on the same projection. The sequence of transformations will be then:\n",
    "\n",
    "Adult-F03 : P32-F16 : P16-F12 : P8-F10 : P4-F06 : P2-F25 : P0-F02.\n",
    "\n",
    "For each transformation, the steps are as follows:\n",
    "\n",
    "1. The original vertices are those of Adult-F03_to_Adult-F03,\n",
    "2. P32-F16_to_Adult-F03 provides the vertices of Adult-F03 projected onto P32-F16,\n",
    "3. I obtain triangle coordinates of the vertices in P32-F16_to_Adult-F03 in the topology P32-F16. Each vertex of coordinates x,y,z is mapped into coordinates i,s,t, where i is the index of a triangle in the topology P32-F16 and s,t are barycentric coordinates,\n",
    "4. I obtain triangle coordinates of the vertices of barycentric coordinates i,s,t in the topology P32-F16 and transform them into barycentric coordinates i,s,t in P16-F12 by using P32-F16_to_P16-F12.\n",
    "\n",
    "For each a_to_b mesh, there is also a a_to_b.sphere version. In the .sphere version the topology is that of a and the geometry that of b. Not really the geometry (because it's a sphere), but the position of the vertices in the sphere, which is such that they move the landmarks from a to match those of b.\n",
    "\n",
    "The computation of barycentri coordinates should be easier using the spherical meshes. The steps in that case would be:\n",
    "\n",
    "1. have vertices Adult-F03.sphere: move them to native space (trivial)\n",
    "2. Use Adult-F03_to_Adult-F01.sphere to move the vertices of Adult-F03 to Adult-F01\n",
    "3. have vertices Adult-F01.sphere: move them to native space,\n",
    "4. Get barycentric coordinates of these vertices relative to the original Adult-F01 sphere topology, that in Adult-F01_to_Adult-F01.sphere,\n",
    "5. Obtain vertices in P32-F16 by using Adult-F01_to_P32-F16.sphere and the previous barycentric coordinates\n",
    "6. have vertices P32-F16.sphere: move them to native space,\n",
    "7. Compute barycentric coordinates relative to the original P32-F16 sphere topology\n",
    "8. Obtain vertices in P16-F12 by using P32-F16_to_P16-F12 and the previous barycentric coordinates\n",
    "9. have vertices P16-F12.sphere: move them to native space,\n",
    "etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[{'np': 2258, 'nt': 4512, 'p': array([[ 1.049709, -1.543577, -3.904279],\n",
       "          [ 1.431044, -0.990687, -3.92231 ],\n",
       "          [ 0.369462, -1.43425 , -3.795803],\n",
       "          ..., \n",
       "          [-3.181885,  0.868488,  3.433966],\n",
       "          [-1.825432, -1.406168,  3.563831],\n",
       "          [-1.335176, -1.507696,  3.524995]]), 't': array([[  3.00000000e+00,   7.57000000e+02,   7.54000000e+02],\n",
       "          [  1.00000000e+00,   7.60000000e+02,   7.55000000e+02],\n",
       "          [  1.00000000e+00,   7.55000000e+02,   7.56000000e+02],\n",
       "          ..., \n",
       "          [  7.48000000e+02,   2.24200000e+03,   2.25500000e+03],\n",
       "          [  7.49000000e+02,   2.25700000e+03,   2.25600000e+03],\n",
       "          [  7.41000000e+02,   2.25000000e+03,   2.25700000e+03]])},\n",
       "  {'np': 2258,\n",
       "   'nt': 4512,\n",
       "   'p': array([[  1.65794200e+00,  -1.41670500e+00,  -3.91746500e+00],\n",
       "          [  2.14248500e+00,  -9.51966000e-01,  -3.88367300e+00],\n",
       "          [  9.37912000e-01,  -1.29562800e+00,  -3.91286700e+00],\n",
       "          ..., \n",
       "          [ -3.88156900e+00,  -3.20200000e-03,   2.51854800e+00],\n",
       "          [ -1.78606800e+00,  -6.81974000e-01,   3.14065700e+00],\n",
       "          [ -1.42045200e+00,  -8.05918000e-01,   3.16781500e+00]]),\n",
       "   't': array([[  3.00000000e+00,   7.57000000e+02,   7.54000000e+02],\n",
       "          [  1.00000000e+00,   7.60000000e+02,   7.55000000e+02],\n",
       "          [  1.00000000e+00,   7.55000000e+02,   7.56000000e+02],\n",
       "          ..., \n",
       "          [  7.48000000e+02,   2.24200000e+03,   2.25500000e+03],\n",
       "          [  7.49000000e+02,   2.25700000e+03,   2.25600000e+03],\n",
       "          [  7.41000000e+02,   2.25000000e+03,   2.25700000e+03]])}],\n",
       " [{'np': 2240, 'nt': 4476, 'p': array([[ 1.01047 , -0.989234, -3.580185],\n",
       "          [ 1.748947,  0.896132, -3.718374],\n",
       "          [ 1.813444,  1.435999, -3.682173],\n",
       "          ..., \n",
       "          [ 3.756866,  1.16471 , -1.616886],\n",
       "          [ 3.504362,  1.050558, -1.71101 ],\n",
       "          [ 3.477273,  0.836235, -1.812509]]), 't': array([[  8.00000000e+00,   7.49000000e+02,   7.48000000e+02],\n",
       "          [  2.00000000e+00,   7.59000000e+02,   7.49000000e+02],\n",
       "          [  3.00000000e+00,   7.70000000e+02,   7.50000000e+02],\n",
       "          ..., \n",
       "          [  2.91000000e+02,   1.15800000e+03,   2.23700000e+03],\n",
       "          [  2.52000000e+02,   2.23700000e+03,   2.23800000e+03],\n",
       "          [  7.44000000e+02,   2.23800000e+03,   2.23900000e+03]])},\n",
       "  {'np': 2240, 'nt': 4476, 'p': array([[ 0.343504, -0.557171, -3.460056],\n",
       "          [ 1.523337,  1.007124, -3.691451],\n",
       "          [ 1.849639,  1.593041, -3.669582],\n",
       "          ..., \n",
       "          [ 3.765165,  1.239948, -1.823711],\n",
       "          [ 3.450772,  1.058145, -1.799291],\n",
       "          [ 3.340108,  0.836795, -1.828906]]), 't': array([[  8.00000000e+00,   7.49000000e+02,   7.48000000e+02],\n",
       "          [  2.00000000e+00,   7.59000000e+02,   7.49000000e+02],\n",
       "          [  3.00000000e+00,   7.70000000e+02,   7.50000000e+02],\n",
       "          ..., \n",
       "          [  2.91000000e+02,   1.15800000e+03,   2.23700000e+03],\n",
       "          [  2.52000000e+02,   2.23700000e+03,   2.23800000e+03],\n",
       "          [  7.44000000e+02,   2.23800000e+03,   2.23900000e+03]])},\n",
       "  {'np': 2240, 'nt': 4476, 'p': array([[ 1.106579, -0.939492, -3.374319],\n",
       "          [ 2.095376,  1.471113, -3.462966],\n",
       "          [ 2.24635 ,  2.014304, -3.392013],\n",
       "          ..., \n",
       "          [ 3.759423,  1.1349  , -1.129631],\n",
       "          [ 3.621883,  1.084423, -1.378771],\n",
       "          [ 3.581671,  0.895369, -1.438303]]), 't': array([[  8.00000000e+00,   7.49000000e+02,   7.48000000e+02],\n",
       "          [  2.00000000e+00,   7.59000000e+02,   7.49000000e+02],\n",
       "          [  3.00000000e+00,   7.70000000e+02,   7.50000000e+02],\n",
       "          ..., \n",
       "          [  2.91000000e+02,   1.15800000e+03,   2.23700000e+03],\n",
       "          [  2.52000000e+02,   2.23700000e+03,   2.23800000e+03],\n",
       "          [  7.44000000e+02,   2.23800000e+03,   2.23900000e+03]])}],\n",
       " [{'np': 2972, 'nt': 5940, 'p': array([[ 2.490973, -1.572496, -4.322743],\n",
       "          [ 2.017202, -1.196542, -4.330479],\n",
       "          [ 1.789447,  1.219857, -4.34371 ],\n",
       "          ..., \n",
       "          [ 3.010911,  1.444337, -1.261308],\n",
       "          [ 2.918555,  1.317347, -1.322652],\n",
       "          [ 2.908904,  1.348094, -1.243366]]), 't': array([[  0.00000000e+00,   9.95000000e+02,   9.92000000e+02],\n",
       "          [  1.00000000e+00,   9.92000000e+02,   9.93000000e+02],\n",
       "          [  6.00000000e+00,   1.00400000e+03,   9.94000000e+02],\n",
       "          ..., \n",
       "          [  3.92000000e+02,   2.97000000e+03,   2.96900000e+03],\n",
       "          [  3.92000000e+02,   1.77600000e+03,   2.97000000e+03],\n",
       "          [  4.19000000e+02,   1.77800000e+03,   2.97100000e+03]])},\n",
       "  {'np': 2972, 'nt': 5940, 'p': array([[ 3.060615, -1.376928, -3.663542],\n",
       "          [ 2.465373, -1.060402, -3.745088],\n",
       "          [ 2.451793,  1.18512 , -3.631171],\n",
       "          ..., \n",
       "          [ 3.326127,  1.105506, -1.166584],\n",
       "          [ 3.26657 ,  1.06536 , -1.153912],\n",
       "          [ 3.310065,  1.044096, -1.183877]]), 't': array([[  0.00000000e+00,   9.95000000e+02,   9.92000000e+02],\n",
       "          [  1.00000000e+00,   9.92000000e+02,   9.93000000e+02],\n",
       "          [  6.00000000e+00,   1.00400000e+03,   9.94000000e+02],\n",
       "          ..., \n",
       "          [  3.92000000e+02,   2.97000000e+03,   2.96900000e+03],\n",
       "          [  3.92000000e+02,   1.77600000e+03,   2.97000000e+03],\n",
       "          [  4.19000000e+02,   1.77800000e+03,   2.97100000e+03]])},\n",
       "  {'np': 2972, 'nt': 5940, 'p': array([[ 3.028245, -1.998072, -3.971659],\n",
       "          [ 2.592375, -1.639717, -4.066964],\n",
       "          [ 2.161827,  1.259957, -4.100427],\n",
       "          ..., \n",
       "          [ 3.353641,  1.165337, -1.376097],\n",
       "          [ 3.269119,  1.168262, -1.360866],\n",
       "          [ 3.315286,  1.114618, -1.381168]]), 't': array([[  0.00000000e+00,   9.95000000e+02,   9.92000000e+02],\n",
       "          [  1.00000000e+00,   9.92000000e+02,   9.93000000e+02],\n",
       "          [  6.00000000e+00,   1.00400000e+03,   9.94000000e+02],\n",
       "          ..., \n",
       "          [  3.92000000e+02,   2.97000000e+03,   2.96900000e+03],\n",
       "          [  3.92000000e+02,   1.77600000e+03,   2.97000000e+03],\n",
       "          [  4.19000000e+02,   1.77800000e+03,   2.97100000e+03]])},\n",
       "  {'np': 2972, 'nt': 5940, 'p': array([[ 2.539167, -1.550919, -4.069567],\n",
       "          [ 2.05482 , -1.211228, -4.089096],\n",
       "          [ 1.698653,  1.51475 , -4.09295 ],\n",
       "          ..., \n",
       "          [ 2.961355,  1.110287, -1.304364],\n",
       "          [ 2.901679,  1.067329, -1.349773],\n",
       "          [ 2.936682,  1.04634 , -1.306279]]), 't': array([[  0.00000000e+00,   9.95000000e+02,   9.92000000e+02],\n",
       "          [  1.00000000e+00,   9.92000000e+02,   9.93000000e+02],\n",
       "          [  6.00000000e+00,   1.00400000e+03,   9.94000000e+02],\n",
       "          ..., \n",
       "          [  3.92000000e+02,   2.97000000e+03,   2.96900000e+03],\n",
       "          [  3.92000000e+02,   1.77600000e+03,   2.97000000e+03],\n",
       "          [  4.19000000e+02,   1.77800000e+03,   2.97100000e+03]])}],\n",
       " [{'np': 5141, 'nt': 10278, 'p': array([[ 2.653107, -1.697588, -5.474947],\n",
       "          [ 2.286556,  1.648852, -5.414738],\n",
       "          [ 3.387123,  1.781849, -5.4173  ],\n",
       "          ..., \n",
       "          [ 4.022955, -2.384107, -0.386038],\n",
       "          [ 3.605758, -3.783285,  0.623482],\n",
       "          [ 3.671371, -3.498383,  0.250342]]), 't': array([[    7.,  1727.,  1715.],\n",
       "          [    8.,  1717.,  1716.],\n",
       "          [    8.,  1718.,  1717.],\n",
       "          ..., \n",
       "          [ 1707.,  3298.,  5138.],\n",
       "          [  928.,  3482.,  5139.],\n",
       "          [  915.,  3482.,  5140.]])},\n",
       "  {'np': 5141, 'nt': 10278, 'p': array([[ 2.002544, -1.581959, -5.116515],\n",
       "          [ 1.937191,  1.53782 , -5.244763],\n",
       "          [ 2.991627,  1.499699, -5.206117],\n",
       "          ..., \n",
       "          [ 3.6104  , -1.89664 , -0.968377],\n",
       "          [ 3.308476, -2.373057, -0.790042],\n",
       "          [ 3.285579, -2.243369, -0.850298]]), 't': array([[    7.,  1727.,  1715.],\n",
       "          [    8.,  1717.,  1716.],\n",
       "          [    8.,  1718.,  1717.],\n",
       "          ..., \n",
       "          [ 1707.,  3298.,  5138.],\n",
       "          [  928.,  3482.,  5139.],\n",
       "          [  915.,  3482.,  5140.]])},\n",
       "  {'np': 5141, 'nt': 10278, 'p': array([[ 1.712895, -1.318834, -4.82521 ],\n",
       "          [ 1.826236,  1.166997, -4.786522],\n",
       "          [ 2.788039,  1.197504, -4.750321],\n",
       "          ..., \n",
       "          [ 3.642941, -2.082922, -0.611755],\n",
       "          [ 3.510339, -3.044082, -0.432727],\n",
       "          [ 3.479362, -2.846353, -0.465475]]), 't': array([[    7.,  1727.,  1715.],\n",
       "          [    8.,  1717.,  1716.],\n",
       "          [    8.,  1718.,  1717.],\n",
       "          ..., \n",
       "          [ 1707.,  3298.,  5138.],\n",
       "          [  928.,  3482.,  5139.],\n",
       "          [  915.,  3482.,  5140.]])},\n",
       "  {'np': 5141, 'nt': 10278, 'p': array([[ 1.331485, -1.512003, -4.703401],\n",
       "          [ 2.064113,  1.03932 , -4.713346],\n",
       "          [ 3.00273 ,  1.058487, -4.572922],\n",
       "          ..., \n",
       "          [ 4.054825, -2.180355, -0.418839],\n",
       "          [ 4.288304, -2.916071, -0.068526],\n",
       "          [ 4.172038, -2.764451, -0.121981]]), 't': array([[    7.,  1727.,  1715.],\n",
       "          [    8.,  1717.,  1716.],\n",
       "          [    8.,  1718.,  1717.],\n",
       "          ..., \n",
       "          [ 1707.,  3298.,  5138.],\n",
       "          [  928.,  3482.,  5139.],\n",
       "          [  915.,  3482.,  5140.]])}],\n",
       " [{'np': 8480, 'nt': 16956, 'p': array([[ 3.614849,  1.153943, -6.286922],\n",
       "          [ 2.930058,  1.602057, -6.396674],\n",
       "          [ 2.852801, -1.859572, -6.268044],\n",
       "          ..., \n",
       "          [ 2.965268, -0.676884, -0.330299],\n",
       "          [ 4.457123, -2.958657, -0.071977],\n",
       "          [ 4.516288, -2.662227, -0.192053]]), 't': array([[  1.00000000e+00,   2.83700000e+03,   2.82800000e+03],\n",
       "          [  1.00000000e+01,   2.83900000e+03,   2.82900000e+03],\n",
       "          [  4.00000000e+00,   2.83200000e+03,   2.83000000e+03],\n",
       "          ..., \n",
       "          [  1.39300000e+03,   5.53000000e+03,   8.47700000e+03],\n",
       "          [  1.41200000e+03,   5.80100000e+03,   8.47800000e+03],\n",
       "          [  1.48700000e+03,   8.47800000e+03,   8.47900000e+03]])},\n",
       "  {'np': 8480, 'nt': 16956, 'p': array([[ 3.19886 ,  0.926241, -6.047628],\n",
       "          [ 2.732341,  1.264967, -6.166062],\n",
       "          [ 2.904816, -2.15678 , -6.117323],\n",
       "          ..., \n",
       "          [ 3.367139, -0.377611, -0.472506],\n",
       "          [ 4.228868, -2.71571 , -0.262187],\n",
       "          [ 4.205086, -2.410999, -0.324713]]), 't': array([[  1.00000000e+00,   2.83700000e+03,   2.82800000e+03],\n",
       "          [  1.00000000e+01,   2.83900000e+03,   2.82900000e+03],\n",
       "          [  4.00000000e+00,   2.83200000e+03,   2.83000000e+03],\n",
       "          ..., \n",
       "          [  1.39300000e+03,   5.53000000e+03,   8.47700000e+03],\n",
       "          [  1.41200000e+03,   5.80100000e+03,   8.47800000e+03],\n",
       "          [  1.48700000e+03,   8.47800000e+03,   8.47900000e+03]])},\n",
       "  {'np': 8480,\n",
       "   'nt': 16956,\n",
       "   'p': array([[  4.58110100e+00,   1.12685300e+00,  -5.56584500e+00],\n",
       "          [  3.90411900e+00,   1.61516700e+00,  -5.69716500e+00],\n",
       "          [  4.78498000e+00,  -3.24227000e+00,  -5.38397100e+00],\n",
       "          ..., \n",
       "          [  3.33482000e+00,  -2.94402000e-01,  -1.89613000e-01],\n",
       "          [  3.97321100e+00,  -3.07547700e+00,   8.60860000e-02],\n",
       "          [  3.97288400e+00,  -2.65926800e+00,  -4.30300000e-03]]),\n",
       "   't': array([[  1.00000000e+00,   2.83700000e+03,   2.82800000e+03],\n",
       "          [  1.00000000e+01,   2.83900000e+03,   2.82900000e+03],\n",
       "          [  4.00000000e+00,   2.83200000e+03,   2.83000000e+03],\n",
       "          ..., \n",
       "          [  1.39300000e+03,   5.53000000e+03,   8.47700000e+03],\n",
       "          [  1.41200000e+03,   5.80100000e+03,   8.47800000e+03],\n",
       "          [  1.48700000e+03,   8.47800000e+03,   8.47900000e+03]])}],\n",
       " [{'np': 14855, 'nt': 29706, 'p': array([[ 4.678298, -2.011525, -7.405587],\n",
       "          [ 3.963161, -2.261656, -7.396012],\n",
       "          [ 3.644455, -1.577236, -7.339177],\n",
       "          ..., \n",
       "          [-1.459111, -7.14569 , -3.335274],\n",
       "          [-7.550813,  3.31196 ,  4.228912],\n",
       "          [-7.892125,  3.453977,  4.185431]]), 't': array([[  0.00000000e+00,   4.95900000e+03,   4.95300000e+03],\n",
       "          [  2.00000000e+00,   4.95500000e+03,   4.95400000e+03],\n",
       "          [  3.00000000e+00,   4.95300000e+03,   4.95500000e+03],\n",
       "          ..., \n",
       "          [  1.20400000e+03,   7.34100000e+03,   1.48520000e+04],\n",
       "          [  4.09700000e+03,   1.48540000e+04,   1.48530000e+04],\n",
       "          [  4.09700000e+03,   1.28690000e+04,   1.48540000e+04]])},\n",
       "  {'np': 14855, 'nt': 29706, 'p': array([[ 5.330573, -2.726277, -7.604253],\n",
       "          [ 4.345645, -2.989469, -7.610884],\n",
       "          [ 3.804988, -2.222224, -7.684813],\n",
       "          ..., \n",
       "          [-3.418753, -8.287624, -2.810532],\n",
       "          [-7.090554,  3.749119,  4.081769],\n",
       "          [-7.451814,  3.987896,  4.219059]]), 't': array([[  0.00000000e+00,   4.95900000e+03,   4.95300000e+03],\n",
       "          [  2.00000000e+00,   4.95500000e+03,   4.95400000e+03],\n",
       "          [  3.00000000e+00,   4.95300000e+03,   4.95500000e+03],\n",
       "          ..., \n",
       "          [  1.20400000e+03,   7.34100000e+03,   1.48520000e+04],\n",
       "          [  4.09700000e+03,   1.48540000e+04,   1.48530000e+04],\n",
       "          [  4.09700000e+03,   1.28690000e+04,   1.48540000e+04]])},\n",
       "  {'np': 14855, 'nt': 29706, 'p': array([[ 4.927196, -2.029512, -7.243867],\n",
       "          [ 4.005247, -2.305599, -7.298337],\n",
       "          [ 3.44169 , -1.536723, -7.281633],\n",
       "          ..., \n",
       "          [-2.779789, -7.339947, -3.091337],\n",
       "          [-6.819797,  3.571472,  4.270683],\n",
       "          [-7.146214,  3.839266,  4.288119]]), 't': array([[  0.00000000e+00,   4.95900000e+03,   4.95300000e+03],\n",
       "          [  2.00000000e+00,   4.95500000e+03,   4.95400000e+03],\n",
       "          [  3.00000000e+00,   4.95300000e+03,   4.95500000e+03],\n",
       "          ..., \n",
       "          [  1.20400000e+03,   7.34100000e+03,   1.48520000e+04],\n",
       "          [  4.09700000e+03,   1.48540000e+04,   1.48530000e+04],\n",
       "          [  4.09700000e+03,   1.28690000e+04,   1.48540000e+04]])},\n",
       "  {'np': 14855, 'nt': 29706, 'p': array([[ 4.824162, -2.123909, -8.085468],\n",
       "          [ 3.903848, -2.455881, -8.011341],\n",
       "          [ 3.475203, -1.62065 , -8.022239],\n",
       "          ..., \n",
       "          [-2.916347, -7.349038, -3.552241],\n",
       "          [-6.768689,  3.384064,  4.881851],\n",
       "          [-7.089322,  3.565151,  4.924067]]), 't': array([[  0.00000000e+00,   4.95900000e+03,   4.95300000e+03],\n",
       "          [  2.00000000e+00,   4.95500000e+03,   4.95400000e+03],\n",
       "          [  3.00000000e+00,   4.95300000e+03,   4.95500000e+03],\n",
       "          ..., \n",
       "          [  1.20400000e+03,   7.34100000e+03,   1.48520000e+04],\n",
       "          [  4.09700000e+03,   1.48540000e+04,   1.48530000e+04],\n",
       "          [  4.09700000e+03,   1.28690000e+04,   1.48540000e+04]])}],\n",
       " [{'np': 22154, 'nt': 44304, 'p': array([[-5.921483,  0.039153, -7.031293],\n",
       "          [-6.74802 ,  0.212333, -7.074325],\n",
       "          [ 3.693253,  0.327191, -6.835725],\n",
       "          ..., \n",
       "          [ 1.759177,  3.248919, -3.150549],\n",
       "          [ 1.502927,  3.528251, -3.199427],\n",
       "          [ 1.311012,  3.511701, -2.855483]]), 't': array([[  1.90000000e+01,   7.45700000e+03,   7.38600000e+03],\n",
       "          [  3.00000000e+00,   7.46900000e+03,   7.38700000e+03],\n",
       "          [  3.00000000e+00,   7.38700000e+03,   7.38800000e+03],\n",
       "          ..., \n",
       "          [  7.38500000e+03,   2.21520000e+04,   2.21510000e+04],\n",
       "          [  7.38500000e+03,   2.21530000e+04,   2.21520000e+04],\n",
       "          [  2.57200000e+03,   1.24910000e+04,   2.21530000e+04]])},\n",
       "  {'np': 22154, 'nt': 44304, 'p': array([[-6.370638,  0.602401, -6.729753],\n",
       "          [-7.142716,  0.605243, -6.727528],\n",
       "          [ 3.7775  ,  0.148195, -6.441604],\n",
       "          ..., \n",
       "          [ 2.28715 ,  3.941809, -4.8693  ],\n",
       "          [ 2.198723,  4.39492 , -4.894017],\n",
       "          [ 1.860036,  4.316628, -4.782752]]), 't': array([[  1.90000000e+01,   7.45700000e+03,   7.38600000e+03],\n",
       "          [  3.00000000e+00,   7.46900000e+03,   7.38700000e+03],\n",
       "          [  3.00000000e+00,   7.38700000e+03,   7.38800000e+03],\n",
       "          ..., \n",
       "          [  7.38500000e+03,   2.21520000e+04,   2.21510000e+04],\n",
       "          [  7.38500000e+03,   2.21530000e+04,   2.21520000e+04],\n",
       "          [  2.57200000e+03,   1.24910000e+04,   2.21530000e+04]])},\n",
       "  {'np': 22154, 'nt': 44304, 'p': array([[-5.51176 ,  1.365883, -7.191366],\n",
       "          [-6.246181,  1.45124 , -7.382205],\n",
       "          [ 3.1572  ,  0.030886, -5.749983],\n",
       "          ..., \n",
       "          [ 1.876025,  4.044248, -3.836338],\n",
       "          [ 1.51799 ,  4.409166, -3.940802],\n",
       "          [ 1.221372,  4.462838, -3.631302]]), 't': array([[  1.90000000e+01,   7.45700000e+03,   7.38600000e+03],\n",
       "          [  3.00000000e+00,   7.46900000e+03,   7.38700000e+03],\n",
       "          [  3.00000000e+00,   7.38700000e+03,   7.38800000e+03],\n",
       "          ..., \n",
       "          [  7.38500000e+03,   2.21520000e+04,   2.21510000e+04],\n",
       "          [  7.38500000e+03,   2.21530000e+04,   2.21520000e+04],\n",
       "          [  2.57200000e+03,   1.24910000e+04,   2.21530000e+04]])}]]"
      ]
     },
     "execution_count": 198,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g_f01=geom(\"Adult-F03\",\"Adult-F01\",\"Adult-F03_to_Adult-F01\");\n",
    "g_f16=topo(g1,\"P32-F16\",\"Adult-F01_to_P32-F16\");\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20\n"
     ]
    }
   ],
   "source": [
    "a=[1,2,3];\n",
    "b=[2,3,4];\n",
    "print numpy.dot(a,b);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Regression\n",
    "\n",
    "Running the surface resampling in python was too slow, I moved to meshgeometry. It turned out that all the pipeline could be done using only the 'resample' function in meshgeometry (code sbn_4D.sh). I wrote a script that produced ~6000 surfaces for the different ages, all with the topology of Adult-F01, using all different paths to go from Adult-F01 to the 2 P0 surfaces. Now, for each vertex in Adult-F01 I have a series of positions describing the transition from P0 to Adult.\n",
    "\n",
    "I want to use a non-linear regression to fit the trajectory of each vertex. I will use an algorithm similar to LOWESS: at each time point I will fit a regression line based on a weighted combination of the neighbouring data. I'll try the algorithm in simulated data first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 363,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEACAYAAAC08h1NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAG+dJREFUeJzt3Xt4HNV9//H32pKQZF2MbLB8Ayc2TsAJYMLFJNAI0v5q\n",
       "8muAQFMCPIQCTaBtEtIkxSY0xeXJrW6bpPlBLk9oKElpSSAtPxIg2AbcJIC5BNsQQGBztQGD8UWW\n",
       "5Jtsnf5xjqtls7ua9ezOmbPzeT3PPJpZjWc+trVfzZ45cw6IiIiIiIiIiIiIiIiIiIiIiIiIiEjg\n",
       "5gO9wBpgQZn9jgP2AGclEUpERGpjLLAWmAE0AquAw0vsdy/wc+DspMKJiMjvGhPzzx+PLfwvAkPA\n",
       "zcAZRfb7FHArsDHm+UREJKa4hX8qsC5ve717rXCfM4DvuG0T85wiIhJD3MIfpYh/E1jo9s25RURE\n",
       "PGmI+edfAabnbU/HXvXnew+2CQhgInAatlno9oL91gIzY+YREcma54BZSZ6wwZ10BtBE6Zu7+9xA\n",
       "6V49oTcBLfIdIIZFvgPEtMh3gJgW+Q4Q0yLfAWJa5DtATBXXzrhX/HuATwJ3Y3vu/AvwNHCp+/73\n",
       "Yh5fRESqLG7hB7jLLflKFfyLqnA+ERGJIe7NXRmx3HeAGJb7DhDTct8BYlruO0BMy30HiGm57wBZ\n",
       "Fnobv4iIDxXXTl3xi4hkjAq/iEjGqPCLSJbksM8TdfkOIpba+EWkllqgYwk074IDdkH7T7GDS4Yu\n",
       "6NoZdHgRSbu2b8CHtsMuA9sNvH8Qmr/oO1UVBF07gw4vImk34VG424Bxy48NTFjiO1UVqFePiEhx\n",
       "e9bA0iG7boClu2HXGq+RRFf8IlJTk2HcOjimD961DdrWAAf6DlUFFdfONA2RvG/YZhGRWhkHnAwM\n",
       "A78EdvqNUxVB105d8YuIVE5t/CIiUp4Kv4hIxqjwi4hkjAq/iEjGqPCLiGSMCr+ISMao8IuIZIwK\n",
       "v4hIxqjwi4hkjAq/iEjGNPgOIJJhLcC52IHC7gVW7udxZgPdwJPApupEEylvPtALrAEWFPn++cBq\n",
       "4HHgfuDIEsfRWD2SJS3Q/gT0DMBf7IL2QeDMyg/T/k/QsR3evRVa+rEDkEm2JF47xwJrgRnYKcxW\n",
       "AYcX7HMi0OnW5wMrShxLhV+y5ONw6iAMu0lBfmmgbUOFxzgZugdgkzvGnQZa36xJWkmzxAdpOx5b\n",
       "+F8EhoCbgTMK9nkQ6HPrDwHTYp5TpB50wZFNI6PpHg4MdVR4jMPg9xiZN3w+sPNAoLlKGaVOxS38\n",
       "U4F1edvr3WulXALcGfOcIvXgXvjBbnsttAX47C5ovrfCY/wW7snBK27zZqDlNepjjHmpobg3dyv5\n",
       "iHEKcDHwvjL7LMpbX+4WkXr0CPRfBB/4Nuxug9Z7oO/8Co/xMOy4Bmb9HXTthq1DsP1DNUkradLj\n",
       "Fm/mAb/I276S4jd4j8Q2Cc0qcyy18Yvsn4OAI7C9hCR7Eq+dDcBz2Ju7TRS/uXsItujPG+VYKvwi\n",
       "IpXzUjtPA57BFvcr3WuXugXgemzf4pVuebjEcVT4RUQqF3TtDDq8iIgnmnNXRETKU+EXEckYFX4R\n",
       "KTQb2r4P428GPug7jNQ3tfGL+DcLWrbB1XvhOwYmDELuXN+hpKyga2fQ4UXqQ9NiuGKvHfvHGFhm\n",
       "YPyzvlNJWbq5KyJxjGmC9ry60A72GR2RmtAVv4h/x0HbIPy7gXsMzB6A5oW+Q0lZQdfOoMOL1JFT\n",
       "oWsFdD0BzZ9lZAhRSaega2fQ4UVEPKm4dmrqRZHKtQJvB95wi0TXDlwEuQPBLAUe8B1I/NIVv4Tg\n",
       "WGjdDNO2QfNOaFX7d3Rt0LYWTt8OV+6F8YOQ+6jvUHUg6NoZdHjJhByM2wC3uK6O6w10DQLH+Q4W\n",
       "iMvgtMGRrqIPGmjTJ6b41J1TpIZaYddEONttTgVOHQbmeMwUkk6Yndc19O3A0DhvaSQVdMUvaZeD\n",
       "lk1wl7ti3Wjg4AHgvb6DBWIutG23D4W9bOCsHdB5q+9QdSDo2hl0eMmMk+2QBnO2Qvt2aP+q70CB\n",
       "+RB0vAwtfdBxM6Ar/viCrp1Bh5dM6QJOAmb6DiJC4LUz6PAiIp7o5q6IiJSnwi8ikjEq/CIiGaPC\n",
       "LyKSMSr8IiIZo8IvMmIscIDvECIhmA/0AmuABSX2+Zb7/mpgbol91J1TfMlB6zXQMARj90DHEtzU\n",
       "UyIBSLx2jgXWAjOARmAVcHjBPh8E7nTrJwArShxLhV+SdjBwLHARvG0AXjGw08Cf7ISOH/kOJxJR\n",
       "4rXzROAXedsL3ZLvu8A5edu9wKQix1LhlwQ1XQwtO2BWH7Tuhq+bkVEjHzfQud53QpGIEn+Aayqw\n",
       "Lm97vXtttH2mxTyvSByHQMO1sLIZ1nTAxxrh14y8fx42kHvVZ0CRWoo7A1fU3zSFc3aW+nOL8taX\n",
       "u0Wk2mbBO3fDO1rs5mJg2jCcsAO6jR09cscnvCaUCEwDcCAwATt+Upfb3reMz1s6CpZdkJuRdOIq\n",
       "6XHLfotb+F8BpudtT8de0ZfbZ5p7rZhFMfOIRPEc9DbZ/gaHYVsfd+2ERy7Ffgq+l9I/o1Iz5gDs\n",
       "fZfC5SBgYsHXCdgb8FuBTcBmt2zJW9YBTwB9ecs2oN8toVrOWy+Kr670AHEL/6PYd84M4FVsW/65\n",
       "BfvcDnwSuBmYh/2Pej3meUXieAn2Xg5HfQum74J1jTB0LvZnVarKNGDv6U0Guossk9xyMHYu443Y\n",
       "+rARO5/xvq/PAG+67U3u61bIDSf4l6kbhU0w++M04JvYHj7/AnwVuNR973vu67XYbp+DwEXAY0WO\n",
       "Y6qURySqbuBQ4DlsUZHIzBjslfdUt0xxS/76ZGzzy5vABuC1vK+vu2VD3vpWyKmTR+Uqrp1pKrQq\n",
       "/BJXDvspdsh3kLCZHLZJ5RBsM+2+r9PylinY5pL12E/7r4x8HdoAl70XVs2F53th6zXYq3SpDRV+\n",
       "yarGC2HMt2GoGdoehW2noybFEkwO27TyNrfMyFsOxRb6ncDLblmXt6x3X1+F3M7ix2+7Dt7+p/C5\n",
       "Vrh/N9y0AQbnAAO1+htlnAq/ZNIJMP5e+HUrvAP46z1wwwroO9l3MH/MGOxV+qyCZSZ2lvOdwAt5\n",
       "y0vAi+7ry5Db35ufjTB2O2xssB1rAN43AA9cBGh+3dqouHbGvbkrkgbvg/MaYI7bvKYBrj3Ba6LE\n",
       "mHbgnQXLbGyR34zturTWfX0Iez/jecj11SjQGCAHzXkvtYFqTaroP0PqwevwyG7Y22T7GDwGNG+p\n",
       "r5YF0wIcAbwbeBf2t9wcbFv8s8BT2H6pP8H2gFkLOQ//ALkzoXMrHH4gfHYM7BiGX+0GliWfRUKg\n",
       "u/myvxqh41dwZD+cOwDjBrFjRAXKTAJzGpgrwfwYzNNgdoB5HMxN7vXTwcwEM9Z32hFjPwbdg3CL\n",
       "gR8YaBuGcfdjP31I7QRdO4MOL941AB8G/gzb0B8IMwnMH4G5GszPwLwCZhOYe8D8I5gLwBwJpsl3\n",
       "0tF1rYa78sY8+icDHf/qO1UGVFw71dQj9WIP8F++Q5RnmoFjsA8yzgOOBzqxD0I+Avwr8CngpYL+\n",
       "7Adh2+2fB7YnGLhSBvbmbe7F/gIQKU0/IFJnzEFgznRX7g+AGQDzKJhr3ZX8Ya5rZRktn4fmnTBt\n",
       "G7Rswf6ySKncOTBxEH5o4DoDrYPYYa+ltoKunUGHFwEzGcxHwXwXzFNg+sDcBeYqMKeAGVfhAY+B\n",
       "AwdhnWs6+amB1jdId7fnM6DrFzD+Nuz8G1J7QdfOoMNLFpkuMGeDuc7dgN0M5jYwfwXmmCrceL0A\n",
       "ztw20mZuDDQOYUeXFNkn6NoZdHjJAtMI5mQwXwLzMJht7or+81Uq9IXm2aaTja7oLzWuuSfNV/yS\n",
       "vKBrZ9DhpV6ZqWD+DMxPwWwB8xswXwHT44YRrrG2r0D7dpizFVr6iTkOu9SloGtn0OGlXpgxYI4F\n",
       "cw2Yx1zzzc1gPma7XnoxEzgJOwa9SKGga2fQ4SVkpgnMH4L5jutH3wvmH1yzjro8S9oFXTuDDi+h\n",
       "MS2uq+WP3FX9g2CuABPQw18iQOC1M+jwEgLTDObDYP4DzFYw94L5SzBTfCcTiSHo2hl0eEkr0+Ca\n",
       "cW50N2fvA3OZx/b6pLVD503Q8RqMX4361tejoGtn0OElTUwOzPFgvgXmdTArwHzaPmCVNR13wZ/s\n",
       "hKcN3GRcz6BDfaeSqgq6dgYdXtLAHALmC+7m7BowfwsmyyNDNsCYPbAj7wGwPx4ALvEdTKoq6NoZ\n",
       "dHjxxbSAOQ/MUuyolt8Fc+LoY+BkQg4ad8GLrugPGzixH/io72BSVUHXzqDDS9LMXDdUwiYwd7sx\n",
       "cppH/3NZ07wApg7C3xs4awe09QItvlNJVQVdO4MOL5HkgA8BnwFOrfyPm3Ywl2JHuHzRNeUcUt2I\n",
       "del0aP5nyF2BmwdR6krQtTPo8DKqHLTfADP74bKd0D0A4/4u2h81R7smnC1g/hPM/BqMiyMSqsRr\n",
       "ZxewFDvn5xJgfJF9pgP3AU8CvwU+XeJYKvz17UiYMAj9rr15g7HjzDOx+O6mGTtm/YNgXgbzxTL9\n",
       "7RuAM4GLgcNqkl4kvRKvnYuBK9z6AuBrRfbpBo52623YiaAPL7KfCn99OxXmbuUtQwxP6sfOLJXH\n",
       "HALmq64b5t3YuWXLDZvQCO33w7v74SMDbvKPP6zh30MkbRKvnb3Avgdhut32aG4DPlDkdRX++jYB\n",
       "WvrgVgPbDVw3DONeA5pcv/tTXDPOJjDfADN71CNaF8AJA7Anb+jitldr+RcRSZnEa+eWvPVcwXYx\n",
       "M4CXKH6DSYW//h0P7S/ZvuUdT8PRR2GHPH4cO2PVn4Op9ObjX8Ond498ithioHFnTdKLpFNNJltf\n",
       "ir2aL3RVkZOXC9AG3ApcDgyU2GdR3vpyt0j9eBj6DwUzFfgL7M/Ww8DngGUFE4xHdT/8aAj+vBFm\n",
       "AVcPQesK6KtibJFU6cHzvAy9jPxSmEzppp5G4G5sN75SdMVf98yxYG5yo2H+PzBVuhHbcAk07YAx\n",
       "e6FjBXBwdY7LEdD5KLRtgs57AA3mJmnk5ebuAre+kOI3d3PAD4FvjHIsFf66ZMaCOQvMr8C8BOZz\n",
       "YIr1/oorBzRV8XjjoXUTfHsYnjewYAjanwHUjVTSxkt3zmX8bnfOKcAdbv0kYBhYBax0y/wix1Lh\n",
       "ryum3Q2M9rzrkvmRUXrnpM3vv7UX0rCBAwex96lE0iTo2hl0eNnHTAezGMybYG6x4+YE6QSY2g+7\n",
       "XeHfbKB5FyWfOxDxJujaGXR4McdjJzjZDObrYGb4ThTTGOhYAicOwpcMzB6Atut8hxIpIujaGXT4\n",
       "bDJjwZwN5teu/f6zYDp9p6qiBuDj0LgYOAd7H0EkbYKunUGHzxbTCeavwLwA5oEA2+9F6knQtTPo\n",
       "8NlgDsPOarXZNetoGj8R/4KunUGHr19mDHbO2p+DeQPMV8BM851KRP5X0LUz6PD1x3S67pi9YFaB\n",
       "uQSMJvAQSZ+ga2fQ4euHmQvme3nNOSehaQxF0izo2hl0+LCZcWAuBrPCjX3/N2Am+04lIpEEXTuD\n",
       "Dh8ek3N97/dd3f9/MP8XzWwlEpqga2fQ4cNhpoC5AsyTYNaC+YIbLVNEwhR07Qw6fLqZ8WAuArMM\n",
       "O2/t9WBOVtu9SF0IunYGHT59TCeY810TTh+Y28D8sXrmiNSdoGtn0OHTwUwFcxmYO8FsA3M7mI/V\n",
       "2TAKIvJWQdfOoMP7YZrA/J57qGoldr7afwNzDph23+lEJBEV1840tfEa0pUnhUwDcAzwfuzUaycD\n",
       "a7Czm90FPAi5Pd7iiYgPFdfONBVaFf63MDlgGnAscDxwIvAe4AXgv91yH+Q2eYsoImmgwh8ekwMm\n",
       "ALOBI4A5wJHAUcBe4BHgN8ADwEOQ2+opqIikkwp/upixQCd28u+DsRPTTwOmAodip/Gb6XZeAzwF\n",
       "PA08jp2qcgPkdO9DRMpR4ffDfBt7xd4KtAEd2ILfDvQDG4HX3bIeeAV40S0vAJtU4EVkP6nw+2F6\n",
       "sLM1bXdLH7AN2Aq5vR6DiUj9C7h2qjuniMj+qLh2jqlFChERSS8VfhGRjIlT+LuApcCzwBJgfJl9\n",
       "xwIrgZ/FOJ+IiFRBnMK/EFv4ZwP3uO1SLsd2VVQ7vohIwHqBSW69220XMw1YBpxC+St+/VIQEalc\n",
       "orVzS956rmA73y3AXOz4Mir8IiLVVXHtbBjl+0uxV/OFripy4mIn/yPgDWz7fk+EPIvy1pe7RURE\n",
       "RvQQrZ7WRC8jvxQmU7yp5yvAOuzTqa8Bg8APSxxPV/wiIpVLtHYuBha49YXA10bZX009IiLVl2jt\n",
       "7MLetC3szjkFuKPI/u8Hbi9zPBV+EZHKBV07gw4vIuKJhmwQEZHyVPhFRDJGhV9EJGNU+EVEMkaF\n",
       "X0QkY1T4RUQyRoVfRCRjVPhFRDJGhV9EJGNU+EVEMkaFX+LKQfPfwLjNMG4LtH7ZviYiMjqN1ROk\n",
       "pk/ArAF42sCzBuYMwAGf8Z1KJEOCrp1Bh8+uCcvgxwaMW243MOEB36lEMkSDtEnShjbCmuGR7TUG\n",
       "9r7pL4+IhERX/GE6DFq2wsW74NLd0LINmOM7lEiGVFw703QTzpCuPBLddOAc7P/fLcCLXtOIZEvQ\n",
       "tVNX/CIilVMbv4iIlKfCLyKSMSr8IiIZo8IvIpIxKvwiIhmjwi8ikjFxCn8XsBR4FlgCjC+x33jg\n",
       "VuBp4ClgXoxzioiIR4uBK9z6AuBrJfa7EbjYrTcAnSX2Uz9+EZHKJVo7e4FJbr3bbRfqBJ6PeDwV\n",
       "fhGRyiVaO7fkrecKtvc5GngIuAF4DPg+0FrieCr8IiKVq7h2Nozy/aXYq/lCVxU5cbGTNwDHAJ8E\n",
       "HgG+CSwE/rbE+RblrS93i4iIjOhxixe9jPxSmEzxpp5u4IW87ZOAn5c4nq74RUQql+hYPbcDF7r1\n",
       "C4HbiuyzAVgHzHbbvw88GeOcIiLiURewjN/tzjkFuCNvv6OwzTyrgf9EvXpERKop6NoZdHgREU80\n",
       "LLOIiJSnwi8ikjEq/CIiGaPCLyKSMSr84csBY32HEBHZH+rVU7HGS+GAQRizBzrvAw70nUhEEhd0\n",
       "7Qw6vAfvhwmD8JSBnQYu2QWdd/oOJSKJC7p2Bh3egy/Cwr1gjF1eNdDc7zuUiCRO/fgz5A14aOfI\n",
       "//kqoGmTz0AiIpXSFX9lmqF9JRzXD+cNQusgdiwkEcmWimtnrhYp9pMhXXlCcABwBnacpOXYcZNE\n",
       "JFuCrp264hcRqZza+EVEpDwVfhGRjFHhFxHJGBV+EZGMUeEXEckYFX4RkYxR4RcRyRgVfhGRjFHh\n",
       "FxHJmAbfAWS/HQB8HBpnwND9wH95ziMiUjEN2RBdA3Q8CKduhy8ZOHQAWr/sO5SIeJFo7ewClmIH\n",
       "BluCHSismCuBJ4EngH/HXqkWo8If3f+B2f2wx43Fv8FAw26g2XcwEUlcomP1LMQW/tnAPW670Azg\n",
       "48AxwLuxc8N+NMY5xWqHacMjU+0eBIw1QIvHTCKSAb3AJLfe7bYLdQHPYOeCbQB+Rukx43XFH90k\n",
       "aOmDG4bhOQN/uQvaf0PAQ7OKyH5LtHZuyVvPFWzn+wTQD7wB/KjM8VT4KzMXOldD25vQeQcw0Xcg\n",
       "EfGi4to5Wq+epdir+UJXFTlxsZPPBD6DbfLpA24BzgduKnG+RXnry90ixa2EvqN8hxCRxPW4xYte\n",
       "Rn4pTKZ4U885wPV52xcA15U4nq74RUQql+jN3duBC936hcBtRfbpBeZhbzrmsO37T8U4p4iIeNQF\n",
       "LON3u3NOAe7I2+8KRrpz3gg0ljiervhFRCoXdO0MOryIiCeac1dERMpT4RcRyRgVfhGRjFHhFxHJ\n",
       "GBV+EZGMUeEXEckYFX4RkYzRDFz1JQecgR0Cew3wE2DYayIRkTL0AFdsbf8MbxuAhXvhXQPQ/hM0\n",
       "VLNIvQu6dgYdPgW6oWUnbHazcm03MHEQe/UvIvVLT+5m2HjoHLJz3oAdF697iNJTYoqIeKcr/nga\n",
       "Ydw6WLwXXjdw/TC0bgI6fQcTkZoKunYGHT4l3g6dj8IBg9DxJGrmEcmCoGtn0OFFRDxRG7+IiJSn\n",
       "wi8ikjEq/CIiGaPCLyKSMSr8IiIZo8IvIpIxKvwiIhmjwi8ikjFxCv9HgCeBvcAxZfabD/Rihwle\n",
       "EON8IiLi2TuB2cB9lC78Y4G1wAygEVgFHF5i39Cf3O3xHSCGHt8BYurxHSCmHt8BYurxHSCmHt8B\n",
       "Ykr0yd1e4NlR9jkeW/hfBIaAm7EThdSjHt8BYujxHSCmHt8BYurxHSCmHt8BYurxHSBptW7jnwqs\n",
       "y9te714TERFPRpt6cSnQXeT1LwA/i3D80JtvRETqTjWm5bsP+BzwWJHvzQMWYW/wAlyJnQP274vs\n",
       "uxaYWYU8IiJZ8hwwK+mT3ge8p8T3GrChZgBNlL+5KyIiKfdhbPv9DmADcJd7fQpwR95+pwHPYK/o\n",
       "r0wyoIiIiIiIeNSFvXn8LLCE4pOCT8c2JT0J/Bb4dGLpiovyMNq33PdXA3MTyhXVaPnPx+Z+HLgf\n",
       "ODK5aJFEfRjwOGAPcFYSoSoQJX8PsBL78748kVTRjZZ/IvALbJPub4E/TSzZ6H4AvA48UWafNL93\n",
       "R8uf9vfu/1oMXOHWFwBfK7JPN3C0W2/DNhn5ukcQ5WG0DwJ3uvUTgBVJhYsgSv4TGZmcfT7h5d+3\n",
       "373Az4GzkwoXQZT847EXOdPc9sSkwkUQJf8i4KtufSKwidF7DiblZGwxL1U40/zehdHzV/Te9TlW\n",
       "z+nAjW79RuDMIvtswP6AAQwAT2PvIfgQ5WG0/L/TQ9g38qSE8o0mSv4HgT63/hAjBSgNoj4M+Cng\n",
       "VmBjYsmiiZL/POCn2OddAN5MKlwEUfK/BnS49Q5s4d+TUL7R/ArYUub7aX7vwuj5K3rv+iz8k7Af\n",
       "XXBfR/tHnoH9jfdQDTOVE+VhtGL7pKV4Vvow3SWMXAGlQdR//zOA77jtND1HEiX/Ydgm0PuAR4EL\n",
       "kokWSZT83wfmAK9imx0uTyZaVaT5vVupUd+7tf4YVuoBsKsKtg3l36Rt2Ku4y7FX/j5ELSKFz0ak\n",
       "pfhUkuMU4GLgfTXKsj+i5P8msNDtm6M6z6lUS5T8jdhxrz4AtGKv4lZg2519i5L/C9hP6D3YZ3KW\n",
       "AkcB/bWLVVVpfe9WItJ7t9aF/w/KfO917C+FDcBk4I0S+zViP/7+G3BbVdNV5hXszeZ9pjPykbzU\n",
       "PtPca2kQJT/Ym0Lfx7YTlvtombQo+d+DbYIA28Z8GrZZ4vaapxtdlPzrsM07O9zyS2zhTEPhj5L/\n",
       "vcCX3fpzwAvAO7CfXtIuze/dqNL63n2LxYz0DFhI8Zu7OeCHwDeSClVGlIfR8m8QzSNdN4ii5D8E\n",
       "2447L9Fk0VT6MOANpKtXT5T87wSWYW+ktmJv5B2RXMSyouT/OnC1W5+E/cXQlVC+KGYQ7eZu2t67\n",
       "+8ygdP40v3ffogv7Q17YnTP/AbCTsEM8rMJ2cVvJyPAPPhR7GO1St+xzrfv+asrPU+DDaPmvx96Q\n",
       "2/dv/XDSAUcR5d9/n7QVfoiW//PYnj1P4L/7cqHR8k/EjuG1Gpv/vKQDlvEf2HsPu7GfrC4mrPfu\n",
       "aPnT/t4VEREREREREREREREREREREREREREREREREcmG/wEiULGWkR6zxwAAAABJRU5ErkJggg==\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x107d0f190>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "N=20;\n",
    "x=numpy.random.rand(N);\n",
    "y=numpy.random.rand(N);\n",
    "g=0.5;\n",
    "y=numpy.add((1-g)*y,g*numpy.log(x));\n",
    "pylab.scatter(x,y);\n",
    "\n",
    "c=0.9;\n",
    "N=100;\n",
    "a=numpy.zeros(N);\n",
    "b=numpy.zeros(N);\n",
    "for tt in range(0,N):\n",
    "    a[tt]=tt/float(N);\n",
    "    f=0;\n",
    "    wsum=0;\n",
    "    for i in range(len(x)):\n",
    "        w=(1-c) + c/(numpy.power(0.1*(x[i]-a[tt]),2)+0.001);\n",
    "        f+=y[i]*w;\n",
    "        wsum+=w;\n",
    "    b[tt]=f/float(wsum);\n",
    "pylab.plot(a,b);\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}