<jupyter_start><jupyter_text>Surface Determination Workflow<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
import os
import h5py
import sys
from tomo_encoders import Patches, DataFile
import tensorflow as tf
import time, glob
import open3d as o3d

sys.path.append('/data02/MyArchive/aisteer_3Dencoders/TomoEncoders/scratchpad/surface_determination/trainer')
from tomo_encoders.neural_nets.surface_segmenter import SurfaceSegmenter
from tomo_encoders.misc.feature_maps_vis import view_midplanes
from tomo_encoders.misc.voxel_processing import normalize_volume_gpu
from tomo_encoders.misc.viewer import view_midplanes
from tqdm import tqdm
gt_path = '/data02/MyArchive/tomo_datasets/ZEISS_try2/GT_VOLS/Sample2'
ct_path = '/data02/MyArchive/tomo_datasets/ZEISS_try2/Sample2_CT'
psize = (16,16,16)
size_lab = '16x16'
ds_gt = DataFile(gt_path, tiff = True, VERBOSITY = 0)
vol_gt = ds_gt.read_full()
# t00 = time.time()
# p = Patches(vol_gt.shape, initialize_by = 'regular-grid', patch_size = psize)
# vols = p.extract(vol_gt, psize)
# ystd = np.std(vols, axis = (1,2,3))
# idxs = np.where(ystd > 0.0)[0]
# vols = vols[idxs,...]
# p = p.select_by_indices(idxs)
# p.dump(f'/data02/MyArchive/tomo_datasets/ZEISS_try2/sample2_edge_patches_{size_lab}.hdf5')
# print('time taken: {time.time() - t00}')
p = Patches(None, initialize_by = 'file', fpath = f'/data02/MyArchive/tomo_datasets/ZEISS_try2/sample2_edge_patches_{size_lab}.hdf5')
y_true = p.extract(vol_gt, psize)
vol_out = np.zeros(p.vol_shape)
vols = np.ones((len(p),) + psize)
p.fill_patches_in_volume(vols, vol_out)
fig, ax = plt.subplots(1,3, figsize = (16,8))
view_midplanes(vol_out, ax = ax)<jupyter_output><empty_output><jupyter_text>Measure IoU in datasets<jupyter_code>from tomo_encoders.misc.img_stats import calc_SNR, calc_jac_acc, calc_dice_coeff
from tomo_encoders import DataFile
import pandas as pd
flist = glob.glob('/data02/MyArchive/tomo_datasets/ZEISS_try2/Sample2_SEG/*')
ct_path = '/data02/MyArchive/tomo_datasets/ZEISS_try2/Sample2_CT'
cols = ['X', 'CT-algo', 'SURF-model', 'IoU', 'SNR', 'ftag']
df = pd.DataFrame(columns = cols)
for ii, fpath in enumerate(flist):
    ftag = os.path.split(fpath)[-1]
    r1, r2, r3 = ftag.split('_')
    
    if 'FDK' in r1:
        ct_algo = 'FDK'
    else:
        ct_algo = 'DR'
    
    X = int(r1.split(ct_algo)[-1].split('X')[0])
    surf_model = r3.split('a0')[-1]
    
    
    # IoU
    ds = DataFile(fpath, tiff = True, VERBOSITY = 0)
    vol_seg = ds.read_full()
    y_seg = p.extract(vol_seg, psize)
    iou = calc_jac_acc(y_true, y_seg)

    # SNR
    ct_fname = os.path.join(ct_path, f'{ct_algo}{X}X')
    ds_ct = DataFile(ct_fname, tiff = True, VERBOSITY = 0)
    vol_ct = ds_ct.read_full()
    snr = calc_SNR(vol_ct[150:-150,300:-300,300:-300], vol_gt[150:-150,300:-300,300:-300])
    # Dice
#     calc_dice_coeff(y_true, y_seg)    
    
    new_item = dict(zip(cols, [X,ct_algo,surf_model,iou,snr, ftag]))
    print(f'{new_item}')
    df = df.append(new_item,ignore_index = True)
    
    
    df.to_csv('/home/atekawade/Dropbox/Arg/transfers/ZEISS_v2/stats.csv', index = False)

    
df
# pc = np.asarray(np.where(vol_dr != vol_fdk)).T
# pcd = o3d.geometry.PointCloud()
# pcd.points = o3d.utility.Vector3dVector(pc[:,::-1])
# o3d.visualization.draw_geometries([pcd])
# from skimage import measure
# verts, faces, normals, values = measure.marching_cubes(vol_fdk[::4,::4,::4], 0)
# mesh = o3d.geometry.TriangleMesh()
# mesh.vertices = o3d.utility.Vector3dVector(verts)
# mesh.triangles = o3d.utility.Vector3iVector(faces)
# mesh.compute_vertex_normals()
# o3d.visualization.draw_geometries([mesh])<jupyter_output><empty_output>