AFNI_ANTs_affine_transforms.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.io import loadmat, savemat\n",
    "from numpy import loadtxt, around, hstack, vstack, zeros, float64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'AffineTransform_double_3_3': array([[ 9.99771194e-01],\n",
       "        [-4.22685205e-05],\n",
       "        [-1.66161657e-05],\n",
       "        [-4.15907894e-06],\n",
       "        [ 8.65957305e-01],\n",
       "        [-4.99806282e-01],\n",
       "        [-7.38437417e-05],\n",
       "        [ 4.99559786e-01],\n",
       "        [ 8.65365078e-01],\n",
       "        [-1.17959414e-03],\n",
       "        [ 8.09175497e+00],\n",
       "        [ 4.96808506e+00]]), 'fixed': array([[0.],\n",
       "        [0.],\n",
       "        [0.]])}"
      ]
     },
     "execution_count": 86,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loadmat('sub-001_desc-deobliquedrotatedrealignedants_T1w0GenericAffine.mat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'AffineTransform_double_3_3': array([[ 9.99771194e-01],\n",
       "        [-4.22685205e-05],\n",
       "        [-1.66161657e-05],\n",
       "        [-4.15907894e-06],\n",
       "        [ 8.65957305e-01],\n",
       "        [-4.99806282e-01],\n",
       "        [-7.38437417e-05],\n",
       "        [ 4.99559786e-01],\n",
       "        [ 8.65365078e-01],\n",
       "        [-1.17959414e-03],\n",
       "        [ 8.09175497e+00],\n",
       "        [ 4.96808506e+00]]), 'fixed': array([[0.],\n",
       "        [0.],\n",
       "        [0.]])}"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loadmat('sub-001_desc-deobliquedrotatedrealignedants_T1w0GenericAffine.mat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "orig_ants_mat = loadmat('sub-001_desc-deobliquedrotatedrealignedants_T1w0GenericAffine.mat')['AffineTransform_double_3_3']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "orig_afni_mat = loadtxt('sub-001_desc-deobliquedrotatedrealigned_T1w.aff12.1D')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  , -0.  , -0.  , -0.  ],\n",
       "       [ 0.87, -0.5 , -0.  ,  0.5 ],\n",
       "       [ 0.87, -0.  ,  8.09,  4.97]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "around(orig_ants_mat, 2).reshape(3,4, order='C')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  , -0.  , -0.  , -0.  ],\n",
       "       [-0.  ,  0.87,  0.5 ,  8.09],\n",
       "       [-0.  , -0.5 ,  0.87,  4.97]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "around(orig_ants_mat, 2).reshape(3,4, order='F')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  , -0.  , -0.  , -0.  ],\n",
       "       [ 0.87, -0.5 , -0.  ,  0.5 ],\n",
       "       [ 0.87, -0.  ,  8.09,  4.97]])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "around(orig_ants_mat, 2).reshape(3,4, order='A')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  , -0.  ,  0.  ,  0.  ],\n",
       "       [-0.  ,  0.87, -0.5 ,  8.08],\n",
       "       [ 0.  ,  0.5 ,  0.87,  4.96]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "around(orig_afni_mat, 2).reshape(3,4, order='C')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  ,  0.  , -0.5 ,  0.5 ],\n",
       "       [-0.  , -0.  ,  8.08,  0.87],\n",
       "       [ 0.  ,  0.87,  0.  ,  4.96]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "around(orig_afni_mat, 2).reshape(3,4, order='F')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  , -0.  ,  0.  ,  0.  ],\n",
       "       [-0.  ,  0.87, -0.5 ,  8.08],\n",
       "       [ 0.  ,  0.5 ,  0.87,  4.96]])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "around(orig_afni_mat, 2).reshape(3,4, order='C')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  , -0.  , -0.  , -0.  ],\n",
       "       [ 0.87, -0.5 , -0.  ,  0.5 ],\n",
       "       [ 0.87, -0.  ,  8.09,  4.97]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "around(orig_ants_mat, 2).reshape(3,4, order='C')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  , -0.  , -0.  ],\n",
       "       [-0.  ,  0.87, -0.5 ],\n",
       "       [-0.  ,  0.5 ,  0.87]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "around(orig_ants_mat[:9], 2).reshape(3,3, order='C')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "mat=[reshape(afftransform(1:9),[3,3])',afftransform(10:12)];\n",
    "\n",
    "m_Translation=mat(:,4);\n",
    "mat=[mat;[0,0,0,1]];\n",
    "\n",
    "\n",
    "for i=1:3\n",
    "    m_Offset(i) = m_Translation(i) + m_Center(i);\n",
    "    for j=1:3\n",
    "       m_Offset(i) = m_Offset(i)-(mat(i,j) * m_Center(j));  % (i,j) should remain the same since in C indexing rows before cols, too.\n",
    "    end\n",
    "end\n",
    "\n",
    "mat(1:3,4)=m_Offset;\n",
    "mat=inv(mat);\n",
    "\n",
    "% convert RAS to LPS (ITK uses LPS)\n",
    "mat=mat.*...\n",
    "    [1  1 -1 -1\n",
    "     1  1 -1 -1\n",
    "    -1 -1  1  1\n",
    "     1  1  1  1];\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.  , -0.  , -0.  , -0.  ],\n",
       "       [-0.  ,  0.87, -0.5 ,  8.09],\n",
       "       [-0.  ,  0.5 ,  0.87,  4.97],\n",
       "       [ 0.  ,  0.  ,  0.  ,  1.  ]])"
      ]
     },
     "execution_count": 140,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def read_ants_affine(input_file, debug=False):\n",
    "    \n",
    "    def calculateOffset(affine_transform, translations, fixed):\n",
    "        offset = zeros(3)\n",
    "        for i in range(3):\n",
    "            offset[i] = translation[i] + fixed[i]\n",
    "            for j in range(3):\n",
    "                offset[i] = offset[i] - (affine_transform[i,j] * fixed[j])\n",
    "        return offset.reshape(3,-1)\n",
    "\n",
    "    mat_content = loadmat(input_file)\n",
    "    if debug:\n",
    "        print(mat_content)\n",
    "    matching_affine_keys = [k for k in mat_content.keys() if k.startswith('AffineTransform')]\n",
    "    assert len(matching_affine_keys) == 1\n",
    "    matching_affine_key = matching_affine_keys[0]\n",
    "    affine_transform = mat_content[matching_affine_key][:9].reshape(3,3, order='C')\n",
    "    translation = mat_content[matching_affine_key][9:]\n",
    "    fixed = mat_content['fixed']\n",
    "    offset = calculateOffset(affine_transform, translations, fixed)\n",
    "    combined = vstack((hstack((affine_transform, offset)),[0,0,0,1]))\n",
    "    return combined\n",
    "    \n",
    "around(read_ants_affine('sub-001_desc-deobliquedrotatedrealignedants_T1w0GenericAffine.mat'), 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_afni_affine(input_file, debug = False):\n",
    "    orig_afni_mat = loadtxt(input_file)\n",
    "    if debug:\n",
    "        print(orig_afni_mat)\n",
    "\n",
    "    combined = vstack((orig_afni_mat.reshape(3,4, order='C'), [0,0,0,1]))\n",
    "    return combined"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.  , 8.08, 4.96])"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "around(read_afni_affine('sub-001_desc-deobliquedrotatedrealigned_T1w.aff12.1D'), 2)[:3,3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 9.99771194e-01, -4.22685205e-05, -1.66161657e-05,\n",
       "        -1.17959414e-03],\n",
       "       [-4.15907894e-06,  8.65957305e-01, -4.99806282e-01,\n",
       "         8.09175497e+00],\n",
       "       [-7.38437417e-05,  4.99559786e-01,  8.65365078e-01,\n",
       "         4.96808506e+00]])"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hstack((affine_transform, calculateOffset(affine_transform, translations, fixed)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_ants_affine(affine, out_file, debug=False):\n",
    "    out_dict = {}\n",
    "    out_dict['AffineTransform_double_3_3'] = hstack((affine[:3,:3].reshape(1, -1), affine[:3,3].reshape(1, -1))).reshape(-1,1).astype(float64)\n",
    "    out_dict['fixed'] = zeros((3,1))\n",
    "    if debug:\n",
    "        print(out_dict)\n",
    "        \n",
    "    savemat(out_file, out_dict, format='4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.999991   -0.00110987 -0.00419122  0.0665791   0.00107965  0.999973\n",
      " -0.00720623 -0.0969192   0.00419911  0.00720163  0.999965    1.06055   ]\n",
      "{'AffineTransform_double_3_3': array([[ 0.999991  ],\n",
      "       [-0.00110987],\n",
      "       [-0.00419122],\n",
      "       [ 0.00107965],\n",
      "       [ 0.999973  ],\n",
      "       [-0.00720623],\n",
      "       [ 0.00419911],\n",
      "       [ 0.00720163],\n",
      "       [ 0.999965  ],\n",
      "       [ 0.0665791 ],\n",
      "       [-0.0969192 ],\n",
      "       [ 1.06055   ]]), 'fixed': array([[0.],\n",
      "       [0.],\n",
      "       [0.]])}\n"
     ]
    }
   ],
   "source": [
    "write_ants_affine(read_afni_affine('C:/stanford_gdrive/data/transformation_test_data/vol0001tovol0020.aff12.1D', True), \n",
    "                  \"test.mat\", True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [],
   "source": [
    "savemat?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(238, 12)"
      ]
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loadtxt(\"C:/stanford_gdrive/data/transformation_test_data/volred.aff12.1D\").shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}