jackman/opencv-test
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity
Files
06260589fd9ab91908cf7ea04c3eb0dc79e6a2f1
opencv-test/2016-06-15-opencv-RGBmasks-xcorr.ipynb

664 lines
135 KiB
Plaintext
Raw Normal View History

py nb
2018-11-16 16:24:21 -08:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os, sys, subprocess, datetime, time\n",
"from timeit import default_timer as timer\n",
"import cv2\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import h5py\n",
"import pandas as pd\n",
"from skimage.draw import polygon\n",
"from skimage import data\n",
"from skimage.util import random_noise\n",
"import tifffile\n",
"#from skimage.feature import match_template\n",
"#import scipy.io as spio\n",
"#from scipy import ndimage as ndi\n",
"#from scipy import signal as signal\n",
"#from joblib import Parallel, delayed\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def importFilelist(filepath):\n",
" #import text file of filenames as pandas dataframe\n",
" files = pd.read_table(filepath,names=['filename','matfilename', 'age.g'],sep=' ',header=None)\n",
" return files\n",
"\n",
"\n",
"def writeData(filepath,data):\n",
" #write results dataframe by appending to file\n",
" #filepath = 'dXcorrObs.txt'\n",
" if os.path.isfile(filepath):\n",
" writeHeader=False\n",
" else:\n",
" writeHeader=True\n",
"\n",
" with open(filepath, 'a') as f:\n",
" data.to_csv(f, index=False, sep='\\t', header=writeHeader) #make if logic to not append with header if file already exists\n",
"\n",
"\n",
"\n",
"def getA3binaryArray(f,region):\n",
" #read in domainData into A3 binary array\n",
" sz = np.int64(region['domainData']['CC/ImageSize'][:,0])\n",
" A3 = np.zeros(sz,dtype='uint8')\n",
" for i in np.arange(region['domainData']['CC/PixelIdxList'].shape[0]):\n",
" pxind = np.array(f[region['domainData']['CC/PixelIdxList'][i,0]], dtype='int64')\n",
" A3.T.flat[pxind[0,:]] = 255\n",
" return A3\n",
"\n",
"\n",
"def getHemisphereMasks(f,region,A3):\n",
" #read region['coords'] and region['name'] to get hemisphere masks\n",
" sz = A3.shape\n",
" leftMask = np.zeros((sz[0],sz[1]),dtype=bool)\n",
" rightMask = np.zeros((sz[0],sz[1]),dtype=bool)\n",
" bothMasks = np.zeros((sz[0],sz[1]),dtype=bool)\n",
" for i in np.arange(region['name'].len()): #this is length 33\n",
" currentString = f[region['name'][i,0]][...].flatten().tostring().decode('UTF-16')\n",
" if (currentString == 'cortex.L') or (currentString == 'cortex.R'):\n",
" #print(i,currentString)\n",
" x = np.array(f[region['coords'][i,0]],dtype='int64')[0,:]\n",
" y = np.array(f[region['coords'][i,0]],dtype='int64')[1,:]\n",
" x = np.append(x,x[0]) #close the polygon coords\n",
" y = np.append(y,y[0]) #close the polygon coords\n",
" rr, cc = polygon(y,x,bothMasks.shape)\n",
" if currentString == 'cortex.L':\n",
" x_L = x\n",
" y_L = y\n",
" leftMask[rr,cc] = True\n",
" elif currentString == 'cortex.R':\n",
" x_R = x\n",
" y_R = y\n",
" rightMask[rr,cc] = True\n",
" bothMasks[rr,cc] = True\n",
" return leftMask, rightMask\n",
"\n",
"\n",
"\n",
"def getHemisphereCoords(f,region):\n",
" for i in np.arange(region['name'].len()): #this is length 33\n",
" currentString = f[region['name'][i,0]][...].flatten().tostring().decode('UTF-16')\n",
" if (currentString == 'cortex.L') or (currentString == 'cortex.R'):\n",
" #print(i,currentString)\n",
" x = np.array(f[region['coords'][i,0]],dtype='int32')[0,:]\n",
" y = np.array(f[region['coords'][i,0]],dtype='int32')[1,:]\n",
" x = np.append(x,x[0]) #close the polygon coords\n",
" y = np.append(y,y[0]) #close the polygon coords\n",
" if currentString == 'cortex.L':\n",
" xyPtsL = np.stack((x,y),axis=-1)\n",
" elif currentString == 'cortex.R':\n",
" xyPtsR = np.stack((x,y),axis=-1)\n",
" return xyPtsL, xyPtsR\n",
"\n",
"\n",
"def stridemask3d(mask,z):\n",
" #make a 2D image mask into a 3D repeated image mask using numpy stride_tricks to fake the 3rd dimension\n",
" mask3d = np.lib.stride_tricks.as_strided(\n",
" mask, # input array\n",
" (mask.shape[0], mask.shape[1],z), # output dimensions\n",
" (mask.strides[0], mask.strides[1],0) # stride length in bytes\n",
" )\n",
" return mask3d\n",
"\n",
"\n",
"def cvblur(image,template):\n",
" #fast gaussian blur with cv2\n",
" image = cv2.GaussianBlur(image,(0,0),3)\n",
" template = cv2.GaussianBlur(template,(0,0),3)\n",
" return image, template\n",
"\n",
"\n",
"def cvcorr(image,template):\n",
" #faster xcorr with cv2\n",
" image = image.astype('float32')\n",
" template = template.astype('float32') #cv2.filter2D needs a float as input\n",
" c = cv2.filter2D(image,-1,template)\n",
" corrValue = c[c.shape[0]/2,c.shape[1]/2]\n",
" return corrValue, c\n",
"\n",
"\n",
"def getA3hemipshereArrays(A3,leftMask,rightMask):\n",
" #calculate lateral column shift for a flipped right hemisphere for a common coord system\n",
" rightMask_flip = np.fliplr(rightMask)\n",
" IND_L = np.nonzero(leftMask) #dim2 of the tuple is the col indices\n",
" IND_Rflip = np.nonzero(rightMask_flip) #dim2 of the tuple is the col indices\n",
" dx = np.amax(IND_L[1]) - np.amax(IND_Rflip[1]) #number of pixesl to shift the flipped right hemisphere column indices by\n",
"\n",
" #get 3d mask arrays, new binary data arrays, and a flipped R hemi array\n",
" leftMask3d = stridemask3d(leftMask,A3.shape[2])\n",
" rightMask3d = stridemask3d(rightMask,A3.shape[2])\n",
"\n",
" A3left = np.logical_and(A3, leftMask3d)\n",
" A3left = A3left.view(np.uint8) #cast to uint8 inplace and multiply by 255 inplace\n",
" A3left *= 255\n",
" #np.place(A3left, A3left>0, 255)\n",
" \n",
" A3right = np.logical_and(A3, rightMask3d)\n",
" A3right = np.fliplr(A3right)\n",
" \n",
" #make a new shifted right hemisphere array\n",
" IND_R3d = np.nonzero(A3right) #dim2 of the tuple is the col indices\n",
" INDnew = IND_R3d[1] + dx #add shift to column indices\n",
" #A3right_shift = np.zeros_like(A3)\n",
" #A3right = np.zeros_like(A3)\n",
" A3right = A3right.view(np.uint8) #cast to uint8 inplace and multiply by 0 inplace\n",
" #np.putmask(A3right, A3right>0, A3right*0)\n",
" #np.place(A3right, A3right>0, 0)\n",
" A3right *= 0\n",
" A3right[IND_R3d[0],INDnew,IND_R3d[2]] = 255\n",
" \n",
" return A3left, A3right\n",
"\n",
"def getCorr(A3left,A3right,i):\n",
" image = A3left[:,:,i]\n",
" template = A3right[:,:,i]\n",
" image,template = cvblur(image,template)\n",
" autoCorr = cvcorr(image,image)\n",
" obsCorr = cvcorr(image,template)\n",
" fr = i + 1\n",
" return fr, autoCorr, obsCorr\n",
"\n",
"\n",
"def cvmatch(image,template):\n",
" #faster xcorr with cv2\n",
" c = cv2.matchTemplate(image,template,cv2.TM_CCORR_NORMED)\n",
" corrValue = c[c.shape[0]/2,c.shape[1]/2]\n",
" return corrValue\n",
"\n",
"def cvmatchN(image,template):\n",
" #faster xcorr with cv2\n",
" c = cv2.matchTemplate(image,template,cv2.TM_CCOEFF_NORMED)\n",
" corrValue = c[c.shape[0]/2,c.shape[1]/2]\n",
" return corrValue\n",
"\n",
"\n",
"def playMovie(A3,newMinMax=False):\n",
" #play movie in opencv\n",
" cv2.startWindowThread()\n",
" cv2.namedWindow(\"raw\", cv2.WINDOW_NORMAL) #Create a resizable window\n",
" i = 0\n",
" \n",
" if newMinMax:\n",
" #im = cv2.normalize(im, None, newMinMax[0], newMinMax[1], cv2.NORM_MINMAX)\n",
" #im = cv2.normalize(im, None, newMinMax[0], newMinMax[1], cv2.NORM_L2)\n",
" #im = cv2.equalizeHist(im)\n",
" #im,cdf = histeq(im)\n",
" cv2.normalize(A3,A3,newMinMax[0],newMinMax[1])\n",
" \n",
" while True:\n",
" #im = np.uint8(A3[:,:,i] * 255)\n",
" im = A3[:,:,i]\n",
" im = cv2.applyColorMap(im, cv2.COLORMAP_HOT)\n",
" im = cv2.GaussianBlur(im,(0,0),3)\n",
" #th,bw = cv2.threshold(im,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)\n",
" #bw = cv2.adaptiveThreshold(im,255,cv2.ADAPTIVE_THRESH_MEAN_C,cv2.THRESH_BINARY,5,0)\n",
" cv2.imshow('raw',im)\n",
" k = cv2.waitKey(10) \n",
" if k == 27: #if esc is pressed\n",
" break\n",
" if k == ord('b'):\n",
" i -= 1000\n",
" elif k == ord('f'):\n",
" i += 1000\n",
" else:\n",
" i += 1\n",
" if (i > (A3.shape[2]-1)) or (i < 0) :\n",
" i = 0\n",
"\n",
" cv2.destroyAllWindows()\n",
"\n",
"\n",
" \n",
" \n",
"def histeq(im, nbr_bins=256):\n",
" \"\"\"histogram equalization of a grayscale image\"\"\"\n",
" #get image histo\n",
" imhist,bins = np.histogram(im.flatten(),nbr_bins,normed=True)\n",
" cdf = imhist.cumsum() #cumulative distribution function\n",
" cdf = 255 * cdf / cdf[-1] #normalize\n",
" #use linear interp of cdf to find new pixel values\n",
" im2 = np.interp(im.flatten(),bins[:-1],cdf)\n",
" return im2.reshape(im.shape), cdf\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def playMovieRGB(A3,pts=False,filename=False):\n",
" #play movie in opencv\n",
" #input is a uint8 numpy rgb array with A3.ndim = 4 and shape (x,y,z,3)\n",
" #pts should be a list of numpy.ndarrays 'list((xyPts1,xyPts2))', since cv2.polylines() can be used to draw multiple lines if a list of arrays is passed to it.\n",
" cv2.startWindowThread()\n",
" cv2.namedWindow(\"raw\", cv2.WINDOW_NORMAL) #Create a resizable window\n",
" i = 0\n",
" toggleNext = True\n",
" tf = True\n",
" #if isinstance(pts, (list, np.ndarray)):\n",
" if isinstance(pts, (list)):\n",
" drawCoords = True\n",
" else:\n",
" drawCoords = False\n",
" \n",
" if filename == False:\n",
" import datetime\n",
" filename = datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d-%H%M%S')\n",
" \n",
" while True:\n",
" im = A3[:,:,i,:]\n",
" img = cv2.cvtColor(im, cv2.COLOR_RGB2BGR) #convert numpy RGB input to BGR colorspace for cv2\n",
" img = cv2.GaussianBlur(img,(0,0),3)\n",
" corrValue = cvmatch(img[:,:,2],img[:,:,1])\n",
" cv2.putText(img, str(i), (5,25), cv2.FONT_HERSHEY_SIMPLEX, 1.0, (155,155,155)) #draw frame text\n",
" cv2.putText(img, '{0:.4f}'.format(corrValue), (img.shape[1]-150,25), cv2.FONT_HERSHEY_SIMPLEX, 1.0, (155,155,155)) #draw frame text\n",
" if drawCoords:\n",
" #cv2.polylines(img, pts, False,(255, 255, 255), thickness=1, lineType=cv2.CV_AA) #lineType for OpenCV 2.4.13\n",
" cv2.polylines(img, pts, False,(255, 255, 255), thickness=1, lineType=cv2.LINE_AA) #lineType for OpenCV 3.x\n",
" cv2.imshow('raw',img)\n",
" k = cv2.waitKey(10) \n",
" if k == 27: #if esc is pressed\n",
" break\n",
" elif (k == ord(' ')) and (toggleNext == True): #if space is pressed\n",
" tf = False\n",
" elif (k == ord(' ')) and (toggleNext == False): #if space is pressed\n",
" tf = True\n",
" toggleNext = tf #toggle the switch\n",
" if k == ord('b') and toggleNext:\n",
" i -= 100\n",
" elif k == ord('f') and toggleNext:\n",
" i += 100\n",
" elif k == ord('>') and (toggleNext == False):\n",
" i += 1\n",
" elif k == ord('<') and (toggleNext == False):\n",
" i -= 1\n",
" elif k == ord('s') and (toggleNext == False):\n",
" #params = list()\n",
" #params.append(cv.CV_IMWRITE_PNG_COMPRESSION)\n",
" #params.append(8)\n",
" #cv2.imwrite('image{0}.png'.format(i), img, params)\n",
" cv2.imwrite(filename + '-fr' + str(i) + '.png', img)\n",
" elif toggleNext:\n",
" i += 1\n",
" \n",
" if (i > (A3.shape[2]-1)) or (i < 0) :\n",
" i = 0\n",
"\n",
" cv2.destroyAllWindows()\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.7707633420359343\n",
"0.030280801933258772\n",
"10.70011685800273\n",
"2.7615016559138894\n",
"(540, 640, 3000, 3)\n"
]
}
],
"source": [
"#matfilename = '/Users/ackman/Data/2photon/120518i/2015-03-13-114620_sigma3_thresh1/120518_07_20150321-153601_d2r.mat'\n",
"#matfilename = '/Users/ackman/Data/2photon/140219/2015-03-13-114620_sigma3_thresh1/140219_07_20150313-122433_d2r.mat'\n",
"#matfilename = '/Users/ackman/Data/2photon/140328/2015-03-13-114620_sigma3_thresh1/140328_05_20150313-125415_d2r.mat'\n",
"matfilename = '/Users/ackman/Data/2photon/140331/2015-03-13-114620_sigma3_thresh1/140331_05_20150314-090750_d2r.mat'\n",
"#matfilename = '/Users/ackman/Data/2photon/140509/2015-03-13-114620_sigma3_thresh1/140509_22_20150313-193050_d2r.mat'\n",
"#matfilename = '/Users/ackman/Data/2photon/141125/2015-03-13-114620_sigma3_thresh1/141125_05_20150313-140653_d2r.mat'\n",
"\n",
"f = h5py.File(matfilename,'r') \n",
"region = f.get('region')\n",
"t0=timer()\n",
"A3 = getA3binaryArray(f,region); print(timer()-t0)\n",
"t0=timer()\n",
"leftMask, rightMask = getHemisphereMasks(f,region,A3); print(timer()-t0)\n",
"t0=timer()\n",
"A3left, A3right_shift = getA3hemipshereArrays(A3,leftMask,rightMask); print(timer()-t0)\n",
"\n",
"t0 = timer()\n",
"A4 = np.zeros((A3left.shape[0],A3left.shape[1],A3left.shape[2],3),dtype='uint8')\n",
"A4[:,:,:,0] = A3left\n",
"A4[:,:,:,1] = A3right_shift\n",
"print(timer() - t0)\n",
"print(A4.shape)\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fn = os.path.splitext(os.path.basename(matfilename))[0]\n",
"xyPtsL,xyPtsR = getHemisphereCoords(f,region)\n",
"#pts = np.stack((xyPtsL,xyPtsR),axis=2)\n",
"playMovieRGB(A4,list((xyPtsL,xyPtsR)),fn)"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Variables already gone...\n"
]
}
],
"source": [
"try:\n",
" del(region, A3, leftMask, rightMask, A3left, A3right_shift, A4, images, image)\n",
"except:\n",
" print('Variables already gone...')\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/ackman/anaconda/envs/py27/lib/python2.7/site-packages/tifffile/tifffile.py:1974: UserWarning: tags are not ordered by code\n",
" warnings.warn(\"tags are not ordered by code\")\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"13.4423279762\n",
"(540, 640, 3000)\n",
"float64\n",
"(345600, 3000)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAS4AAAEACAYAAAAN5psFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXm4pldZ5vt7v/3teazaVZWqSiWpSiqVgUwkkABGTQgQ\nZoIoArYKtLaibZ/LoWmlz5G0esTrnO62bYduFYTuVpoDiMgoM4QEDSEDCUkqY83zuIfa8/7e88da\nv72eb5MKtsecgsv9Xte+9t7f945rPet+7ud+nrXeqq5rVraVbWVb2b6XtsaZvoGVbWVb2Va2/9Vt\nBbhWtpVtZfue21aAa2Vb2Va277ltBbhWtpVtZfue21aAa2Vb2Va277ltBbhWtpVtZfue254x4Kqq\n6qVVVW2vqurRqqr+zTN1nZVtZVvZ/ult1TNRx1VVVQN4FLgJ2A/cBbyhruvt/+gXW9lWtpXtn9z2\nTDGua4HH6rreVdf1PPAB4DXP0LVWtpVtZfsntj1TwHU2sCf8vzd/trKtbCvbyvb/eVsR51e2lW1l\n+57bms/QefcB54b/N+XPlraqqlYmSa5sK9vK9rRbXdfVU33+TAHXXcDWqqrOAw4AbwDe+O27vZCk\n33fm//3dAVQkQriYP+/Mn3cCXfnWG0ALGAL6gG5gIP/dAYzk83RDoyq7D4Rb8JJz+fC58F0DmLkV\n+m9Np+nN37fyd3X+6cm/5/Jl63B8F7CQP6vy74V83Sqcw8/78t+z+W/yfdW5KTzvfDj+6K0wemu6\nL5uqzr9n8+8FyvcVMJN/+xm5SWfy+bvz97PhOebzfvO5LU6FLjgfmADGgHNIvX4+8Gro3TbG9B8M\nw1dvhb5b03NNhrbzuZr5ej35fjtC2/k9pO722qcopjGd72shf1flY1r5Zzafq5Wfp5W/nwt9MnYr\nrLq19NdcOP9Uvv5ibp/ZfB7bzOvaR7P5fup8bw2K3ZCv35F/9+frVbTbS0e4b03+xK3Qc2t5xvn8\n+Vzex89iG3qdGcrwWqT9HLaVNux9Lubv7P+lzQ8mKEY5kW9kGvgvwJvzRefziaZDwy5QDHkx/O/2\nC8svuLQ9I8BV1/ViVVX/Evgs6YneU9f1w9++Z0f+sQVblBHj7QlW9prA5Yjrzb8dXb0ky14M58qg\n1Zt/PDWkdtRo/Q3FMDXIQUond+bjhsIxrfDTnz+H1Gferv3TnX932GCkfuwKTeOt1/kcYvFi3tcm\n0tA0cAHGQdJFsaO+/LubBA4n8jPY9Iv5py/ffzPv5/nm8/1E4FwA1uVrT+e/j5AA5SjwBzC9MFza\nsz9fo5cyuJqh/bvz530koHAQdS27h1mKaQhGgn50MF5rkgIEMxRg8Dx9+T4mKP2rg/K+tRnb1mvr\nWz23+/TkY7tDP8yRHGeVr9FBMe35fMw4sCrcv6BD3qcjn2OCYo91vo77NfIzCezihe3uNcn3CAXo\ndcCex/7RDlp12BkSQRA158MNeZGe/NkkxXNGLzxFO1LrUU+/PVOMi7qu/wa46Dvv2UG6eS3UFvPh\nuigjs5eCAo64nvAjwGnFjXJoF2VARrajwWpw8/lw+8H2FxQ6wukdeHokmYKA16IMfAduFc6hV5vP\n546srcm3MyzP05v3maTYisAs9ktWBduhcG7vbTTfUyt8ZvPZFbKIKQrDicAsMT5AGkxzlME/Rhmo\nHi8Q2Beyv57QD03gZGh/QXme5EDigF8+CHXYXrcO+wv4+jbP6cC2L6MjE5R0EIKd/bKQ7528r+xb\ngFmgOCRtxPbuC/d5igJuPuNTbfaTbO5U/j1FYVtQbEnGp9PTufqZ7RaxIjIuQcu2awCtfJGqKqxs\nyWOuIVVA2aEOihmSYUqN9Tp69hiueKHTb88YcP39tgsoIznCv3RDJtUZPtOtxDiqj2SJA6RR1SgD\ntofUebKcAcpgtI0mSW1nWOVAaQG9NxQWJkDFjtebGoqJu/EYQXKRMugFMkFGA4iAajgUvxfEZGE1\nMHRDCQdkQVBAOzovDVvQlUF5bvdfSwKPGBLKxqCE1oMkZjUU2rCH0lUyhy6AG9LxhkQj+Rr9+ecE\nBfTdp4cyKG1T2V6T9kijh8KUoPTvIgUc9XWT+UfQcWD335DOO0ZxeDIpwcJ94z1EBzeQjxungJPX\ntm90jB0koO/O3w1QxrXhvOC5mM/fd0P63pBdx2Xb++zigazNe43yRYt2O1CyaFLQYXHZ327d+R6X\nOqs377A2P8S1JCA7QjKUcYoxOAgXKCHNTD4+Epin3p6RAtS/z5bE+d+l3KQtJx2KTEuLHch/95As\nHWCY9NCyrQHoqNLhPRRw8dR6HNlsDH9q0iDRSDzWcCLuJxnUM8sWekjGIkD5eNMUz+o5vR/BSdan\nEcvCBkkDyf0MK/X22sB0vo8uykD2O59TcG3SHprpyW1qAUTfMUEZ8DIIZYs+kn2eIIFR1I/WUEBd\nBipb7MnPVZFseoLCWCMoOIgcxDKvqDDImuZJYHmKAsLzFPCwbQUOAdi+89nVtqD0+6l83ciayfdx\ngtTnneHc9r9MWdtarmF1hM+jftrMz6Q9Qnv4pyZqSG14vNz+Zil97bECfnR2Eehkz7abhKmZzxHb\nbhZo1Lk9BJ7xfGMaySLJS6l1GNqcoqCl+/r7bf+/i/P/C1tUUaG4Cv+XWw/kzxRo5Nr+5O+qqgze\nrnB6T+PfCp6xM6POpEHNhXNo9HaywKNX7cv/91JAZIDUhzK+GP67qWMJEg56xfVTeT89n2K9xgvF\nwAwRvGfDA8OJmIxQZ1EHMjztyvfdoj1EiHqLhLhJkRTXUrSgJgnEIlm2jWpKKLgm39cqUt5ZMKgo\noNbK9xmZnIzYwSvj7KQ4B7UsQVxmebrEhICmXzR87CcB7QDtjND76g3tFQX6KHTLrjXvyBQrEmiT\n22E87APF6XgfftdJsS3D7g5KqGlfRmero+qiJKMim4rJGyhjQIYXh6rnjeHjkoEY/ZykCKDeqAdO\nkjpvluKFNOzv6lDRlExUWLXOmGJRw4JkJX0Uzu7f2Zr800FVUbyPIKOAvkgaHFJ09xH47PwITjIo\nBWAHjcKqLCkOYCjRsOCoHjVHYR+ea44SCgqsA5RMnIMsJmXUgxyIamuGtmozkc1FjWaCBECydUK3\nqKPExILZOwG+P99vF7CeogEqewzmZ2iQSpO7gNWkgWdfyEivIKkI3wTuILEZwT2qBAJQBCSVBAe7\n2pGAMUVhqwv5+9UUBqSGLMC2KFqiIGYSbTHsR+gbdbIu2s1aMiLBEGAiEE1SQEKQkTlpj2a+jQTG\nKJJS1PHGw/W1A8V2P4cSZhq6+gzaTZSevHeBezGcB0KHRMHOhos3XYUTKYBqcFG/eOrtDAOX1qTr\ns+XscUUUR4f/mzseKd9VJCNycHgaL6NICaW9mpRsnaFB7CgBxjAw0nqdip8PhWMXw+95ShZpMe8H\nBa/FXkFNwxJg1FUMjwQrgSRKhAuUcElvb/MKDNpEbOZOUsRthC7gyiC8vr5kkOIvZJ4S31FSZpF8\nnXGKuK/tmt2dJTEMQ5ADJDmkmwQOE/m+JkN7O1AU8e0L2zmG8z67bev928+zlIqZ5an/mNFzW54Y\ni2F2FMXjQBdEom3I6nQEjmvv0RDsFIUp2texfAQKmzYk7aUkR2TfstCYnBL8fHaBT/1WGxCEY3Ki\nKxzr+SKg1X7pF6tzw9mgxpp6QTUNvbodevrtDANXdFsG6BGkuii9NUxxXatYKouwMTspoqhGIOOK\nmo56VIv2zrC2SCakYGrY5zkEOY1AYFOzMrsnePTlc1mKoJalGD2c/zbUF4Q8Xz+F/Tm4ZDNQUv1G\ny9FeBB7ZlTqYGpjMUx0lDj7vP2YtDftGSQAs4FG6g7WkXPJ0bqetpCTT0byv5zlBCaXV8Sby/4dI\nYeNUvk5NYhWWA9nGU5TwK2ZSZVZVuJ5+0bIAQUcR3P6N+qNhYSPso2OMLNR2c1ALqjFk15T9zu/7\nwnO1wr6Tub2WA2NkSTJ
"text/plain": [
"<matplotlib.figure.Figure at 0x11e3af810>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t0 = timer()\n",
"#fn = '/Users/ackman/Desktop/140509_21.tif'\n",
"fn = '/Users/ackman/Data/2photon/120518i/120518_07_fr1-300.tif'\n",
"with tifffile.TiffFile(fn) as tif:\n",
" A = tif.asarray()\n",
"del(tif)\n",
"A = np.rollaxis(A, 0, 3)\n",
"A = np.float64(A)\n",
"print(timer()-t0)\n",
"print(A.shape)\n",
"print(A.dtype)\n",
"sz = A.shape\n",
"\n",
"plt.imshow(A[:,:,0])\n",
"A = np.reshape(A, (A.size/sz[2], sz[2]))\n",
"print(A.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This takes just 0.32 s to read in array at its ** default uint16.**\n",
"And takes 0.78 s to read in array and cast as float64.\n",
"\n",
"Matlab takes 2.11 s to read in same array into a float64.\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"13.2841491699\n"
]
}
],
"source": [
"t0=timer()\n",
"Amean = np.mean(A,axis=1)\n",
"Amean2D = np.lib.stride_tricks.as_strided(Amean, (Amean.shape[0],A.shape[1]),(Amean.strides[0],0))\n",
"A /= Amean2D\n",
"A -= 1\n",
"print(timer()-t0)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"A += abs(np.amin(A))\n",
"A /= np.amax(A)\n",
"A = np.reshape(A,sz)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Without stride tricks, the above F/F0-1 sequence was slower than matlab (2s vs 1s). With stride tricks we get 0.44 s, a big improvement. "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"55.4744589329\n",
"22.1628820896\n",
"75.272135973\n"
]
}
],
"source": [
"t0=timer()\n",
"meanA = np.mean(A) \n",
"stdA = np.std(A)\n",
"newMin = meanA - 2*stdA\n",
"newMax = meanA + 5*stdA\n",
"#playMovie(A2,(newMin,newMax))\n",
"print(timer()-t0)\n",
"\n",
"t0=timer()\n",
"A -= newMin\n",
"A *= 255.0/(newMax-newMin)\n",
"print(timer()-t0)\n",
"\n",
"t0=timer()\n",
"#np.clip(A, 0, 255, out=A)\n",
"A.clip(0, 255) #faster than the documented np.clip(A,min,max,out=A) which should be inplace. Possible version/memory leak issue\n",
"print(timer()-t0)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A2[logical]: 1.28217601776\n",
"\n",
"np.clip: 1.14894986153\n"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"playMovie(np.uint8(A))"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"del(A,Amean,Amean2D)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Variable Type Data/Info\n",
"---------------------------------------------------------------\n",
"cv2 module <module 'cv2' from '/User<...>.7/site-packages/cv2.so'>\n",
"cvblur function <function cvblur at 0x11e318488>\n",
"cvcorr function <function cvcorr at 0x11e318de8>\n",
"cvmatch function <function cvmatch at 0x11e318140>\n",
"cvmatchN function <function cvmatchN at 0x11e318050>\n",
"data module <module 'skimage.data' fr<...>image/data/__init__.pyc'>\n",
"datetime module <module 'datetime' from '<...>lib-dynload/datetime.so'>\n",
"f File <HDF5 file \"140331_05_201<...>090750_d2r.mat\" (mode r)>\n",
"fn str 140331_05_20150314-090750_d2r\n",
"getA3binaryArray function <function getA3binaryArray at 0x11e21f578>\n",
"getA3hemipshereArrays function <function getA3hemipshereArrays at 0x11e318f50>\n",
"getCorr function <function getCorr at 0x11e318230>\n",
"getHemisphereCoords function <function getHemisphereCoords at 0x11e318500>\n",
"getHemisphereMasks function <function getHemisphereMasks at 0x11e21f0c8>\n",
"h5py module <module 'h5py' from '/Use<...>kages/h5py/__init__.pyc'>\n",
"histeq function <function histeq at 0x11b996ed8>\n",
"importFilelist function <function importFilelist at 0x11e21fc80>\n",
"matfilename str /Users/ackman/Data/2photo<...>5_20150314-090750_d2r.mat\n",
"meanA float64 0.133782654165\n",
"newMax float64 0.226439697864\n",
"newMin float64 0.0967198366862\n",
"np module <module 'numpy' from '/Us<...>ages/numpy/__init__.pyc'>\n",
"os module <module 'os' from '/Users<...>27/lib/python2.7/os.pyc'>\n",
"pd module <module 'pandas' from '/U<...>ges/pandas/__init__.pyc'>\n",
"playMovie function <function playMovie at 0x11b996f50>\n",
"playMovieRGB function <function playMovieRGB at 0x11b996e60>\n",
"plt module <module 'matplotlib.pyplo<...>s/matplotlib/pyplot.pyc'>\n",
"polygon builtin_function_or_method <built-in function polygon>\n",
"pts ndarray 20x2x2: 80 elems, type `int32`, 320 bytes\n",
"random_noise function <function random_noise at 0x11643d7d0>\n",
"stdA float64 0.0185314087397\n",
"stridemask3d function <function stridemask3d at 0x11e318ed8>\n",
"subprocess module <module 'subprocess' from<...>ython2.7/subprocess.pyc'>\n",
"sys module <module 'sys' (built-in)>\n",
"sz tuple n=3\n",
"t0 float 1466189975.43\n",
"tifffile module <module 'tifffile' from '<...>s/tifffile/__init__.pyc'>\n",
"time module <module 'time' from '/Use<...>2.7/lib-dynload/time.so'>\n",
"timer builtin_function_or_method <built-in function time>\n",
"writeData function <function writeData at 0x11e21f398>\n",
"xyPtsL ndarray 20x2: 40 elems, type `int32`, 160 bytes\n",
"xyPtsR ndarray 20x2: 40 elems, type `int32`, 160 bytes\n"
]
}
],
"source": [
"%whos"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"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.5.1"
}
},
"nbformat": 4,
"nbformat_minor": 0
}
Reference in New Issue Copy Permalink