{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"coregistration.ipynb","provenance":[],"toc_visible":true,"mount_file_id":"1LrAu38-K_ZjOKSZkDLmlfTOjR38stV47","authorship_tag":"ABX9TyPiZkcNQ57aRU7WOSPIv5jW"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"markdown","metadata":{"id":"gyd8tZLembn_"},"source":["# Coregistration"]},{"cell_type":"markdown","metadata":{"id":"On_kC6p_mjxm"},"source":["- Three different images need to be coregistered: EM, 2p structural stack, and 2p video.\n","- 2p structural image is used as reference. EM and 2p video are coregistered to 2p structural image"]},{"cell_type":"code","metadata":{"id":"MY4vuyJCmYNL","executionInfo":{"status":"ok","timestamp":1616646974032,"user_tz":-540,"elapsed":704,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}}},"source":["import numpy as np\n","import pandas as pd\n","\n","from pydrive.auth import GoogleAuth\n","from pydrive.drive import GoogleDrive\n","from google.colab import auth\n","from oauth2client.client import GoogleCredentials"],"execution_count":9,"outputs":[]},{"cell_type":"code","metadata":{"id":"VT2yLztirRA4","executionInfo":{"status":"ok","timestamp":1616647141182,"user_tz":-540,"elapsed":33769,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}}},"source":["auth.authenticate_user()"],"execution_count":12,"outputs":[]},{"cell_type":"code","metadata":{"id":"DVYF9fe3q-X_","executionInfo":{"status":"ok","timestamp":1616647930689,"user_tz":-540,"elapsed":1190,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}}},"source":["gauth = GoogleAuth()\n","gauth.credentials = GoogleCredentials.get_application_default()\n","drive = GoogleDrive(gauth)"],"execution_count":17,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"J7vPpTV2ntk-"},"source":["## EM to 2p structural stack\n","- Affine transformation computed from manually annotated correspondence points.\n","- From $\\mu m$ to $\\mu m$."]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":419},"id":"R3hYBJRFu-7v","executionInfo":{"status":"ok","timestamp":1616648817874,"user_tz":-540,"elapsed":1162,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}},"outputId":"bafe79d3-674a-4ebb-c1a8-b366289f9a9f"},"source":["downloaded = drive.CreateFile({'id': '1E-8BIJ-5gSkGPMs1D_7VUz-05BHB8Tkd'})\n","downloaded.GetContentFile('correspondence_points.csv')\n","correspondence_df = pd.read_csv('correspondence_points.csv')\n","correspondence_df"],"execution_count":29,"outputs":[{"output_type":"execute_result","data":{"text/html":["
\n","\n","\n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n"," \n","
em_x_nmem_y_nmem_z_nm2p_x_um2p_y_um2p_z_um
0343611.4492250417.957914958.552610185.15625228.90625182.0
1355013.9059185993.06282590.900730194.53125225.78125113.0
2286345.4157182402.22803763.086764139.06250229.68750107.0
3289383.5417193791.09773976.920274140.62500231.25000118.0
4442196.8375225395.32903789.988336266.40625221.87500155.0
.....................
507304441.2743183007.608885893.699400146.09375153.90625126.0
508240674.3707204383.271684512.51684094.53125164.84375147.0
509238628.4303275948.970386287.09017089.06250175.78125221.0
510264332.4059258343.367386435.445200110.15625166.40625200.0
511250600.1024253139.636985827.925000100.00000168.75000195.0
\n","

512 rows × 6 columns

\n","
"],"text/plain":[" em_x_nm em_y_nm em_z_nm 2p_x_um 2p_y_um 2p_z_um\n","0 343611.4492 250417.9579 14958.552610 185.15625 228.90625 182.0\n","1 355013.9059 185993.0628 2590.900730 194.53125 225.78125 113.0\n","2 286345.4157 182402.2280 3763.086764 139.06250 229.68750 107.0\n","3 289383.5417 193791.0977 3976.920274 140.62500 231.25000 118.0\n","4 442196.8375 225395.3290 3789.988336 266.40625 221.87500 155.0\n",".. ... ... ... ... ... ...\n","507 304441.2743 183007.6088 85893.699400 146.09375 153.90625 126.0\n","508 240674.3707 204383.2716 84512.516840 94.53125 164.84375 147.0\n","509 238628.4303 275948.9703 86287.090170 89.06250 175.78125 221.0\n","510 264332.4059 258343.3673 86435.445200 110.15625 166.40625 200.0\n","511 250600.1024 253139.6369 85827.925000 100.00000 168.75000 195.0\n","\n","[512 rows x 6 columns]"]},"metadata":{"tags":[]},"execution_count":29}]},{"cell_type":"code","metadata":{"id":"btYuuOlSn7ax","executionInfo":{"status":"ok","timestamp":1616648845780,"user_tz":-540,"elapsed":686,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}}},"source":["n = correspondence_df.shape[0]\n","\n","em_coords = np.zeros((n,3))\n","twop_coords = np.zeros((n,3))\n","i = 0\n","for i in range(n):\n"," \n"," em_coords[i,:] = np.array([correspondence_df[\"em_x_nm\"][i],\n"," correspondence_df[\"em_y_nm\"][i],\n"," correspondence_df[\"em_z_nm\"][i]])/1000 # Divide by 1000 to convert unit to um.\n"," twop_coords[i,:] = np.array([correspondence_df[\"2p_x_um\"][i],\n"," correspondence_df[\"2p_y_um\"][i],\n"," correspondence_df[\"2p_z_um\"][i]])"],"execution_count":31,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"Yv2U42XMyF1o"},"source":["### Compute affine transformation\n","- This affine transformation transforms coordinates in EM space to 2p structural space. Both in physical dimension with unit $\\mu m$."]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"FKof7x_jxPui","executionInfo":{"status":"ok","timestamp":1616650062182,"user_tz":-540,"elapsed":841,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}},"outputId":"7d3c8c99-77ac-4441-bfcd-cac6931783aa"},"source":["X = np.concatenate((em_coords[:,:],np.ones((n,1))),axis=1)\n","\n","trans = np.dot(np.linalg.inv(np.dot(X.T,X)), np.dot(X.T,twop_coords[:,:]))\n","\n","A = trans[:3,:]\n","b = trans[3,:]\n","\n","print(A, b)"],"execution_count":36,"outputs":[{"output_type":"stream","text":["[[ 0.81958006 -0.08841949 0.01176497]\n"," [-0.02841662 0.16955195 1.01670694]\n"," [-0.10903189 -0.89056016 0.2186504 ]] [-89.93552264 227.50461563 -79.87499339]\n"],"name":"stdout"}]},{"cell_type":"code","metadata":{"id":"CsMHinFL6BUr","executionInfo":{"status":"ok","timestamp":1616650999631,"user_tz":-540,"elapsed":959,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}}},"source":["def em2struct(p):\n","\n"," return np.dot(p, A) + b"],"execution_count":38,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"vrjrldlfzXQB"},"source":["### Example transformation"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"fV-9Nlp-zWXM","executionInfo":{"status":"ok","timestamp":1616651009486,"user_tz":-540,"elapsed":912,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}},"outputId":"e38d1171-6482-4c1e-88ba-b83b304dcf31"},"source":["em_sample = np.array([[343611, 250417, 14959],\n"," [355014, 185993, 2591],\n"," [286345, 182402, 3763]])\n","em_sample = em_sample/1000 # Convert from nm to um.\n","\n","twop_sample = em2struct(em_sample)\n","for i in range(em_sample.shape[0]):\n"," print(\"{} um in EM --> {} um in 2P\".format(em_sample[i,:], twop_sample[i,:]))"],"execution_count":39,"outputs":[{"output_type":"stream","text":["[343.611 250.417 14.959] um in EM --> [182.93418791 226.25950707 182.03907437] um in 2P\n","[355.014 185.993 2.591] um in EM --> [195.45907833 225.34249304 113.96863426] um in 2P\n","[286.345 182.402 3.763] um in EM --> [139.15359385 229.76157338 109.76600892] um in 2P\n"],"name":"stdout"}]},{"cell_type":"markdown","metadata":{"id":"YGWrJ8sD2OeL"},"source":["## 2P video to 2P structural stack\n","- Translation that gives maximum correlation.\n","- From coordinate in video (256 x 256) to coordinate in structural stack (512 x 512)"]},{"cell_type":"markdown","metadata":{"id":"iOWWyluV4duP"},"source":["### Transformation\n","- 9 (scans) x 3 (slices) array for each translation.\n","- vshift_i : video shift in x dimension.\n","- vshift_j : video shift in y dimension.\n","- sshift_i : structure stack shift in x dimension.\n","- sshift_j : sturcture stack shift in y dimension."]},{"cell_type":"code","metadata":{"id":"BRG9l-r_2CHT","executionInfo":{"status":"ok","timestamp":1616650531429,"user_tz":-540,"elapsed":779,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}}},"source":["vshift_i = np.array([[0,0,-1],\n"," [-6,-6,-6],\n"," [-5,-5,-6],\n"," [-6,-6,-6],\n"," [-5,-5,-5],\n"," [-3,-4,-4],\n"," [-3,-3,-3],\n"," [-3,-3,-4]])\n","\n","vshift_j = np.array([[2,2,2],\n"," [0,0,-1],\n"," [0,-1,0],\n"," [-1,-1,-1],\n"," [-2,-2,-2],\n"," [-3,-3,-3],\n"," [-1,-1,-1],\n"," [-2,-2,-2]])\n","\n","sshift_i = np.array([[1,1,0],\n"," [0,0,0],\n"," [1,1,0],\n"," [0,0,0],\n"," [1,1,0],\n"," [1,0,0],\n"," [0,0,0],\n"," [0,0,0]])\n","\n","sshift_j = np.array([[1,1,1],\n"," [1,1,0],\n"," [1,0,1],\n"," [0,0,0],\n"," [0,1,0],\n"," [-1,0,0],\n"," [0,0,0],\n"," [0,0,-1]])"],"execution_count":37,"outputs":[]},{"cell_type":"code","metadata":{"id":"8iBXgqwH586x","executionInfo":{"status":"ok","timestamp":1616651331251,"user_tz":-540,"elapsed":926,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}}},"source":["def vid2struct(p, scan_id, slice_idx):\n","\n"," scan_id -= 1\n"," slice_idx -= 1\n"," \n"," dxv = vshift_i[scan_id, slice_idx]\n"," dxs = sshift_i[scan_id, slice_idx]\n"," dyv = vshift_j[scan_id, slice_idx]\n"," dys = sshift_j[scan_id, slice_idx]\n","\n"," x = (p[:,0] - dxv)*2 + dxs\n"," y = (p[:,1] - dyv)*2 + dys\n","\n"," return np.concatenate((x.reshape(-1,1), y.reshape(-1,1)), axis=1)"],"execution_count":41,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"jDRikr2n4bDA"},"source":["### Example transformation"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"W5zSGOEZ4W8p","executionInfo":{"status":"ok","timestamp":1616651368040,"user_tz":-540,"elapsed":824,"user":{"displayName":"J. Alexander Bae","photoUrl":"","userId":"06779053848267069857"}},"outputId":"3dc07da3-89bb-4c33-a8b5-aee9a3718b84"},"source":["scan_id = 2 # 1~9 : 9 scans\n","slice_idx = 1 # 1~3 : 3 slices\n","\n","vid_sample = np.array([[100,50],\n"," [150,50],\n"," [150,100]])\n","\n","twop_sample = vid2struct(vid_sample, scan_id, slice_idx)\n","for i in range(vid_sample.shape[0]):\n"," print(\"{} in video --> {} in 2P\".format(vid_sample[i,:], twop_sample[i,:]))"],"execution_count":43,"outputs":[{"output_type":"stream","text":["[100 50] in video --> [212 101] in 2P\n","[150 50] in video --> [312 101] in 2P\n","[150 100] in video --> [312 201] in 2P\n"],"name":"stdout"}]}]}