Dataset¶
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.
- Importamos las librerÃas necesarias
In [1]:
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 de una pagina web en tiempo real.
In [2]:
#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)
Out[2]:
In [3]:
# como el archivo es un comprimido, lo procedemos a descomprimir
decompress = "unzip -P '12345' datos.zip"
os.system(decompress)
Out[3]:
- Extraemos algunos parámetros dentro de las imagenes Dicom con pydicom, y lo guardamos en un DataFrame
In [4]:
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])
In [5]:
df.head()
Out[5]:
In [6]:
# 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])
- Convertimos nuestros datos a formato .root
In [7]:
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
- Utilizamos las funciones de Root en python para hacer histogramas de los datos.
In [8]:
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
Out[8]:
In [9]:
# Creamos una carpeta para guardar las imagenes creadas
crear_directorio = "mkdir I_Presentacion"
os.system(crear_directorio)
Out[9]:
In [10]:
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")
In [11]:
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")
In [12]:
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")
In [13]:
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")
- Trabajaremos con la matriz de pixeles de una imagen de un paciente del conjunto de Datos.
In [14]:
ds = pydicom.read_file("datos/paciente16-A/00000054-A") # Lenemos la imagen
matriz = ds.pixel_array # Extraemos la matriz de pixeles
- Vemos algunas propiedades de la imagen y la abrimos con Matplotlib.
In [15]:
print(type(matriz))
print(matriz.shape)
plt.imshow(matriz, cmap=plt.cm.bone)
Out[15]:
- Extraemos los coeficiente de atenuación de rayos X en una lÃnea (1D) de la imagen, para guardarla en un archivo root y proceder a graficarla
In [16]:
#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
In [17]:
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")
- Manipulamos la matriz de pixeles para poder graficar en relieve la imagen de CT utilizando root.
In [18]:
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)
In [19]:
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
In [20]:
Houns = ROOT.TFile.Open("Pixel.root") # Leemos el archivo root con los valores de atenuacion
tree = Houns.Get("tree")
tree.GetEntries()
Out[20]:
In [21]:
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()
- Gráfica en relibe
In [22]:
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")
In [ ]: