Commit acf1ce36 authored by Olga Stamati's avatar Olga Stamati
Browse files

new spam-pixelSearchPropagate script plus helper functions

parent da8375db
Pipeline #60318 passed with stages
in 13 minutes and 25 seconds
#!/usr/bin/env python
"""
This script performs a point-by-point "pixel search" which computes a correlation coefficient
of an imagette extracted in im1 to a brute-force search in a given search range in z, y, x in image 2.
Imagettes in im1 can either be defined with a nodeSpacing and a halfWindowSize or a labelled image.
Copyright (C) 2020 SPAM Contributors
This program is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the Free
Software Foundation, either version 3 of the License, or (at your option)
any later version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
"""
import os
os.environ['OPENBLAS_NUM_THREADS'] = '1'
import spam.DIC
import spam.mesh
import spam.helpers
import numpy
import progressbar
import argparse
import tifffile
from scipy.spatial import KDTree
# Define argument parser object
parser = argparse.ArgumentParser(description="spam-pixelSearchPropagate "+spam.helpers.optionsParser.GLPv3descriptionHeader +\
"This script performs a pixel search from im1 to im2 propagating the motion from the top guiding point\n",
formatter_class=argparse.RawTextHelpFormatter)
# Parse arguments with external helper function
args = spam.helpers.optionsParser.pixelSearchPropagate(parser)
print("Current Settings:")
argsDict = vars(args)
for key in sorted(argsDict):
print("\t{}: {}".format(key, argsDict[key]))
# Fill in search range
searchRange = numpy.array([args.SEARCH_RANGE[0], args.SEARCH_RANGE[1], args.SEARCH_RANGE[2], args.SEARCH_RANGE[3], args.SEARCH_RANGE[4], args.SEARCH_RANGE[5]])
pixelSearchCCmin = args.CC_MIN
weightingDistance = args.DIST
# Load reference image
im1 = tifffile.imread(args.im1.name)
# Detect unpadded 2D image first:
#if len(im1.shape) == 2: im1 = im1[numpy.newaxis, ...]
#if im1.shape[0] == 1: twoD = True
#else: twoD = False
# Load deformed image
im2 = tifffile.imread(args.im2.name)
#if len(im2.shape) == 2: im2 = im2[numpy.newaxis, ...]
if args.MASK1:
im1mask = tifffile.imread(args.MASK1.name) != 0
#if len(im1mask.shape) == 2:
# im1mask = im1mask[numpy.newaxis, ...]
else:
im1mask = None
if args.MASK2:
im2mask = tifffile.imread(args.MASK2.name) != 0
#if len(im2mask.shape) == 2:
# im2mask = im2mask[numpy.newaxis, ...]
else:
im2mask = None
print("\n\tRanking points")
gp = numpy.genfromtxt(args.gp.name)
#gp = numpy.abs(gp[:, -1:1:-1]) #swap x-z axes and give abs coordinates
guidingPoints = spam.mesh.rankPoints(gp, neighbourRadius=args.RADIUS)
numberOfPoints = guidingPoints.shape[0]
# Possible to do it on a structured grid, too!
# (provided that the fist point is well chosen...)
#gp, nodesDim = spam.DIC.makeGrid(im1.shape, args.NS)
#gp[0] = [26, 25, 245] #overwrite the first point (this can easily be an input actually)
#guidingPoints = spam.mesh.rankPoints(gp, neighbourRadius=args.RADIUS)
#numberOfPoints = guidingPoints.shape[0]
# Initialise arrays
PhiField = numpy.zeros((numberOfPoints, 4, 4))
for point in range(numberOfPoints):
PhiField[point] = numpy.eye(4)
pixelSearchCC = numpy.zeros((numberOfPoints), dtype=float)
# Returns compatible with register()
error = numpy.zeros((numberOfPoints), dtype=float)
returnStatus = numpy.zeros((numberOfPoints), dtype=int)
deltaPhiNorm = numpy.ones((numberOfPoints), dtype=int)
iterations = numpy.zeros((numberOfPoints), dtype=int)
print("\n\tStarting sequential Pixel Search")
widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfPoints)
pbar.start()
# Step 1: simple PS for first point
imagetteReturnsTop = spam.DIC.getImagettes(im1, guidingPoints[0], args.HWS, PhiField[0].copy(), im2, searchRange.copy(), applyF='no')
# If getImagettes was successful (size check and mask coverage check)
if imagetteReturnsTop['returnStatus'] == 1:
PSreturnsTop = spam.DIC.correlate.pixelSearch(imagetteReturnsTop['imagette1'],
imagetteReturnsTop['imagette2'],
imagette1mask = imagetteReturnsTop['imagette1mask'],
imagette2mask = imagetteReturnsTop['imagette2mask'],
returnError = True)
PhiField[0, 0:3, -1] += imagetteReturnsTop['pixelSearchOffset'] + PSreturnsTop[0]
pixelSearchCC[0] = PSreturnsTop[1]
error[0] = PSreturnsTop[2]
# Failed to extract imagettes or something
else:
print("That's a bit of a problem I'm afraid!")
# Step 2: Loop sequentially over the guiding points list
# 2.1: create the tree of the coordinates to find easily neighbours
treeCoord = KDTree(guidingPoints)
for point in range(1, numberOfPoints):
indices = []
radius = args.RADIUS
# 2.2: Extract good neighbours
# double the radius until it finds at least 1 point in the vicinity
while len(indices) < 1:
indices = numpy.array(treeCoord.query_ball_point(guidingPoints[point], radius))
# Discard current point and points with low CC from indices
indices = numpy.delete(indices, numpy.where(numpy.logical_or(indices == point, pixelSearchCC[indices] < pixelSearchCCmin))[0])
radius *= 2
# 2.3: Estimate initial displacement
# by a gaussian weighting of extracted good neighbours
distances = numpy.linalg.norm(guidingPoints[point] - guidingPoints[indices], axis=1)
weights = numpy.exp(-distances**2/weightingDistance**2)
initialDisplacement = numpy.sum(PhiField[indices, 0:3, -1]*weights[:, numpy.newaxis], axis=0)/weights.sum()
# 2.4: Call PS around the estimated position
PhiField[point, 0:3, -1] += initialDisplacement
imagetteReturns = spam.DIC.getImagettes(im1, guidingPoints[point], args.HWS, PhiField[point].copy(), im2, searchRange.copy(), applyF='no')
if imagetteReturns['returnStatus'] == 1:
PSreturns = spam.DIC.pixelSearch(imagetteReturns['imagette1'],
imagetteReturns['imagette2'],
imagette1mask = imagetteReturns['imagette1mask'],
imagette2mask = imagetteReturns['imagette2mask'],
returnError = True)
PhiField[point, 0:3, -1] += imagetteReturns['pixelSearchOffset'] + PSreturns[0]
pixelSearchCC[point] = PSreturns[1]
error[point] = PSreturns[2]
returnStatus[point] = imagetteReturns['returnStatus']
widgets[0] = progressbar.FormatLabel(" CC={:0>7.5f} ".format(PSreturns[1]))
pbar.update(point)
else:
PhiField[point, 0:3, -1] = [numpy.nan]*3
error[point] = numpy.inf
returnStatus[point] = imagetteReturns['returnStatus']
print("")
if args.TSV:
# Make one big array for writing:
# First the node number,
# 3 node positions,
# F[0:3,0:3]
# Pixel-search CC
TSVheader = "NodeNumber\tZpos\tYpos\tXpos\tFzz\tFzy\tFzx\tZdisp\tFyz\tFyy\tFyx\tYdisp\tFxz\tFxy\tFxx\tXdisp\tpixelSearchCC\treturnStatus\terror\tdeltaPhiNorm\titerations"
outMatrix = numpy.array([numpy.array(range(numberOfPoints)),
guidingPoints[:, 0], guidingPoints[:, 1], guidingPoints[:, 2],
PhiField[:, 0, 0], PhiField[:, 0, 1], PhiField[:, 0, 2], PhiField[:, 0, 3],
PhiField[:, 1, 0], PhiField[:, 1, 1], PhiField[:, 1, 2], PhiField[:, 1, 3],
PhiField[:, 2, 0], PhiField[:, 2, 1], PhiField[:, 2, 2], PhiField[:, 2, 3],
pixelSearchCC,
returnStatus,
error,
deltaPhiNorm,
iterations]).T
numpy.savetxt(args.OUT_DIR+"/"+args.PREFIX+".tsv",
outMatrix,
fmt='%.7f',
delimiter='\t',
newline='\n',
comments='',
header=TSVheader)
# Collect data into VTK output
if args.VTK:
# boring nans overwriting
disp = PhiField[:, 0:3, -1]
disp[numpy.logical_not(numpy.isfinite(disp))] = 0
magDisp = numpy.linalg.norm(disp, axis=1)
magDisp[numpy.logical_not(numpy.isfinite(magDisp))] = 0
VTKglyphDict = {'displacements': disp,
'mag(displacements)': magDisp ,
'pixelSearchCC': pixelSearchCC}
spam.helpers.writeGlyphsVTK(guidingPoints, VTKglyphDict, fileName=args.OUT_DIR + "/" + args.PREFIX + ".vtk")
......@@ -147,6 +147,7 @@ TTK_modules_pb = ["mesh", "label", "filters", "measurements", "kalisphera", "DIC
scripts = ['scripts/spam-mmr',
'scripts/spam-reg',
'scripts/spam-pixelSearch',
'scripts/spam-pixelSearchPropagate',
'scripts/spam-gdic',
"scripts/spam-ddic",
"scripts/spam-ldic",
......
......@@ -2192,3 +2192,188 @@ def pixelSearch(parser):
args.PREFIX += "-pixelSearch"
return args
def pixelSearchPropagate(parser):
parser.add_argument('im1',
metavar='im1',
type=argparse.FileType('r'),
help="Greyscale image of reference state for correlation")
parser.add_argument('im2',
metavar='im2',
type=argparse.FileType('r'),
help="Greyscale image of deformed state for correlation")
parser.add_argument('gp',
metavar='gp',
type=argparse.FileType('r'),
help='Path to file containing the guiding points')
parser.add_argument('-mf1',
'--maskFile1',
dest='MASK1',
default=None,
type=argparse.FileType('r'),
help="Path to tiff file containing the mask of image 1 -- masks zones not to correlate, which should be == 0")
parser.add_argument('-mf2',
'--maskFile2',
dest='MASK2',
default=None,
type=argparse.FileType('r'),
help="Path to tiff file containing the mask of image 2 -- masks zones not to correlate, which should be == 0")
parser.add_argument('-mc',
'--mask-coverage',
type=float,
default=0.5,
dest='MASK_COVERAGE',
help="In case a mask is defined, tolerance for a subvolume's pixels to be masked before it is skipped with RS=-5. Default = 0.5")
parser.add_argument('-sr',
'--search-range',
nargs=6,
type=int,
default=[-3, 3, -3, 3, -3, 3],
dest='SEARCH_RANGE',
help='Z- Z+ Y- Y+ X- X+ ranges (in pixels) for the pixel search. Default = +-3px')
# Default: window size equal in all three directions
parser.add_argument('-hws',
'--half-window-size',
nargs=1,
type=int,
default=10,
dest='HWS',
help="Half correlation window size (in pixels), measured each side of the node pixel (assumed equal in all 3 directions -- see -hws3 for different setting). Default = 10 px")
# Possible: window size different in all three directions
parser.add_argument('-hws3',
'--half-window-size-3',
nargs=3,
type=int,
default=[10, 10, 10],
dest='HWS',
help="Half correlation window size (in pixels), measured each side of the node pixel (different in 3 directions). Default = 10, 10, 10px")
#parser.add_argument('-ns3',
#'--node-spacing-3',
#nargs=3,
#type=int,
#default=None,
#dest='NS',
#help="Node spacing in pixels (different in 3 directions). Default = 10, 10, 10px")
parser.add_argument('-nr',
'--neighbourhood-radius',
type=float,
default=None,
dest='RADIUS',
help="Radius (in pixels) inside which to select neighbours. Default = mean(hws)+mean(sr)")
parser.add_argument('-gwd',
'--gaussian-weighting-distance',
type=float,
default=None,
dest='DIST',
help="Distance (in pixels) over which the neighbour's distance is weighted. Default = sum(hws)+sum(sr)")
parser.add_argument('-cct',
'--CC-threshold',
type=float,
default=0.9,
dest='CC_MIN',
help="Pixel search correlation coefficient threshold BELOW which the point is considered badly correlated. Default = 0.9")
parser.add_argument('-glt',
'--grey-low-threshold',
type=float,
default=-numpy.inf,
dest='GREY_LOW_THRESH',
help="Grey threshold on mean of reference imagette BELOW which the correlation is not performed. Default = -infinity")
parser.add_argument('-ght',
'--grey-high-threshold',
type=float,
default=numpy.inf,
dest='GREY_HIGH_THRESH',
help="Grey threshold on mean of reference imagette ABOVE which the correlation is not performed. Default = infinity")
parser.add_argument('-od',
'--out-dir',
type=str,
default=None,
dest='OUT_DIR',
help="Output directory, default is the dirname of im1 file")
parser.add_argument('-pre',
'--prefix',
type=str,
default=None,
dest='PREFIX',
help='Prefix for output files (without extension). Default is basename of im1 and im2 files')
parser.add_argument('-vtk',
'--VTKout',
action="store_true",
dest='VTK',
help='Activate VTK output format. Default = False')
parser.add_argument('-notsv',
'--noTSVout',
action="store_false",
dest='TSV',
help='Disactivate TSV output format. Default = False')
args = parser.parse_args()
# 2019-04-05 EA: 2D image detection approved by Christophe Golke, update for shape 2019-08-29
#tiff = tifffile.TiffFile(args.im1.name)
#if len(tiff.pages) == 1 and len(tiff.series[0].shape) == 2:
# twoD = True
#else:
# twoD = False
#tiff.close()
# If we have no out dir specified, deliver on our default promise -- this can't be done inline before since parser.parse_args() has not been run at that stage.
if args.OUT_DIR is None:
args.OUT_DIR = os.path.dirname(args.im1.name)
# However if we have no dir, notice this and make it the current directory.
if args.OUT_DIR == "":
args.OUT_DIR = "./"
else:
# Check existence of output directory
try:
if args.OUT_DIR:
os.makedirs(args.OUT_DIR)
else:
args.DIR_out = os.path.dirname(args.im1.name)
except OSError:
if not os.path.isdir(args.OUT_DIR):
raise
# Output file name prefix
if args.PREFIX is None:
args.PREFIX = os.path.splitext(os.path.basename(args.im1.name))[0] + "-" + os.path.splitext(os.path.basename(args.im2.name))[0] + "-pixelSearchPropagate"
else:
args.PREFIX += "-pixelSearchPropagate"
# Catch 3D options
if len(args.HWS) == 1:
args.HWS = [args.HWS[0]]*3
# Fill radius and dist if not given
if args.RADIUS is None:
args.RADIUS = numpy.mean(args.HWS) + (args.SEARCH_RANGE[1]+args.SEARCH_RANGE[3]+args.SEARCH_RANGE[5])/3.0
if args.DIST is None:
args.DIST = numpy.sum(args.HWS) + numpy.sum(args.SEARCH_RANGE[1]+args.SEARCH_RANGE[3]+args.SEARCH_RANGE[5])
## Catch and overwrite 2D options
#if twoD:
#args.HWS[0] = 0
#args.SEARCH_RANGE[0] = 0
#args.SEARCH_RANGE[1] = 0
return args
......@@ -1061,3 +1061,69 @@ def create8Patche(lengths, lcar1, lcar2, origin=[0., 0., 0.], verbose=False, gms
return
def rankPoints(points, neighbourRadius=20):
'''
This function ranks an array of points around the top point
Parameters
----------
points: numpy array N x 3
Coordinates (zyx) of the points
neighbourRadius: float, optional
Distance from the current point to include neighbours
If no neighbour is found, then the closest point is taken
Default: 20
Returns
-------
pointsRanked: numpy array N x 3
Coordinates (zyx) of the ranked points
Note
-----
Based on https://gricad-gitlab.univ-grenoble-alpes.fr/DVC/pt4d
'''
from scipy.spatial import KDTree
import progressbar
# counters
p = pR = 0
# create the ranked array, first ranked point is the top point of the input array
pointsRanked = numpy.zeros_like(points)
pointsRanked[0] = points[0]
# remove ranked point from input array
points = numpy.delete(points, 0, 0)
# Create progressbar
numberOfPoints = pointsRanked.shape[0]
widgets = [progressbar.FormatLabel(''), ' ', progressbar.Bar(), ' ', progressbar.AdaptiveETA()]
pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfPoints)
pbar.start()
while points.shape[0] >= 1:
# Try to see how many points can be found around the current point based on the given radius
treeCoord = KDTree(points)
indRadius = numpy.array(treeCoord.query_ball_point(pointsRanked[pR], neighbourRadius))
indRadius = indRadius[numpy.argsort(indRadius)]
# if no points inside the given radius, just find the closest point
if len(indRadius)<1:
distance, ind = treeCoord.query(pointsRanked[pR], k=1)
indRadius = numpy.array([ind])
# fill in the ranked array with the point(s) found
pointsRanked[p+1:p+1+len(indRadius)] = points[indRadius]
# remove ranked point(s) from input array
points = numpy.delete(points, indRadius, 0)
# update counters
p += len(indRadius)
pR += 1
widgets[0] = progressbar.FormatLabel("{:.1f}%%".format((pR/numberOfPoints)*100))
pbar.update(pR)
return pointsRanked
Markdown is supported
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