Vous avez reçu un message "Your GitLab account has been locked ..." ? Pas d'inquiétude : lisez cet article https://docs.gricad-pages.univ-grenoble-alpes.fr/help/unlock/

Commit f94ea529 authored by Edward Andò's avatar Edward Andò
Browse files

Drop opencv dep

parent bef23b7c
-r requirements.txt
coverage
sphinx>=2.0
sphinx-gallery>=0.1.13
......
numpy==1.18.5
scipy==1.4.1
opencv-python==3.4.10.35
SimpleITK==1.2.4
mpi4py==3.0.3
......
......@@ -399,9 +399,48 @@ if mpiRank == boss or not mpi:
# 2020-07-05 try applying F to im1 this is expected to help with pixel searching
PhiNoDisp = PhiField[label].copy()
PhiNoDisp[0:3,-1] = 0.0
imagette1def = spam.DIC.applyPhi(imagette1, PhiNoDisp)
pixelSearchReturns = spam.DIC.correlate.pixelSearch(imagette1def,
# Deform imagette1, using relative COM as point of application
imagette1def = spam.DIC.applyPhi(imagette1, PhiNoDisp, PhiPoint=getLabel['centreOfMassREL'])
imagette1defCrop = imagette1def[args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN,
args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN, # This is to remove edge artefacts of applyPhi
args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN]
imagette1crop = imagette1[args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN,
args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN, # This is to remove edge artefacts of applyPhi
args.GRID_POINT_MARGIN:-args.GRID_POINT_MARGIN]
if args.MASK:
imagette1toCorrelate = imagette1def.copy()
maskette1def = spam.DIC.applyPhi(maskette1, PhiNoDisp, PhiPoint=getLabel['centreOfMassREL'], interpolationOrder=0)
imagette1toCorrelate[maskette1def == 0] = numpy.nan
else:
# add crop??
imagette1toCorrelate = imagette1def.copy()
if args.DEBUG:
import matplotlib.pyplot as plt
print("t before pixel search:\n", spam.deformation.decomposePhi(PhiField[label])['t'])
print("r before pixel search:\n", spam.deformation.decomposePhi(PhiField[label])['r'])
plt.subplot(1,5,1)
plt.title("ref label")
plt.imshow(maskette1[maskette1.shape[0]//2])
plt.subplot(1,5,2)
plt.title("ref image")
plt.imshow(imagette1[imagette1.shape[0]//2])
plt.subplot(1,5,3)
plt.title("ref image deformed (no disp)")
plt.imshow(imagette1def[imagette1def.shape[0]//2])
plt.subplot(1,5,4)
plt.title("ref image deformed (no disp) + mask")
plt.imshow(imagette1toCorrelate[imagette1toCorrelate.shape[0]//2])
plt.subplot(1,5,5)
plt.title("def image search range")
plt.imshow(imagette2[imagette2.shape[0]//2])
plt.show()
pixelSearchReturns = spam.DIC.correlate.pixelSearch(imagette1toCorrelate,
imagette2,
searchCentre = searchCentre,
searchRange = searchRangeForThisLabel)
......@@ -409,9 +448,14 @@ if mpiRank == boss or not mpi:
#print(pixelSearchReturns['transformation']['t'])
#print(pixelSearchReturns['cc'], "\n\n")
# 2020-01-21 EA: Add displacements back in
PhiField[label, 0:3, 3] = pixelSearchReturns['transformation']['t'] + labelDisplacementInt
PSCC[label] = pixelSearchReturns['cc']
if args.DEBUG:
print("Phi after pixel search:\n", PhiField[label])
if args.GRID_POINT:
labelDisplacementInt = numpy.round(PhiField[label][0:3, -1])
......
......@@ -74,8 +74,8 @@ class testAll(unittest.TestCase):
self.assertEqual(numpy.add(stack, -snow).sum(), 0)
# test all cases
spam.helpers.stackToArray("tmp_", stack=range(100), createMask=True, verbose=True)
spam.helpers.stackToArray("tmp_", stack=range(100), createMask=True, erosion=True)
spam.helpers.stackToArray("tmp_", stack=range(100), verbose=True)
spam.helpers.stackToArray("tmp_", stack=range(100), erosion=True)
def test_shift(self):
# create slices from snow.tif, restack
......
......@@ -356,10 +356,10 @@ class testAll(unittest.TestCase):
# put 0 in the middle
centres -= numpy.mean(centres, axis=0)
rMax = numpy.amax(radii)
# pad box size
boxSizeDEM = boxSizeDEM + 5 * rMax
# add half box size to centres
centres += numpy.array(boxSizeDEM)/2.0
boxSizePx = (boxSizeDEM / pixelSize).astype(int)
......@@ -379,7 +379,7 @@ class testAll(unittest.TestCase):
tifffile.imsave(testFolder + "Lab0.tif", labIm0.astype(spam.label.labelType))
#test of rigid translation and rotation
#Create Phi and Apply (25 px displacement on Y-axis, and 5 degree rotation along Z axis)
#Create Phi and Apply (2 px displacement on X-axis, and 5 degree rotation along Z axis)
translationStep1 = [0, 0, 2]
rotationStep1 = [5, 0, 0]
transformation = {'t': translationStep1, 'r': rotationStep1}
......@@ -401,7 +401,6 @@ class testAll(unittest.TestCase):
tifffile.imsave(testFolder + "Step1.tif", boxDeformed.astype('<f4'))
#test of rigid translation and rotation
#Create Phi and Apply (25 px displacement on Y-axis, and 5 degree rotation along Z axis)
translationStep2 = [0, 0, 18]
rotationStep2 = [0, 0, 0]
transformation = {'t': translationStep2, 'r': rotationStep2}
......@@ -471,7 +470,8 @@ class testAll(unittest.TestCase):
testFolder + "Lab0.tif",
testFolder + "Step1.tif",
"-od", testFolder + "",
"-pf", testFolder + "Step0-Step1-discreteDVC.tsv"])
"-pf", testFolder + "Step0-Step1-discreteDVC.tsv",
"-pfd"])
self.assertEqual(exitCode, 0)
#Check number of converged grains
......@@ -579,7 +579,8 @@ class testAll(unittest.TestCase):
testFolder + "Lab0.tif",
testFolder + "Step2.tif",
"-od", testFolder + "",
"-pf", testFolder + "Step0-Step2-discreteDVC.tsv"])
"-pf", testFolder + "Step0-Step2-discreteDVC.tsv",
"-pfd"])
self.assertEqual(exitCode, 0)
# Check TSV for local results
......@@ -606,7 +607,7 @@ class testAll(unittest.TestCase):
self.assertAlmostEqual(rotationStep2[2], averageRot[2], places=0)
#######################################################
### 4. Check -skp
### 4. Check -skp
#######################################################
# Manually modify the status of the last label
TSV['PhiField'][-1][:] = 0
......@@ -622,7 +623,8 @@ class testAll(unittest.TestCase):
TSV2 = spam.helpers.readCorrelationTSV(testFolder + "Step0-Step2-discreteDVC.tsv")
# Compare "deltaPhiNorm" for all but the last one (they should be the same since they were skipped)
for i in range(1, TSV2['PhiField'].shape[0]-1):
self.assertEqual(TSV['deltaPhiNorm'][i], TSV2['deltaPhiNorm'][i])
for u in range(3):
self.assertAlmostEqual(TSV['PhiField'][i][u,-1], TSV2['PhiField'][i][u,-1], places=1)
def test_discreteStrain(self):
# make sure it runs the help without error
......
......@@ -32,7 +32,7 @@ import progressbar
#numpy.set_printoptions(precision=3, suppress=True)
def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=None, maxIterations=25, deltaPhiMin=0.001, updateGradient=False, interpolationOrder=1, interpolator='python', verbose=False, imShowProgress=None, imShowProgressNewFig=False):
def register(im1, im2, im1mask=None, PhiInit=None, PhiRigid=False, PhiInitBinRatio=1.0, margin=None, maxIterations=25, deltaPhiMin=0.001, updateGradient=False, interpolationOrder=1, interpolator='python', verbose=False, imShowProgress=None, imShowProgressNewFig=False):
"""
Perform subpixel image correlation between im1 and im2.
......@@ -65,6 +65,10 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
Initial deformation to apply to im1.
Default = numpy.eye(4), `i.e.`, no transformation
PhiRigid : bool, optional
Run a rigid correlation? Only the rigid part of your PhiInit will be kept.
Default = False
PhiInitBinRatio : float, optional
Change translations in PhiInit, if it's been calculated on a differently-binned image. Default = 1
......@@ -138,10 +142,6 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
Since we're using classical DIC (without multimodal registration) there is no need to update the gradient image,
so this is calculated once at the beginning.
"""
# History:
# 2018-03-20 EA: Trying to make margin a list of 3 values as z, y, x margins in order to try
# clean-ish 2D compatilibility
# Explicitly set input images to floats
im1 = im1.astype('<f4')
im2 = im2.astype('<f4')
......@@ -230,8 +230,8 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
if imShowProgress is not None:
import matplotlib.pyplot as plt
# Plot ranges for signed residual
vmin = -im1crop.max() / 2
vmax = im1crop.max() / 2
vmin = -im1crop.max()
vmax = im1crop.max()
if not imShowProgressNewFig:
if imShowProgress == "Z" or imShowProgress == "z":
plt.subplot(1,3,1)
......@@ -284,6 +284,10 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
# Copying into different variable
# Apply binning on displacement
Phi = PhiInit.copy()
# If we're in rigid mode, keep only translations and rotations for this guess
# If you don't do this it goes mad (i.e., rigid updates to non-rigid guess don't seem to work)
if PhiRigid: Phi = spam.deformation.computeRigidPhi(Phi.copy())
Phi[0:3, -1] *= PhiInitBinRatio
# invert PhiInit to apply it to im2
......@@ -316,21 +320,6 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
#if verbose: print("done ")
return imGradZ, imGradY, imGradX
# Calculate gradient of the non-moving im1
# 2017-11-06 EA: This is going in a try, since it's possible to have an error like this:
# "Shape of array too small to calculate a numerical gradient, "
# ValueError: Shape of array too small to calculate a numerical gradient, at least (edge_order + 1) elements are required.
#try:
#if updateGradient:
#imGradZ, imGradY, imGradX = computeGradient(im2def, twoD)
#else:
#imGradZ, imGradY, imGradX = computeGradient(im1, twoD)
## If gradient calculation failed, set singular to true, means early exit
#except ValueError:
## Override max iteration and also set singular
#maxIterations = 0
#singular = True
# Apply stationary im1 mask
if im1mask is not None:
im1crop[im1mask[crop1] == False] = numpy.nan
......@@ -352,6 +341,7 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
# pbar = progressbar.ProgressBar(widgets=widgets, maxval=numberOfNodes)
pbar.start()
# --- Start Iterations ---
while iterations < maxIterations and deltaPhiNorm > deltaPhiMin:
errorPrev = error
......@@ -409,23 +399,39 @@ def register(im1, im2, im1mask=None, PhiInit=None, PhiInitBinRatio=1.0, margin=N
singular = True
break
deltaPhiNorm = numpy.linalg.norm(deltaPhi)
# If we're doing a rigid registration...
if PhiRigid:
# Add padding zeros
deltaPhi = numpy.hstack([deltaPhi, numpy.zeros(4)]).reshape((4, 4))
# Add padding zeros
deltaPhi = numpy.hstack([deltaPhi, numpy.zeros(4)]).reshape((4, 4))
deltaPhiPlusI = (numpy.eye(4) + deltaPhi)
# Keep only rigid part of deltaPhi
deltaPhiPlusIrigid = spam.deformation.computeRigidPhi(deltaPhiPlusI.copy())
# Apply Delta Phi correction to Phi In Roux X-N paper equation number 11
#Phi = numpy.dot((numpy.eye(4) + deltaPhi), Phi)
Phi = numpy.dot(Phi, (numpy.eye(4) + deltaPhi))
# Subtract I from the rigid dPhi+1, and compute norm only on first 3 rows
# ...basically recompute deltaPhiNorm only on rigid part
deltaPhiNorm = numpy.linalg.norm((deltaPhiPlusIrigid-numpy.eye(4))[0:3].ravel())
# Apply Delta Phi correction to Phi In Roux X-N paper equation number 11
Phi = numpy.dot(Phi, deltaPhiPlusIrigid)
else:
# The general, non-rigid case
deltaPhiNorm = numpy.linalg.norm(deltaPhi)
# Add padding zeros
deltaPhi = numpy.hstack([deltaPhi, numpy.zeros(4)]).reshape((4, 4))
# Update Phi
Phi = numpy.dot(Phi, (numpy.eye(4) + deltaPhi))
# Solve for delta Phi
try: PhiInv = numpy.linalg.inv(Phi.copy())
except numpy.linalg.linalg.LinAlgError:
singular = True
break
# reset im1def as emtpy matrix for deformed image
if interpolator == 'C': im2def = spam.DIC.applyPhi(im2, Phi=PhiInv, interpolationOrder=interpolationOrder)
if interpolator == 'C': im2def = spam.DIC.applyPhi(im2, Phi=PhiInv, interpolationOrder=interpolationOrder)
elif interpolator == 'python': im2def = spam.DIC.applyPhiPython(im2, Phi=PhiInv, interpolationOrder=interpolationOrder)
# Error calculation
......
......@@ -44,7 +44,7 @@ def applyPhi(im, Phi=None, PhiPoint=None, interpolationOrder=1):
i.e., the centre of the image
interpolationOrder : int, optional
Order of image interpolation to use. This value is passed directly to ``scipy.ndimage.map_coordinates`` as "order".
Order of image interpolation to use, options are either 0 (strict nearest neighbour) or 1 (trilinear interpolation)
Default = 1
Returns
......
......@@ -22,7 +22,7 @@ import numpy
import tifffile
def stackToArray(prefix, sufix='.tif', stack=range(10), digits='05d', createMask=False, erosion=False, verbose=False):
def stackToArray(prefix, sufix='.tif', stack=range(10), digits='05d', erosion=False, verbose=False):
"""
Convert of stack of 2D sequential tif images into a 3D array.
......@@ -40,9 +40,6 @@ def stackToArray(prefix, sufix='.tif', stack=range(10), digits='05d', createMask
digits : string, default='05d'
The format (number of digits) of the numbers (add leading zeros).
createMask : bool, default=False
Create a mask while stacking the slices. Has been tested for cylinders.
erosion : bool, default=None
Apply an erosion of 1px to the mask in order to avoid border noise.
......@@ -83,8 +80,8 @@ def stackToArray(prefix, sufix='.tif', stack=range(10), digits='05d', createMask
print('\tStack slice number {}/{} ({s:{d}})'.format(i + 1, len(stack), s=s, d=digits))
im[i, :, :] = tifffile.imread(slice_name)
# we fill the mask
if createMask:
mask[i, :, :] = _mask2D(im[i, :, :], erosion=erosion)
#if createMask:
#mask[i, :, :] = _mask2D(im[i, :, :], erosion=erosion)
return im, mask
......@@ -460,37 +457,37 @@ def _binarisation(im, threshold=0.0, boolean=False, op='>', mask=None):
# private functions
def _mask2D(im, erosion=False, structure=None, ):
"""get contour of 2D image.
"""
import cv2
from scipy import ndimage
# step 2: convert into uint8 if not the case
if im.dtype != 'uint8':
# actually it rescales the image but it doesn't really amtter
im = rescaleToInteger(im, nBytes=1)
# Step 3: ...
blur = cv2.GaussianBlur(im, (5, 5), 0)
_, thresh = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
_, contours, _ = cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)
largest = 0
biggest = []
for contour in contours:
area = cv2.contourArea(contour)
if largest < area:
largest = area
biggest = contour
mask = numpy.zeros(im.shape, dtype='<u1')
cv2.drawContours(mask, [biggest], 0, 1, -1)
# Step 4: apply erosion of the mask (which corresponds to an erosion of the specimen)
if erosion:
mask = ndimage.morphology.binary_erosion(
mask, structure=structure).astype(mask.dtype)
return mask
#def _mask2D(im, erosion=False, structure=None, ):
#"""
#get contour of 2D image.
#"""
#import cv2
#from scipy import ndimage
## step 2: convert into uint8 if not the case
#if im.dtype != 'uint8':
## actually it rescales the image but it doesn't really amtter
#im = rescaleToInteger(im, nBytes=1)
## Step 3: ...
#blur = cv2.GaussianBlur(im, (5, 5), 0)
#_, thresh = cv2.threshold(blur, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
#_, contours, _ = cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)
#largest = 0
#biggest = []
#for contour in contours:
#area = cv2.contourArea(contour)
#if largest < area:
#largest = area
#biggest = contour
#mask = numpy.zeros(im.shape, dtype='<u1')
#cv2.drawContours(mask, [biggest], 0, 1, -1)
## Step 4: apply erosion of the mask (which corresponds to an erosion of the specimen)
#if erosion:
#mask = ndimage.morphology.binary_erosion(
#mask, structure=structure).astype(mask.dtype)
#return mask
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