import matplotlib.pyplot as plt
import numpy
[docs]
def gauss(x,a,mu,sigma, offset):
return a*numpy.exp(-1*(x-mu)**2/(2*sigma**2))+offset
[docs]
def gaussianFitFunction(x1, y1):
import scipy.optimize
import warnings
from scipy.optimize import OptimizeWarning
warnings.filterwarnings("ignore", category=OptimizeWarning)
# probably there is a smarter way to guess A instead of 1...
mu1 = sum(x1 * y1) / sum(y1)
sigma1 = numpy.sqrt(sum(y1 * (x1 - mu1)**2) / sum(y1))
popt1, pcov1 = scipy.optimize.curve_fit(gauss, x1, y1, p0=[max(y1),mu1, sigma1, 0], maxfev=1000)
a1 = popt1[0]
mu1 = popt1[1]
s1 = numpy.abs(popt1[2])
offset1 = popt1[3]
return a1, mu1, s1, offset1
[docs]
def findHistogramPeaks(image, valley1=None, valley2=None, mask=None, greyRange=None, bins=256, gaussianFit=False):
"""
This function finds the peaks of the 2 or 3 phases of a greylevel image.
Parameters
-----------
image: 3D numpy array
The greyscale image.
valley1 : float (default=None)
A arbitrary greylevel value between peak of phase 1 and phase 2
If None then uses Otsu thresholding method to determine the valley.
valley2 : float (default=None)
A arbitrary greylevel value between peak of phase 2 and phase 3
If None, only looks for 2 phases.
mask: 3D boolean numpy array, optional
Ignore masked grey values
greyRange : list (default=None)
Range from which the histograms are computed.
If None it is guessed from the image type.
bins : int (default=256)
Number of bins used for computing the histogram.
gaussianFit : bool (default=False)
Finds the peaks with a Gaussian fit.
Returns
--------
1D numpy array:
The peaks of the phases of the image (size 2 or 3).
1D numpy array:
Bins of the histogram
1D numpy array:
Counts of the histogram
"""
import spam.plotting.greyLevelHistogram
import skimage.filters
if valley2 and valley1 and valley1>=valley2:
raise ValueError(f"valley2 > valley1 in findHistogramPeaks")
# determine valley1 if not set
valley1 = skimage.filters.threshold_otsu(image) if not valley1 else float(valley1)
if mask is not None:
image = image.astype('<f4')
image[mask==0] = numpy.nan
reshist = spam.plotting.greyLevelHistogram.plotGreyLevelHistogram(image[numpy.isfinite(image)], greyRange=greyRange, bins=bins)
else:
reshist = spam.plotting.greyLevelHistogram.plotGreyLevelHistogram(image, greyRange=greyRange, bins=bins)
# rehist = [midBins, binLimits, counts]
totalCounts = numpy.array(reshist[2])
totalBins = numpy.array(reshist[1])
totalMidBins = numpy.array(reshist[0])
if valley2 is None:
greyPhase1= totalMidBins[totalMidBins<=valley1]
countsPhase1 = totalCounts[0:greyPhase1.shape[0]]
peakPhase1 = greyPhase1[numpy.argmax(countsPhase1)]
greyPhase2= totalMidBins[totalMidBins>=valley1]
countsPhase2 = totalCounts[(bins-greyPhase2.shape[0]):bins]
peakPhase2 = greyPhase2[numpy.argmax(countsPhase2)]
if gaussianFit:
try: #Lets try to make the fitting
#Compute midpoint between two peaks
midPointGrey = (greyPhase1[numpy.argmax(countsPhase1)] + greyPhase2[numpy.argmax(countsPhase2)] ) / 2
midPointArg = numpy.argmin(numpy.abs(totalMidBins - midPointGrey))
midPointCount = totalCounts[midPointArg]
#Compute index of peaks on total
indexPeak1 = numpy.where(totalCounts == countsPhase1.max())
indexPeak2 = numpy.where(totalCounts == countsPhase2.max())
#Compute midpoint in Y for subsets
midPointCount1 = midPointCount + 0.25*(countsPhase1.max() - midPointCount)
midPointCount2 = midPointCount + 0.25*(countsPhase2.max() - midPointCount)
#Get new subsets index
startIndex1 = numpy.argmin(numpy.abs( totalCounts[0:indexPeak1[0][0]]- midPointCount1))
stopIndex1 = numpy.argmin(numpy.abs( totalCounts[indexPeak1[0][0]:midPointArg]- midPointCount1)) + indexPeak1[0][0] +1
startIndex2 = numpy.argmin(numpy.abs( totalCounts[midPointArg:indexPeak2[0][0]]- midPointCount2)) + midPointArg
stopIndex2 = numpy.argmin(numpy.abs( totalCounts[indexPeak2[0][0]:]- midPointCount2)) + indexPeak2[0][0] +1
a1, peakPhase1, sigmaPhase1, offset1 = gaussianFitFunction(totalMidBins[startIndex1:stopIndex1], totalCounts[startIndex1:stopIndex1])
a2, peakPhase2, sigmaPhase2, offset2 = gaussianFitFunction(totalMidBins[startIndex2:stopIndex2], totalCounts[startIndex2:stopIndex2])
sigmas = numpy.array([sigmaPhase1,sigmaPhase2])
except TypeError:
print("spam.helpers.findHistogramPeaks: Scipy optimisation for Gaussian Fit did not work")
pass
peakPhases = numpy.array([peakPhase1, peakPhase2])
else:
greyPhase1= totalMidBins[totalMidBins<=valley1]
countsPhase1 = totalCounts[0:greyPhase1.shape[0]]
peakPhase1 = greyPhase1[countsPhase1==countsPhase1.max()][0]
if gaussianFit:
a1, peakPhase1, sigmaPhase1, offset1 = gaussianFitFunction(greyPhase1, countsPhase1)
greyPhase3= totalMidBins[totalMidBins>=valley2]
countsPhase3 = totalCounts[(bins-greyPhase3.shape[0]):bins]
peakPhase3 = greyPhase3[countsPhase3==countsPhase3.max()][0]
if gaussianFit:
a3, peakPhase3, sigmaPhase3, offset3 = gaussianFitFunction(greyPhase3, countsPhase3)
greyPhase2 = totalMidBins[(totalMidBins>valley1)&(totalMidBins<valley2)]
countsPhase2 = totalCounts[greyPhase1.shape[0]:(bins-greyPhase3.shape[0])]
peakPhase2 = greyPhase2[countsPhase2==countsPhase2.max()][0]
if gaussianFit:
a2, peakPhase2, sigmaPhase2, offset2 = gaussianFitFunction(greyPhase2, countsPhase2)
peakPhases = numpy.array([peakPhase1,peakPhase2,peakPhase3])
if gaussianFit:
sigmas = numpy.array([sigmaPhase1,sigmaPhase2,sigmaPhase3])
return peakPhases, totalBins, totalCounts
[docs]
def histogramNorm(im, twoPeaks, peaksNormed=[0.25, 0.75], cropGreyvalues=[-numpy.inf, numpy.inf]):
"""
This function normalise the histogram in order to range beween 0 and 1, presenting two peaks at p1 and p2 (p1<p2)
Parameters
-----------
im : 3D numpy array
The image to normalise
twoPeaks : list of two floats
First and second peak of the original histogram
peaksNormed : list of two floats (default=[0.25, 0.75]
The desired level for the first and second peak of the normalized histogram.
cropGreyvalues : list of two floats, optional
The limits on the generated normalised values.
Default = [-numpy.inf, numpy.inf]
Returns
--------
3D numpy array:
Normalised image
"""
if len(twoPeaks) != 2:
raise ValueError("spam.helpers.histogramNorm: The number of peaks should be 2")
peak1 = twoPeaks[0]
peak2 = twoPeaks[1]
p1 = peaksNormed[0]
p2 = peaksNormed[1]
if p1 > p2:
raise ValueError("spam.helpers.histogramNorm: p1 should be less than p2")
if peak1 > peak2:
raise ValueError("spam.helpers.histogramNorm: peak1 should be less than peak2")
m = (p2 - p1) / (peak2 - peak1)
b = p1 - m * peak1
imNormed = m * im + b
imNormed[imNormed<cropGreyvalues[0]]=cropGreyvalues[0]
imNormed[imNormed>cropGreyvalues[1]]=cropGreyvalues[1]
return imNormed
