Merge pull request #1 from CUQI-DTU/dev

Dev

EITLib/KTCRegularization_NLOpt.py

@@ -0,0 +1,118 @@
+#%%
+import numpy as np
+import matplotlib.pyplot as plt
+
+class SMPrior:
+    def __init__(self, ginv, corrlength, var, mean, covariancetype=None):
+        self.corrlength = corrlength
+        self.mean = mean
+        self.c = 1e-9  # default value
+        if covariancetype is not None:
+            self.covariancetype = covariancetype
+        else:
+            self.covariancetype = 'Squared Distance'  # default
+        self.compute_L(ginv, corrlength, var)
+
+    def compute_L(self, g, corrlength, var):
+        ng = g.shape[0]
+        a = var - self.c
+        b = np.sqrt(-corrlength**2 / (2 * np.log(0.01)))
+        Gamma_pr = np.zeros((ng, ng))
+
+        for ii in range(ng):
+            for jj in range(ii, ng):
+                dist_ij = np.linalg.norm(g[ii, :] - g[jj, :])
+                if self.covariancetype == 'Squared Distance':
+                    gamma_ij = a * np.exp(-dist_ij**2 / (2 * b**2))
+                elif self.covariancetype == 'Ornstein-Uhlenbeck':
+                    gamma_ij = a * np.exp(-dist_ij / corrlength)
+                else:
+                    raise ValueError('Unrecognized prior covariance type')
+                if ii == jj:
+                    gamma_ij = gamma_ij + self.c
+                Gamma_pr[ii, jj] = gamma_ij
+                Gamma_pr[jj, ii] = gamma_ij
+
+
+        self.L = np.linalg.cholesky(np.linalg.inv(Gamma_pr)).T
+
+    def draw_samples(self, nsamples):
+        samples = self.mean + np.linalg.solve(self.L, np.random.randn(self.L.shape[0], nsamples))
+        return samples
+
+    def eval_fun(self, args):
+        sigma = args[0]
+        res = 0.5 * np.linalg.norm(self.L @ (sigma - self.mean))**2
+        return res
+
+    def evaluate_target_external(self, x, compute_grad=False):
+        x = x.reshape((-1,1))
+        # print("x.shape: ", x.shape)
+        # print("self.mean.shape: ", self.mean.shape)
+        if compute_grad:
+            grad = self.L.T @ self.L @ (x - self.mean)
+        else:
+            grad = None
+
+        return self.eval_fun(x), grad
+
+
+    def compute_hess_and_grad(self, args, nparam):
+        sigma = args[0]
+        Hess = self.L.T @ self.L
+        grad = Hess @ (sigma - self.mean)
+
+        if nparam > len(sigma):
+            Hess = np.block([[Hess, np.zeros((len(sigma), nparam - len(sigma)))],
+                             [np.zeros((nparam - len(sigma), len(sigma))), np.zeros((nparam - len(sigma), nparam - len(sigma)))]])
+            grad = np.concatenate([grad, np.zeros(nparam - len(sigma))])
+
+        return Hess, grad
+
+
+if __name__ ==  '__main__':
+    from utils import EITFenics, create_disk_mesh
+    from dolfin import *
+    import pickle
+    L = 32
+    F = 25
+    n = 300
+    radius = 0.115
+    mesh = create_disk_mesh(radius, n, F)
+    myeit = EITFenics(mesh, L, background_conductivity=0.8)
+    H = FunctionSpace(myeit.mesh, 'CG', 1)
+
+    plot(myeit.mesh)
+
+    v2d = vertex_to_dof_map(H)
+    d2v = dof_to_vertex_map(H)
+
+    sigma0 = 0.8*np.ones((myeit.mesh.num_vertices(), 1)) #linearization point
+    corrlength =  radius#* 0.115 #used in the prior
+    var_sigma = 0.05 ** 2 #prior variance
+    mean_sigma = sigma0
+    smprior = SMPrior(myeit.mesh.coordinates()[d2v], corrlength, var_sigma, mean_sigma)
+
+
+    sample = smprior.draw_samples(1)
+    fun = Function(H)
+    fun.vector().set_local(sample)
+    im = plot(fun)
+    plt.colorbar(im)
+
+    mesh_file =XDMFFile('mesh_file_'+str(L)+'_'+str(n)+'.xdmf')
+    mesh_file.write(myeit.mesh)
+    mesh_file.close()
+    #%%
+    file = open('smprior_'+str(L)+'_'+str(n)+'.p', 'wb')
+    pickle.dump(smprior, file)
+    file.close()
+
+
+
+
+
+
+
+
+# %%

EITLib/KTCScoring.py

@@ -0,0 +1,129 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import scipy.signal as sps
+
+def Otsu(image, nvals, figno):
+    # binary Otsu's method for finding the segmentation level for sigma
+    histogramCounts, x = np.histogram(image.ravel(), nvals)
+    # plt.figure(figno)
+    # plt.clf()
+    # plt.hist(image.ravel(), 256)
+    # plt.hold(True)
+
+    total = np.sum(histogramCounts)
+    top = 256
+    sumB = 0
+    wB = 0
+    maximum = 0.0
+    sum1 = np.dot(np.arange(top), histogramCounts)
+    for ii in range(1, top):
+        wF = total - wB
+        if wB > 0 and wF > 0:
+            mF = (sum1 - sumB) / wF
+            val = wB * wF * (((sumB / wB) - mF) ** 2)
+            if val >= maximum:
+                level = ii
+                maximum = val
+        wB = wB + histogramCounts[ii]
+        sumB = sumB + (ii - 1) * histogramCounts[ii]
+
+    # plt.plot([x[level]] * 2, [0, np.max(histogramCounts)], linewidth=2, color='r')
+    # plt.title('histogram of image pixels')
+    # plt.gcf().set_size_inches(9, 5)
+    # plt.show()
+
+    return level, x
+
+def Otsu2(image, nvals, figno):
+    # three class Otsu's method to find the semgentation point of sigma
+    histogramCounts, tx = np.histogram(image.ravel(), nvals)
+    x = (tx[0:-1] + tx[1:])/2
+    # plt.figure(figno)
+    # plt.clf()
+    # plt.stairs(histogramCounts, tx)
+    # plt.hold(True)
+
+    #total = np.sum(histogramCounts)
+    #top = 256
+    maximum = 0.0
+    muT = np.dot(np.arange(1, nvals+1), histogramCounts) / np.sum(histogramCounts)
+    for ii in range(1, nvals):
+        for jj in range(1, ii):
+            w1 = np.sum(histogramCounts[:jj])
+            w2 = np.sum(histogramCounts[jj:ii])
+            w3 = np.sum(histogramCounts[ii:])
+            if w1 > 0 and w2 > 0 and w3 > 0:
+                mu1 = np.dot(np.arange(1, jj+1), histogramCounts[:jj]) / w1
+                mu2 = np.dot(np.arange(jj+1, ii+1), histogramCounts[jj:ii]) / w2
+                mu3 = np.dot(np.arange(ii+1, nvals+1), histogramCounts[ii:]) / w3
+
+                val = w1 * ((mu1 - muT) ** 2) + w2 * ((mu2 - muT) ** 2) + w3 * ((mu3 - muT) ** 2)
+                if val >= maximum:
+                    level = [jj-1, ii-1]
+                    maximum = val
+
+    # plt.plot([x[level[0]]] * 2, [0, np.max(histogramCounts)], linewidth=2, color='r')
+    # plt.plot([x[level[1]]] * 2, [0, np.max(histogramCounts)], linewidth=2, color='r')
+    # plt.title('histogram of image pixels')
+    # plt.gcf().set_size_inches(9, 5)
+    # plt.show()
+
+    return level, x
+
+def scoringFunction(groundtruth, reconstruction):
+
+    if (np.any(groundtruth.shape != np.array([256, 256]))):
+        raise Exception('The shape of the given ground truth is not 256 x 256!')
+
+    if (np.any(reconstruction.shape != np.array([256, 256]))):
+        return 0
+
+    truth_c = np.zeros(groundtruth.shape)
+    truth_c[np.abs(groundtruth - 2) < 0.1] = 1
+    reco_c = np.zeros(reconstruction.shape)
+    reco_c[np.abs(reconstruction - 2) < 0.1] = 1
+
+    score_c = ssim(truth_c, reco_c)
+
+    truth_d = np.zeros(groundtruth.shape)
+    truth_d[np.abs(groundtruth - 1) < 0.1] = 1
+    reco_d = np.zeros(reconstruction.shape)
+    reco_d[np.abs(reconstruction - 1) < 0.1] = 1
+
+    score_d = ssim(truth_d, reco_d)
+
+    score = 0.5*(score_c + score_d)
+
+    return score
+
+def ssim(truth, reco):
+
+    c1 = 1e-4
+    c2 = 9e-4
+    r = 80
+
+    ws = np.ceil(2*r)
+    wr = np.arange(-ws, ws+1)
+    X, Y = np.meshgrid(wr, wr)
+    ker = (1/np.sqrt(2*np.pi)) * np.exp(-0.5 * np.divide((np.square(X) + np.square(Y)), r**2))
+    correction = sps.convolve2d(np.ones(truth.shape), ker, mode='same')
+
+    gt = np.divide(sps.convolve2d(truth, ker, mode='same'), correction)
+    gr = np.divide(sps.convolve2d(reco, ker, mode='same'), correction)
+
+    mu_t2 = np.square(gt)
+    mu_r2 = np.square(gr)
+    mu_t_mu_r = np.multiply(gt, gr)
+
+    sigma_t2 = np.divide(sps.convolve2d(np.square(truth), ker, mode='same'), correction) - mu_t2
+    sigma_r2 = np.divide(sps.convolve2d(np.square(reco), ker, mode='same'), correction) - mu_r2
+    sigma_tr = np.divide(sps.convolve2d(np.multiply(truth, reco), ker, mode='same'), correction) - mu_t_mu_r;
+
+    num = np.multiply((2*mu_t_mu_r + c1), (2*sigma_tr + c2))
+    den = np.multiply((mu_t2 + mu_r2 + c1), (sigma_t2 + sigma_r2 + c2))
+    ssimimage = np.divide(num, den)
+
+    score = np.mean(ssimimage)
+
+    return score
+