FAQ | This is a LIVE service | Changelog

Skip to content
Snippets Groups Projects
Unverified Commit 0bad96be authored by Eduardo Gonzalez Solares's avatar Eduardo Gonzalez Solares
Browse files

Style format changes

parent d19b748c
No related branches found
No related tags found
1 merge request!3Dask
Pipeline #163279 failed
Stage: test
Hide whitespace changes
Inline Side-by-side
stpt_pipeline/chi_functions.py
+ 194
− 132
View file @ 0bad96be
import numpy as np
from scipy.signal import medfilt2d
#
def MAD(x):
return np.median(np.abs(x-np.median(x)))
return np.median(np.abs(x - np.median(x)))
#
def read_mosaicifile_stpt(filename):
#
# Reads the mosaic.txt file that sets up the stage and
# Reads the mosaic.txt file that sets up the stage and
# returns the expected positions in microns
#
file_p=open(filename,'r')
# delta_x/_y are the positions in microns,
file_p = open(filename, 'r')
# delta_x/_y are the positions in microns,
# label_x/_y are the names of each position,
# that normally indicate the column/row position
# but are not used later
delta_x=[]
label_x=[]
delta_y=[]
label_y=[]
delta_x = []
label_x = []
delta_y = []
label_y = []
for this_line in file_p.readlines():
if this_line.find('XPos')>-1:
temp=this_line.split(':')
if this_line.find('XPos') > -1:
temp = this_line.split(':')
delta_x.append(float(temp[1]))
label_x.append(temp[0])
if this_line.find('YPos')>-1:
temp=this_line.split(':')
if this_line.find('YPos') > -1:
temp = this_line.split(':')
delta_y.append(float(temp[1]))
label_y.append(temp[0])
#
file_p.close()
#
return np.array(delta_x),np.array(delta_y),np.array(label_x),np.array(label_y)
return np.array(delta_x), np.array(delta_y), np.array(label_x), np.array(label_y)
#
def prepare_image_conf(img,conf,ORIENTATION='X',apply_filter=False,
DOUBLE_MEDIAN=False,IMG_STD=-1):
def prepare_image_conf(
img, conf, ORIENTATION='X', apply_filter=False, DOUBLE_MEDIAN=False, IMG_STD=-1
):
#
# Generates properly aligned images, as the crossmatching code works
# assuming a fixed relative orientation of the images.
#
# Returns the rotated image, a rotated filtered image to be
# Returns the rotated image, a rotated filtered image to be
# used in the crossmatching plus a confidence map. If apply_filter
# is switched on, a median substracted image (either 1D or 2D median)
# is returned as fimg_filt, otherwise a copy of the original image
#
N_SIGMA_CLIP=1.5
N_SIGMA_CLIP = 1.5
#
if ORIENTATION=='X':
img_fin=np.rot90(img,1)
conf_fin=np.rot90(conf,1)
if ORIENTATION == 'X':
img_fin = np.rot90(img, 1)
conf_fin = np.rot90(conf, 1)
else:
img_fin=img.copy()
conf_fin=conf.copy()
img_fin = img.copy()
conf_fin = conf.copy()
#
if apply_filter:
if DOUBLE_MEDIAN:
img_filt=img_fin-\
0.5*(medfilt2d(img_fin,(51,1))+medfilt2d(img_fin,(1,51)))
img_filt = img_fin - 0.5 * (
medfilt2d(img_fin, (51, 1)) + medfilt2d(img_fin, (1, 51))
)
else:
img_filt=img_fin-medfilt2d(img_fin,(51,1))
img_filt = img_fin - medfilt2d(img_fin, (51, 1))
else:
img_filt=img_fin.copy()
img_filt = img_fin.copy()
#
mask_fin=np.ones_like(img_fin)
mask_fin = np.ones_like(img_fin)
#
# if STD=-1 the std for the confidence mask is computed over the image,
# otherwise the valued passed as STD is used. By default the code uses
# a STD computed over the whole cube, that is mucho more robust that
# a STD computed over the whole cube, that is mucho more robust that
# local estimations
#
if IMG_STD==-1:
mask_fin[np.abs(img_filt)<=N_SIGMA_CLIP*1.48*MAD(img_filt)]=0.
if IMG_STD == -1:
mask_fin[np.abs(img_filt) <= N_SIGMA_CLIP * 1.48 * MAD(img_filt)] = 0.0
else:
mask_fin[img_fin<=IMG_STD]=0.
mask_fin[img_fin <= IMG_STD] = 0.0
# input conf map
mask_fin[conf_fin<0.7]=0.
mask_fin[conf_fin < 0.7] = 0.0
#
return img_fin,img_filt,mask_fin
return img_fin, img_filt, mask_fin
#
#
def find_overlap_conf(img_ref,conf_ref,img_obj,conf_obj,ORIENTATION='X',\
produce_image=False,blind_start=True,init_desp=[-50,1858],\
DOUBLE_MEDIAN=False,return_chi=False,IMG_STD=-1):
def find_overlap_conf(
img_ref,
conf_ref,
img_obj,
conf_obj,
ORIENTATION='X',
produce_image=False,
blind_start=True,
init_desp=[-50, 1858],
DOUBLE_MEDIAN=False,
return_chi=False,
IMG_STD=-1,
):
#
# This calculates the displacementes that minimize the squared difference
# between images.
#
REFINE_CHI=True
REFINE_CHI = True
#
STEPS = 30 # The chi surface will be of STEPSxSTEPS size
#
STEPS=30 # The chi surface will be of STEPSxSTEPS size
#
# DELTA is the step that dx and dy are increased at each
# point of the evaliations of chi.
#
# If an initial estimate is provided, the (dx,dy) space
# of possible displacements to explore is narrower
# If an initial estimate is provided, the (dx,dy) space
# of possible displacements to explore is narrower
#
if blind_start:
DELTA=[16,8,1]
DELTA = [16, 8, 1]
else:
DELTA=[8,2,1]
dx=init_desp[0]
dy=init_desp[1]
DELTA = [8, 2, 1]
dx = init_desp[0]
dy = init_desp[1]
#
# By default the code assumes that img0 is the reference image and that
# By default the code assumes that img0 is the reference image and that
# img1 is aligned along the X direction (in the python matrix coordinate
# frame). If this is not the case, ORIENTATION is set to Y and the images
# rotated inside prepare_image
#
img0,img0_filt,mask0=prepare_image_conf(img_ref,conf_ref,apply_filter=False,\
ORIENTATION=ORIENTATION,DOUBLE_MEDIAN=DOUBLE_MEDIAN,IMG_STD=IMG_STD)
img1,img1_filt,mask1=prepare_image_conf(img_obj,conf_obj,apply_filter=False,\
ORIENTATION=ORIENTATION,DOUBLE_MEDIAN=DOUBLE_MEDIAN,IMG_STD=IMG_STD)
DET_SIZE=img0.shape[1]
img0, img0_filt, mask0 = prepare_image_conf(
img_ref,
conf_ref,
apply_filter=False,
ORIENTATION=ORIENTATION,
DOUBLE_MEDIAN=DOUBLE_MEDIAN,
IMG_STD=IMG_STD,
)
img1, img1_filt, mask1 = prepare_image_conf(
img_obj,
conf_obj,
apply_filter=False,
ORIENTATION=ORIENTATION,
DOUBLE_MEDIAN=DOUBLE_MEDIAN,
IMG_STD=IMG_STD,
)
DET_SIZE = img0.shape[1]
#
#
for i_delta in range(len(DELTA)):
print('Iteration {0:1d}'.format(i_delta+1))
# the first iteration sets up the displacements to be
print('Iteration {0:1d}'.format(i_delta + 1))
# the first iteration sets up the displacements to be
# evaluated
if (i_delta==0)&blind_start:
if (i_delta == 0) & blind_start:
#
# add 10 so that at least there are 10px to compare overlap
# desp_large runs along the direction where the displacement is
# desp_large runs along the direction where the displacement is
# the largest, so typically will have values between 1700 and 2100,
# while desp_short runs in the other, and moves between -200 and 200.
#
desp_large=DET_SIZE-10-np.arange(STEPS+1)*DELTA[i_delta]-1.
desp_short=np.arange(STEPS+1)*DELTA[i_delta]-\
(STEPS-1)*DELTA[i_delta]*0.5
desp_large = DET_SIZE - 10 - np.arange(STEPS + 1) * DELTA[i_delta] - 1.0
desp_short = (
np.arange(STEPS + 1) * DELTA[i_delta]
- (STEPS - 1) * DELTA[i_delta] * 0.5
)
else:
# in case there is a prior estimation of the displacement
# that minimizes chi, the vectors are centered at that point
desp_short=np.arange(STEPS)*DELTA[i_delta]-\
(STEPS-1)*DELTA[i_delta]*0.5+dx
desp_short = (
np.arange(STEPS) * DELTA[i_delta]
- (STEPS - 1) * DELTA[i_delta] * 0.5
+ dx
)
# this just makes sure that we don't run out of image
if (STEPS-1)*DELTA[i_delta]*0.5+dy>=DET_SIZE:
delta_cor=(STEPS-1)*DELTA[i_delta]*0.5+dy-DET_SIZE
if (STEPS - 1) * DELTA[i_delta] * 0.5 + dy >= DET_SIZE:
delta_cor = (STEPS - 1) * DELTA[i_delta] * 0.5 + dy - DET_SIZE
else:
delta_cor=0
delta_cor = 0
#
desp_large=np.arange(STEPS)*DELTA[i_delta]-\
(STEPS-1)*DELTA[i_delta]*0.5+dy-delta_cor
# because these will be used as indexes, need to be
desp_large = (
np.arange(STEPS) * DELTA[i_delta]
- (STEPS - 1) * DELTA[i_delta] * 0.5
+ dy
- delta_cor
)
# because these will be used as indexes, need to be
# cast as ints.
desp_y=desp_large.astype(int)
desp_x=desp_short.astype(int)
desp_y = desp_large.astype(int)
desp_x = desp_short.astype(int)
#
chi=np.zeros((len(desp_x),len(desp_y)))
npx=np.ones((len(desp_x),len(desp_y)))
chi = np.zeros((len(desp_x), len(desp_y)))
npx = np.ones((len(desp_x), len(desp_y)))
for i in range(len(desp_y)):
# a slice of the reference image
t_ref=img0[:,desp_y[i]:]
err_ref=np.sqrt(img0[:,desp_y[i]:]) # assuming poisson errors
mask_ref=mask0[:,desp_y[i]:]
t_ref = img0[:, desp_y[i] :]
err_ref = np.sqrt(img0[:, desp_y[i] :]) # assuming poisson errors
mask_ref = mask0[:, desp_y[i] :]
# same for the other image
t1=img1[:,0:-desp_y[i]]
maskt=mask1[:,0:-desp_y[i]]
err_1=np.sqrt(img1[:,0:-desp_y[i]])
t1 = img1[:, 0 : -desp_y[i]]
maskt = mask1[:, 0 : -desp_y[i]]
err_1 = np.sqrt(img1[:, 0 : -desp_y[i]])
#
for j in range(len(desp_x)):
if desp_x[j]<0:
if desp_x[j] < 0:
# this is the difference of the overlap of the reference image
# and the displaced image. Depending on whether the small displacemente
# is positive, negative or zero the indexes need to be generated,
# so there's on branch of the if
temp=t1[np.abs(desp_x[j]):,]*maskt[np.abs(desp_x[j]):,:]-\
t_ref[0:desp_x[j],]*mask_ref[0:desp_x[j],:]
# combined mask
mask_final=maskt[np.abs(desp_x[j]):,:]+mask_ref[0:desp_x[j],:]
temp = (
t1[np.abs(desp_x[j]) :,] * maskt[np.abs(desp_x[j]) :, :]
- t_ref[0 : desp_x[j],] * mask_ref[0 : desp_x[j], :]
)
# combined mask
mask_final = (
maskt[np.abs(desp_x[j]) :, :] + mask_ref[0 : desp_x[j], :]
)
# and erros added in cuadrature
div_t=np.sqrt(err_ref[0:desp_x[j],:]**2+err_1[np.abs(desp_x[j]):,]**2)
#
elif desp_x[j]==0:
div_t = np.sqrt(
err_ref[0 : desp_x[j], :] ** 2
+ err_1[np.abs(desp_x[j]) :,] ** 2
)
#
elif desp_x[j] == 0:
#
temp=t1*maskt-t_ref*mask_ref
temp = t1 * maskt - t_ref * mask_ref
#
div_t=np.sqrt(err_ref**2+err_1**2)
div_t = np.sqrt(err_ref ** 2 + err_1 ** 2)
#
mask_final=maskt+mask_ref
mask_final = maskt + mask_ref
else:
#
temp=t1[0:-desp_x[j],]*maskt[0:-desp_x[j],]-\
t_ref[desp_x[j]:,]*mask_ref[desp_x[j]:,]
temp = (
t1[0 : -desp_x[j],] * maskt[0 : -desp_x[j],]
- t_ref[desp_x[j] :,] * mask_ref[desp_x[j] :,]
)
#
mask_final=maskt[0:-desp_x[j],]+mask_ref[desp_x[j]:,]
mask_final = maskt[0 : -desp_x[j],] + mask_ref[desp_x[j] :,]
#
div_t=np.sqrt(err_ref[np.abs(desp_x[j]):,]**2+err_1[0:-desp_x[j],]**2)
div_t = np.sqrt(
err_ref[np.abs(desp_x[j]) :,] ** 2 + err_1[0 : -desp_x[j],] ** 2
)
#
# if there are not enough good pixels we set chi to 10e7
if mask_final.sum()>0:
chi[j,i]=np.mean((temp[mask_final>0]/div_t[mask_final>0])**2)
#
npx[j,i]=np.sum(mask_final)
if mask_final.sum() > 0:
chi[j, i] = np.mean(
(temp[mask_final > 0] / div_t[mask_final > 0]) ** 2
)
#
npx[j, i] = np.sum(mask_final)
else:
chi[j,i]=1e7
chi[j, i] = 1e7
#
if REFINE_CHI:
# this generates an estimation of the "background" of the chi surface,
# that sometimes is heplful, mostly for cosmetic reasons, to have a flat
# chi with a very clear minumum
temp_bg0=np.array(list(np.median(chi,0))*chi.shape[1]).reshape(chi.shape)
temp_bg1=np.array(list(np.median(chi,1))*chi.shape[0]).reshape(chi.shape)
chi-=temp_bg0
temp_bg0 = np.array(list(np.median(chi, 0)) * chi.shape[1]).reshape(
chi.shape
)
temp_bg1 = np.array(list(np.median(chi, 1)) * chi.shape[0]).reshape(
chi.shape
)
chi -= temp_bg0
#
# now finfing the minimum
i_x,i_y=np.where(chi==chi.min())
# now finfing the minimum
i_x, i_y = np.where(chi == chi.min())
#
try:
# in case of many nans, this fails, default to middle
dx=desp_x[i_x][0]
dy=desp_y[i_y][0]
dx = desp_x[i_x][0]
dy = desp_y[i_y][0]
except:
dx=int(np.median(desp_x))
dy=int(np.median(desp_y))
dx = int(np.median(desp_x))
dy = int(np.median(desp_y))
#
print('Found min at dx={0:4d},dy={1:4d}'.format(dx,dy))
print('Found min at dx={0:4d},dy={1:4d}'.format(dx, dy))
#
# now on to generate the relevant outputs
#
......@@ -223,19 +280,21 @@ def find_overlap_conf(img_ref,conf_ref,img_obj,conf_obj,ORIENTATION='X',\
# This produces a mosaiced image, good to check
# if the procedure worked properly, but takes time.
#
big_picture,overlap=overlap_images(img0,img1,dx,dy)
big_picture, overlap = overlap_images(img0, img1, dx, dy)
#
if return_chi:
return dx,dy,big_picture,overlap,chi
return dx, dy, big_picture, overlap, chi
else:
return dx,dy,big_picture,overlap
return dx, dy, big_picture, overlap
else:
if return_chi:
return dx,dy,chi
return dx, dy, chi
else:
return dx,dy
return dx, dy
#
def overlap_images(img0,img1,dx,dy):
def overlap_images(img0, img1, dx, dy):
#
# Produces a mosaiced image displacing img1
# by dx,dy and overlapping it to img0. The
......@@ -246,26 +305,29 @@ def overlap_images(img0,img1,dx,dy):
# and overlap keeps track of how many real pixels went into
# each of the pixels of big_picture
#
dx=int(dx)
dy=int(dy)
big_picture=np.zeros((img0.shape[0]+np.abs(dx),\
img0.shape[1]+np.abs(dy)))
if dx<0:
delta_x=np.abs(dx)
dx = int(dx)
dy = int(dy)
big_picture = np.zeros((img0.shape[0] + np.abs(dx), img0.shape[1] + np.abs(dy)))
if dx < 0:
delta_x = np.abs(dx)
else:
delta_x=0
if dy<0:
delta_y=np.abs(dy)
delta_x = 0
if dy < 0:
delta_y = np.abs(dy)
else:
delta_y=0
overlap=np.zeros_like(big_picture)
big_picture[delta_x:delta_x+img0.shape[0],\
delta_y:delta_y+img0.shape[1]]=img0
overlap[delta_x:delta_x+img0.shape[0],\
delta_y:delta_y+img0.shape[1]]+=1
big_picture[delta_x+dx:delta_x+img0.shape[0]+dx,\
delta_y+dy:delta_y+dy+img0.shape[1]]+=img1
overlap[delta_x+dx:delta_x+img0.shape[0]+dx,\
delta_y+dy:delta_y+dy+img0.shape[1]]+=1
delta_y = 0
overlap = np.zeros_like(big_picture)
big_picture[
delta_x : delta_x + img0.shape[0], delta_y : delta_y + img0.shape[1]
] = img0
overlap[delta_x : delta_x + img0.shape[0], delta_y : delta_y + img0.shape[1]] += 1
big_picture[
delta_x + dx : delta_x + img0.shape[0] + dx,
delta_y + dy : delta_y + dy + img0.shape[1],
] += img1
overlap[
delta_x + dx : delta_x + img0.shape[0] + dx,
delta_y + dy : delta_y + dy + img0.shape[1],
] += 1
#
return big_picture,overlap
return big_picture, overlap
stpt_pipeline/mosaic_functions.py
+ 53
− 35
View file @ 0bad96be
......@@ -3,80 +3,98 @@ try:
except:
from imaxt_image.external import tifffile as tf
import numpy as np
from os import listdir,walk,mkdir
from os import listdir, walk, mkdir
from scipy.misc import imresize
#
def find_delta(dx,dy):
def find_delta(dx, dy):
#
# In theory each position in the mosaic stage is offset
# by a fixed value, but there's some jitter in the microscope
# so this finds this delta. It is used later to get the column
# and row position of each single image
#
r=np.sqrt((dx.astype(float)-dx[0])**2+(dy.astype(float)-dy[0])**2)
r[0]=r.max() # avoiding minimum at itself
r = np.sqrt((dx.astype(float) - dx[0]) ** 2 + (dy.astype(float) - dy[0]) ** 2)
r[0] = r.max() # avoiding minimum at itself
# first candidate
dx_1=np.abs(dx[r.argmin()]-dx[0])
dy_1=np.abs(dy[r.argmin()]-dy[0])
dx_1 = np.abs(dx[r.argmin()] - dx[0])
dy_1 = np.abs(dy[r.argmin()] - dy[0])
# second candidate
r[r.argmin()]=r.max()
dx_2=np.abs(dx[r.argmin()]-dx[0])
dy_2=np.abs(dy[r.argmin()]-dy[0])
r[r.argmin()] = r.max()
dx_2 = np.abs(dx[r.argmin()] - dx[0])
dy_2 = np.abs(dy[r.argmin()] - dy[0])
#
return np.max([dx_1,dx_2]),np.max([dy_1,dy_2])
return np.max([dx_1, dx_2]), np.max([dy_1, dy_2])
#
#
def get_img_list(root_dir,mosaic_size=56):
def get_img_list(root_dir, mosaic_size=56):
#
# Scans the directory, finds the images
# and returns the img names, the channel number,
# the corresponding optical slice and the
# and returns the img names, the channel number,
# the corresponding optical slice and the
# running image number. This requires the naming
# convention to be constant, so it is likely to
# crash periodically when the STPT machine changes
# convention to be constant, so it is likely to
# crash periodically when the STPT machine changes
# naming conventions.
#
imgs=[]
img_num=[]
optical_slice=[]
channel=[]
imgs = []
img_num = []
optical_slice = []
channel = []
for this_file in listdir(root_dir):
if (this_file.find('tif')>-1)&(this_file.find('imaxt')==-1)&(this_file[0]!='.'):
if (
(this_file.find('tif') > -1)
& (this_file.find('imaxt') == -1)
& (this_file[0] != '.')
):
imgs.append(this_file)
temp=this_file.split('.')[0].split('-')[-1].split('_')
temp = this_file.split('.')[0].split('-')[-1].split('_')
img_num.append(float(temp[0]))
channel.append(float(temp[1]))
optical_slice.append(int(img_num[-1]/mosaic_size)+1)
imgs=np.array(imgs)
i_ord=np.argsort(img_num)
imgs=imgs[i_ord]
img_num=np.array(img_num)[i_ord]
channel=np.array(channel)[i_ord]
optical_slice=np.array(optical_slice)[i_ord]
optical_slice.append(int(img_num[-1] / mosaic_size) + 1)
imgs = np.array(imgs)
i_ord = np.argsort(img_num)
imgs = imgs[i_ord]
img_num = np.array(img_num)[i_ord]
channel = np.array(channel)[i_ord]
optical_slice = np.array(optical_slice)[i_ord]
#
return imgs,img_num,channel,optical_slice
return imgs, img_num, channel, optical_slice
#
def get_mosaic_file(root_dir):
# just locates the name of the stage mosaic text file
# inside a provided dir
for this_file in listdir(root_dir):
if (this_file.find('osaic')>-1)&(this_file.find('txt')>-1):
if (this_file.find('osaic') > -1) & (this_file.find('txt') > -1):
return this_file
#
return ''
#
def get_img_cube(root_dir,imgs,img_num,channel,optical_slice,channel_to_use=4,slice_to_use=1):
def get_img_cube(
root_dir, imgs, img_num, channel, optical_slice, channel_to_use=4, slice_to_use=1
):
#
# Reads all the tiff files for a given channel and slice, and stores them into a cube.
# Reads all the tiff files for a given channel and slice, and stores them into a cube.
#
ii=np.where((channel==channel_to_use)&((optical_slice-np.min(optical_slice)+1)==slice_to_use))[0]
im_t=[tf.imread(root_dir+t).astype('float32') for t in imgs[ii]]
ii = np.where(
(channel == channel_to_use)
& ((optical_slice - np.min(optical_slice) + 1) == slice_to_use)
)[0]
im_t = [tf.imread(root_dir + t).astype('float32') for t in imgs[ii]]
#
return np.array(im_t)
#
#
def imread(img_name):
# this is just for convenience in case we change tiff library
im_t=tf.imread(img_name).astype('float32')
im_t = tf.imread(img_name).astype('float32')
#
return im_t
stpt_pipeline/stpt_displacement.py
+ 169
− 128
View file @ 0bad96be
import matplotlib
matplotlib.use('Agg')
try:
import tifffile as tf
......@@ -7,124 +8,137 @@ except:
from matplotlib import pyplot as pl
import numpy as np
from os import listdir
from chi_functions import read_mosaicifile_stpt,find_overlap_conf,MAD
from mosaic_functions import get_mosaic_file,get_img_cube,get_img_list,find_delta
from chi_functions import read_mosaicifile_stpt, find_overlap_conf, MAD
from mosaic_functions import get_mosaic_file, get_img_cube, get_img_list, find_delta
from scipy.ndimage.filters import gaussian_filter
from scipy.ndimage import geometric_transform
#
##################################################################################
#
##################################################################################
#
# these two functions are only used to filter and defringe the flat
def med_box(y,half_box=2):
def med_box(y, half_box=2):
#
ym=[]
ym = []
for i in range(len(y)):
# too close to cero
i_min = (i-half_box) if (i-half_box) >=0 else 0
i_min = (i - half_box) if (i - half_box) >= 0 else 0
# too close to end
i_min = i_min if i_min+2*half_box <= len(y)-1 else len(y)-1-2*half_box
ym.append(np.median(y[i_min:i_min+2*half_box]))
i_min = (
i_min if i_min + 2 * half_box <= len(y) - 1 else len(y) - 1 - 2 * half_box
)
ym.append(np.median(y[i_min : i_min + 2 * half_box]))
#
return np.array(ym)
#
def defringe(img):
#
fr_img=img.copy()
fr_img = img.copy()
for i in range(fr_img.shape[1]):
if i<5:
t=np.median(img[:,0:10],1)
elif i>fr_img.shape[1]-5:
t=np.median(img[:,-10:],1)
if i < 5:
t = np.median(img[:, 0:10], 1)
elif i > fr_img.shape[1] - 5:
t = np.median(img[:, -10:], 1)
else:
t=np.median(img[:,i-5:i+5],1)
t = np.median(img[:, i - 5 : i + 5], 1)
#
fr_img[:,i]=img[:,i]-med_box(t,5)
fr_img[:, i] = img[:, i] - med_box(t, 5)
#
return fr_img
#
# get_coords is the function that feeds geometric_transform
# get_coords is the function that feeds geometric_transform
# in order to correct for the optical distortion of the detector.
#
def get_coords(coords,cof,center_x,max_x,direct=True):
def get_coords(coords, cof, center_x, max_x, direct=True):
#
max_desp=cof[0]*coords[1]**2+cof[1]*coords[1]+cof[2]
dy_cof=max_desp/(max_x-center_x)**2
max_desp = cof[0] * coords[1] ** 2 + cof[1] * coords[1] + cof[2]
dy_cof = max_desp / (max_x - center_x) ** 2
if direct:
sign=(coords[0]-center_x)/np.abs(coords[0]-center_x)
sign = (coords[0] - center_x) / np.abs(coords[0] - center_x)
if np.isnan(sign):
sign=1.0
xi=np.abs(coords[0]-center_x)
return (center_x+sign*(xi+dy_cof*xi**2),coords[1])
sign = 1.0
xi = np.abs(coords[0] - center_x)
return (center_x + sign * (xi + dy_cof * xi ** 2), coords[1])
else:
xi=np.abs(coords[0]-center_x-cof[2])
sign=(coords[0]-center_x-cof[2])/np.abs(coords[0]-center_x-cof[2])
xi = np.abs(coords[0] - center_x - cof[2])
sign = (coords[0] - center_x - cof[2]) / np.abs(coords[0] - center_x - cof[2])
if np.isnan(sign):
sign=1.0
return (center_x+sign*(np.sqrt(1+4*dy_cof*xi)-1)/(2*dy_cof),coords[1])
sign = 1.0
return (
center_x + sign * (np.sqrt(1 + 4 * dy_cof * xi) - 1) / (2 * dy_cof),
coords[1],
)
#
##################################################################################
#
# x/y_min and max are the lower/upper boundaries of the useful section of the
# detector.
#
# x/y_min and max are the lower/upper boundaries of the useful section of the
# detector.
#
# Norm val is a normalization value for the stored values, ideally the gain
#
# do_flat activates flat correction prior to crossmatching. defringe does that to
# do_flat activates flat correction prior to crossmatching. defringe does that to
# the flat before dividing by it, but it is time consuming and not very
# useful at the moment
#
# cof_dist are the coefficients of the optical distortion, as measured from
# the images themselves.
#
x_min=2
x_max=2080
y_min=80
y_max=1990
norm_val=10000.
do_flat=True
defring=False
channel_to_use=4
cof_dist=np.array([3.93716645e-05, -7.37696218e-02, 2.52457306e+01])/2.
x_min = 2
x_max = 2080
y_min = 80
y_max = 1990
norm_val = 10000.0
do_flat = True
defring = False
channel_to_use = 4
cof_dist = np.array([3.93716645e-05, -7.37696218e-02, 2.52457306e01]) / 2.0
#
# this is the directory where all the subdirectories (one per physical slice)
# are stored. Ideally this should be passed as an argument
#
root_dir='../20190201_tumour_opticalsections_extraoverlap/'
root_dir = '../20190201_tumour_opticalsections_extraoverlap/'
#
#
##################################################################################
#
# This function transform the raw images into the ones used for crossmatching
# and mosaicing
magic_function = \
lambda x: np.flipud(x/nflat)[x_min:x_max,y_min:y_max]/norm_val
magic_function = lambda x: np.flipud(x / nflat)[x_min:x_max, y_min:y_max] / norm_val
#
# Read flat
#
if do_flat:
fl=np.load('flat2.npy')
fl = np.load('flat2.npy')
if defringe:
fr_img=defringe(fl[:,:,channel_to_use-1])
nflat=(fl[:,:,channel_to_use-1]-fr_img)/np.median(fl[:,:,channel_to_use-1])
fr_img = defringe(fl[:, :, channel_to_use - 1])
nflat = (fl[:, :, channel_to_use - 1] - fr_img) / np.median(
fl[:, :, channel_to_use - 1]
)
else:
nflat=(fl[:,:,channel_to_use-1])/np.median(fl[:,:,channel_to_use-1])
nflat[nflat<0.5]=1.0
nflat = (fl[:, :, channel_to_use - 1]) / np.median(fl[:, :, channel_to_use - 1])
nflat[nflat < 0.5] = 1.0
else:
nflat=1.0
nflat = 1.0
#
# This lists all the subdirectories. Once there is an
# This lists all the subdirectories. Once there is an
# standarized naming convention this will have to be edited,
# as at the moment looks for dir names with the '4t1' string
# on them, as all the experiments had this
#
dirs=[]
dirs = []
for this_file in listdir(root_dir):
if this_file.find('4t1')>-1:
if this_file.find('.txt')>-1:
if this_file.find('4t1') > -1:
if this_file.find('.txt') > -1:
continue
dirs.append(root_dir+'/'+this_file+'/')
dirs.append(root_dir + '/' + this_file + '/')
dirs.sort()
#
# Now we loop over all subirectories
......@@ -133,128 +147,155 @@ for this_dir in dirs:
#
print(this_dir)
# read mosaic file, get delta in microns
mosaic_file=get_mosaic_file(this_dir)
dx,dy,lx,ly=read_mosaicifile_stpt(this_dir+mosaic_file)
delta_x,delta_y=find_delta(dx,dy)
mosaic_file = get_mosaic_file(this_dir)
dx, dy, lx, ly = read_mosaicifile_stpt(this_dir + mosaic_file)
delta_x, delta_y = find_delta(dx, dy)
# dx0,dy0 are just the coordinates of each image
# in columns/rows
dx0=np.round((dx-dx.min())/delta_x).astype(int)
dy0=np.round((dy-dy.min())/delta_y).astype(int)
dx0 = np.round((dx - dx.min()) / delta_x).astype(int)
dy0 = np.round((dy - dy.min()) / delta_y).astype(int)
# get image names and select the first optical slice
imgs,img_num,channel,optical_slice=get_img_list(this_dir,mosaic_size=len(dx))
optical_slice=optical_slice-optical_slice.min()+1
imgs, img_num, channel, optical_slice = get_img_list(this_dir, mosaic_size=len(dx))
optical_slice = optical_slice - optical_slice.min() + 1
# read image cube
img_cube=get_img_cube(\
this_dir,imgs,img_num,channel,optical_slice,channel_to_use=channel_to_use,\
slice_to_use=optical_slice.min()\
)
img_cube = get_img_cube(
this_dir,
imgs,
img_num,
channel,
optical_slice,
channel_to_use=channel_to_use,
slice_to_use=optical_slice.min(),
)
#
# std_img is the per pixel std of the cube,
# shapes is the reference detector size
#
std_img=magic_function(np.std(img_cube,0))
shapes=magic_function(img_cube[0,...]).shape
std_img = magic_function(np.std(img_cube, 0))
shapes = magic_function(img_cube[0, ...]).shape
#
print('Found {0:d} images, cube STD {1:.3f}'.\
format(img_cube.shape[0],np.mean(std_img)))
print(
'Found {0:d} images, cube STD {1:.3f}'.format(
img_cube.shape[0], np.mean(std_img)
)
)
# Default confidence map, as the array is square but
# the geometrical transform to correct distortion
# leaves some pixels empty, so these are set to zero in
# leaves some pixels empty, so these are set to zero in
# this conf. map. As all images go through the same initial
# transformation, only on conf map is needed
general_conf=np.ones_like(magic_function(img_cube[0,...]))
dist_conf=geometric_transform(general_conf,\
get_coords,output_shape=(int(shapes[0]),\
shapes[1]),extra_arguments=(cof_dist,\
shapes[0]*0.5,shapes[0]),\
extra_keywords={'direct':True},mode='constant',cval=0.0,order=1)
general_conf = np.ones_like(magic_function(img_cube[0, ...]))
dist_conf = geometric_transform(
general_conf,
get_coords,
output_shape=(int(shapes[0]), shapes[1]),
extra_arguments=(cof_dist, shapes[0] * 0.5, shapes[0]),
extra_keywords={'direct': True},
mode='constant',
cval=0.0,
order=1,
)
#
# These are the matrixes where the displacements that come out
# of cross-correlation are stored. dx_mat[i,j] is the optimal
# displacement between image i and j in the x direction. If these
# images do not overlap, it is set to -9999
#
dx_mat=np.zeros((len(dx),len(dx)))-9999
dy_mat=np.zeros((len(dx),len(dx)))-9999
dx_mat = np.zeros((len(dx), len(dx))) - 9999
dy_mat = np.zeros((len(dx), len(dx))) - 9999
#
for i in range(len(dx0)):
# pick up if crash by reloading already calculated
# displacements
try:
dx_mat=np.load(this_dir+'desp_dist_x.npy')
dy_mat=np.load(this_dir+'desp_dist_y.npy')
dx_mat = np.load(this_dir + 'desp_dist_x.npy')
dy_mat = np.load(this_dir + 'desp_dist_y.npy')
except:
print('No saved displacements found')
#
# This is more or less the distance between images in detectors,
# so we only crossmatch images withn 1 detector size of the
# so we only crossmatch images withn 1 detector size of the
# reference image
#
r=(dx0-dx0[i])**2+(dy0-dy0[i])**2
r[i]=100
i_t=np.where(r<=1)[0]
r = (dx0 - dx0[i]) ** 2 + (dy0 - dy0[i]) ** 2
r[i] = 100
i_t = np.where(r <= 1)[0]
#
for this_img in i_t:
if dx_mat[i,this_img]!=-9999:
if dx_mat[i, this_img] != -9999:
# already done
print('({0:2d},{1:2d}) already done'.format(i,this_img))
print('({0:2d},{1:2d}) already done'.format(i, this_img))
continue
print('Finding shifts for ({0:2d},{1:2d})'.format(i,this_img))
print('Finding shifts for ({0:2d},{1:2d})'.format(i, this_img))
#
# Find relative orientation. As the crossmatching code always assumes
# that the images are aligned along the X axis and the reference is to
# that the images are aligned along the X axis and the reference is to
# the left, we need to juggle things if that's not the case
#
ind_ref=i
ind_obj=this_img
if np.abs(dx0[i]-dx0[this_img])>\
np.abs(dy0[i]-dy0[this_img]):
ORIENTATION='X'
ind_ref = i
ind_obj = this_img
if np.abs(dx0[i] - dx0[this_img]) > np.abs(dy0[i] - dy0[this_img]):
ORIENTATION = 'X'
# Relative position
if dx0[i]>dx0[this_img]:
ind_ref=this_img
ind_obj=i
continue
if dx0[i] > dx0[this_img]:
ind_ref = this_img
ind_obj = i
continue
else:
ORIENTATION='Y'
ORIENTATION = 'Y'
# Relative position
if dy0[i]>dy0[this_img]:
ind_ref=this_img
ind_obj=i
if dy0[i] > dy0[this_img]:
ind_ref = this_img
ind_obj = i
continue
#
# Transformed images
#
new_ref=geometric_transform(magic_function(img_cube[ind_ref,...]),\
get_coords,output_shape=(int(shapes[0]),\
shapes[1]),extra_arguments=(cof_dist,\
shapes[0]*0.5,shapes[0]),\
extra_keywords={'direct':True},mode='constant',cval=0.0,order=1)
new_obj=geometric_transform(magic_function(img_cube[ind_obj,...]),\
get_coords,output_shape=(int(shapes[0]),\
shapes[1]),extra_arguments=(cof_dist,\
shapes[0]*0.5,shapes[0]),\
extra_keywords={'direct':True},mode='constant',cval=0.0,order=1)
new_ref = geometric_transform(
magic_function(img_cube[ind_ref, ...]),
get_coords,
output_shape=(int(shapes[0]), shapes[1]),
extra_arguments=(cof_dist, shapes[0] * 0.5, shapes[0]),
extra_keywords={'direct': True},
mode='constant',
cval=0.0,
order=1,
)
new_obj = geometric_transform(
magic_function(img_cube[ind_obj, ...]),
get_coords,
output_shape=(int(shapes[0]), shapes[1]),
extra_arguments=(cof_dist, shapes[0] * 0.5, shapes[0]),
extra_keywords={'direct': True},
mode='constant',
cval=0.0,
order=1,
)
# finding shift
dx,dy,mer=find_overlap_conf(\
new_ref,dist_conf,\
new_obj,dist_conf,\
ORIENTATION=ORIENTATION,produce_image=False,\
return_chi=True,DOUBLE_MEDIAN=False,IMG_STD=np.mean(std_img)
)
dx, dy, mer = find_overlap_conf(
new_ref,
dist_conf,
new_obj,
dist_conf,
ORIENTATION=ORIENTATION,
produce_image=False,
return_chi=True,
DOUBLE_MEDIAN=False,
IMG_STD=np.mean(std_img),
)
# undoing the index shifts
if ORIENTATION=='Y':
dx_mat[ind_ref,ind_obj]=dx
dy_mat[ind_ref,ind_obj]=dy
if ORIENTATION == 'Y':
dx_mat[ind_ref, ind_obj] = dx
dy_mat[ind_ref, ind_obj] = dy
#
dx_mat[ind_obj,ind_ref]=-dx
dy_mat[ind_obj,ind_ref]=-dy
if ORIENTATION=='X':
dx_mat[ind_ref,ind_obj]=dy
dy_mat[ind_ref,ind_obj]=-dx
dx_mat[ind_obj, ind_ref] = -dx
dy_mat[ind_obj, ind_ref] = -dy
if ORIENTATION == 'X':
dx_mat[ind_ref, ind_obj] = dy
dy_mat[ind_ref, ind_obj] = -dx
#
dx_mat[ind_obj,ind_ref]=-dy
dy_mat[ind_obj,ind_ref]=dx
dx_mat[ind_obj, ind_ref] = -dy
dy_mat[ind_obj, ind_ref] = dx
# storing
np.save(this_dir+'desp_dist_x',dx_mat)
np.save(this_dir+'desp_dist_y',dy_mat)
np.save(this_dir + 'desp_dist_x', dx_mat)
np.save(this_dir + 'desp_dist_y', dy_mat)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment