Coverage for /usr/local/lib/python3.8/site-packages/spam/helpers/histogramTools.py: 97%

1import matplotlib.pyplot as plt 

2import numpy 

3 

4 

5def findHistogramPeaks(image, valley1=20000, valley2=[], phases=2, gaussianFit=False, returnSigma=False, mask=True, showGraph=False, greyRange=[0,65535], bins=256, saveFigPath=None): 

6 """ 

7 This function finds the peaks of the phases of a greylevel image (16bit or 8bit). The peaks are in greylevel units. Minimum number of Phases is 2 and Maximum is 3. 

8 

9 Parameters 

10 ----------- 

11 

12 image : 3D numpy array 

13 

14 valley1 : float 

15 An initial guess (arbitrary) greylevel value between peak phase 1 and peak phase 2, where peak phase 1 < valley1 < peak phase 2. 

16 

17 valley2 : float 

18 An initial guess (arbitrary) greylevel value between peak phase 2 and peak phase 3, where peak phase 2 < valley2 < peak phase 3. 

19 If the image has 2 phases, value 2 is []. 

20 

21 phases : int (2 or 3) 

22 The phases that exist in the greylevel image. 

23 Easy way to determine by looking how many peaks the histogram has. 

24 

25 gaussianFit : bool, optional 

26 Should the peaks be fitted with a Gaussian, or just the bin with the max value returned? 

27 Default = True 

28 

29 returnSigma : bool, optional 

30 Return a list of standard deviations from the Gaussian Fit? 

31 Requires gaussianFit=True 

32 Default= False 

33 

34 mask: bool 

35 If image is masked with mask set to 0. 

36 

37 showGraph: bool, optional 

38 If True, showing two matplotlib graphs - one corresponding to the histogram as returned in the spam.plotting.greyLevelHistogram and one corresponding to the histogram with the peaks of the phases. 

39 Default = False 

40 

41 greyRange : list, optional 

42 If a 16-bit greylevel image --> greyRange = [0,65535] 

43 If a 8-bit greylevel image --> greyRange = [0,255] 

44 Default = [0,65535] 

45 

46 bins : int, optional 

47 Default = 256 

48 

49 saveFigPath : string, optional 

50 Path to save figure to. 

51 Default = None 

52 

53 Returns 

54 -------- 

55 The peaks of the phases of the image. 

56 """ 

57 import scipy.optimize 

58 import spam.plotting.greyLevelHistogram 

59 

60 def gauss(x,a,mu,sigma, offset): 

61 return a*numpy.exp(-1*(x-mu)**2/(2*sigma**2))+offset 

62 

63 def gaussianFitFunction(x1, y1): 

64 # probably there is a smarter way to guess A instead of 1... 

65 mu1 = sum(x1 * y1) / sum(y1) 

66 sigma1 = numpy.sqrt(sum(y1 * (x1 - mu1)**2) / sum(y1)) 

67 popt1,pcov1 = scipy.optimize.curve_fit(gauss, x1, y1, p0=[max(y1),mu1, sigma1, 0],maxfev=1000) 

68 a1 = popt1[0] 

69 mu1 = popt1[1] 

70 s1 = numpy.abs(popt1[2]) 

71 offset1 = popt1[3] 

72 return a1, mu1, s1, offset1 

73 

74 if phases == 1 or phases>=4: 

75 print("spam.helpers.histogramTools.findHistogramPeaks: Need to give me correct number of phases, remember always 2 or 3") 

76 return 

77 

78 if phases==3: 

79 if valley1>=valley2: 

80 print("spam.helpers.histogramTools.findHistogramPeaks: Need to give me the correct values, where always valley1 < valley2") 

81 return 

82 

83 if mask==True: 

84 image=image.astype(float) 

85 image[image==0]=numpy.nan 

86 reshist = spam.plotting.greyLevelHistogram.plotGreyLevelHistogram(image[numpy.isfinite(image)], greyRange=greyRange, bins=bins) 

87 else: 

88 reshist = spam.plotting.greyLevelHistogram.plotGreyLevelHistogram(image, greyRange=greyRange, bins=bins) 

89 

90 

91 totalCounts = numpy.array(reshist[1]) 

92 totalgreylevel = numpy.array(reshist[0]) 

93 

94 if phases==2: 

95 

96 greyPhase1= totalgreylevel[totalgreylevel<=valley1] 

97 countsPhase1 = totalCounts[0:greyPhase1.shape[0]] 

98 peakPhase1 = greyPhase1[numpy.argmax(countsPhase1)] 

99 # peakPhase1 = greyPhase1[countsPhase1==countsPhase1.max()] 

100 greyPhase2= totalgreylevel[totalgreylevel>=valley1] 

101 countsPhase2 = totalCounts[(bins-greyPhase2.shape[0]):bins] 

102 peakPhase2 = greyPhase2[numpy.argmax(countsPhase2)] 

103 # peakPhase2 = greyPhase2[countsPhase2==countsPhase2.max()] 

104 

105 if gaussianFit: 

106 try: #Lets try to make the fitting 

107 #Compute midpoint between two peaks 

108 midPointGrey = (greyPhase1[numpy.argmax(countsPhase1)] + greyPhase2[numpy.argmax(countsPhase2)] ) / 2 

109 midPointArg = numpy.argmin(numpy.abs(totalgreylevel - midPointGrey)) 

110 midPointCount = totalCounts[midPointArg] 

111 #Compute index of peaks on total 

112 indexPeak1 = numpy.where(totalCounts == countsPhase1.max()) 

113 indexPeak2 = numpy.where(totalCounts == countsPhase2.max()) 

114 #Compute midpoint in Y for subsets 

115 midPointCount1 = midPointCount + 0.25*(countsPhase1.max() - midPointCount) 

116 midPointCount2 = midPointCount + 0.25*(countsPhase2.max() - midPointCount) 

117 #Get new subsets index 

118 startIndex1 = numpy.argmin(numpy.abs( totalCounts[0:indexPeak1[0][0]]- midPointCount1)) 

119 stopIndex1 = numpy.argmin(numpy.abs( totalCounts[indexPeak1[0][0]:midPointArg]- midPointCount1)) + indexPeak1[0][0] +1 

120 startIndex2 = numpy.argmin(numpy.abs( totalCounts[midPointArg:indexPeak2[0][0]]- midPointCount2)) + midPointArg 

121 stopIndex2 = numpy.argmin(numpy.abs( totalCounts[indexPeak2[0][0]:]- midPointCount2)) + indexPeak2[0][0] +1 

122 a1, peakPhase1, sigmaPhase1, offset1 = gaussianFitFunction(totalgreylevel[startIndex1:stopIndex1], totalCounts[startIndex1:stopIndex1]) 

123 a2, peakPhase2, sigmaPhase2, offset2 = gaussianFitFunction(totalgreylevel[startIndex2:stopIndex2], totalCounts[startIndex2:stopIndex2]) 

124 sigmas = numpy.array([sigmaPhase1,sigmaPhase2]) 

125 except: 

126 print("spam.helpers.histogramTools.findHistogramPeaks: There is not enough data to make the fitting, aborting Gaussian Fit") 

127 gaussianFit = False 

128 

129 peakPhases = numpy.array([peakPhase1,peakPhase2]) 

130 

131 

132 if showGraph == True: 

133 fig = plt.figure(1) 

134 plt.plot(totalgreylevel,totalCounts,label="Histogram") 

135 

136 plt.plot(peakPhase1,countsPhase1[countsPhase1==countsPhase1.max()],"ro",label="Peak Phase 1 @ {:5.0f}".format(float(peakPhase1))) 

137 plt.plot(peakPhase2,countsPhase2[countsPhase2==countsPhase2.max()],"yo",label="Peak Phase 2 @ {:5.0f}".format(float(peakPhase2))) 

138 plt.xlabel("Greylevel") 

139 plt.ylabel("Counts") 

140 

141 if gaussianFit: 

142 plt.plot(totalgreylevel[startIndex1:stopIndex1], gauss(totalgreylevel[startIndex1:stopIndex1], a1, peakPhase1, sigmaPhase1, offset1), label='Fit 1') 

143 plt.plot(totalgreylevel[startIndex2:stopIndex2], gauss(totalgreylevel[startIndex2:stopIndex2], a2, peakPhase2, sigmaPhase2, offset2), label='Fit 2') 

144 

145 plt.legend() 

146 fig.tight_layout() 

147 if saveFigPath is not None: 

148 plt.savefig(saveFigPath) 

149 else: 

150 plt.show() 

151 

152 plt.close() 

153 

154 

155 if phases==3: 

156 #if gaussianFit: 

157 # print("spam.helpers.findHistogramPeaks(): Three-peak fitting not yet implemented") 

158 greyPhase1= totalgreylevel[totalgreylevel<=valley1] 

159 countsPhase1 = totalCounts[0:greyPhase1.shape[0]] 

160 peakPhase1 = greyPhase1[countsPhase1==countsPhase1.max()] 

161 if gaussianFit: 

162 a1, peakPhase1, sigmaPhase1, offset1 = gaussianFitFunction(greyPhase1, countsPhase1) 

163 

164 greyPhase3= totalgreylevel[totalgreylevel>=valley2] 

165 countsPhase3 = totalCounts[(bins-greyPhase3.shape[0]):bins] 

166 peakPhase3 = greyPhase3[countsPhase3==countsPhase3.max()] 

167 if gaussianFit: 

168 a3, peakPhase3, sigmaPhase3, offset3 = gaussianFitFunction(greyPhase3, countsPhase3) 

169 

170 greyPhase2 = totalgreylevel[(totalgreylevel>valley1)&(totalgreylevel<valley2)] 

171 countsPhase2 = totalCounts[greyPhase1.shape[0]:(bins-greyPhase3.shape[0])] 

172 peakPhase2 = greyPhase2[countsPhase2==countsPhase2.max()] 

173 

174 if gaussianFit: 

175 a2, peakPhase2, sigmaPhase2, offset2 = gaussianFitFunction(greyPhase2, countsPhase2) 

176 

177 #print("spam.helpers.findHistogramPeaks: The peak of the Phase1 is at {:5.0f} of Greylevel".format(float(peakPhase1))) 

178 #print("spam.helpers.findHistogramPeaks: The peak of the Phase2 is at {:5.0f} of Greylevel".format(float(peakPhase2))) 

179 #print("spam.helpers.findHistogramPeaks: The peak of the Phase3 is at {:5.0f} of Greylevel".format(float(peakPhase3))) 

180 

181 peakPhases = numpy.array([peakPhase1,peakPhase2,peakPhase3]) 

182 if gaussianFit: 

183 sigmas = numpy.array([sigmaPhase1,sigmaPhase2,sigmaPhase3]) 

184 

185 if showGraph == True: 

186 fig = plt.figure() 

187 plt.plot(totalgreylevel,totalCounts,label="Histogram") 

188 plt.plot(peakPhase1,countsPhase1[countsPhase1==countsPhase1.max()],"ro",label="Peak Phase 1 @ {:5.0f}".format(float(peakPhase1))) 

189 plt.plot(peakPhase2,countsPhase2[countsPhase2==countsPhase2.max()],"yo",label="Peak Phase 2 @ {:5.0f}".format(float(peakPhase2))) 

190 plt.plot(peakPhase3,countsPhase3[countsPhase3==countsPhase3.max()],"bo",label="Peak Phase 3 @ {:5.0f}".format(float(peakPhase3))) 

191 plt.xlabel("Greylevel") 

192 plt.ylabel("Counts") 

193 if gaussianFit: 

194 plt.plot(greyPhase1, gauss(greyPhase1, a1, peakPhase1, sigmaPhase1), label='Fit 1') 

195 plt.plot(greyPhase2, gauss(greyPhase2, a2, peakPhase2, sigmaPhase2), label='Fit 2') 

196 plt.plot(greyPhase3, gauss(greyPhase3, a3, peakPhase3, sigmaPhase3), label='Fit 3') 

197 plt.legend() 

198 fig.tight_layout() 

199 plt.show() 

200 

201 if returnSigma and not gaussianFit: 

202 print("spam.helpers.histogramTools.findHistogramPeaks: cannot return sigma if fitGaussian is not activated") 

203 if returnSigma and gaussianFit: 

204 return peakPhases, sigmas 

205 else: 

206 return peakPhases 

207 

208 

209 

210def histogramNorm(im_Or, twoPeaks, peaksNormed=[0.25, 0.75], cropGreyvalues=[-numpy.inf, numpy.inf]): 

211 """ 

212 This function normalise the histogram in order to range beween 0 and 1, presenting two peaks at p1 and p2 (p1<p2) 

213 

214 

215 Parameters 

216 ----------- 

217 

218 im_Or : 3D numpy array 

219 

220 twoPeaks : list of two floats 

221 First and second peak of the original histogram 

222 

223 peaksNormed : list of two floats, optional 

224 The desired level for the first and second peak of the normalized histogram. 

225 Default = [0.25, 0.75] 

226 

227 cropGreyvalues : list of two floats, optional 

228 The limits on the generated normalised values. 

229 Default = [-numpy.inf, numpy.inf] 

230 

231 Returns 

232 -------- 

233 im_Norm : 3D numpy array. Image with grey range between [0,1] and peaks at 0.25 and 0.75 

234 

235 

236 """ 

237 if len(twoPeaks) != 2: 

238 print("spam.helpers.histogramTools.histogramTools.histogramNorm: The number of peaks should be 2") 

239 return 

240 peak1 = twoPeaks[0] 

241 peak2 = twoPeaks[1] 

242 p1 = peaksNormed[0] 

243 p2 = peaksNormed[1] 

244 if p1 > p2: 

245 print("spam.helpers.histogramTools.histogramTools.histogramNorm: p1 should be less than p2") 

246 return 

247 if peak1 > peak2: 

248 print("spam.helpers.histogramTools.histogramTools.histogramNorm: peak1 should be less than peak2") 

249 return 

250 #if peaksNormed is None: 

251 # peaksNormed = [0.25, 0.75] 

252 

253 m = (p2-p1)/(peak2 - peak1) 

254 b = p1 - m*peak1 

255 im_Norm = m*im_Or + b 

256 im_Norm[im_Norm<cropGreyvalues[0]]=cropGreyvalues[0] 

257 im_Norm[im_Norm>cropGreyvalues[1]]=cropGreyvalues[1] 

258 

259 return im_Norm 

260