In [2]:
%pylab inline
import io3d
import skelet3d
from pathlib import Path
import copy
import glob
import time
import pandas as pd
# import sed
In [3]:
def vessel_skeleton_extraction(pth1):
datap = io3d.read(pth1, datap=True)
from skimage.filters import threshold_otsu
threshold = threshold_otsu(datap["data3d"].ravel()[::1000])
imthr = datap["data3d"] > threshold
# imshow(imthr[int(imthr.shape[0]/2),:,:])
volume_data = imthr
skelet = skelet3d.skelet3d(volume_data)
return skelet, volume_data, datap["voxelsize_mm"]
def skeleton_analysis(skelet, volume_data, voxelsize_mm):
skan = skelet3d.SkeletonAnalyser(skelet, volume_data=volume_data, voxelsize_mm=voxelsize_mm)
stats = skan.skeleton_analysis()
df = stats_as_dataframe(skan)
return df
def extract_df(df):
dfs = df[["id",
"nodeA_ZYX 0",
"nodeA_ZYX 1",
"nodeA_ZYX 2",
"nodeB_ZYX 0",
"nodeB_ZYX 1",
"nodeB_ZYX 2",
"nodeA_ZYX_mm 0",
"nodeA_ZYX_mm 1",
"nodeA_ZYX_mm 2",
"nodeB_ZYX_mm 0",
"nodeB_ZYX_mm 1",
"nodeB_ZYX_mm 2",
"radius_mm",
"connectedEdgesA 0",
"connectedEdgesA 1",
"connectedEdgesA 2",
"connectedEdgesB 0",
"connectedEdgesB 1",
"connectedEdgesB 2",
]]
dfs.to_csv(pth1 + ".csv")
return dfs
In [4]:
def stats_as_dataframe(self):
import pandas as pd
import exsu.dili
if self.stats is None:
msg = "Run skeleton_analyser before stats_as_dataframe()"
logger.error(msg)
raise RuntimeError(msg)
# import imtools.dili
df = pd.DataFrame()
for stats_key in self.stats:
one_edge = copy.copy(self.stats[stats_key])
k = "orderedPoints_mm_X"
if k in one_edge:
one_edge[k] = str(one_edge[k])
k = "orderedPoints_mm_Y"
if k in one_edge:
one_edge[k] = str(one_edge[k])
k = "orderedPoints_mm_Z"
if k in one_edge:
one_edge[k] = str(one_edge[k])
k = "orderedPoints_mm"
if k in one_edge:
one_edge[k] = str(one_edge[k])
# one_edge[]
one_dct = exsu.dili.flatten_dict_join_keys(one_edge, simplify_iterables=True)
# df_one = pd.DataFrame(one_dct)
df_one = pd.DataFrame([list(one_dct.values())], columns=list(one_dct.keys()))
df = df.append(df_one, ignore_index=True)
return df
In [5]:
# pth = Path(io3d.datasets.join_path("medical/processed/porcine_liver_ct_raw", get_root=True))
pths = glob.glob("C:/Users/Jirik/Downloads/porcine_liver_ct_raw/*.mhd")
# pth1 = pth / "P01_MakroCT_HEAD_5_0_H31S_0004.mhd"
pth1 = pths[-1]
# pth1 = pths[2]
pths
Out[5]:
In [6]:
# pth1 = pth / "P01_MakroCT_HEAD_5_0_H31S_0004.mhd"
# pth1 = "C:/Users/Jirik/Downloads/porcine_liver_ct_raw/P01_MakroCT_HEAD_5_0_H31S_0004.mhd"
# pth1 = r"C:/Users/Jirik/Downloads/porcine_liver_ct_raw/P01_a_MikroCT-nejhrubsi_rozliseni_DICOM_liver-1st-important_Macro_pixel-size53.0585um.mhd"
# pth1 = Path("g:\Můj disk\data\medical\processed\porcine_liver_ct_raw\P01_MakroCT_HEAD_5_0_H31S_0004.mhd")
print(pth1)
print(Path(pth1).exists())
t0 = time.time()
datap = io3d.read(pth1, datap=True)
print(time.time()-t0)
In [7]:
import scipy.stats
data3d = datap["data3d"]
dsc = scipy.stats.describe(data3d.ravel()[::100])
# print(f"mn={np.min(datap['data3d'])}, mx={np.max(datap['data3d'])} ")
dsc
Out[7]:
In [8]:
from skimage.filters import threshold_otsu, threshold_multiotsu
# classes = 5 if Path(pth1).name == 'P01_MakroCT_HEAD_5_0_H31S_0004.mhd'else 2
# thresholds = threshold_multiotsu(datap["data3d"].ravel()[::1000], classes=classes)
# threshold = thresholds[-1]
# classes
if Path(pth1).name in ('P01_MakroCT_HEAD_5_0_H31S_0004.mhd', 'P01_MakroCT_po_rozrezani_HEAD_0_6_H20S_0003.mhd'):
threshold = 0
data3d = data3d[:,:400,:] # cut the table
else:
threshold = threshold_otsu(data3d.ravel()[::1000])
In [9]:
imthr = data3d > threshold
imshow(imthr[int(imthr.shape[0]/2),:,:])
Out[9]:
In [10]:
volume_data = imthr
t0 = time.time()
skelet = skelet3d.skelet3d(volume_data)
print(time.time()-t0)
# skan = skelet3d.SkeletonAnalyser(skelet, volume_data=volume_data, voxelsize_mm=datap["voxelsize_mm"])
t0 = time.time()
skan = skelet3d.SkeletonAnalyser(
skelet,
# volume_data=volume_data,
voxelsize_mm=datap["voxelsize_mm"]
)
print(time.time()-t0)
# skan.
In [11]:
skan.sklabel.dtype
Out[11]:
In [12]:
# sklabel_nodes = skan.sklabel.copy()
# sklabel_nodes[sklabel_nodes > 0] = 0
nz = np.nonzero(skan.sklabel < 0)
nz
Out[12]:
In [13]:
voxelsize_mm = datap["voxelsize_mm"]
dfs = pd.DataFrame({
"node_0_px":nz[0],
"node_1_px":nz[1],
"node_2_px":nz[2],
"node_0_mm":nz[0]*voxelsize_mm[0],
"node_1_mm":nz[1]*voxelsize_mm[1],
"node_2_mm":nz[2]*voxelsize_mm[2],
})
dfs.to_csv(pth1 + ".nodes.csv")
# df["Z_px"]
In [14]:
# io3d.write({"data3d"})
np.savez(pth1 + ".npz", sklabel=skan.sklabel, imthr=imthr)
In [15]:
# dir(skan)#.__connection_analysis
# skan._SkeletonAnalyser__connection_analysis(3)
In [ ]:
t0 = time.time()
stats = skan.skeleton_analysis()
print(time.time()-t0)
In [12]:
df = stats_as_dataframe(skan)
In [32]:
df.keys()
Out[32]:
In [25]:
datap["voxelsize_mm"]
Out[25]:
In [26]:
dfs = df[["id",
"nodeA_ZYX 0",
"nodeA_ZYX 1",
"nodeA_ZYX 2",
"nodeB_ZYX 0",
"nodeB_ZYX 1",
"nodeB_ZYX 2",
"nodeA_ZYX_mm 0",
"nodeA_ZYX_mm 1",
"nodeA_ZYX_mm 2",
"nodeB_ZYX_mm 0",
"nodeB_ZYX_mm 1",
"nodeB_ZYX_mm 2",
"radius_mm",
"connectedEdgesA 0",
"connectedEdgesA 1",
"connectedEdgesA 2",
"connectedEdgesB 0",
"connectedEdgesB 1",
"connectedEdgesB 2",
]]
dfs
Out[26]:
In [28]:
dfs.to_csv(pth1 + ".csv")
In [36]:
df[["id", "lengthEstimation", "radius_mm", "tortuosity", 'phiAc', "phiBa", "phiBb", 'connectedEdgesA 0', 'connectedEdgesA 1',
'connectedEdgesB 0', 'connectedEdgesB 1',
'connectedEdgesA 2',
'connectedEdgesA 3',
'connectedEdgesB 2', 'connectedEdgesB 3']]
Out[36]:
In [45]:
In [8]:
import glob
pths = glob.glob("C:/Users/Jirik/Downloads/porcine_liver_ct_raw/*.mhd")
# pth1 = pth / "P01_MakroCT_HEAD_5_0_H31S_0004.mhd"
# pth1
pths
Out[8]:
In [11]:
# for pth1 in pths:
# vesse
outputs = vessel_skeleton_extraction(pths[0])
# df = keleton_analysis(*outputs)
# dfs = extract_df(df)
In [12]:
outputs
Out[12]: