Skip to content

Instantly share code, notes, and snippets.

@wbwakeman

wbwakeman/gist:b8719e0ad5b1b6486482974f410ab086

Last active January 13, 2022 17:45
Show Gist options
  • Save wbwakeman/b8719e0ad5b1b6486482974f410ab086 to your computer and use it in GitHub Desktop.
Save wbwakeman/b8719e0ad5b1b6486482974f410ab086 to your computer and use it in GitHub Desktop.
Download ZIP
sine_dewarp_grid
Raw
gistfile1.txt
import sys
import os
import h5py
import scipy.ndimage.measurements as measurements
import numpy as np
#import Queue
import math
from PIL import Image
from pylab import *
#Get stack.h5
#Generate 8-bit grayscale maximum intensity projection image, maybe with this code:
# import imageio
# import h5py
# with h5py.File('/path/input_image_XYTZ.h5', 'r') as f:
# x = f['data'][:]
# x = (x - x.min()) / (x.max() - x.min())
# x *= 255
# x = x.astype('uint8')
# max_proj = x.max(axis=0)
# imageio.imsave('/path/mip_image.png', max_proj)
# print("Next step (rig 3): /shared/bioapps/infoapps/lims2_modules/lib/python/run_python.sh ../sine_dewarp_gridtest_last.py /path/mip_image.png 3 /path/output_image.png ")
#View 8bit.png warped image and use that as input to:
#/allen/aibs/technology/waynew/ophys/sine_dewarp_gridtest_last.py gridimagein 2pnumber gridimageout
#output is the four numbers, and an image of the dewarped input
#If image looks good, done; else manual tweak
#Record history in that file.
# Add these to solve no X server problem:
# do this before importing pylab or pyplot
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
#2p3 history
'''
date Txx Txy Tyx Tyy
date aL aR bL bR
08/12/2016 160.0 112.0 110.0 110.0
10/25/2018 160.0 160.0 110.0 110.0
01/29/2019 160.0 160.0 85.0 90.0
01/13/2022 160.0 160.0 100.0 100.0 <-try this setting
'''
#before 10/25/2018
'''
EQUIPMENT_PARAMETER_FILES = {
"CAM2P.3": {
"aL": "160.0",
"aR": "112.0",
"bL": "110.0",
"bR": "110.0"
} ,
"CAM2P.4": {
"aL": "160.0",
"aR": "160.0",
"bL": "95.0",
"bR": "130.0"
} ,
"CAM2P.5": {
"aL": "160.0",
"aR": "160.0",
"bL": "80.0",
"bR": "150.0"
}
}
'''
#after 10/25/2018
#after 01/29/2019 this is not used
EQUIPMENT_PARAMETER_FILES = {
"CAM2P.3": {
"aL": "160.0",
"aR": "160.0",
#"bL": "110.0",
#"bR": "110.0"
"bL": "85.0",
"bR": "100.0"
} ,
"CAM2P.4": {
"aL": "160.0",
"aR": "150.0",
"bL": "85.0",
"bR": "100.0"
} ,
"CAM2P.5": {
"aL": "160.0",
"aR": "160.0",
"bL": "80.0",
"bR": "120.0"
}
}
x = np.zeros( 256, np.int )
#xindex = np.zeros( 256, np.int )
xindexL = np.zeros( 256, np.int )
xindexLB = np.zeros( 256, np.float ) # between pixels
xindexR = np.zeros( 256, np.int )
xindexRB = np.zeros( 256, np.float ) # between pixels
modevalue = 0.0
bgfactorL = 0.0
bgfactorR = 0.0
L = 0
R = 0
maxlevel = 1.0
# default 2p4 parameters
aL = 160.0
aR = 160.0
bL = 95.0
bR = 110.0
def set_params(system):
global aL
global aR
global bL
global bR
if system == 3:
aL = 160.0
aR = 160.0
#bL = 110.0 #10252018
#bR = 110.0 #10252018
#bL = 85.0 #01292019
#bR = 90.0 #01292019
bL = 100.0 #01132022
bR = 100.0 #01132022
elif system == 4:
aL = 160.0
aR = 150.0
bL = 85.0
bR = 100.0
elif system == 5:
#aL = 160.0 # before 2/8/2019
aL = 150.0
#aR = 160.0 # before 2/8/2019
aR = 150.0
#bL = 80.0 # before 2/8/2019
bL = 85.0
#bR = 120.0 # before 2/8/2019
bR = 120.0
else:
aL = 160.0
aR = 160.0
bL = 80.0
bR = 110.0
def mode(data):
counts={}
for x in data.flatten():
counts[x]=counts.get(x,0) +1
maxcount = max(counts.values())
modelist=[]
for x in counts:
if counts[x]==maxcount:
modelist.append(x)
return modelist, maxcount
def create_xtable() :
#assume flat area 160 - 382
#Left side
for j in range (0,int(aL)):
x[j] = j
#xindexL[j] = j - int((80.0) * (1.0 - math.sin((j/(160.0 * 3.0) + 1.0/6.0) * 3.14159265) ) +0.5) # ok
xindexL[j] = j - int((bL) * (1.0 - math.sin((j/(aL * 3.0) + 1.0/6.0) * 3.14159265) ) +0.5) # ok
#xindexLB[j] = (j+0.5) - ((80.0) * (1.0 - math.sin(((j+0.5)/(160.0 * 3.0) + 1.0/6.0) * 3.14159265) ) ) # ok
xindexLB[j] = (j+0.5) - ((bL) * (1.0 - math.sin(((j+0.5)/(aL * 3.0) + 1.0/6.0) * 3.14159265) ) ) # ok
#Right side
for j in range (0,int(aR)):
x[j] = j
xindexR[j] = j - int((bR) * (1.0 - math.sin((j/(aR * 3.0) + 1.0/6.0) * 3.14159265) ) +0.5) # default
xindexRB[j] = (j+0.5) - ((bR) * (1.0 - math.sin(((j+0.5)/(aR * 3.0) + 1.0/6.0) * 3.14159265) ) ) # default
print (256-j, j, xindexL[j], xindexLB[j], xindexR[j], xindexRB[j])
# compute the blobal variables
global modevalue
global bgfactorL
global bgfactorR
global L
global R
modelist, maxcount = mode(img)
modevalue = modelist[0]
j = 0
while xindexLB[j] < 0.0 or xindexL[j] < 0:
j = j + 1
else:
bgfactorL = xindexLB[j+1] - xindexLB[j]
L=j
j = 0
while xindexRB[j] < 0.0 or xindexR[j] < 0:
j = j + 1
else:
bgfactorR = xindexRB[j+1] - xindexRB[j]
R=j
print ('modevalue=', modevalue)
print ('bgfactorL=', bgfactorL)
print ('bgfactorR=', bgfactorR)
print ('L R ', L, R)
def xdewarp(img, imgout) :
for i in range(0,1):
#Center
#for j in range(140,400): #near flat area
# imgout[i,j] = img[i,j]
imgout[:,(int(aL)):(512-(int(aR)))] = img[:,(int(aL)):(512-(int(aR)))]
sum=np.zeros(512, np.float)
#Left side
for j in range(0,int(aL)):
sum[:] = 0.0 # reset
if xindexL[j] >= 0 :
if xindexLB[j-1] >= 0.0 :
s = int(math.floor(xindexLB[j-1]))
e = int(math.floor(xindexLB[j]))
#print j, s, xindexLB[j-1], (s+1 - xindexLB[j-1]), e, xindexLB[j], (xindexLB[j] - e), img[i,xindexL[j]], sum
sum[:] = (s+1 - xindexLB[j-1]) * img[:,s] + (xindexLB[j] - e) * img[:,e]
if (e-s) > 1 : # have a middle pixel
sum[:] = sum[:] + img[:, s+1]
mask = (sum > maxlevel) # saturated? for max image==1.0
sum[mask] = maxlevel
imgout[:,j] = sum
#if i==250:
# print j, s, xindexLB[j-1], (s+1 - xindexLB[j-1]), e, xindexLB[j], (xindexLB[j] - e), img[i,xindexL[j]], sum
else:
imgout[:,j] = img[:,xindexL[j]] * bgfactorL
else:
imgout[:,j] = modevalue * bgfactorL
# Right side
for j in range(0,(int(aR))):
sum[:] = 0.0
if xindexR[j] >= 0 :
if xindexRB[j-1] >= 0.0 :
s = int(math.floor(xindexRB[j-1]))
e = int(math.floor(xindexRB[j]))
sum[:] = (s+1 - xindexRB[j-1]) * img[:, 511-s] + (xindexRB[j] - e) * img[:,511-e]
if (e-s) > 1 : # have a middle pixel
sum[:] = sum[:] + img[:, 511-(s+1)]
mask = (sum > maxlevel) # saturated? for max image==1.0
sum[mask] = maxlevel
imgout[:,511-j] = sum
#if i==250:
# print j, s, xindexRB[j-1], (s+1 - xindexRB[j-1]), e, xindexRB[j], (xindexRB[j] - e), img[i,511-xindexR[j]], sum
else:
imgout[:,511-j] = img[:,511-xindexR[j]] * bgfactorR
else:
imgout[:,511-j] = modevalue * bgfactorR
if __name__ == '__main__':
'''
#For testing experiment movie
if len(sys.argv) < 7:
print ("Reguire 6 arguments: input_dir input_file_name output_file_name output_file_name system(3/4/5) FOVwidth")
exit(1)
input_dir = sys.argv[1]
input_file_name = sys.argv[2]
output_dir = sys.argv[3]
output_file_name = sys.argv[4]
system = sys.argv[5]
print ("2p." +sys.argv[5])
FOVwidth = int(sys.argv[6])
print FOVwidth
if FOVwidth == 512:
print ('output fixed width FOV')
else:
print ('output variable width FOV')
#print "input_dir: " + input_dir
#print "output_h5_file: " + output_file_name
'''
'''
data_file = h5py.File(os.path.join(input_dir, input_file_name), "r")
stack = data_file["data"]
print ('read', stack.shape)
img = stack[0]
stackout =np.zeros(stack.shape, np.uint16 )
'''
#------ Select a grid image from a Scientifica system -------
#2p.3 grid data
#system=3
#gridimg = "522609838_maxInt_a13a.png" # float
#gridimg = "522437527_maxInt_a13a.png"
#gridimg = "nofilter2p3_2pmode_1_512.png"
#gridimg = "520351219_maxInt_a13a.png"
#2p.4 grid data
#system=4
#gridimg = "524046446_maxInt_a13a.png"
#gridimg = "20160706_2p4_grid_locA_1p2x_XYT.tif" # funny image
#gridimg = "20160706_2p4_grid_locB_1p2x_XYT.png"
#gridimg = "max_20160707_2p4_GridImages_2P.png" # from Kate
#gridimg = "med_20160720_2p4_GridImages_2P.png"
#2p.5 grid data
#system=5
#gridimg = "525802500_maxInt_a13a.png"
#gridimg = "527051153_max_downsample_4Hz_0.png"
#gridimg = "2p5_grid_locC_1p2x_8bit.png"
if len(sys.argv) < 4:
print ("Reguire 3 arguments: gridimagein 2pnumber gridimageout")
exit(1)
gridimg = sys.argv[1]
system = int(sys.argv[2])
img = plt.imread(gridimg)
print ('input gridimg: ' + gridimg)
print ('read grid image shape:', img.shape)
#imgin =np.zeros(img.shape )
#imgout =np.zeros(img.shape, dtype=np.uint8 )
imgout =np.zeros(img.shape )
set_params(system)
print ("2p." +str(system))
print ('aL bL', aL, bL)
print ('aR bR', aR, bR)
create_xtable()
maxlevel= 1.0 #grid 8bit png
xdewarp(img, imgout )
'''
maxlevel= 65535 #16bit
for n in range (0,stack.shape[0]):
imgin = stack[n]
xdewarp(imgin, imgout )
stackout[n] = imgout.astype(np.uint16) # just for test
#stack[n] = imgout.astype(np.uint16) # just for test copying back - scary!
print 'frame done', n
print 'last frame inimg max=', np.amax(imgin)
print 'last frame imgout max=', np.amax(imgout)
'''
#print imgout.shape
#plt.imsave('imgout.png', imgout[:,:], cmap=plt.cm.gray) #rgb version
# to PIL
image = Image.fromarray((imgout[:,:]*255).astype(np.uint8)) # for 8bit grid png
#image = Image.fromarray((imgout[:,L:512-R]*255).astype(np.uint8)) # for 8bit grid png
#image.save('imgout_unit8.png') # grayscale version
#image = Image.fromarray((imgout/256).astype(np.uint8)) # for 16bit movie w/ 2 sides
#image = Image.fromarray((imgout[:,L:512-R]/256).astype(np.uint8)) # for 16bit movie
#image.save('gridimgout_uint8.png') # grayscale version dewarped grid
#print "output: gridimgout_uint8.png"
image.save(sys.argv[3]) # grayscale version dewarped grid
print ("output: " + sys.argv[3])
'''
f = h5py.File(os.path.join(output_dir,output_file_name), 'w')
if FOVwidth == 512:
f["data"] = stackout # 512 x 512 keep 2 sides
#f["data"] = stack # 512 x 512 keep 2 sides
else:
f["data"] = stackout[:,:,L:512-R] # remove 2 sides
#f["data"] = stack[:,:,L:512-R] # remove 2 sides
f.close()
'''
exit(0)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment