Commit 7b83a76b authored by Edward Andò's avatar Edward Andò
Browse files

Small cosmetic changes to globalDVC, without BC we now do numpy.solve()

parent 3eac8fc9
Pipeline #46383 passed with stages
in 29 minutes and 18 seconds
......@@ -203,9 +203,9 @@ if args.STRAIN:
print("spam-gdic: Computing strains", end='...')
Ffield = spam.deformation.FfieldBagi(inputMesh['points'], inputMesh['tetra'], displacements)
if (not args.SMALL_STRAINS):
components = ['vol', 'dev', 'U']
components = ['vol', 'dev']
else:
components = ['volss', 'devss', 'e']
components = ['volss', 'devss']
decomposedFfield = spam.deformation.decomposeFfield(Ffield, components)
for component in components:
......
......@@ -713,6 +713,49 @@ def pixelSearch(imagette1, imagette2, searchCentre=None, searchRange={'zRange':
def globalCorrelation(im1, im2, mesh, convergenceCriterion=0.01, debugFiles=False, maxIterations=20, initialDisplacements=None, boundaryConditions=None, prefix="./"):
"""
Global DVC (works only with 3D images).
Parameters
----------
im1 : 3D numpy array
Reference image in which the mesh is defined
im2 : 3D numpy array
Deformed image, should be same size as im1
mesh : dictionary containing mesh elements:
-
convergenceCriterion : float
Convergence criterion for change in displacements in px
Default = 0.01
debugFiles : bool
Write temporary results to file for debugging?
Default = False
maxIterations : int
Number of iterations to stop after if convergence has not been reached
Default = 20
initialDisplacements : Nx3 numpy array of floats
Initial guess for nodal displacements, must be coherent with input mesh
Default = None
boundaryConditions : Nx3 numpy array of floats
Nodal boundary conditions
Default = None
prefix : string
Output file prefix for debugFiles
Default = "./"
Returns
-------
displacements : Nx3 numpy array of floats
Nodal displacements
"""
import spam.helpers
print("gdic: convergenceCriterion = {}".format(convergenceCriterion))
print("gdic: maxIterations = {}".format(maxIterations))
......@@ -801,16 +844,16 @@ def globalCorrelation(im1, im2, mesh, convergenceCriterion=0.01, debugFiles=Fals
else:
K11 = globalMatrix.copy()
print("\tgdic: Inversing global matrix")
# print globalMatrix, "\n\n\n\n\n"
# + 0.0*numpy.eye( globalMatrix.shape[0] ) )
# K11 = globalMatrix.copy()
if applyBC:
print("\tgdic: Inversing global matrix...")
tmp = K11 - numpy.dot(K12, numpy.dot(K22i, K21))
globalMatrixInv = numpy.linalg.inv(tmp)
else:
globalMatrixInv = numpy.linalg.inv(K11)
# print globalMatrixInv
print("\tgdic: NOT Inversing global matrix...")
#globalMatrixInv = numpy.linalg.inv(K11)
#def errorCalc(im1, im2, im2ref, meshPaddingSlice):
#errorInitial = numpy.sqrt(numpy.square(im2ref[meshPaddingSlice] - im1[meshPaddingSlice]).sum())
......@@ -846,7 +889,9 @@ def globalCorrelation(im1, im2, mesh, convergenceCriterion=0.01, debugFiles=Fals
displacements[nbc] += dx
else:
dx = numpy.dot(globalMatrixInv, globalVector).astype('<f8')
# Just try to solve directly
#dx = numpy.dot(globalMatrixInv, globalVector).astype('<f8')
dx = numpy.linalg.solve(K11, globalVector).astype('<f8')
displacements += dx
dxNorm = numpy.linalg.norm(dx)
......
......@@ -1045,7 +1045,7 @@ def gdicParser(parser):
type=numpy.float,
default=None,
dest='MESH_CUBOID',
help='7 floats parameters to mesh a cuboid. first 3 are position of the bottom corner (Z, Y X), next 3 are the lenghts of the cube (Z, Y X) and the last one is the average length between two nodes.')
help='7 floats parameters to mesh a cuboid. first 3 are position of the bottom corner (Z, Y, X), next 3 are the lenghts of the cube (Z, Y, X) and the last one is the average length between two nodes.')
args = parser.parse_args()
......
Supports Markdown
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