Coverage for /usr/local/lib/python3.10/site-packages/spam/DIC/ddic.py: 100%

1#!/usr/bin/env python 

2 

3# Copyright (C) 2020 SPAM Contributors 

4# 

5# This program is free software: you can redistribute it and/or modify it 

6# under the terms of the GNU General Public License as published by the Free 

7# Software Foundation, either version 3 of the License, or (at your option) 

8# any later version. 

9# 

10# This program is distributed in the hope that it will be useful, but WITHOUT 

11# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 

12# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for 

13# more details. 

14# 

15# You should have received a copy of the GNU General Public License along with 

16# this program. If not, see <http://www.gnu.org/licenses/>. 

17 

18import multiprocessing 

19import os 

20 

21import numpy 

22import progressbar 

23import spam.deformation 

24import spam.DIC 

25import spam.helpers 

26 

27os.environ["OPENBLAS_NUM_THREADS"] = "1" 

28numpy.seterr(all="ignore") 

29 

30 

31def ddic( 

32 im1, 

33 im2, 

34 lab1, 

35 labelsToCorrelate=None, 

36 PhiField=None, 

37 boundingBoxes=None, 

38 centresOfMass=None, 

39 processes=None, 

40 labelDilate=1, 

41 margin=5, 

42 maskOthers=True, 

43 volThreshold=100, 

44 multiScaleBin=1, 

45 updateGrad=False, 

46 correlateRigid=True, 

47 maxIter=50, 

48 deltaPhiMin=0.001, 

49 interpolationOrder=1, 

50 debug=False, 

51 twoD=False, 

52): 

53 numberOfLabels = lab1.max() 

54 

55 if labelsToCorrelate is None: 

56 labelsToCorrelate = numpy.arange(1, numberOfLabels) 

57 

58 PhiFieldOut = numpy.zeros([numberOfLabels, 4, 4]) 

59 if PhiField is None: 

60 for i in range(numberOfLabels): 

61 PhiFieldOut[i] = numpy.eye(4) 

62 else: 

63 PhiFieldOut = PhiField.copy() 

64 

65 if boundingBoxes is None: 

66 boundingBoxes = spam.label.boundingBoxes(lab1) 

67 

68 if centresOfMass is None: 

69 centresOfMass = spam.label.centresOfMass(lab1) 

70 

71 if processes is None: 

72 processes = multiprocessing.cpu_count() 

73 

74 if debug: 

75 print("spam.DIC.ddic(): I was passed debug=True, so setting #processes to 1") 

76 processes = 1 

77 

78 numberOfLabels = (lab1.max() + 1).astype("u4") 

79 PSCC = numpy.zeros((numberOfLabels), dtype="<f4") 

80 error = numpy.zeros((numberOfLabels), dtype="<f4") 

81 iterations = numpy.zeros((numberOfLabels), dtype="<u2") 

82 returnStatus = numpy.zeros((numberOfLabels), dtype="<i2") 

83 deltaPhiNorm = numpy.zeros((numberOfLabels), dtype="<f4") 

84 labelDilateList = numpy.zeros((numberOfLabels), dtype="<u2") 

85 

86 global _multiprocessingCorrelateOneLabel 

87 

88 def _multiprocessingCorrelateOneLabel(label): 

89 # label, labelDilateCurrent = q.get() 

90 

91 # WARNING HACK BAD FIXME 

92 labelDilateCurrent = labelDilate 

93 initialDisplacement = numpy.round(PhiFieldOut[label][0:3, 3]).astype(int) 

94 

95 if debug: 

96 print("\n\n\nWorking on label:", label, "\n") 

97 if debug: 

98 print("Position (ZYX):", centresOfMass[label]) 

99 

100 imagetteReturns = spam.label.getImagettesLabelled( 

101 lab1, 

102 label, 

103 PhiFieldOut[label], 

104 im1, 

105 im2, 

106 [0, 0, 0, 0, 0, 0], # Search range, don't worry about it 

107 boundingBoxes, 

108 centresOfMass, 

109 margin=margin, 

110 labelDilate=labelDilateCurrent, 

111 maskOtherLabels=maskOthers, 

112 applyF="no", 

113 volumeThreshold=volThreshold, 

114 ) 

115 

116 if twoD: 

117 imagetteReturns["imagette1"] = imagetteReturns["imagette1"][int(imagetteReturns["imagette1"].shape[0] - 1) // 2, :, :] 

118 imagetteReturns["imagette2"] = imagetteReturns["imagette2"][int(imagetteReturns["imagette2"].shape[0] - 1) // 2, :, :] 

119 imagetteReturns["imagette1mask"] = imagetteReturns["imagette1mask"][int(imagetteReturns["imagette1mask"].shape[0] - 1) // 2, :, :] 

120 

121 badPhi = numpy.eye(4) 

122 badPhi[0:3, 3] = numpy.nan 

123 

124 # In case the label is missing or the Phi is duff 

125 if imagetteReturns["returnStatus"] != 1 or not numpy.all(numpy.isfinite(PhiFieldOut[label])): 

126 return label, badPhi, -7, numpy.inf, 0, numpy.inf, labelDilateCurrent 

127 

128 else: 

129 # Remove int() part of displacement since it's already used to extract imagette2 

130 PhiTemp = PhiFieldOut[label].copy() 

131 PhiTemp[0:3, -1] -= initialDisplacement 

132 if debug: 

133 print("Starting lk iterations with Phi - int(disp):\n", PhiTemp) 

134 if debug: 

135 print("\nStarting lk iterations with int(disp):\n", initialDisplacement) 

136 

137 registerReturns = spam.DIC.registerMultiscale( 

138 imagetteReturns["imagette1"], 

139 imagetteReturns["imagette2"], 

140 multiScaleBin, 

141 im1mask=imagetteReturns["imagette1mask"], 

142 margin=1, 

143 PhiInit=PhiTemp, 

144 PhiRigid=correlateRigid, 

145 updateGradient=updateGrad, 

146 maxIterations=maxIter, 

147 deltaPhiMin=deltaPhiMin, 

148 interpolationOrder=interpolationOrder, 

149 verbose=debug, 

150 imShowProgress=debug, 

151 ) 

152 goodPhi = registerReturns["Phi"] 

153 goodPhi[0:3, -1] += initialDisplacement 

154 return ( 

155 label, 

156 goodPhi, 

157 registerReturns["returnStatus"], 

158 registerReturns["error"], 

159 registerReturns["iterations"], 

160 registerReturns["deltaPhiNorm"], 

161 labelDilateCurrent, 

162 ) 

163 

164 print(f"\n\tStarting discrete dic for {len(labelsToCorrelate)} labels (with {processes} process{'es' if processes > 1 else ''})") 

165 widgets = [ 

166 progressbar.FormatLabel(""), 

167 " ", 

168 progressbar.Bar(), 

169 " ", 

170 progressbar.AdaptiveETA(), 

171 ] 

172 pbar = progressbar.ProgressBar(widgets=widgets, maxval=len(labelsToCorrelate)) 

173 pbar.start() 

174 

175 finishedLabels = 0 

176 

177 with multiprocessing.Pool(processes=processes) as pool: 

178 for returns in pool.imap_unordered(_multiprocessingCorrelateOneLabel, labelsToCorrelate): 

179 finishedLabels += 1 

180 

181 # Update progres bar if point is not skipped 

182 if returns[2] > 0: 

183 widgets[0] = progressbar.FormatLabel(" it={:0>3d} dPhiNorm={:0>6.4f} rs={:+1d} ".format(returns[4], returns[5], returns[2])) 

184 pbar.update(finishedLabels) 

185 label = returns[0] 

186 PhiFieldOut[label] = returns[1] 

187 returnStatus[label] = returns[2] 

188 error[label] = returns[3] 

189 iterations[label] = returns[4] 

190 deltaPhiNorm[label] = returns[5] 

191 labelDilateList[label] = returns[6] 

192 

193 pbar.finish() 

194 

195 print("\n") 

196 

197 return PhiFieldOut, returnStatus, error, iterations, deltaPhiNorm, labelDilateList, PSCC