Skip to content

GitLab

Commit 5c3198e9 authored by Edward Andò's avatar Edward Andò
Browse files

implemented -pfd for pixelSearch and avoided dangerous condition when not in -pfd

parent 03328023
Hide whitespace changes
Inline Side-by-side
scripts/spam-pixelSearch
View file @ 5c3198e9
......@@ -150,42 +150,50 @@ if args.PHIFILE is not None:
# If the read Phi-file contains multiple lines it's an F field!
else:
# Read the coordinates and values of the input F field
fieldCoords = PhiFromFile["fieldCoords"]
fieldValues = PhiFromFile["PhiField"]
if args.PHIFILE_DIRECT:
print("spam-pixelSearch: Assuming loaded PhiFile is coherent with the current run.")
PhiField = PhiFromFile["PhiField"]
# Create the k-d tree of the coordinates of the input F field
from scipy.spatial import KDTree
nNeighbours = args.NEIGHBOURS
fieldTree = KDTree(fieldCoords)
else:
# We don't trust this completely, re-interpolate it onto the grid
# Read the coordinates and values of the input F field
fieldCoords = PhiFromFile["fieldCoords"]
fieldValues = PhiFromFile["PhiField"]
# Loop over each point of the current grid
for node in range(nodePositions.shape[0]):
# if verbose: print "\nWorking on node {} {:0.2f}%".format( node, (node/float(numberofPoints))*100)
# Create the k-d tree of the coordinates of the input F field
from scipy.spatial import KDTree
nNeighbours = args.NEIGHBOURS
fieldTree = KDTree(fieldCoords)
# Loop over each point of the current grid
for node in range(nodePositions.shape[0]):
# if verbose: print "\nWorking on node {} {:0.2f}%".format( node, (node/float(numberofPoints))*100)
# reset F to zero, since we'll be doing an additive interpolation at the bottom here
PhiField[node][0:3, 0:3] = 0
# Calculate the distance of the point of the current grid to the points of the input F field
distance, indices = fieldTree.query(nodePositions[node], k=nNeighbours)
# Calculate the distance of the point of the current grid to the points of the input F field
distance, indices = fieldTree.query(nodePositions[node], k=nNeighbours)
# Check if we've hit the same point in the two grids
if numpy.any(distance == 0):
# Check if we've hit the same point in the two grids
if numpy.any(distance == 0):
# Copy F of that point to the F in the current grid point
PhiField[node] = fieldValues[indices][numpy.where(distance == 0)].copy()
# Copy F of that point to the F in the current grid point
PhiField[node] = fieldValues[indices][numpy.where(distance == 0)].copy()
# If not, consider the closest neighbours
else:
# If not, consider the closest neighbours
else:
# Compute the "Inverse Distance Weighting" since the closest points should have the major influence
weightSumInv = sum(1/distance)
# Compute the "Inverse Distance Weighting" since the closest points should have the major influence
weightSumInv = sum(1/distance)
# Loop over each neighbour
for neighbour in range(nNeighbours):
# Loop over each neighbour
for neighbour in range(nNeighbours):
# Calculate it's weight
weightInv = (1/distance[neighbour])/float(weightSumInv)
# Calculate it's weight
weightInv = (1/distance[neighbour])/float(weightSumInv)
# Fill the F of the current grid point with the weighted F components of the ith nearest neighbour in the input F field
PhiField[node][:-1] += fieldValues[indices[neighbour]][:-1]*weightInv
# Fill the F of the current grid point with the weighted F components of the ith nearest neighbour in the input F field
PhiField[node][:-1] += fieldValues[indices[neighbour]][:-1]*weightInv
### This vvv is the pixelSearchOnGrid function from grid.py
......
tests/test_scripts.py
View file @ 5c3198e9
......@@ -89,7 +89,7 @@ class testAll(unittest.TestCase):
pass
def _test_gridDICtools(self):
def test_gridDICtools(self):
# Make test directory
if not os.path.isdir(testFolder):
os.mkdir(testFolder)
......@@ -671,7 +671,7 @@ class testAll(unittest.TestCase):
#### 3. Rerun ddic -- Step0 -> Step1 with prev result
########################################################
### 3a ... pixel search reading prev pixel search and +- 2 search around that
exitCode = subprocess.call(["spam-pixelSearch", "-pf", testFolder + "balls-2b-pixelSearch.tsv",
exitCode = subprocess.call(["spam-pixelSearch", "-pf", testFolder + "balls-2b-pixelSearch.tsv", "-pfd"
"-sr",
"0", "0", "0", "0", "0", "0",
testFolder + "Step0.tif", testFolder + "Step1.tif",
......
tools/helpers/optionsParser.py
View file @ 5c3198e9
......@@ -1962,6 +1962,13 @@ def pixelSearch(parser):
dest='PHIFILE_BIN_RATIO',
help="Ratio of binning level between loaded Phi file and current calculation. If the input Phi file has been obtained on a 500x500x500 image and now the calculation is on 1000x1000x1000, this should be 2. Default = 1")
parser.add_argument('-pfd',
'--phiFile-direct',
action="store_true",
#default=1,
dest='PHIFILE_DIRECT',
help="Trust the Phi file completely? This option ignores and overrides -pfni and requires same nodes in same positions. Default = False")
parser.add_argument('-pfni',
'--neighbours-for-phi-field-interpolation',
type=int,
......@@ -2101,32 +2108,9 @@ def pixelSearch(parser):
if not os.path.isdir(args.OUT_DIR):
raise
# Catch interdependent node spacing and correlation window sizes
if args.NS is None:
print("\nUsing default node spacing: "),
if args.HWS is None:
print("2x default half window size"),
args.HWS = [10]
print("({}) which is".format(args.HWS[0])),
args.NS = [args.HWS[0] * 2]
else:
print("2x user-set half window size"),
if len(args.HWS) == 1:
print("({}) which is".format(args.HWS[0])),
args.NS = [int(args.HWS[0] * 2)]
elif len(args.HWS) == 3:
print("({} -- selecting smallest) which is".format(args.HWS)),
args.NS = [int(min(args.HWS) * 2)]
print(args.NS)
# Catch 3D options
if len(args.NS) == 1:
args.NS = [args.NS[0], args.NS[0], args.NS[0]]
if len(args.HWS) == 1:
args.HWS = [args.HWS[0], args.HWS[0], args.HWS[0]]
if args.LAB1 is not None:
# We have a labelled image and so no nodeSpacing or halfWindowSize
print("I have been passed a labelled image and so I am disactivating node spacing and half-window size and mask and setting mask coverage to 0")
args.HWS = None
args.NS = None
......@@ -2136,12 +2120,39 @@ def pixelSearch(parser):
if twoD:
args.SEARCH_RANGE[0] = 0
args.SEARCH_RANGE[1] = 0
# Catch and overwrite 2D options
elif twoD:
args.NS[0] = 1
args.HWS[0] = 0
args.SEARCH_RANGE[0] = 0
args.SEARCH_RANGE[1] = 0
else:
# We are in grid, with a nodeSpacing and halfWindowSize
# Catch interdependent node spacing and correlation window sizes
if args.NS is None:
print("\nUsing default node spacing: "),
if args.HWS is None:
print("2x default half window size"),
args.HWS = [10]
print("({}) which is".format(args.HWS[0])),
args.NS = [args.HWS[0] * 2]
else:
print("2x user-set half window size"),
if len(args.HWS) == 1:
print("({}) which is".format(args.HWS[0])),
args.NS = [int(args.HWS[0] * 2)]
elif len(args.HWS) == 3:
print("({} -- selecting smallest) which is".format(args.HWS)),
args.NS = [int(min(args.HWS) * 2)]
print(args.NS)
# Catch 3D options
if len(args.NS) == 1:
args.NS = [args.NS[0], args.NS[0], args.NS[0]]
if len(args.HWS) == 1:
args.HWS = [args.HWS[0], args.HWS[0], args.HWS[0]]
# Catch and overwrite 2D options
if twoD:
args.NS[0] = 1
args.HWS[0] = 0
args.SEARCH_RANGE[0] = 0
args.SEARCH_RANGE[1] = 0
# Output file name prefix
if args.PREFIX is None:
......
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