'''
License
    This file is part of foam-extend.

    foam-extend is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by the
    Free Software Foundation, either version 3 of the License, or (at your
    option) any later version.

    foam-extend is distributed in the hope that it will be useful, but
    WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Author
    Matthias Rauter matthias.rauter@uibk.ac.at
'''

import numpy as np
import shapefile 
 
#pdi = self.GetPolyDataInput()
#block = pdi.GetBlock(0).GetBlock(0)


# find source
polyname = list(GetSources().keys())[list(GetSources().values()).index(GetActiveSource())][0]
contour1 = FindSource(polyname)
cdata = servermanager.Fetch(contour1) 

block = cdata.GetBlock(0).GetBlock(0)

#print block
#print block.GetNumberOfPoints()
#print block.GetNumberOfCells() 




import PyQt4.QtGui
polyname = str(PyQt4.QtGui.QInputDialog.getText(None, 'Enter polygon name', '*.shp', text=polyname)[0])

p = []

for i in range(block.GetNumberOfPoints()):
    p.append(block.GetPoint(i))

p = np.array(p) 

z = []
for i in range(block.GetNumberOfCells()):
    #print block.GetCell(i).GetPointIds().GetId(0), block.GetCell(i).GetPointIds().GetId(1)
    z.append([block.GetCell(i).GetPointIds().GetId(0), block.GetCell(i).GetPointIds().GetId(1)])

 
polys = []
closed = []
pointsleft = True

while(pointsleft):
    poly = []
    #print 'new polygon'
    i = 0
    while i < len(z):
        if z[i][0] > -1:
            poly.append(z[i][0])
            poly.append(z[i][1])
            #print '    found beginning'
            del z[i]
            if i == len(z)-1:
                pointsleft = False
            break
        i = i+1
    reverseSearch = True
    if len(poly) > 0:
        i = 0
        while i < len(z):
            if z[i][0] == poly[-1]:
                poly.append(z[i][1])
                del z[i]
                #print '        found point, ', len(z)
                i = 0
            elif z[i][1] == poly[-1]:
                poly.append(z[i][0])
                del z[i]
                #print '        found point, ', len(z)
                i = 0
            else:
                i = i+1 
            if poly[-1] == poly[0]:
                break
            if poly[-1] != poly[0] and i == len(z) and reverseSearch:
                reverseSearch = False
                poly = list(reversed(poly))
                i = 0
    polys.append(poly)
    closed.append(poly[-1] == poly[0])
    if len(z) < 1:
        break

bp = []
ep = []
closedpolys = []

for i in range(len(polys)):
    bp.append(p[polys[i][0]])
    ep.append(p[polys[i][-1]]) 
i = 0
while i < len(polys):
    if closed[i]:
        closedpolys.append(polys[i])
        i = i+1
    else:
        mindist = float('inf')
        mindistpoly = -1
        reverse = False
        for j in range(i+1, len(polys)):
            if not closed[j]:
                d = np.linalg.norm(np.array(ep[i])-np.array(bp[j]))
                if d < mindist:
                    mindist = d
                    mindistpoly = j
                    reverse = False
                d = np.linalg.norm(np.array(ep[i])-np.array(ep[j]))
                if d < mindist:
                    mindist = d
                    mindistpoly = j
                    reverse = True
        if mindistpoly == -1:
            polys[i].append(polys[i][0])
            closedpolys.append(polys[i])
            closed[i] = True
            i = i+1
        else:
            if reverse:
                polys[mindistpoly] = list(reversed(polys[mindistpoly]))
            polys[i] = polys[i] + polys[mindistpoly]
            ep[i] = p[polys[i][-1]]
            del polys[mindistpoly]
            del bp[mindistpoly]
            del ep[mindistpoly]
            del closed[mindistpoly]


polys = closedpolys
#for poly in polys:
#    plt.plot(p[poly,0], p[poly,1])
#plt.axis('equal')
#plt.show()

polys2 = []
for poly in polys:
    polys2.append([[p[poly[i],0], p[poly[i],1]] for i in range(len(poly))])


w = shapefile.Writer()
w.line(parts=polys2)

w.field('FIRST_FLD','C','40')
w.field('SECOND_FLD','C','40')
w.record('First','Polygon')

w.save(polyname)

