977 lines
30 KiB
Python
Executable File
977 lines
30 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
|
|
from __future__ import division
|
|
|
|
import argparse, pdb
|
|
import numpy as np
|
|
import os
|
|
|
|
from scipy.sparse import csr_matrix, coo_matrix
|
|
from scipy.sparse.csgraph import minimum_spanning_tree, depth_first_tree
|
|
from scipy.spatial import Delaunay, ConvexHull
|
|
|
|
import pulp
|
|
import timeit as T
|
|
|
|
class Vertex(object):
|
|
'''
|
|
Defines vertex.
|
|
'''
|
|
|
|
def __init__(self, x=None, y=None, phase=None, compNumber=None, index=None):
|
|
self.x = x
|
|
self.y = y
|
|
self.phase = phase
|
|
self.compNumber = compNumber
|
|
self.index = index
|
|
self.pts = None
|
|
self.sigma = None
|
|
self.source = None
|
|
self.dist = None
|
|
self.n = None
|
|
|
|
def __str__(self):
|
|
ostr = 'Location: (%d, %d)'%(self.y, self.x)
|
|
#ostr += '\nComponent: (%d, %d)'%(self.source, self.compNumber)
|
|
return ostr
|
|
|
|
def __eq__(self, other):
|
|
# if other is not None:
|
|
try:
|
|
return (self.x == other.x) and (self.y == other.y)
|
|
except:
|
|
pass
|
|
return None
|
|
|
|
def __hash__(self):
|
|
return hash((self.x,self.y))
|
|
|
|
def updatePhase(self, n):
|
|
self.n = n
|
|
|
|
def getIndex(self):
|
|
return self.index
|
|
|
|
def getPhase(self):
|
|
return self.n
|
|
|
|
def getUnwrappedPhase(self):
|
|
return self.phase - 2 * self.n * np.pi
|
|
|
|
class Edge(object):
|
|
'''
|
|
Defines edge of Delaunay Triangulation.
|
|
'''
|
|
|
|
def __init__(self, source, dest, dist=1, name=None):
|
|
self.__CONST = 10000
|
|
self.__MAXFLOW = 100
|
|
self.src = source
|
|
self.dst = dest
|
|
self.triIdx = None
|
|
self.adjtriIdx = None
|
|
self.flow = None
|
|
self.cost = self.__computeCost()
|
|
self.__dist__ = dist
|
|
|
|
if name is None:
|
|
name = "E(%d,%d)"%(source.getIndex(), dest.getIndex())
|
|
|
|
# Using PuLP to define the variable
|
|
self.__var__ = pulp.LpVariable(name, 0, 1, pulp.LpContinuous)
|
|
|
|
def isBarrier(self):
|
|
if None not in (self.src.compNumber, self.dst.compNumber):
|
|
return (self.src.compNumber == self.dst.compNumber)
|
|
else:
|
|
return False
|
|
|
|
def isCorner(self):
|
|
return self.adjtriIdx == None
|
|
|
|
def __computeCost(self):
|
|
if self.isBarrier():
|
|
return self.__CONST
|
|
else:
|
|
return 1
|
|
|
|
def isUnwrapped(self):
|
|
return None not in (self.src.getPhase(), self.dst.getPhase())
|
|
|
|
def diff(self):
|
|
return int(np.round((self.dst.phase - self.src.phase)/(2*np.pi)))
|
|
|
|
def updateTri(self, index):
|
|
if self.triIdx is not None:
|
|
self.triIdx.append(index)
|
|
else:
|
|
self.triIdx = [index]
|
|
|
|
def updateAdj(self, index):
|
|
if self.adjtriIdx is not None:
|
|
self.adjtriIdx.append(index)
|
|
else:
|
|
self.adjtriIdx = [index]
|
|
|
|
def updateFlow(self, flow):
|
|
if self.isBarrier():
|
|
# Check if the solver's solution for high cost node is zero
|
|
if (flow != 0):
|
|
raise ValueError("Solver Solution Incorrect")
|
|
self.flow = flow
|
|
|
|
def unwrap(self, rEdge):
|
|
if self.src.getPhase() is not None:
|
|
self.dst.updatePhase(self.src.getPhase() + (self.flow - rEdge.flow) + self.diff())
|
|
else:
|
|
return None
|
|
|
|
def plot(self, plt, c='g'):
|
|
plt.plot([self.src.x, self.dst.x], [self.src.y, self.dst.y], c=c)
|
|
|
|
if self.flow != 0 and self.isBarrier():
|
|
plt.plot([self.src.x, self.dst.x], [self.src.y, self.dst.y], c=c)
|
|
|
|
def getFlow(self, neutralNode):
|
|
if self.adjtriIdx == None:
|
|
return 'a %d %d 0 %d %d\n' % (neutralNode, self.triIdx[0], self.__MAXFLOW, self.cost)
|
|
else:
|
|
return 'a %d %d 0 %d %d\n' % (self.adjtriIdx[0], self.triIdx[0], self.__MAXFLOW, self.cost)
|
|
|
|
def getCost(self):
|
|
return self.cost
|
|
|
|
def getNeutralWeight(self):
|
|
# Number of times the edge was used with negative weight
|
|
if self.adjtriIdx is None:
|
|
numReverseLoop = 0
|
|
else:
|
|
numReverseLoop = len(self.adjtriIdx)
|
|
|
|
# Number of times the edge was used with positive weight
|
|
if self.triIdx is None:
|
|
numLoop = 0
|
|
else:
|
|
numLoop = len(self.triIdx)
|
|
|
|
return numReverseLoop - numLoop
|
|
|
|
def getLPVar(self):
|
|
return self.__var__
|
|
|
|
def sign(self):
|
|
if self.src.compNumber < self.dst.compNumber:
|
|
return 1
|
|
elif self.src.compNumber > self.dst.compNumber:
|
|
return -1
|
|
else:
|
|
return 0
|
|
|
|
def __str__(self):
|
|
ostr = 'Edge between : \n'
|
|
ostr += str(self.src) + '\n'
|
|
ostr += str(self.dst) + '\n'
|
|
ostr += 'with Residue %f'%(self.flow)
|
|
return ostr
|
|
|
|
def __eq__(self, other):
|
|
return (self.src == other.src) and (self.dst == other.dst)
|
|
|
|
def __neg__(self):
|
|
return Edge(self.dst, self.src)
|
|
|
|
# Defines each delaunay traingle
|
|
class Loop(object):
|
|
'''
|
|
Collection of edges in loop - Expects the vertices to be in sequence
|
|
'''
|
|
def __init__(self, vertices=None, edges=None, index=None):
|
|
self.edges = None
|
|
self.index = None
|
|
self.residue = None
|
|
|
|
# Reverse edges that are used during unwrapping - Done contribute
|
|
# in the residue computation
|
|
self.__edges = None
|
|
self.__center = None
|
|
|
|
def_vertices = True
|
|
if 'any' not in dir(vertices) and vertices == None:
|
|
def_vertices = False
|
|
|
|
def_edges = True
|
|
if 'any' not in dir(edges) and edges == None:
|
|
def_edges = False
|
|
|
|
def_index = True
|
|
if 'any' not in dir(index) and index == None:
|
|
def_index = False
|
|
|
|
# initializes only when all the vertices are availables
|
|
# if None not in (vertices, edges, index):
|
|
if def_vertices and def_edges and def_index:
|
|
self.index = index
|
|
self.edges = self.__getEdges(vertices, edges)
|
|
self.__edges = self.__getReverseEdges(vertices, edges)
|
|
self.residue = self.computeResidue()
|
|
self.__updateEdges()
|
|
self.__center = self.__getCenter(vertices)
|
|
|
|
# Edges traversing the vertices in a sequence
|
|
@staticmethod
|
|
def __getReverseEdges(vertices, edges):
|
|
rSeqEdges = []
|
|
for vx, vy in zip(vertices[1:] + [vertices[0]], vertices):
|
|
rSeqEdges.append(edges[vx, vy])
|
|
return rSeqEdges
|
|
|
|
# Edges traversing the vertices in a sequence
|
|
@staticmethod
|
|
def __getEdges(vertices, edges):
|
|
seqEdges = []
|
|
for vx, vy in zip(vertices, vertices[1:] + [vertices[0]]):
|
|
seqEdges.append(edges[vx, vy])
|
|
return seqEdges
|
|
|
|
# Returns a string in the RelaxIV format
|
|
def getNodeSupply(self):
|
|
return "n %d %d\n"%(self.index, self.residue)
|
|
|
|
# Returns a string in the RelaxIV format for flows
|
|
def getFlowConstraint(self, neutralNode):
|
|
edgeFlow = []
|
|
for edge in self.edges:
|
|
edgeFlow.append(edge.getFlow(neutralNode))
|
|
return edgeFlow
|
|
|
|
def getLPFlowConstraint(self):
|
|
edgeFlow = []
|
|
lpConstraint = []
|
|
for edge, rEdge in zip(self.edges, self.__edges):
|
|
lpConstraint.extend([edge.getLPVar(), -(rEdge.getLPVar())])
|
|
return (pulp.lpSum(lpConstraint) == -self.residue)
|
|
|
|
# Updates the cost of the edges
|
|
def updateEdgeFlow(self, flow):
|
|
for cost, edge in zip(flow, self.edges):
|
|
edge.updateFlow(cost)
|
|
|
|
# Computes the Residue in the triangle
|
|
def computeResidue(self):
|
|
isBarrier = map(lambda x: x.isBarrier(), self.edges)
|
|
if any(isBarrier):
|
|
return 0
|
|
else:
|
|
residue = 0
|
|
for edge in self.edges:
|
|
residue = residue + edge.diff()
|
|
return residue
|
|
|
|
# Each edge keeps a list of triangles it is part of
|
|
def __updateEdges(self):
|
|
for edge in self.edges:
|
|
edge.updateTri(self.index)
|
|
for edge in self.__edges:
|
|
edge.updateAdj(self.index)
|
|
|
|
def unwrap(self):
|
|
for edge, redge in zip(self.edges, self.__edges):
|
|
edge.unwrap(redge)
|
|
|
|
# Return None if not corner; else the rEdge
|
|
def Corner(self):
|
|
cornerEdge = []
|
|
for edge, redge in zip(self.edges, self.__edges):
|
|
if edge.isCorner():
|
|
cornerEdge.append(redge)
|
|
return cornerEdge
|
|
|
|
# test function
|
|
@staticmethod
|
|
def __getCenter(vertices):
|
|
center = (0,0)
|
|
for v in vertices:
|
|
center = center + (v.x, v.y)
|
|
return (center[0]/len(vertices), center[1]/len(vertices))
|
|
|
|
def isUnwrapped(self):
|
|
unWrapped = map(lambda x: x.isUnwrapped(), self.edges)
|
|
return all(unWrapped)
|
|
|
|
def printEdges(self):
|
|
for v in self.edges:
|
|
print(v)
|
|
|
|
def printFlow(self):
|
|
flow = []
|
|
for edge in self.edges:
|
|
flow.append(edge.flow)
|
|
print(flow)
|
|
|
|
def plot(self, ax):
|
|
if self.residue != 0:
|
|
ax.plot(self.__center[0], self.__center[1], '*', c='r')
|
|
|
|
# Packs all the traingles together
|
|
class PhaseUnwrap(object):
|
|
def __init__(self, x=None, y=None, phase=None, compNum=None, redArcs=0):
|
|
# Expects a list of ve:tices
|
|
self.loops = None
|
|
self.neutralResidue = None
|
|
self.neutralNodeIdx = None
|
|
self.__unwrapSeq = None
|
|
self.__unwrapped = None
|
|
self.__cornerEdges = []
|
|
|
|
# used only for plotting, and finally returning in sequence
|
|
# Dont use these variables for computation
|
|
self.__x = None
|
|
self.__y = None
|
|
self.__delauneyTri = None
|
|
self.__vertices = None
|
|
self.__edges = None
|
|
self.__compNum = None
|
|
self.__adjMat = None
|
|
self.__spanningTree = None
|
|
self.__CSRspanTree = None
|
|
self.__redArcs = None
|
|
|
|
# Edges used for unwrapping
|
|
self.__unwrapEdges = None
|
|
|
|
# Using PuLP as an interface for defining the LP problem
|
|
self.__prob__ = pulp.LpProblem("Unwrapping as LP optimization problem", pulp.LpMinimize)
|
|
|
|
if compNum is None:
|
|
compNum = [None]*len(x)
|
|
else:
|
|
self.__compNum = compNum
|
|
|
|
def_x = True
|
|
if 'any' not in dir(x) and x == None:
|
|
def_x = False
|
|
|
|
def_y = True
|
|
if 'any' not in dir(y) and y == None:
|
|
def_y = False
|
|
|
|
def_phase = True
|
|
if 'any' not in dir(phase) and phase == None:
|
|
def_phase = False
|
|
|
|
# if None not in (x, y, phase):
|
|
if def_x and def_y and def_phase:
|
|
# Create
|
|
vertices = self.__createVertices(x, y, phase, compNum)
|
|
self.nVertices = len(vertices)
|
|
delauneyTri = self.__createDelaunay(vertices)
|
|
self.__adjMat = self.__getAdjMat(delauneyTri.vertices)
|
|
self.__spanningTree = self.__getSpanningTree()
|
|
edges = self.__createEdges(vertices, redArcs)
|
|
|
|
if (redArcs >= 0):
|
|
self.loops = self.__createLoop(vertices, edges)
|
|
else:
|
|
self.loops = self.__createTriangulation(delauneyTri, vertices, edges)
|
|
self.neutralResidue, self.neutralNodeIdx = self.__computeNeutralResidue()
|
|
|
|
# Saving some variables for plotting
|
|
self.__redArcs = redArcs
|
|
self.__x = x
|
|
self.__y = y
|
|
self.__vertices = vertices
|
|
self.__delauneyTri = delauneyTri
|
|
self.__edges = edges
|
|
self.__unwrapEdges = []
|
|
|
|
def __getSpanningTree(self):
|
|
'''
|
|
Computes spanning tree from adjcency matrix
|
|
'''
|
|
from scipy.sparse import csr_matrix
|
|
from scipy.sparse.csgraph import minimum_spanning_tree
|
|
|
|
# Spanning Tree
|
|
spanningTree = minimum_spanning_tree(csr_matrix(self.__adjMat))
|
|
spanningTree = spanningTree.toarray().astype(int)
|
|
|
|
# Converting into a bi-directional graph
|
|
spanningTree = np.logical_or(spanningTree, spanningTree.T).astype(int)
|
|
return spanningTree
|
|
|
|
@staticmethod
|
|
def __getTriSeq(va, vb, vc):
|
|
'''
|
|
Get Sequence of triangle points
|
|
'''
|
|
line = lambda va, vb, vc: (((vc.y - va.y)*(vb.x - va.x))) - ((vc.x - va.x)*(vb.y - va.y))
|
|
|
|
# Line equation through pt0 and pt1
|
|
# Test for pt3 - Does it lie to the left or to the right ?
|
|
pos3 = line(va, vb, vc)
|
|
if(pos3 > 0):
|
|
# left
|
|
return (va, vc, vb)
|
|
else: # right
|
|
return (va, vb, vc)
|
|
|
|
# Create Delaunay Triangulation.
|
|
def __createDelaunay(self, vertices, ratio=1.0):
|
|
pts = np.zeros((len(vertices), 2))
|
|
|
|
for index, point in enumerate(vertices):
|
|
pts[index,:] = [point.x, ratio*point.y]
|
|
|
|
return Delaunay(pts)
|
|
|
|
def __createVertices(self, x, y, phase, compNum):
|
|
vertices = []
|
|
for i, (cx, cy, cphase, cNum) in enumerate(zip(x, y, phase, compNum)):
|
|
vertex = Vertex(cx, cy, cphase, cNum, i)
|
|
vertices.append(vertex)
|
|
return vertices
|
|
|
|
# Edges indexed with the vertices
|
|
def __createEdges(self, vertices, redArcs):
|
|
edges = {}
|
|
def add_edge(i,j,dist):
|
|
if (i,j) not in edges:
|
|
cand = Edge(i,j,dist,name="E(%d,%d)"%(i.getIndex(),j.getIndex()))
|
|
edges[i, j] = cand
|
|
if (j,i) not in edges:
|
|
cand = Edge(j,i,dist,name="E(%d,%d)"%(j.getIndex(),i.getIndex()))
|
|
edges[j, i] = cand
|
|
return
|
|
|
|
# Mask used to keep track of already created edges for efficiency
|
|
mask = np.zeros((self.nVertices, self.nVertices))
|
|
|
|
# Loop count depending on redArcs or MCF
|
|
if redArcs <= 0:
|
|
loopCount = 0
|
|
else:
|
|
loopCount = redArcs
|
|
|
|
# Edges are added using the adj Matrix
|
|
distMat = self.__adjMat.copy()
|
|
dist = 1
|
|
while (loopCount >= 0):
|
|
# Remove self loops
|
|
np.fill_diagonal(distMat, 0)
|
|
|
|
# Add edges with incrementing length
|
|
Va, Vb = np.where(np.logical_and((distMat > 0), np.invert(mask.astype(bool))))
|
|
for (va, vb) in zip(Va, Vb):
|
|
add_edge(vertices[va], vertices[vb], dist)
|
|
|
|
# Update mask
|
|
mask = np.logical_or(self.__adjMat.astype(bool), mask)
|
|
|
|
# Generating adj matrix with redundant arcs
|
|
distMat = np.dot(distMat, self.__adjMat.T)
|
|
loopCount = loopCount-1
|
|
|
|
return edges
|
|
|
|
def __getAdjMat(self, vertices):
|
|
'''
|
|
Get adjacency matrix using the vertices
|
|
'''
|
|
adjMat = np.zeros((self.nVertices, self.nVertices))
|
|
for ia, ib, ic in vertices:
|
|
adjMat[ia, ib] = adjMat[ib, ic] = adjMat[ic, ia] = 1
|
|
adjMat[ib, ia] = adjMat[ic, ib] = adjMat[ia, ic] = 1
|
|
|
|
return adjMat
|
|
|
|
# Inefficient - need to be replaced
|
|
def __getSequence(self, idxA, idxB):
|
|
from scipy.sparse import csr_matrix
|
|
from scipy.sparse.csgraph import depth_first_order
|
|
|
|
def traverseToRoot(nodeSeq, pred):
|
|
# scipy uses -9999
|
|
__SCIPY_END = -9999
|
|
seqV = [idxB]
|
|
parent = pred[idxB]
|
|
while parent != __SCIPY_END:
|
|
seqV = [parent] + seqV
|
|
parent = pred[parent]
|
|
return seqV
|
|
|
|
if self.__CSRspanTree is None:
|
|
self.__CSRspanTree = csr_matrix(self.__spanningTree)
|
|
|
|
(nodeSeq, pred) = depth_first_order(self.__CSRspanTree, i_start=idxA, \
|
|
directed=False, \
|
|
return_predecessors=True)
|
|
|
|
# Traverse through predecessors to the root node
|
|
seqV = traverseToRoot(nodeSeq, pred)
|
|
if (seqV[0] != idxA):
|
|
raise ValueError("Traversal Incorrect")
|
|
else:
|
|
return seqV
|
|
|
|
def __createLoop(self, vertices, edges):
|
|
# Creating the traingle object
|
|
# triangles.append(Loop([va, vb, vc], edges, i+1))
|
|
loops = []
|
|
zeroVertex = Vertex(0,0)
|
|
index = 0
|
|
loopExist = {}
|
|
for (va, vb), edge in edges.items():
|
|
if self.__spanningTree[va.index, vb.index]:
|
|
# Edge belongs to spanning tree
|
|
continue
|
|
else:
|
|
# Traverse through the spanning tree
|
|
seqV = self.__getSequence(va.index, vb.index)
|
|
seqV = [vertices[i] for i in seqV]
|
|
orientV = self.__getTriSeq(seqV[0], seqV[1], seqV[-1])
|
|
|
|
# Reverse SeqV if needed to get proper loop direction
|
|
if orientV[1] != seqV[1]:
|
|
seqV = seqV[::-1]
|
|
|
|
# Get a uniform direction of loop
|
|
if (seqV[0], seqV[-1]) in loopExist:
|
|
# Loop already added
|
|
continue
|
|
else:
|
|
loopExist[(seqV[0], seqV[-1])] = 1
|
|
|
|
# Create Loop
|
|
loops.append(Loop(seqV, edges, index+1))
|
|
index = index + 1
|
|
|
|
return loops
|
|
|
|
def __createTriangulation(self, delauneyTri, vertices, edges):
|
|
'''
|
|
Creates the Triangle residule for MCF formulation
|
|
'''
|
|
|
|
# Generate the points in a sequence
|
|
def getSeq(va, vb, vc):
|
|
line = lambda va, vb, vc: (((vc.y - va.y)*(vb.x - va.x))) - ((vc.x - va.x)*(vb.y - va.y))
|
|
|
|
# Line equation through pt0 and pt1
|
|
# Test for pt3 - Does it lie to the left or to the right ?
|
|
pos3 = line(va, vb, vc)
|
|
if(pos3 > 0):
|
|
# left
|
|
return (va, vc, vb)
|
|
else: # right
|
|
return (va, vb, vc)
|
|
|
|
loop = []
|
|
for i, triIdx in enumerate(delauneyTri.simplices):
|
|
va, vb, vc = self.__indexList(vertices, triIdx)
|
|
va, vb, vc = getSeq(va, vb, vc)
|
|
loop.append(Loop([va, vb, vc], edges, i+1))
|
|
return loop
|
|
|
|
def __unwrapNode__(self, idx, nodes, predecessors):
|
|
'''
|
|
Unwraps a specific node while traversing to the root
|
|
'''
|
|
if self.__unwrapped[predecessors[idx]] == 0:
|
|
# Go unwrap the predecessor before unwrapping current idx
|
|
self.__unwrapNode__(predecessors[idx], nodes, predecessors)
|
|
|
|
# Unwrap Node
|
|
srcV = self.__vertices[predecessors[idx]]
|
|
destV = self.__vertices[idx]
|
|
edge = self.__edges[srcV, destV]
|
|
rEdge = self.__edges[destV, srcV]
|
|
self.__unwrapEdges.append([edge, rEdge])
|
|
edge.unwrap(rEdge)
|
|
|
|
# Update Flag
|
|
self.__unwrapped[idx] = 1
|
|
return
|
|
|
|
def unwrapLP(self):
|
|
from scipy.sparse.csgraph import breadth_first_order as DFS
|
|
|
|
# Depth First search to get sequence of nodes
|
|
(nodes, predecessor) = DFS(self.__spanningTree, 0)
|
|
|
|
# Init Vertex
|
|
self.__unwrapped = np.zeros((len(self.__vertices)))
|
|
self.__vertices[nodes[0]].updatePhase(0)
|
|
self.__unwrapped[0] = 1
|
|
|
|
# Loop until there is a node unwrapped
|
|
while (0 in self.__unwrapped):
|
|
idx = next((i for i, x in enumerate(self.__unwrapped) if not(x)), None)
|
|
self.__unwrapNode__(idx, nodes, predecessor)
|
|
|
|
# Returns unwrapped vertices n values
|
|
nValue = []
|
|
for vertex in self.__vertices:
|
|
nValue.append(vertex.getPhase())
|
|
|
|
return nValue
|
|
|
|
# Unwrapping the triangles
|
|
# Traingle unwrap
|
|
def __unwrap__(self):
|
|
# Start unwrapping with the root
|
|
initTriangle = self.loops[self.__unwrapSeq[0]]
|
|
initTriangle.edges[0].src.updatePhase(0)
|
|
|
|
for triIdx in self.__unwrapSeq:
|
|
tri = self.loops[triIdx]
|
|
tri.unwrap()
|
|
|
|
# Returns unwrapped vertices n values
|
|
nValue = []
|
|
for vertex in self.__vertices:
|
|
nValue.append(vertex.getPhase())
|
|
|
|
return nValue
|
|
|
|
# Unwrap triangle wise
|
|
def __unwrapSequence(self, neighbors):
|
|
# Generate sequence
|
|
nodeSequence = [0]
|
|
for i in range(len(self.loops)):
|
|
|
|
# Finding adjacent nodes to current triangle
|
|
cNode = nodeSequence[i]
|
|
adjlist = np.array(neighbors[cNode])
|
|
adjNode = adjlist[np.where(adjlist != -1)[0]]
|
|
|
|
# adding list of new nodes by carefully removing already existing nodes
|
|
newNodes = list(set(adjNode) - set(nodeSequence))
|
|
nodeSequence = nodeSequence + newNodes
|
|
|
|
return nodeSequence
|
|
|
|
# Balances the residue in the entire network
|
|
def __computeNeutralResidue(self):
|
|
neutralResidue = 0
|
|
for t in self.loops:
|
|
neutralResidue = neutralResidue - t.residue
|
|
return (neutralResidue, len(self.loops) + 1)
|
|
|
|
def __nodeSupply(self):
|
|
nodeSupply = []
|
|
for i, loop in enumerate(self.loops):
|
|
nodeSupply.append(loop.getNodeSupply())
|
|
nodeSupply.append("n %d %d\n"%(self.neutralNodeIdx, self.neutralResidue))
|
|
return nodeSupply
|
|
|
|
def __flowConstraints(self):
|
|
edgeFlow = []
|
|
for loop in self.loops:
|
|
edgeFlow.extend(loop.getFlowConstraint(self.neutralNodeIdx))
|
|
|
|
# Edges between neutral and corner trianglea
|
|
for loop in self.loops:
|
|
cornerEdge = loop.Corner()
|
|
if cornerEdge != []:
|
|
for rEdge in cornerEdge:
|
|
# Dirty - To be cleaned later
|
|
rEdge.updateTri(self.neutralNodeIdx)
|
|
|
|
# Edge flow constraints for corner edges
|
|
edgeFlow.append(rEdge.getFlow(self.neutralNodeIdx))
|
|
self.__cornerEdges.append(rEdge)
|
|
|
|
return edgeFlow
|
|
|
|
# Phase unwrap using relax IV
|
|
@staticmethod
|
|
def __createRelaxInput(nodes, edges, fileName):
|
|
# Input file for Relax IV
|
|
f = open(fileName, "w")
|
|
|
|
# Prepending nodes and edges
|
|
f.write("p min %d %d\n"%(len(nodes), len(edges)))
|
|
|
|
# Write node
|
|
for line in nodes:
|
|
f.write(line)
|
|
|
|
# Write edge
|
|
for line in edges:
|
|
f.write(line)
|
|
|
|
# Done writting network
|
|
f.close()
|
|
|
|
@staticmethod
|
|
def __MCFRelaxIV(edgeLen, fileName="network.dmx"):
|
|
''' Uses MCF from RelaxIV '''
|
|
try:
|
|
from . import unwcomp
|
|
except:
|
|
raise Exception("MCF requires RelaxIV solver - Please drop the RelaxIV software \
|
|
into the src folder and re-make")
|
|
return unwcomp.relaxIVwrapper_Py(fileName)
|
|
|
|
def solve(self, solver, filename="network.dmx"):
|
|
# Choses LP or Min cost using redArcs
|
|
if (self.__redArcs == -1):
|
|
self.__solveMinCost__(filename)
|
|
else:
|
|
self.__solveEdgeCost__(solver)
|
|
|
|
def __solveEdgeCost__(self, solver, fileName="network.dmx"):
|
|
# Add objective
|
|
objective = []
|
|
for v, edge in self.__edges.items():
|
|
objective.append(edge.getCost() * edge.getLPVar())
|
|
self.__prob__ += pulp.lpSum(objective)
|
|
|
|
# Add Constraints
|
|
for loop in self.loops:
|
|
self.__prob__.addConstraint(loop.getLPFlowConstraint())
|
|
|
|
# Solve the objective function
|
|
if solver == 'glpk':
|
|
print('Using GLPK MIP solver')
|
|
MIPsolver = lambda: self.__prob__.solve(pulp.GLPK(msg=0))
|
|
elif solver == 'pulp':
|
|
print('Using PuLP MIP solver')
|
|
MIPsolver = lambda: self.__prob__.solve()
|
|
elif solver == 'gurobi':
|
|
print('Using Gurobi MIP solver')
|
|
MIPsolver = lambda: self.__prob__.solve(pulp.GUROBI_CMD())
|
|
|
|
print('Time Taken (in sec) to solve: %f'%(T.timeit(MIPsolver, number=1)))
|
|
|
|
# Get solution
|
|
for v, edge in self.__edges.items():
|
|
flow = pulp.value(edge.getLPVar())
|
|
edge.updateFlow(flow)
|
|
|
|
def __writeMPS__(self, filename):
|
|
self.__prob__.writeMPS(filename)
|
|
|
|
# Phase Unwrap using picos
|
|
def __solveMinCost__(self, fileName="network.dmx"):
|
|
|
|
# Creating RelaxIV input
|
|
nodes = self.__nodeSupply()
|
|
edgeFlow = self.__flowConstraints()
|
|
|
|
# Running RelaxIV Interface
|
|
self.__createRelaxInput(nodes, edgeFlow, fileName)
|
|
edgeWeight = self.__MCFRelaxIV(len(edgeFlow), fileName)
|
|
|
|
# Creating dictionary for interior edges
|
|
triEdgeCost = zip(*(iter(edgeWeight[:3*len(self.loops)]),) *3)
|
|
for triCost, tri in zip(triEdgeCost, self.loops):
|
|
tri.updateEdgeFlow(triCost)
|
|
|
|
# Dirty - to be changed
|
|
cornerEdgeCost = edgeWeight[3*len(self.loops):]
|
|
for rEdge, cost in zip(self.__cornerEdges, cornerEdgeCost):
|
|
rEdge.updateFlow(cost)
|
|
|
|
@staticmethod
|
|
def __indexList(a,b):
|
|
from operator import itemgetter
|
|
|
|
return itemgetter(*b)(a)
|
|
|
|
# Test function
|
|
def isUnwrapped(self):
|
|
for tri in self.loops:
|
|
if tri.isUnwrapped():
|
|
continue
|
|
else:
|
|
return False
|
|
return True
|
|
|
|
# Plot functions
|
|
def plotResult(self, fileName):
|
|
from mpl_toolkits.mplot3d import Axes3D
|
|
import matplotlib.pyplot as plt
|
|
|
|
fig = plt.figure()
|
|
plt.title('Phase as 3D plot')
|
|
ax = fig.add_subplot(111, projection='3d')
|
|
for v in self.__vertices:
|
|
ax.scatter(v.x, v.y, v.phase, marker='o', label='Wrapped Phase')
|
|
ax.scatter(v.x, v.y, v.getUnwrappedPhase(), marker='^', c='r', label='Unwrapped Phase')
|
|
plt.savefig(fileName)
|
|
|
|
def plotDelaunay(self, fileName):
|
|
import matplotlib.pyplot as plt
|
|
fig = plt.figure()
|
|
ax = fig.add_subplot(1,1,1)
|
|
plt.title('Delaunay Traingulation with Residue and non-zero edge cost')
|
|
ax.triplot(self.__x, self.__y, self.__delauneyTri.simplices.copy(), c='b', linestyle='dotted')
|
|
|
|
cmap = plt.get_cmap('hsv')
|
|
numComponents = np.unique(self.__compNum).shape[0]
|
|
colors = [cmap(i) for i in np.linspace(0, 1, numComponents)]
|
|
|
|
ax.plot(self.__x, self.__y, '.')
|
|
|
|
for v, edge in self.__edges.items():
|
|
edge.plot(ax)
|
|
for tri in self.loops:
|
|
tri.plot(ax)
|
|
|
|
plt.savefig(fileName)
|
|
|
|
def plotSpanningTree(self, fileName):
|
|
import matplotlib.pyplot as plt
|
|
fig = plt.figure()
|
|
ax = fig.add_subplot(1,1,1)
|
|
ax.plot(self.__x, self.__y, '.')
|
|
|
|
for edge, rEdge in self.__unwrapEdges:
|
|
if edge.flow == rEdge.flow:
|
|
edge.plot(ax, c='g')
|
|
else:
|
|
edge.plot(ax, c='r')
|
|
plt.savefig(fileName)
|
|
|
|
# Removes collinear points on the convex hull and removing duplicates
|
|
def filterPoints(x, y):
|
|
def indexPoints(points):
|
|
# Dirty: To be cleaned later
|
|
# Numbering points
|
|
pt2Idx = {}
|
|
for idx, pt in enumerate(points):
|
|
pt2Idx[pt] = idx
|
|
return list(pt2Idx)
|
|
|
|
# Removes duplicate
|
|
points = list(zip(x, y))
|
|
points = np.array(indexPoints(points))
|
|
|
|
return (points[:,0], points[:,1])
|
|
|
|
def wrapValue(x):
|
|
if((x <= np.pi) and (x > -np.pi)):
|
|
return x
|
|
elif(x > np.pi):
|
|
return wrapValue(x - 2*np.pi)
|
|
else:
|
|
return wrapValue(x + 2*np.pi)
|
|
|
|
# Command Line argument parser
|
|
def firstPassCommandLine():
|
|
|
|
# Creating the parser for the input arguments
|
|
parser = argparse.ArgumentParser(description='Phase Unwrapping')
|
|
|
|
# Positional argument - Input XML file
|
|
parser.add_argument('--inputType', choices=['plane', 'sinc', 'sine'],
|
|
help='Type of input to unwrap', default='plane', dest='inputType')
|
|
parser.add_argument('--dim', type=int, default=100,
|
|
help='Dimension of the image (square)', dest='dim')
|
|
parser.add_argument('-c', action='store_true',
|
|
help='Component-wise unwrap test', dest='compTest')
|
|
parser.add_argument('-MCF', action='store_true',
|
|
help='Minimum Cost Flow', dest='mcf')
|
|
parser.add_argument('--redArcs', type=int, default=0,
|
|
help='Redundant Arcs', dest='redArcs')
|
|
parser.add_argument('--solver', choices=['glpk', 'pulp', 'gurobi'],
|
|
help='Type of solver', default='pulp', dest='solver')
|
|
|
|
# Parse input
|
|
args = parser.parse_args()
|
|
return args
|
|
|
|
def main():
|
|
# Parsing the input arguments
|
|
args = firstPassCommandLine()
|
|
inputType = args.inputType
|
|
dim = args.dim
|
|
compTest = args.compTest
|
|
mcf = args.mcf
|
|
|
|
print("Input Type: %s"%(inputType))
|
|
print("Dimension: %d"%(dim))
|
|
|
|
# Random seeding to reapeat random numbers in case
|
|
np.random.seed(100)
|
|
|
|
inputImg = np.empty((dim, dim))
|
|
if inputType == 'plane':
|
|
# z = a*x + b*y
|
|
a = np.random.randint(0,10,size=1)/50
|
|
b = np.random.randint(0,10,size=1)/50
|
|
for i in range(dim):
|
|
for j in range(dim):
|
|
inputImg[i,j] = (a*i) + (b*j)
|
|
|
|
elif inputType == 'sinc':
|
|
mag = np.random.randint(1,10,size=1)
|
|
n = np.random.randint(1,3,size=1)
|
|
for i in range(dim):
|
|
for j in range(dim):
|
|
inputImg[i,j] = mag * np.sinc(i*n*np.pi/100) * np.sinc(i*n*np.pi/100)
|
|
|
|
elif inputType == 'sine':
|
|
mag = np.random.randint(1,10,size=1)
|
|
n = np.random.randint(1,3,size=1)
|
|
for i in range(dim):
|
|
for j in range(dim):
|
|
inputImg[i,j] = mag * np.sin(i*n*np.pi/100) * np.sin(i*n*np.pi/100)
|
|
|
|
if compTest:
|
|
# Component wise unwrap testing
|
|
n1 = np.random.randint(0,10,size=4)
|
|
n2 = np.random.randint(0,10,size=4)
|
|
compImg = np.empty((dim, dim))
|
|
wrapImg = np.empty((dim, dim))
|
|
for i in range(dim):
|
|
for j in range(dim):
|
|
compImg[i,j] = (i > (dim/2)) + 2*(j > (dim/2))
|
|
n = n1[compImg[i,j]] + n2[compImg[i,j]]
|
|
wrapImg[i,j] = inputImg[i,j] + n*(2*np.pi)
|
|
wrapImg = np.array(wrapImg)
|
|
compImg = np.array(compImg)
|
|
else:
|
|
# Wrapping input image
|
|
wrapImg = np.array([[wrapValue(x) for x in row] for row in inputImg])
|
|
compImg = None
|
|
|
|
# Choosing random samples
|
|
xidx = np.random.randint(0,dim,size=400).tolist()
|
|
yidx = np.random.randint(0,dim,size=400).tolist()
|
|
xidx, yidx = filterPoints(xidx, yidx)
|
|
|
|
# Creating the Minimum Cost Flow problem
|
|
if mcf is True:
|
|
# We use redArcs to run Minimum Cost Flow
|
|
redArcs = -1
|
|
print('Relax IV used as the solver for Minimum Cost Flow')
|
|
solver = None
|
|
else:
|
|
redArcs = args.redArcs
|
|
solver = args.solver
|
|
|
|
if compImg is None:
|
|
phaseunwrap = PhaseUnwrap(xidx, yidx, wrapImg[xidx, yidx], redArcs=redArcs)
|
|
else:
|
|
phaseunwrap = PhaseUnwrap(xidx, yidx, wrapImg[xidx, yidx], compImg[xidx, yidx], redArcs=redArcs)
|
|
|
|
# Including the neutral node for min cost flow
|
|
phaseunwrap.solve(solver)
|
|
|
|
# Unwrap
|
|
phaseunwrap.unwrapLP()
|
|
#phaseunwrap.plotDelaunay("final%d.png"%(redArcs))
|
|
#phaseunwrap.plotSpanningTree("spanningTree%d.png"%(redArcs))
|
|
phaseunwrap.plotResult("final%d.png"%(redArcs))
|
|
|
|
#import pdb; pdb.set_trace()
|
|
#fig = plt.figure()
|
|
#ax = fig.add_subplot(111, projection='3d')
|
|
|
|
if __name__ == '__main__':
|
|
|
|
# Main Engine
|
|
main()
|