From 0bad96bee4b05c000b4bba3871b1757851a9c9f9 Mon Sep 17 00:00:00 2001
From: Eduardo Gonzalez Solares <eglez@ast.cam.ac.uk>
Date: Mon, 29 Apr 2019 15:32:51 +0100
Subject: [PATCH] Style format changes
---
stpt_pipeline/chi_functions.py | 326 +++++++++++++++++------------
stpt_pipeline/mosaic_functions.py | 88 ++++----
stpt_pipeline/stpt_displacement.py | 297 +++++++++++++++-----------
3 files changed, 416 insertions(+), 295 deletions(-)
diff --git a/stpt_pipeline/chi_functions.py b/stpt_pipeline/chi_functions.py
index 9326b28..358a0e8 100644
--- a/stpt_pipeline/chi_functions.py
+++ b/stpt_pipeline/chi_functions.py
@@ -1,220 +1,277 @@
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
diff --git a/stpt_pipeline/mosaic_functions.py b/stpt_pipeline/mosaic_functions.py
index bd0454e..2bc541e 100644
--- a/stpt_pipeline/mosaic_functions.py
+++ b/stpt_pipeline/mosaic_functions.py
@@ -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
diff --git a/stpt_pipeline/stpt_displacement.py b/stpt_pipeline/stpt_displacement.py
index 13f3821..e83443f 100644
--- a/stpt_pipeline/stpt_displacement.py
+++ b/stpt_pipeline/stpt_displacement.py
@@ -1,4 +1,5 @@
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)
--
GitLab