In [1]:
drive_path = '/Volumes/Brain2016 1'

%gui qt

import os
from allensdk.core.mouse_connectivity_cache import MouseConnectivityCache
import pandas as pd
import scipy.ndimage as ndi
import numpy as np
import time
import get_connectivity as gc
import vis3DConnect as vc

import pyqtgraph as pg
import pyqtgraph.metaarray as metaarray
from pyqtgraph.Qt import QtGui, QtCore
import pyqtgraph.opengl as pgl
pg.mkQApp()


Out[1]:
<PyQt4.QtGui.QApplication at 0x1162493e0>

In [2]:
# When downloading 3D connectivity data volumes, what resolution do you want (in microns)?  
# Options are: 10, 25, 50, 100
resolution_um=25

# The manifest file is a simple JSON file that keeps track of all of
# the data that has already been downloaded onto the hard drives.
# If you supply a relative path, it is assumed to be relative to your
# current working directory.
manifest_file = os.path.join(drive_path, "MouseConnectivity","manifest.json")

mcc = MouseConnectivityCache(manifest_file=manifest_file, resolution=resolution_um)

template, template_info = mcc.get_template_volume()

ontology = mcc.get_ontology()

In [ ]:

Set up your variables for the density plots


In [7]:
# FIND CONNECTIONS
fid_PRNc_III, fpd_PRNc_III=gc.get_connectivity('PRNc','III',ontology,mcc)
piIDs = list(fpd_PRNc_III.experiment_id)

In [8]:
# FIND CONNECTIONS
fid_SCm_PRNc, fpd_SCm_PRNc=gc.get_connectivity('SCm','PRNc',ontology,mcc)
spIDs = list(fpd_SCm_PRNc.experiment_id)

In [9]:
# GET PROJECTION DENSITIES
pdens1,_ = mcc.get_projection_density(158838128)
pdens2,_= mcc.get_projection_density(175158132)

In [10]:
# GET MASK DATA
PRNc_id = ontology.df[ontology.df.acronym=='PRNc'].id.values[0]
PRNc_mask, _ = mcc.get_structure_mask(PRNc_id)

SCm_id = ontology.df[ontology.df.acronym=='SCm'].id.values[0]
SCm_mask, _ = mcc.get_structure_mask(SCm_id)

In [11]:
# DEFINE THE INJECTION INPUT(S)
lt = [([255,15,0],pdens1),([0,15,255],pdens2)]
# lt = [([0,15,240],pdens2)]

# Note that the dominant RGB color should be 255, and the added color you want at the peak intensity
# (i.e., where the injection occurred) should be set to about 10 or 15.

RUNNING THE VISUALIZATION


In [12]:
# RUN THE ENTIRE VISUALIZATION
view = vc.vis3D(template,lt)


pyqtgraph/functions.py:2067: VisibleDeprecationWarning: converting an array with ndim > 0 to an index will result in an error in the future
  faceTableI[faceTableInds[:,0]] = np.array([triTable[j] for j in faceTableInds])
build brain isosurface 1.91
build injection volume 27.29
rendering 0.00
total run time: 29.19

In [ ]:
# ROTATE AROUND THE IMAGE 
for i in range(360): # determine the number of steps
    view.orbit(-2,0)  # determine the step size in horizontal and vertical dimensions
    time.sleep(0.1) # determine the time between steps
    pg.QtGui.QApplication.processEvents()

In [18]:
# JUST PLOT THE GLASS BRAIN
view = vc.vis3D_glassBrain(template,pad=30,ds_factor=6)

In [11]:
# PLOT SOME PROJECTIONS (ASSUMING AN EXISTING GLASS BRAIN)
view = vc.vis3D_projections(view,lt,ds_factor=6)

In [16]:
# PLOT A STRUCTURE MASK
view = vc.vis3D_structureMask(view,PRNc_mask,maskCol=[0.5,1,0.5],ds_factor=1)
view = vc.vis3D_structureMask(view,SCm_mask,maskCol=[1,0.5,1],ds_factor=1)

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:

MODIFYING THE VISUALIZATION


In [87]:
# REMOVE ALL ITEMS FROM VIEW
for i in range(len(view.items)):
    view.removeItem(view.items[0])  
# view.items

In [103]:
# REMOVE NON-GLASS-BRAIN ITEMS FROM VIEW
for i in range(len(view.items)-1):
    view.removeItem(view.items[1])

In [ ]:


In [ ]:


In [ ]:


In [ ]:

Plotting pathways


In [19]:
import urllib
import json

In [99]:
def get_path( target_voxel, experiment_id ) :
    url = "http://api.brain-map.org/api/v2/data/query.json?criteria=service::mouse_connectivity_target_spatial"
    url = url + "[seed_point$eq%s]" % ','.join([str(s) for s in target_voxel])
    url = url + "[section_data_set$eq%d]" % experiment_id
    
    response = urllib.urlopen(url)
    data = json.loads(response.read())
    data = [s['coord'] for s in data['msg'][0]['path']]
    return data

def vis3D_showPaths(view,paths,ds_factor=2):
    pts = paths[::ds_factor]
    plt = pgl.GLLinePlotItem(pos=pts, color=pg.glColor([255,0,0,255]), width=2, antialias=True)
    plt.translate(-528./2,-320./2,-456./2)
    plt.rotate(-90, 1, 0, 0)
    view.addItem(plt)
    
    pos = np.empty((2, 3))
    color = np.empty((2, 4))
    size = np.empty((2))
    
    pos[0]=paths[0,::]
    pos[1]=paths[len(paths)-1::]
    color[0]=[255,0,0,255]
    color[1]=[255,0,0,255]
    size[0]=10
    size[1]=10
    sp1 = pgl.GLScatterPlotItem(pos=pos, size=size, color=color, pxMode=True)
    sp1.translate(-528./2,-320./2,-456./2)
    sp1.rotate(-90, 1, 0, 0)
    view.addItem(sp1)
    
    view.show()

    return view

In [81]:
target_voxel = [9400,1250,6200]
experiment_id = 309113907
data = get_path( target_voxel, experiment_id )

In [82]:
# unpack it for plotting
x, y, z = zip(*data)

# divide by 25 
nx = np.array(x).astype(float)/25
ny = np.array(y).astype(float)/25
nz = np.array(z).astype(float)/25

xyz = np.vstack([nx,ny,nz]).transpose()

In [ ]:


In [89]:
# JUST PLOT THE GLASS BRAIN
view = vc.vis3D_glassBrain(template,pad=30,ds_factor=6)

In [104]:
# GET MASK DATA
v1_id = ontology.df[ontology.df.acronym=='VISp'].id.values[0]
v1_mask, _ = mcc.get_structure_mask(v1_id)
SCsg_id = ontology.df[ontology.df.acronym=='SCsg'].id.values[0]
SCsg_mask, _ = mcc.get_structure_mask(SCsg_id)

# PLOT A STRUCTURE MASK
view = vc.vis3D_structureMask(view,v1_mask,maskCol=[1,0.5,1],ds_factor=1)
view = vc.vis3D_structureMask(view,SCsg_mask,maskCol=[0.5,1,0.5],ds_factor=1)

In [105]:
# PLOT THE PATHS
view = vis3D_showPaths(view,xyz)

In [ ]:


In [62]:


In [79]:
xyz


Out[79]:
array([[ 376.  ,   52.  ,  248.  ],
       [ 375.52,   55.52,  248.  ],
       [ 374.2 ,   58.64,  248.  ],
       [ 371.52,   60.88,  248.44],
       [ 368.  ,   62.64,  249.32],
       [ 364.  ,   63.52,  250.64],
       [ 360.44,   64.  ,  252.  ],
       [ 357.32,   64.  ,  253.76],
       [ 354.64,   64.44,  256.  ],
       [ 352.  ,   65.76,  258.2 ],
       [ 349.32,   68.44,  260.  ],
       [ 347.08,   71.52,  261.32],
       [ 344.88,   74.64,  262.64],
       [ 342.64,   77.32,  263.52],
       [ 340.  ,   80.44,  264.  ],
       [ 337.76,   83.52,  264.  ],
       [ 336.  ,   86.64,  264.44],
       [ 334.64,   89.32,  265.32],
       [ 332.88,   92.44,  266.64],
       [ 331.08,   96.  ,  267.52],
       [ 328.88,   99.52,  268.  ],
       [ 326.64,  102.64,  268.44],
       [ 324.  ,  104.88,  269.76],
       [ 321.76,  107.08,  272.44],
       [ 320.  ,  108.88,  275.52],
       [ 318.2 ,  110.64,  278.64],
       [ 316.  ,  111.52,  280.88],
       [ 313.32,  112.44,  283.08],
       [ 311.08,  113.76,  284.88],
       [ 309.32,  116.  ,  287.08],
       [ 308.  ,  118.2 ,  289.32],
       [ 306.64,  119.52,  292.44],
       [ 305.32,  120.44,  296.  ],
       [ 304.  ,  121.76,  299.52],
       [ 302.64,  124.  ,  302.64],
       [ 300.88,  126.2 ,  305.32],
       [ 298.64,  127.52,  308.  ],
       [ 296.  ,  128.44,  310.2 ],
       [ 293.32,  129.76,  311.52],
       [ 290.64,  132.  ,  312.44],
       [ 288.  ,  134.2 ,  313.76],
       [ 285.32,  135.52,  316.  ],
       [ 283.08,  136.44,  318.2 ],
       [ 281.32,  137.76,  320.  ],
       [ 280.44,  140.44,  321.32],
       [ 280.  ,  143.52,  323.08],
       [ 280.  ,  146.2 ,  325.32],
       [ 280.  ,  147.52,  328.44],
       [ 280.  ,  148.  ,  332.  ],
       [ 280.  ,  147.52,  336.  ],
       [ 280.  ,  146.2 ,  340.  ],
       [ 280.  ,  144.  ,  343.52],
       [ 280.44,  141.32,  346.2 ],
       [ 281.32,  138.64,  348.  ],
       [ 282.64,  135.52,  349.76],
       [ 283.52,  132.44,  352.  ],
       [ 284.  ,  129.32,  354.2 ],
       [ 284.44,  126.64,  355.52],
       [ 285.32,  123.52,  356.  ],
       [ 286.64,  120.  ,  356.  ],
       [ 287.52,  116.  ,  356.  ],
       [ 288.44,  112.44,  356.  ],
       [ 289.76,  109.32,  356.  ],
       [ 292.  ,  106.64,  356.  ],
       [ 294.2 ,  103.52,  356.  ],
       [ 296.  ,  100.44,  356.  ],
       [ 297.32,   97.32,  356.  ],
       [ 299.08,   94.64,  356.  ],
       [ 300.88,   91.52,  356.  ],
       [ 303.08,   88.44,  356.  ],
       [ 304.88,   85.32,  356.  ],
       [ 307.08,   83.08,  356.  ],
       [ 308.88,   81.32,  356.  ],
       [ 311.08,   80.44,  356.  ],
       [ 313.32,   79.52,  356.  ],
       [ 316.44,   78.2 ,  356.  ],
       [ 319.52,   75.52,  355.52],
       [ 322.64,   72.44,  354.64],
       [ 324.88,   69.32,  352.88],
       [ 327.08,   66.64,  351.08],
       [ 330.64,   63.52,  348.  ],
       [ 336.  ,   60.  ,  344.  ],
       [ 347.6 ,   51.6 ,  311.2 ]])

In [ ]:


In [ ]:
# import sys
# print sys.path
# sys.path.append('somepath')


# list of tuples
# lt = [('r',np.ones((5,5))),('g',np.ones((8,2)))]
# print(lt[0][1])


# render volume
#vol = np.empty(img.shape + (4,), dtype='ubyte')
#vol[:] = img[..., None]
#vol = np.ascontiguousarray(vol.transpose(1, 2, 0, 3))
#vi = pgl.GLVolumeItem(vol)
#self.glView.addItem(vi)
#vi.translate(-vol.shape[0]/2., -vol.shape[1]/2., -vol.shape[2]/2.)


# remove item
# view.removeItem(mesh)


# g = pgl.GLGridItem()
# g.scale(10, 10, 1)
# view.addItem(g)


# # render volume
# vol = np.zeros(pdensity.shape + (4,), dtype='ubyte')
# vol[...,3] = pdensity*255
# vol[...,1] = 255
# #vol = np.ascontiguousarray(vol.transpose(1, 2, 0, 3))
# vi = pgl.GLVolumeItem(vol)
# vi.translate(-vol.shape[0]/8., -vol.shape[1]/8., -vol.shape[2]/8.)
# vi.scale(1./4,1./4,1./4)


# verts_pd, faces_pd = pg.isosurface(ndi.gaussian_filter(pd.astype('float32'), (2, 2, 2)), .05)
# md_d = pgl.MeshData(vertexes=verts_pd, faces=faces_pd)
# mesh_d = pgl.GLMeshItem(meshdata=md_d, smooth=True, color=[1, 0.5, 1, 0.9], shader='balloon')
# mesh_d.setGLOptions('additive')
# mesh_d.translate(-pd.shape[0]/2., -pd.shape[1]/2., -pd.shape[2]/2.)
# view.addItem(mesh_d)

In [ ]:

Define the visualization functions


In [3]:
# def vis3D(brain_array,inj_array,pad = 30,ds_factor=6):
    
#     # set up time variables
#     now = time.time()
#     now_start = now
    
#     view = vis3D_glassBrain(template,pad,ds_factor)
#     print "build brain isosurface %0.2f" % (time.time() - now); now = time.time() 
    
#     view = vis3D_projections(view,inj_array,ds_factor)
#     print "build injection volume %0.2f" % (time.time() - now); now = time.time() 
    
#     view.show()
    
#     print "rendering %0.2f" % (time.time() - now); now = time.time() 
#     print "total run time: %0.2f" % (time.time() - now_start)
    
#     return view

In [4]:
# def vis3D_glassBrain(brain_array,pad,ds_factor):
    
#     # initialize the window
#     view = pgl.GLViewWidget()   
    
#     # downsample the brain image using the ds_factor
#     img = brain_array[::ds_factor,::ds_factor,::ds_factor]
    
#     # do padding of the brain to avoid holes during rendering
#     pad_img = np.zeros((img.shape[0]+pad, img.shape[1]+pad, img.shape[2]+pad), dtype=img.dtype)
#     pad_img[pad/2:pad/2+img.shape[0], pad/2:pad/2+img.shape[1], pad/2:pad/2+img.shape[2]] = img
    
#     # build the brain isosurface
#     verts, faces = pg.isosurface(ndi.gaussian_filter(pad_img.astype('float32'), (1, 1, 1)), 5.0)
#     md = pgl.MeshData(vertexes=verts, faces=faces)
#     mesh = pgl.GLMeshItem(meshdata=md, smooth=True, color=[0.5, 0.5, 0.5, 0.1], shader='balloon')
#     mesh.setGLOptions('additive')
#     mesh.translate(-pad_img.shape[0]/2., -pad_img.shape[1]/2., -pad_img.shape[2]/2.)
#     mesh.rotate(-90, 1, 0, 0)
#     view.addItem(mesh)
#     view.setCameraPosition(distance=200, elevation=20, azimuth=90)
#     view.show()

#     return view

In [5]:
# def vis3D_projections(view,inj_array,ds_factor):
    
#     # render the injection(s) as a volume
#     # inj_array should be a list of tuples, with the first element in the tuple
#     # being the plotting color (a RGB value), and the second element being the 
#     # ND-array of the volumetric data for a given injection
#     vols = np.zeros(inj_array[0][1].shape + (4,), dtype='float32')
#     for inj in range(len(inj_array)):
#         col = inj_array[inj][0]
#         vols[...,0] += col[0] * inj_array[inj][1] # red channel
#         vols[...,1] += col[1] * inj_array[inj][1] # green channel
#         vols[...,2] += col[2] * inj_array[inj][1] # blue channel
#         vols[...,3] += inj_array[inj][1] * 255    # alpha channel

#     # Set alpha and make sure the maximum alpha is 255
#     vols[...,3] *= 5
#     vols[...,3] = np.clip(vols[...,3],0,255)

#     # now add the volume to the view window
#     vi = pgl.GLVolumeItem(vols)
#     vi.translate(-vols.shape[0]/(2.*ds_factor), -vols.shape[1]/(2.*ds_factor), -vols.shape[2]/(2.*ds_factor))
#     vi.scale(1./ds_factor,1./ds_factor,1./ds_factor)
#     vi.setGLOptions('additive')
#     vi.rotate(-90, 1, 0, 0)
#     view.setCameraPosition(distance=200, elevation=20, azimuth=90)
#     view.addItem(vi)
    
#     return view

In [6]:
# def vis3D_structureMask(view,mask,maskCol,ds_factor):

#     # downsample the brain image using the ds_factor
#     img = mask[::ds_factor,::ds_factor,::ds_factor]

#     # build the brain isosurface
#     verts, faces = pg.isosurface(ndi.gaussian_filter(img.astype('float32'), (0.5, 0.5, 0.5)), .5)
#     md = pgl.MeshData(vertexes=verts, faces=faces)
#     meshMask = pgl.GLMeshItem(meshdata=md, smooth=True, color=[maskCol[0], maskCol[1], maskCol[2], 0.2], shader='balloon')
#     meshMask.setGLOptions('additive')
#     meshMask.translate(-img.shape[0]/2., -img.shape[1]/2., -img.shape[2]/2.)
#     meshMask.rotate(-90, 1, 0, 0)
#     view.addItem(meshMask)
#     view.setCameraPosition(distance=200, elevation=20, azimuth=90)
#     view.show()

#     return view