Added PixelRep in Python

Runs PixelRep and GMMfit
Need to fix a few bugs in Ellipse generation

SRC/MATLAB/PixelReplicate.m

@@ -74,6 +74,6 @@ function [bwd bw] = DepthImageFromPoints(pts,im)
 bw = false(size(im));
 pts=round(pts);
 pts=unique(pts,'rows');
-pts = sub2ind(size(bw),pts(:,1),pts(:,2));
+pts = sub2ind(size(bw),pts(:,2),pts(:,1));
 bw(pts)=1;
 bwd=bwdist(~bw);
\ No newline at end of file

SRC/MATLAB/demoPR.m

@@ -66,10 +66,10 @@ objPR = fitgmdist(ptsReplicated, K, 'replicates',5);
 
 % draw the clustering
 [y,x]=find(bw);
-idx = objPR.cluster([x,y]);
+idx = objPR.cluster([y,x]);
 for kk=1:K
     idxK=find(idx==kk);
-    plot(x(idxK),y(idxK),'.','color',cmap(kk,:))
+    plot(x(idxK),y(idxK),'.')
 end
 
 %draw the gaussian isocontours for the triangular elliptical approximation

SRC/Python/GenerateEllipseLib.py

+# GENERATE_ELLIPSE_LIB houses randomly generated as ellipse
+
+# Copyright (c) 2016, Drexel University
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+#
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+#
+# * Neither the name of PixelRep nor the names of its
+#   contributors may be used to endorse or promote products derived from
+#   this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+def GetRandomEllipseImage(K):
+    import numpy as np
+    IMSIZE = 1000
+    mu1 = [500, 500]
+    bw = np.zeros([IMSIZE, IMSIZE], dtype=bool)
+    ptsCombined = np.empty([1,2])
+    pts = []
+    ptsEdge = []
+    bValidEllipse = False
+
+    while ~bValidEllipse:
+        for k in np.arange(0, K, 1, dtype=int):
+            pts_tmp, ptsEdge_tmp = GetRandomEllipsePts(mu1, IMSIZE)
+            ptsCombined = np.append(ptsCombined, pts_tmp, axis=0)
+            pts.append(pts_tmp)
+            ptsEdge.append(ptsEdge_tmp)
+
+        for k in np.arange(1, K, 1, dtype=int):
+            ptsConcat = np.concatenate((pts[k], pts[k-1]))
+            idxK = 1 + np.round((pts[k].shape[0] - 1) * np.random.uniform())
+            idxCombined = 1 + np.round((ptsConcat.shape[0] - 1) * np.random.uniform())
+            offset = ptsConcat[int(idxCombined), 0:] - pts[k][int(idxK), 0:]
+            pts[k] = pts[k] + np.matlib.repmat(offset, pts[k].shape[0], 1)
+            ptsEdge[k] = ptsEdge[k] + np.matlib.repmat(offset, ptsEdge[k].shape[0], 1)
+
+        idxPts = []
+
+        for k in np.arange(0, K, 1, dtype=int):
+            ptsSub2Ind = np.ravel_multi_index([np.array(pts[k][0:,0]), np.array(pts[k][0:,1])], [bw.shape[0], bw.shape[1]])
+            idxPts.append(ptsSub2Ind)
+
+        bNoFullOverlap = GetOverlap(idxPts)
+
+        if bNoFullOverlap:
+            bValidEllipse = True
+
+    ptsCombined = np.concatenate(pts[0:K])
+    bw[np.array([ptsCombined[0:,0]]), np.array([ptsCombined[0:,1]])] = 1
+    return bw
+
+
+def genEllipse(mu, cv, c, bFill):
+    import numpy as np
+    import scipy.ndimage.morphology as bwmorph
+    bw = np.zeros([1000, 1000])
+    theta = np.arange(0, 2 * np.pi, (2 * np.pi)/1000, dtype=float)
+    [eigval, eigvec] = np.linalg.eig(cv)
+    if len([eigval.shape]) == 1:
+        eigval = np.sort(eigval)[::-1]
+        order = np.argsort(eigval)[::-1]
+    else:
+        eigval = np.sort(np.diagonal(eigval))[::-1]
+        order = np.argsort(np.diagonal(eigval))[::-1]
+
+    eigvec = eigvec[0:, order]
+
+    a1 = np.sqrt(eigval[0] * c * eigvec[0:, 0])
+    a2 = np.sqrt(eigval[-1] * c * eigvec[0:, 1])
+
+    x = np.array(mu[0] + a1[0] * np.sin(theta) + a2[0] * np.cos(theta))
+    y = np.array(mu[1] + a1[1] * np.sin(theta) + a2[1] * np.cos(theta))
+    # clamp values to fall between [1, 1000]
+    x = np.maximum(x, np.ones_like(x))
+    y = np.maximum(y, np.ones_like(y))
+    x = np.minimum(x, 1000 * np.ones_like(x))
+    y = np.minimum(y, 1000 * np.ones_like(y))
+
+    bw[np.array(np.round(x), dtype=int), np.array(np.round(y), dtype=int)] = 1
+
+    if bFill:
+        bw = bwmorph.binary_fill_holes(bw)
+
+    idx = np.where(bw==1)
+    pts = np.transpose(np.vstack((idx[0], idx[1])))
+
+    pts = np.array(pts)
+    return pts, bw
+
+
+
+def GetRandomEllipsePts(mu, IMSIZE):
+    import numpy as np
+    from numpy.matlib import repmat
+    import scipy.ndimage.morphology as bwmorph
+    pts = []
+    ptsEdge = []
+
+    RADII_0 = 10
+    RADII_1 = 20
+
+    r1 = RADII_0 + RADII_1 * np.random.rand()
+    r2 = RADII_0 + RADII_1 * np.random.rand()
+
+    cv = np.zeros([2, 2])
+    cv[0,0] = r1**2
+    cv[1,1] = r2**2
+
+    pts = genEllipse(mu, cv, 1, 1)
+    ptsEdge = genEllipse(mu, cv, 1, 0)
+
+    theta = 2 * np.pi * np.random.rand()
+
+    pts = pts - np.matlib.repmat(mu, pts.shape[0], 1)
+    ptsEdge = ptsEdge - np.matlib.repmat(mu, ptsEdge.shape[0], 1)
+
+    rotate = np.zeros([2, 2])
+    rotate[0,0] = np.cos(theta)
+    rotate[0,1] = -np.sin(theta)
+    rotate[1,0] = np.sin(theta)
+    rotate[1,1] = np.cos(theta)
+
+    pts = np.dot(pts,rotate)
+    ptsEdge = np.dot(ptsEdge,rotate)
+
+    pts = pts + np.matlib.repmat(mu, pts.shape[0], 1)
+    ptsEdge = ptsEdge + np.matlib.repmat(mu, ptsEdge.shape[0], 1)
+
+    bw = np.zeros([IMSIZE, IMSIZE])
+    roundPts = np.around(pts[0:,0:])
+    roundPts = roundPts.astype(int)
+    bw[roundPts[0:,0], roundPts[0:,1]] = 1
+    bw = bwmorph.binary_fill_holes(bw)
+
+    idx = np.where(bw==1)
+    pts = np.transpose(np.vstack((idx[0], idx[1])))
+
+    pts = np.array(pts)
+    ptsEdge = np.array(ptsEdge)
+    return pts, ptsEdge
+
+
+def GetOverlap(idxPts):
+    import numpy as np
+    bNoFullOverlap = True
+    K = len(idxPts)
+
+    for k1 in np.arange(0,K,1):
+        for k2 in np.arange(0,K,1):
+            if k1 == k2:
+                continue
+            idxPts_test1 = idxPts[k1]
+            idxPts_test2 = idxPts[k2]
+            diffPts = np.setdiff1d(idxPts_test1, idxPts_test2)
+            if diffPts.size == 0:
+                bNoFullOverlap = False
+                break
+
+    return bNoFullOverlap
\ No newline at end of file

SRC/Python/PixelRepLib.py

+# PIXEL_REP_LIB houses all PixelReplication functions including GMM fitting
+# and cluster drawing
+
+# Copyright (c) 2016, Drexel University
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+#
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+#
+# * Neither the name of PixelRep nor the names of its
+#   contributors may be used to endorse or promote products derived from
+#   this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+# PixelReplicate(im, label)
+# im can be a logical image - if so it should contain only one connected
+#   component
+# im can be a label image - it can contain any number of connected components
+#   and label is used to select the component to replicate
+
+def PixelReplicate(bwim):
+    from numpy.matlib import repmat
+    import numpy as np
+    import cv2
+
+    print("Computing distance transform...")
+    bwd = cv2.distanceTransform(bwim, distanceType=cv2.DIST_L2, maskSize=5)
+    bwd = np.asmatrix(bwd, dtype=float)
+    roundBwd = np.round(bwd)
+    nPtsRep = np.sum(roundBwd, dtype=int)
+
+    ptsRep = np.zeros((nPtsRep.astype(int), 2))
+    r, c = np.where(bwim == 1)
+    pts = np.asarray([r, c], dtype=int)
+    pts = np.transpose(pts)
+    nInsert = int(0)
+    sizePts = pts.shape
+    lengthPts = int(sizePts[0])
+
+    nrepTotal = []
+    print("Replicating Pixels...")
+    for i in np.arange(0, lengthPts, 1, dtype=int):
+        idxPts = np.asmatrix([pts[i,0], pts[i,1]])
+        nrep = np.round(bwd[idxPts[0, 0], idxPts[0, 1]])
+        nrep = int(nrep)
+        nrepTotal = np.append(nrepTotal,nrep)
+        rp = np.matlib.repmat(idxPts, nrep, 1)
+        if (nInsert + nrep) > nPtsRep:
+            ptsRep = np.append(ptsRep, rp, axis=0)
+        else:
+            ptsRep[int(nInsert):int(nInsert + nrep), 0:] = rp
+        nInsert = nInsert + nrep
+
+    return ptsRep
+
+
+def fitGMM(ptsRep, k):
+    from sklearn.mixture import GMM
+    print("Fitting GMM...")
+    gmm0 = GMM(n_components=k, covariance_type='full', n_iter=400, n_init=20)
+    gmmfit = gmm0.fit(ptsRep)
+    return gmmfit
+
+
+def drawClusters(gmmfit, bwim):
+    import numpy as np
+    import matplotlib.pyplot as plt
+
+    r, c = np.where(bwim == 1)
+    pts = np.asarray([r, c], dtype=int)
+    pts = np.transpose(pts)
+
+    print("Clustering...")
+    clusters = gmmfit.predict(pts)
+    print("Clusters Found: ", np.unique(clusters)+1)
+
+    plt.imshow(bwim, cmap='gray')
+    plt.hold(True)
+    print("Plotting...")
+    colors = ['m', 'r', 'c', 'g', 'y']
+    for h in np.unique(clusters):
+        clustIdx = np.where(clusters == h)
+        clustIdx = clustIdx[0]
+        print("Cluster ", h + 1, "Indices: ", clustIdx)
+        plt.scatter(x=pts[clustIdx, 1], y=pts[clustIdx, 0], c=colors[h])
+    plt.show()
+    plt.hold(False)
+    return
\ No newline at end of file

SRC/Python/PixelReplication.py

+
+# Copyright (c) 2016, Drexel University
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+# * Redistributions of source code must retain the above copyright notice, this
+#   list of conditions and the following disclaimer.
+#
+# * Redistributions in binary form must reproduce the above copyright notice,
+#   this list of conditions and the following disclaimer in the documentation
+#   and/or other materials provided with the distribution.
+#
+# * Neither the name of PixelRep nor the names of its
+#   contributors may be used to endorse or promote products derived from
+#   this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
+# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+import cv2
+import PixelRepLib as PR
+import GenerateEllipseLib as GenEllipse
+
+# Generate random ellipse image
+print("Generating Ellipse...")
+K = 5 # number of ellipse to fit
+# bwim = GenEllipse.GetRandomEllipseImage(K)
+
+# You may import an image for demo by uncommenting the following code:
+# im = cv2.imread('filename',cv2.CV_8UC1)
+# th, bwim = cv2.threshold(im, 0, 1, cv2.THRESH_BINARY)
+# plt.imshow(im)
+# plt.show()
+
+im = cv2.imread('U:\Angeline\PixelRep\Ellipse.tif',cv2.CV_8UC1)
+th, bwim = cv2.threshold(im, 0, 1, cv2.THRESH_BINARY)
+
+# IMPORTANT NOTE:
+# bwim only has 1 connected component. if bwim has multiple connected
+# components, process each one separately
+
+ptsRep = PR.PixelReplicate(bwim)
+
+# Fit GMM to PR points
+# NOTE - more replicates will be more accurate, but will take longer.
+
+gmmfit = PR.fitGMM(ptsRep, K)
+
+# Draw clusters
+
+PR.drawClusters(gmmfit, bwim)
+
+# Draw Gaussian isocontours for triangular elliptical approximation (see online methods).
+# These are exact for circles and approximate for ellipticals.
+
+
+
+
+
+