Skip to content
Snippets Groups Projects
Commit 7530d9b3 authored by Angeline Aguinaldo's avatar Angeline Aguinaldo
Browse files

Added python to SRC folder

parent 6a7c3f70
No related branches found
No related tags found
No related merge requests found
Pipeline #
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="jdk" jdkName="Python 3.5.1 (C:\Users\angeline\Anaconda3\python.exe)" jdkType="Python SDK" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
<component name="TestRunnerService">
<option name="PROJECT_TEST_RUNNER" value="Unittests" />
</component>
</module>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectLevelVcsManager" settingsEditedManually="false">
<OptionsSetting value="true" id="Add" />
<OptionsSetting value="true" id="Remove" />
<OptionsSetting value="true" id="Checkout" />
<OptionsSetting value="true" id="Update" />
<OptionsSetting value="true" id="Status" />
<OptionsSetting value="true" id="Edit" />
<ConfirmationsSetting value="0" id="Add" />
<ConfirmationsSetting value="0" id="Remove" />
</component>
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.5.1 (C:\Users\angeline\Anaconda3\python.exe)" project-jdk-type="Python SDK" />
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/Python.iml" filepath="$PROJECT_DIR$/.idea/Python.iml" />
</modules>
</component>
</project>
\ No newline at end of file
This diff is collapsed.
# 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.
# -------------------------------------------------------------------------------
# GetRandomEllipseImage() accepts [a] k-number of clusters (int). This function generates
# k-number of ellipse(s) with randomly initialized radii, rotation, and center of mass.
# The ellipse(s) is/are converted into a binary image. A binary image and the true boundary
# pixels
def GetRandomEllipseImage(K):
import numpy as np
# Initialize parameters of random ellipse
IMSIZE = 1000
mu1 = [500, 500]
bw = np.zeros([IMSIZE, IMSIZE])
ptsCombined = np.empty([1,2])
pts = []
ptsEdge = []
bValidEllipse = False
# Valid Ellipse Creation
while bValidEllipse == False: # validity is dependent on overlap score of ellipse
# Collect points and ellipse(s) edges of generated ellipse(s) centered at mu1
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)
# Offset ellipse(s) by randomly selecting a point of another ellipse as the new center of mass
for k in np.arange(1, K, 1, dtype=int):
ptsConcat = np.concatenate((pts[k], pts[k-1]))
idxK = np.round((pts[k].shape[0] - 1) * np.random.uniform())
idxCombined = 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)
# Flatten 2-D indices to make 1-D in Column-order
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)
# Check overlap of ellipse
bGoodOverlap = GetOverlap(idxPts)
if bGoodOverlap: # good overlap has pixel diff > 400pix
bValidEllipse = True
# Concatenate all the points of the ellipse(s) together
ptsCombined = np.concatenate(pts[0:K])
# Form binary image
bw[np.array([ptsCombined[0:,0]]), np.array([ptsCombined[0:,1]])] = 1
# Concatenate edges points of ellipse(s) together
ptsEdgeCombined = np.concatenate(ptsEdge[0:K])
# Identify crop boundaries
ellipseBounds = np.where(bw == 1)
minx = np.min(ellipseBounds[1]) - 10
maxx = np.max(ellipseBounds[1]) + 10
miny = np.min(ellipseBounds[0]) - 10
maxy = np.max(ellipseBounds[0]) + 10
# Crop boundaries on bw image
bw = bw[miny:maxy, minx:maxx]
bw = np.uint8(bw)
# Normalize edge points according to new image boundaries
ptsEdgeY = ptsEdgeCombined[0:, 0] - miny
ptsEdgeX = ptsEdgeCombined[0:, 1] - minx
ptsEdge = np.transpose(np.vstack((ptsEdgeY, ptsEdgeX)))
return bw, ptsEdge # returns binary image (uint8, array) and true edge points(float, array)
# GenEllipse() accepts [a] mean values (array), [b] covariance matrix (array), [c] isocontour scalar(int),
# and [d] okay to 'fill ellipse' (bool). This function uses principle axes determined by the eigendecomposition
# of the convariance matrix to form one ellipse and generate its binary image. It returns the points and
def GenEllipse(mu, cv, c, bFill):
import numpy as np
import scipy.ndimage.morphology as bwmorph
# Eigen-decomposition of covariance matrix
[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]
# Compute principle axes
a1 = np.sqrt(eigval[0] * c * eigvec[0:, 0])
a2 = np.sqrt(eigval[-1] * c * eigvec[0:, 1])
# General ellipse equation (centered at mu) in polar coordinates
theta = np.arange(0, 2 * np.pi, (2 * np.pi) / 1000, dtype=float) # range of theta
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))
# Form binary image
bw = np.zeros([1000, 1000])
bw[np.array(np.round(x), dtype=int), np.array(np.round(y), dtype=int)] = 1
# Fill image if the 'fill ellipse' bool is true
if bFill:
bw = bwmorph.binary_fill_holes(bw)
# Collect points of generated ellipse
idx = np.where(bw == 1)
pts = np.transpose(np.vstack((idx[0], idx[1])))
return pts, bw # returns points (float, array) and binary image of ellipse (int, array)
# GetRandomEllipsePts() accepts [a] mean and [b] size of the binary image. This generates one ellipse
# using random radii and random rotation. The ellipse is converted to a binary image.
def GetRandomEllipsePts(mu, IMSIZE):
import numpy as np
from numpy.matlib import repmat
import scipy.ndimage.morphology as bwmorph
# Initialize max radii parameters
RADII_0 = 10
RADII_1 = 20
# Generate random radii
r1 = RADII_0 + RADII_1 * np.random.rand()
r2 = RADII_0 + RADII_1 * np.random.rand()
# Populate covariance matrix [r1^2 0; 0 r2^2]
cv = np.zeros([2, 2])
cv[0,0] = r1**2
cv[1,1] = r2**2
# Generate ellipse points and edge points
pts = GenEllipse(mu, cv, 1, 1)
ptsEdge = GenEllipse(mu, cv, 1, 0)
pts = pts[0]
ptsEdge = ptsEdge[0]
# Center points around origin (0, 0)
pts = pts - np.matlib.repmat(mu, pts.shape[0], 1)
ptsEdge = ptsEdge - np.matlib.repmat(mu, ptsEdge.shape[0], 1)
# Rotate according to random theta by matrix multiplication of points with rotation matrix
theta = 2 * np.pi * np.random.rand() # randomly generated theta
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)
# Re-center points around mu
pts = pts + np.matlib.repmat(mu, pts.shape[0], 1)
ptsEdge = ptsEdge + np.matlib.repmat(mu, ptsEdge.shape[0], 1)
# Form binary image
bw = np.zeros([IMSIZE, IMSIZE])
roundPts = np.around(pts[0:,0:]) # must round points to index into binary image
roundPts = roundPts.astype(int) # must be integer (not float)
bw[roundPts[0:,0], roundPts[0:,1]] = 1
bw = bwmorph.binary_fill_holes(bw) # fill holes generated from rounding
# Identify points from newly formed binary image
idx = np.where(bw==1)
pts = np.transpose(np.vstack((idx[0], idx[1])))
pts = np.array(pts)
return pts, ptsEdge # returns ellipse points (int, array) and ellipse edge points (float, array)
# GetOverlap() accepts [a] linear index points (list) where each index of the list are the index points
# for each cluster. This using set difference theory to identify the non-intersecting points of each ellipse with
# the union of the remaining ellipse points. If the length of non-intersecting points is less than 400 pixels,
# the ellipses receive a bad overlap bool. Additionally, if the ellipses do not overlap, the ellipses receives
# a bad overlap bool.
def GetOverlap(idxPts):
import numpy as np
# Initialize good overlap boolean
bGoodOverlap = True
# Identify number of clusters
K = len(idxPts)
# Get intersection of each ellipse and remaining ellipse and evaluate overlap
for k in np.arange(0, K, 1):
kOther = np.setdiff1d(np.arange(0, K, 1), k)
idxPts_other = []
# Concatenate points of ellipses other than current k
for i in kOther:
idxPts_other = np.append(idxPts_other, idxPts[i])
idxPts_test = idxPts[k]
# Calculate intersection
intersectPts = np.intersect1d(idxPts_other, idxPts_test)
# Obtain non-intersecting points
diffPts = idxPts_test.size - intersectPts.size
# Evaluation of overlap
if diffPts < 400:
bGoodOverlap = False
break
if intersectPts.size == 0:
bGoodOverlap = False
break
return bGoodOverlap # returns if ellipses have good/True or bad/False overlap (bool)
# 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() accepts [a] logical image containing only one connected component.
# This function uses the pixel replication algorithm of replicating according the the
# respective depth of the pixel in the Euclidean distance transform.
def PixelReplicate(bwim):
import numpy as np
import cv2
# Computer Euclidean Distance Transform
print("Computing Distance Transform...")
bwd = cv2.distanceTransform(bwim, distanceType=cv2.DIST_L2, maskSize=cv2.DIST_MASK_5)
bwd = np.asmatrix(bwd, dtype=float)
roundBwd = np.round(bwd)
nPtsRep = np.sum(roundBwd, dtype=int) # total number of replicated points is the sum of the
# rounded pixel depths in the distance transform
ptsRep = np.zeros((nPtsRep.astype(int), 2)) # initialize replicated points array
# Extract the points of the ellipse from the binary image
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])
# POINT REPLICATION STEP
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 # returns replicated points(array)
# fitGMM() accepts [a] points replicated (array) and [b] k-number of clusters (int).
# It uses an Expectation Maximum algorithm designed to fit Gaussian mixture models (GMM).
# The GMM fitting initializes 10 times, with search iterations capped at 1000 max iterations
# or min tolerance of 1e-25. The covariance matrix type is "GENERIC" (symmetric and
# positive-definite matrix)
def fitGMM(ptsRep, k):
import numpy as np
import cv2
print("Fitting GMM...")
bestGMM = []
bestLL = float("-inf") # initialize lowest possible log likelihood value (-inf)
# Get unique points from replicated points
coorPtsRep = [tuple(x) for x in ptsRep]
uniqPts = sorted(set(coorPtsRep), key=lambda x: coorPtsRep.index(x))
uniqPts = np.array(uniqPts)
# GAUSSIAN MIXTURE MODEL FITTING
for i in np.arange(0, 10, 1): # Initialization Loop
# Random generation of indices into non-replicated points array
rndIdx = np.random.randint(low=0, high=uniqPts.shape[0], size=(1, 1, k), dtype=int)
rndMeans = np.matrix(uniqPts[np.array(rndIdx), 0:])
# Initialize GMM model
gmm0 = cv2.ml.EM_create()
gmm0.setCovarianceMatrixType(2)
gmm0.setClustersNumber(k)
termCrit = (1000, 1000, 1e-6)
gmm0.setTermCriteria(termCrit)
# Train the GMM
newGMM = gmm0.trainE(samples=ptsRep, means0=rndMeans)
# Calculate Log Likelihood of fit
LL = sum(newGMM[1]) # sum of log likelihoods per sample
print('Likelihood Logarithm of GMM fit (Rep. ', i, ') = ', LL)
# Select better fit using log likelihood comparison
if LL > bestLL: # closer to zero Negative-LL is better
bestGMM = gmm0
bestLL = LL
return bestGMM # return best gmm model (openCV EM class)
# drawClusters() accepts [a] gmm model, [b] the binary image, [c] the true edges of each
# ellipse. This function extracts ellipse points, clusters points using gmm model, and plots
# the clusters with the true boundaries drawn. This returns a bwLabel matrix which is the same
# size as the input image. Each pixel is labeled according which cluster it belongs to.
def drawClusters(gmm0, bwim, ptsEdge):
import numpy as np
import matplotlib.pyplot as plt
# Get ellipse points
r, c = np.where(bwim == 1)
pts = np.asarray([r, c], dtype=int)
pts = np.transpose(pts)
# Use the em::predict2() method to cluster ellipse points
print("Clustering...")
clusters = -1*np.ones((pts.shape[0], 1))
for i in np.arange(0, pts.shape[0], 1):
prob = gmm0.predict2(sample=np.matrix(pts[i, 0:]))
clusters[i] = prob[0][1] # Extract most probable cluster index
clusters = clusters.astype(int)
print("Clusters Found: ", np.unique(clusters)+1)
# Create bwLabel matrix
bwLabel = np.zeros((bwim.shape[0], bwim.shape[1]))
print("Plotting...")
for h in np.arange(0, np.max(clusters)+1, 1):
clustIdx = np.where(clusters == h)
bwLabel[np.array([pts[clustIdx[0], 0]]), np.array([pts[clustIdx[0], 1]])] = h + 1
bwLabel = np.uint8(bwLabel)
# Plot Result
plt.figure()
cmap = plt.cm.hsv
cmap.set_under(color="black")# set zero valued pixels to black
plt.imshow(bwLabel, cmap=cmap, vmin=0.00000001)
plt.hold(True)
# Plots boundaries if not all zero
if sum(np.unique(ptsEdge)) != 0:
plt.scatter(ptsEdge[0:, 1], ptsEdge[0:, 0], s=45, c='w', lw=0.5, alpha=1, edgecolors='w')
# Configure figure display parameters
plt.axis('off')
plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
plt.margins(0)
plt.autoscale(enable=True, axis='both', tight='True')
plt.show()
plt.hold(False)
return bwLabel # returns bwLabel matrix (uint8)
\ No newline at end of file
# PixelRep_Demo is the MAIN demo file. Run this file to simulate Pixel Rep.
# -------------------------------------------------------------------------------
# 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 PIXEL_REPLICATION_LIB as PR
import GENERATE_ELLIPSE_LIB as GenEllipse
# Generate random ellipse image
print("Generating Ellipse...")
K = 2 # number of ellipse to fit
bwim, trueEdges = GenEllipse.GetRandomEllipseImage(K)
# # You may import an image for demo by uncommenting the following code:
# import cv2
# im = cv2.imread('U:\Angeline\PixelRep\Ellipse_k4.tif', cv2.CV_8UC1)
# th, bwim = cv2.threshold(im, 0, 1, cv2.THRESH_BINARY)
# trueEdges = [0, 0] # keep as [0, 0] if true edges are unknown.
# # Include (r, c) array of boundary pixels if known
# # (Remember to modify K number of ellipses)
# 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
bwLabel = PR.drawClusters(gmmfit, bwim, trueEdges)
File added
File added
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment