In [1]:
import sys
sys.path.insert(0, '../../../ndreg/')
sys.path.insert(0,'../code/functions/')
import ndreg
import math
import cv2
import pickle
import subprocess
import matlab.engine
import os
import numpy as np
import SimpleITK as itk
import matplotlib.pyplot as plt
import plosLib as pLib
import connectLib as cLib
import hyperReg as hype
import scipy.io as sio
from cluster import Cluster
from affine import Affine
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from random import randrange
from random import uniform as floatRand
from PIL import Image
Input Space: A cost matrix for two sets of points $P_1$ and $P_2$, where both point sets have identical length $l$
Output Space: An optimal pairing of all $p_1 \in P_1$ with exactly one partner $p_2 \in P_2$ such that $\Sigma \text{Cost}(p_1, p_2)$ is minimized
Algorithm:
In [ ]:
######################################
###THIS IS PSEUDOCODE, WILL NOT RUN###
######################################
def hungarian(costMatrix):
p1CostMatrix = costMatrix.copy()
p2CostMatrix = costMatrix.copy().T
#for every point in the first set
for all p1 in p1CostMatrix:
#find its minimum weighted edge
minVal = min(costMatrix[p1])
#subtract minimum weight from all edges
p1CostMatrix[p1] = p1CostMatrix[p1] - minVal
#for every point in the second set
for all p2 in p2CostMatrix:
#find the minimum weighted edge
minVal = min(p2CostMatrix[p2])
#subtract the minimum weight from all edges
p1CostMatrix[p1] = p1CostMatrix[p1] - minVal
#generate adjacency matrix of only the 0 weight
#after the initial 2 steps
initialMatrix = zeros_like(costMatrix)
for y, x in p1CostMatrix.zero():
initialMatrix[y][x] = 1
for y, x, in p2CostMatrix.zero():
initialMatrix[y][x] = 1
#get the maximal matching after the initial step
matching = minimalMatching(initialMatrix)
#if the initialization solves the problem, we are done
if matching.fullRank():
return matching
#if not, run iterative step until convergence
while not matching.fullRank():
#get the minimum edge that is not yet paired
minRemainingWeight = min(matching.rowsWithoutPivot())
minRemainingP1 = argmin(matching.rowsWithoutPivot())
#subtract that weight from the remaining graph at that edge
initialMatrix[minRemainingP1] -= minRemainingWeight
matching = minimalMatching(initialMatrix)
return matching
In [25]:
#The algorithm can be called with the following
def hungarian(costMat):
#write the matlab argument to disk so the native process can access it
sio.savemat('matlabArg.mat', mdict={'matlabArg':costMat})
#run the matlab process
os.system('matlab -nodisplay -r \"load(\'/home/bstadt/Desktop/lab/background/matlabArg.mat\'); munkres(matlabArg); exit()\"')
#load the results
matlabOut = sio.loadmat('assignment.mat')['assignment']
os.system('rm assignment.mat matlabArg.mat')
#return matlab output as numpy array
return np.array(matlabOut)-1
In [3]:
def loss(cluster1, cluster2):
c1Centroid = cluster1.centroid
c2Centroid = cluster2.centroid
error = math.sqrt(
(c1Centroid[0] - c2Centroid[0])**2 +
(c1Centroid[1] - c2Centroid[1])**2 +
(c1Centroid[2] - c2Centroid[2])**2 +
(cluster1.compactness - cluster2.compactness)**2 +
(cluster1.volume - cluster2.volume)**2
)
return error
#The function for creating the cost matrix
def genCostMatrix(clusterList1, clusterList2):
#check for bipartality
'''
if not len(clusterList1) == len(clusterList2):
print 'Cluster lists must have the same size'
return
'''
costMatrix = np.zeros((len(clusterList1), len(clusterList2)))
for cIdx1 in range(len(clusterList1)):
for cIdx2 in range(len(clusterList2)):
costMatrix[cIdx1][cIdx2] = loss(clusterList1[cIdx1], clusterList2[cIdx2])
return costMatrix
#the function used to force a bipartite graph
def forceBipartite(clusterList1, clusterList2):
l1 = len(clusterList1)
l2 = len(clusterList2)
#if already compliant
if l1 == l2:
return clusterList1, clusterList2
#if there are too many points in l1
elif l1 > l2:
diff = l1 - l2
for i in range(diff):
delIdx = randrange(0, len(clusterList1))
del clusterList1[delIdx]
else:
diff = l2 - l1
for i in range(diff):
delIdx = randrange(0, len(clusterList2))
del clusterList2[delIdx]
return clusterList1, clusterList2
Practically, this algorithm will succeed when there is one or a few very similar assignments that minimize the cost funciton.
Practically, this algorithm will fail when there are multiple differing assignments that minimize the cost function
The following two are two data sets with solutions that can be run to ensure the success of the hungarian algorithm in the simple case
In [33]:
funcData1 = np.identity(2)
funcData2 = np.array([[4., 1., 3.], [2., 0., 5.], [3., 2., 2.]])
print 'Data:'
print funcData1
print 'Optimal Pairing:'
print '[[0, 1], [1, 0]]'
print '\nData:'
print funcData2
print 'Optimal Pairing:'
print '[[0, 1], [1, 0], [2, 2]]'
This data exists in 2 dimensions with 2 features and should converge after the initialization step. It has exactly 1 optimum pairing
This data exists in 2 dimensions with 3 features. It should converge after at least 1 iterative loop, and has exactly 1 optimum pairing
In [5]:
def toDiff(imgA, imgB):
ret = np.empty((imgA.shape[0], imgA.shape[1], 3), dtype=np.uint8)
for y in range(imgA.shape[0]):
for x in range(imgA.shape[1]):
if imgA[y][x] and not imgB[y][x]:
ret[y][x][0] = 255
ret[y][x][1] = 0
ret[y][x][2] = 0
elif not imgA[y][x] and imgB[y][x]:
ret[y][x][0] = 0
ret[y][x][1] = 255
ret[y][x][2] = 0
elif imgA[y][x] and imgB[y][x]:
ret[y][x][0] = 255
ret[y][x][1] = 255
ret[y][x][2] = 0
else:
ret[y][x][0] = 255
ret[y][x][1] = 255
ret[y][x][2] = 255
return ret
def visDiff(sliceA, sliceB):
disp = toDiff(sliceA, sliceB)
return disp
In [35]:
funcTest1 = zip(*hungarian(funcData1))
funcTest2 = zip(*hungarian(funcData2))
print 'Test1: ', funcTest1 == [(1.,), (0.,)]
print '\n\tExpected: ',[(1.,), (0.,)]
print '\n\tActual: ', funcTest1
print '\n'
print 'Test2: ', funcTest2 == [(1.,), (0.,), (2.,)]
print '\n\tExpected: ',[(1.,), (0.,), (2.,)]
print '\n\tActual: ', funcTest2
Functionality testing had perfect results. This means that the algorithm is ready to move on to the validation testing phase
In [5]:
#import the pickled versions of the real data
tp1 = pickle.load(open('../code/tests/synthDat/realDataRaw_t0.synth', 'r'))
tp2 = pickle.load(open('../code/tests/synthDat/realDataRaw_t1.synth', 'r'))
In [ ]:
#cut the data to a reasonable size for testing
tp1TestData = tp1[:7]
tp2TestData = tp2[:7]
In [ ]:
#run the data through the pipeline
tp1PostPipe = cLib.otsuVox(pLib.pipeline(tp1TestData))
tp2PostPipe = cLib.otsuVox(pLib.pipeline(tp2TestData))
#cut out the ill defined sections
tp1PostPipe = tp1PostPipe[1:6]
tp2PostPipe = tp2PostPipe[1:6]
In [ ]:
#Display the data to be used for testing
for i in range(tp1PostPipe.shape[0]):
fig = plt.figure()
plt.title('Time Point 1 Pipe Output at z='+str(i))
plt.imshow(tp1PostPipe[i], cmap='gray')
plt.show()
In [ ]:
#Display the data to be used for testing
for i in range(tp2PostPipe.shape[0]):
fig = plt.figure()
plt.title('Time Point 2 Pipe Output at z='+str(i))
plt.imshow(tp2PostPipe[i], cmap='gray')
plt.show()
In [ ]:
transform = hype.get3DRigid(pitch=0., yaw=.15, roll=0., xT=0., yT=0., zT=0.)
transformVolume = hype.apply3DRigid(tp1PostPipe, transform, True)
In [ ]:
for i in range(tp1PostPipe.shape[0]):
plt.figure()
plt.title('Initial Disperity at z='+str(i))
disp = visDiff(tp1PostPipe[i], transformVolume[i])
plt.imshow(disp)
plt.show()
In [13]:
#get the cluster lists from both the base and the transformation
tp1BaseClusters = cLib.connectedComponents(tp1PostPipe)
tp1TransClusters = cLib.connectedComponents(transformVolume)
In [ ]:
#generate cost matrix
costMat = genCostMatrix(tp1TransClusters, tp1BaseClusters)
In [28]:
assignments = hungarian(costMat)
In [ ]: