# -*- coding: utf-8 -*-
"""
Created on Wed Mar 14 11:18:57 2018

@author: Tom
"""

# -*- coding: utf-8 -*-
"""
Created on Mon Mar 12 15:18:46 2018

@author: Tom
"""
import numpy as np
import porespy as ps
import matplotlib.pyplot as plt
import skimage as sk
import OpenPNM as op

plt.close('all')
np.random.seed(1)

blobs = ps.generators.blobs(100)



def edge_val(arr, value):
    for ax in range(3):
        for side in range(2):
            side *= -1
            for edge in range(2):
                edge *= -1
                arr.swapaxes(0, ax)[side, edge, :] = value
                arr.swapaxes(0, ax)[side, :, edge] = value
    return arr

snow_im = ps.network_extraction.snow(blobs)
regions = snow_im.regions
r_shape = np.shape(regions)
#im = np.zeros([r_shape[0]+2, r_shape[1]+2, r_shape[2]+2], dtype=int)
n_id = regions.max() + 1
print(n_id)

boundary_pores = []
im = np.pad(regions, 1, 'edge')
for ax in range(3):
    for side in range(2):
        side *= -1
        face = im.swapaxes(0, ax)[side, :, :]
        print(ax, side, 'pores', len(np.unique(face)))
        for p_id in np.unique(face):
            face[face == p_id] = n_id
            boundary_pores.append(n_id)
            n_id += 1
bps = np.asarray(boundary_pores)-1

blobs_pad = np.pad(blobs, 1, 'edge')
blobs_pad = edge_val(blobs_pad, False)

plt.figure()
plt.imshow(im[50, :, :])

im[~blobs_pad] = 0
print(np.sum(im==118))
net = ps.network_extraction.extract_pore_network(im, voxel_size=1e-6)
pn = op.Network.GenericNetwork()
pn.update(net)
#pn['pore.boundary']=False
#pn['pore.boundary'][bps]=True

pn.trim(pn.check_network_health()['trim_pores'])
op.Network.tools.find_surface_pores(pn)
op.Utilities.IO.VTK.save(pn, 'extracted_blobs')