El dataset es un conjunto de datos de imágenes de tomografÃa computarizada(CT), tomadas de un Hospital de Ecuador, las imágenes fueron anonimizadas con anterioridad.
import ROOT
import root_pandas
%jsroot on
import pydicom
import sys
import os
import pandas as pd
import numpy as np
from os import listdir
from os.path import isfile, join
import os
from pandas import DataFrame
import matplotlib.pyplot as plt
import array
#Descargamos los datos que esta alojados en la direccion, los datos estan en formato .zip
bashCommand = "wget http://201.159.223.36/experimental/jenn_files/recursos/datos.zip -O datos.zip"
os.system(bashCommand)
# como el archivo es un comprimido, lo procedemos a descomprimir
decompress = "unzip -P '12345' datos.zip"
os.system(decompress)
def ruta_carpeta(carpeta, contenido_carpeta):
# genera una ruta
ruta = carpeta + '/' + contenido_carpeta
return ruta
def return_id(primer_ident, segundo_ident):
# ingresa los numeros de identificacion y retorna listos para leer.
primer = '0x'+ str(primer_ident)
segundo = '0x'+ str(segundo_ident)
return primer, segundo
def return_data(ruta_imagen, primer_ident, segundo_ident):
# se ingresa los numero con los id de la informacion
img = pydicom.dcmread(ruta_imagen)
id1, id2 = return_id(primer_ident, segundo_ident)
nombre_dato = img[id1, id2].name
valor_dato = img[id1, id2].value
return nombre_dato, valor_dato
path = 'datos' # directorio donde se encuentran las imagenes.
vector_1 = []
vector_2 = []
vector_3 = []
vector_4 = []
vector_5 = []
vector_6 = []
vector_7 = []
vector_8 = []
vector_9 = []
vector_10 = []
vector_11 = []
vector_12 = []
vector_13 = []
vector_14 = []
directorio_principal = os.listdir(path) # leemos los archivos dentro del directorio.
for contenido in directorio_principal: # recorremos el directorio.
ruta1 = ruta_carpeta(path , contenido)
imagenes = os.listdir(ruta1)
ruta2 = ruta_carpeta(ruta1 , imagenes[0])
ruta3 = ruta_carpeta(ruta1 , imagenes[1])
###################################################################################
nombre_1, valor_1 = return_data(ruta2, 10, 10) # nombre paciente
vector_1.append(valor_1)
###################################################################################
nombre_2, valor_2 = return_data(ruta2, 8, 60) # modalidad
vector_2.append(valor_2)
###################################################################################
nombre_6, valor_6 = return_data(ruta2, 18, 1150 ) # time exposure
vector_6.append(valor_6)
###################################################################################
nombre_7, valor_7 = return_data(ruta2, 18, 1170 ) # generator power
vector_7.append(valor_7)
###################################################################################
nombre_8, valor_8 = return_data(ruta2, 28, 10 ) # rows
vector_8.append(valor_8)
###################################################################################
nombre_9, valor_9 = return_data(ruta2, 28, 11 ) # columns
vector_9.append(valor_9)
###################################################################################
nombre10, valor_10 = return_data(ruta2, 20, 1041 ) # Slice Location
vector_10.append(valor_10)
nombre_10=nombre10+'_i'
###################################################################################
nombre11, valor_11 = return_data(ruta3, 20, 1041 ) # Thickness
vector_11.append(valor_11)
nombre_11=nombre11+'_f'
###################################################################################
nombre_12, valor_12 = return_data(ruta2, 18, 50 ) # columns
vector_12.append(valor_12)
###################################################################################
nombre_13, valor_13 = return_data(ruta2, 28, 1053 ) # columns
vector_13.append(valor_13)
###################################################################################
nombre_14='Number_of_images'
valor_14= abs((valor_11-valor_10)/valor_12)
vector_14.append(valor_14)
Data = {nombre_1: vector_1,
nombre_2: vector_2,
nombre_6: vector_6,
nombre_7: vector_7,
nombre_8: vector_8,
nombre_9: vector_9,
nombre_10: vector_10,
nombre_11: vector_11,
nombre_12: vector_12,
nombre_13: vector_13,
nombre_14: vector_14
}
df = DataFrame(Data,columns=[ nombre_1,nombre_2, nombre_6, nombre_7, nombre_8, nombre_9,
nombre_10, nombre_11, nombre_12, nombre_13, nombre_14])
df.head()
# Definimos un DataFrame con los valores que se utilizara para el análisis
Data2 = {'Modalidad':vector_2,
'Time':vector_6,
'Power':vector_7,
'Thickness':vector_12,
nombre_14:vector_14}
df2 = DataFrame(Data2,columns=['Modalidad','Thickness','Time','Power',nombre_14])
df2.to_root('Datos_num.root', key='mytree') # Guardamos en formato root, utilizando root_pandas
f2 = ROOT.TFile.Open("Datos_num.root") # Abrimos el file.root
canvas = ROOT.TCanvas("Canvas","a first way to plot a variable",400,300)# Definimos el canvas
tree = f2.Get("mytree") # Definimos el tree
tree.GetEntries() # Obtenemos los datos del file.root
# Creamos una carpeta para guardar las imagenes creadas
crear_directorio = "mkdir I_Presentacion"
os.system(crear_directorio)
hist1 = ROOT.TH1F("variable","Exposure Time; Time; Events ",50,0,2000) # definimos el histograma
for event in tree:
hist1.Fill(tree.Time) # Definimos la variable ha hacer el histograma
print("Done!")
hist1.SetFillColor(2) # Definimos los colores a usarse
hist1.Draw() # Dibujamos el Histograma
canvas.Draw() # Plasmamos en Histograma en el canvas
canvas.SaveAs("I_Presentacion/Histo_Time.png")
hist2 = ROOT.TH1F("variable","Thickness; Thickness; Events ",50,0,10)
for event in tree:
hist2.Fill(tree.Thickness) # Definimos la variable ha hacer el histograma
print("Done!")
hist2.SetFillColor(2) # Definimos los colores a usarse
hist2.Draw() # Dibujamos el Histograma
canvas.Draw() # Plasmamos en Histograma en el canvas
canvas.SaveAs("I_Presentacion/Histo_Thickness.png")
hist3 = ROOT.TH1F("variable","Generator Power; Generator Power; Events ",50,0,80000)
for event in tree:
hist3.Fill(tree.Power) # Definimos la variable ha hacer el histograma
print("Done!")
hist3.SetFillColor(2) # Definimos los colores a usarse
hist3.Draw() # Dibujamos el Histograma
canvas.Draw() # Plasmamos en Histograma en el canvas
canvas.SaveAs("I_Presentacion/Histo_Power.png")
hist2 = ROOT.TH1F("variable","Histograma del Numero de Imagenes; Numero de Imagenes; Events ",50,0,200)
for event in tree:
hist2.Fill(tree.Number_of_images) # Definimos la variable ha hacer el histograma
print("Done!")
hist2.SetFillColor(20) # Definimos los colores a usarse
hist2.Fit('gaus') # ajustamos a una gaussiana
hist2.Draw() # Dibujamos el Histograma
canvas.Draw() # Plasmamos en Histograma en el canvas
canvas.SaveAs("I_Presentacion/Histo_Number_gaus.png")
ds = pydicom.read_file("datos/paciente16-A/00000054-A") # Lenemos la imagen
matriz = ds.pixel_array # Extraemos la matriz de pixeles
print(type(matriz))
print(matriz.shape)
plt.imshow(matriz, cmap=plt.cm.bone)
#extraemos los valores para la atenuacion de rayos X en 1D,
# para la coordenada 300
x = array . array ("d",range (511))
y = array . array ("d",range (511))
for i in range(0,511):
y[i] = matriz[345][i]-1024 # extraemos los valores y los transformamos a hounsfiels
g = ROOT . TGraph ( len (x) ,x,y)
g. Draw ("AL")
g.GetXaxis().SetTitle("X")
g.GetYaxis().SetTitle("Unidades Hounsfield")
g.SetTitle("Unidades Hounsfied en un linea de la Imagen CT (1D)")
canvas.Draw("Grafica")
canvas.SaveAs("I_Presentacion/1D.png")
x = np.arange(0, matriz.shape[0], 1)
y = np.arange(0, matriz.shape[1], 1)
def HU_2D(matriz,x,y):
'''Funcion que extrae cada valor proporcional al
coeficiente de atenuacion de los rayos X con sus respectivos
valores en X y Y,y los retorna en forma de arrays.
'''
Houns_list = []
x_list = []
y_list = []
for i in x:
for j in y:
Houns = matriz[x[i],y[j]] -1024.
Houns_list.append(Houns)
x_list.append(x[i])
y_list.append(y[j])
return np.array(Houns_list),np.array(x_list),np.array(y_list)
H, x, y = HU_2D(matriz,x,y)
Data2 = {'x':x,
'y':y,
'H':H
}
pixel = DataFrame(Data2,columns=['x','y','H'])
pixel.to_root('Pixel.root', key='tree') # guardamos en formato root
Houns = ROOT.TFile.Open("Pixel.root") # Leemos el archivo root con los valores de atenuacion
tree = Houns.Get("tree")
tree.GetEntries()
hist10 = ROOT.TH1F("variable",
"Histograma de la cantidad de materiales en la imagen; Housnfield; Numero de pixel ",
50,-2000,2000)
for event in tree:
hist10.Fill(tree.H)
print("Done!")
hist10.SetFillColor(40)
hist10.Draw()
canvas.SaveAs("I_Presentacion/Histo_houns.png")
canvas.Draw()
import ROOT
Number_of_points = len(H)
graph = ROOT.TGraph2D(Number_of_points)
for i in range(Number_of_points): graph.SetPoint(i, x[i], y[i],H[i] )
c = ROOT.TCanvas("c")
graph.GetXaxis().SetTitle("X")
graph.GetYaxis().SetTitle("Y")
graph.GetZaxis().SetTitle("Unidades Hounsfield")
graph.SetTitle("Grafica en Relieve de una Imagen de CT.")
graph.SetMarkerStyle(20); graph.Draw("PCOL Z")
c.Draw()
c.SaveAs("I_Presentacion/3D.png")