我正在尝试基于像素强度值对2D MR图像的不同区域进行自动图像分割。第一步是在图像的直方图上实现高斯混合模型。
我需要将从该score_samples方法获得的结果高斯绘制到直方图上。我已尝试按照答案(了解高斯混合模型)中的代码进行操作。
score_samples
但是,所得的高斯根本无法匹配直方图。如何使高斯与直方图匹配?
import numpy as np import cv2 import matplotlib.pyplot as plt from sklearn.mixture import GaussianMixture # Read image img = cv2.imread("test.jpg",0) hist = cv2.calcHist([img],[0],None,[256],[0,256]) hist[0] = 0 # Removes background pixels # Fit GMM gmm = GaussianMixture(n_components = 3) gmm = gmm.fit(hist) # Evaluate GMM gmm_x = np.linspace(0,255,256) gmm_y = np.exp(gmm.score_samples(gmm_x.reshape(-1,1))) # Plot histograms and gaussian curves fig, ax = plt.subplots() ax.hist(img.ravel(),255,[1,256]) ax.plot(gmm_x, gmm_y, color="crimson", lw=4, label="GMM") ax.set_ylabel("Frequency") ax.set_xlabel("Pixel Intensity") plt.legend() plt.show()
我还尝试用总和手动构造高斯。
import numpy as np import cv2 import matplotlib.pyplot as plt from sklearn.mixture import GaussianMixture def gauss_function(x, amp, x0, sigma): return amp * np.exp(-(x - x0) ** 2. / (2. * sigma ** 2.)) # Read image img = cv2.imread("test.jpg",0) hist = cv2.calcHist([img],[0],None,[256],[0,256]) hist[0] = 0 # Removes background pixels # Fit GMM gmm = GaussianMixture(n_components = 3) gmm = gmm.fit(hist) # Evaluate GMM gmm_x = np.linspace(0,255,256) gmm_y = np.exp(gmm.score_samples(gmm_x.reshape(-1,1))) # Construct function manually as sum of gaussians gmm_y_sum = np.full_like(gmm_x, fill_value=0, dtype=np.float32) for m, c, w in zip(gmm.means_.ravel(), gmm.covariances_.ravel(), gmm.weights_.ravel()): gauss = gauss_function(x=gmm_x, amp=1, x0=m, sigma=np.sqrt(c)) gmm_y_sum += gauss / np.trapz(gauss, gmm_x) * w # Plot histograms and gaussian curves fig, ax = plt.subplots() ax.hist(img.ravel(),255,[1,256]) ax.plot(gmm_x, gmm_y, color="crimson", lw=4, label="GMM") ax.plot(gmm_x, gmm_y_sum, color="black", lw=4, label="Gauss_sum", linestyle="dashed") ax.set_ylabel("Frequency") ax.set_xlabel("Pixel Intensity") plt.legend() plt.show()
用 ax.hist(img.ravel(),255,[1,256], normed=True)
ax.hist(img.ravel(),255,[1,256], normed=True)
问题在于将直方图而不是像素强度数组传递给GaussianMixture.fit gmm = gmm.fit(hist)。我还发现,n_components = 6视觉上适合此特定直方图的最低要求。
gmm = gmm.fit(hist)
n_components = 6
import numpy as np import cv2 import matplotlib.pyplot as plt from sklearn.mixture import GaussianMixture # Read image img = cv2.imread("test.jpg",0) hist = cv2.calcHist([img],[0],None,[256],[0,256]) hist[0] = 0 # Removes background pixels data = img.ravel() data = data[data != 0] data = data[data != 1] #Removes background pixels (intensities 0 and 1) # Fit GMM gmm = GaussianMixture(n_components = 6) gmm = gmm.fit(X=np.expand_dims(data,1)) # Evaluate GMM gmm_x = np.linspace(0,253,256) gmm_y = np.exp(gmm.score_samples(gmm_x.reshape(-1,1))) # Plot histograms and gaussian curves fig, ax = plt.subplots() ax.hist(img.ravel(),255,[2,256], normed=True) ax.plot(gmm_x, gmm_y, color="crimson", lw=4, label="GMM") ax.set_ylabel("Frequency") ax.set_xlabel("Pixel Intensity") plt.legend() plt.show()