Introdução ao processamento de imagens médicas com Python: CT Lung e segmentação de embarcações sem rótulos

Hora de algum tutorial prático sobre imagens médicas. No entanto, desta vez não usaremos a IA louca, mas algoritmos básicos de processamento de imagens. O objetivo é familiarizar o leitor com conceitos em torno de imagens médicas e tomografia computadorizada especificamente (CT).

É fundamental entender o quão longe alguém pode ir sem Aprendizagem profunda, para entender quando é melhor usá -lo. Novos profissionais tendem a ignorar essa parte, mas a análise de imagem médica ainda é o processamento de imagens em 3D.

Também incluo partes do código para facilitar a compreensão do meu processo de pensamento.

O notebook do Google Colab que acompanha pode ser encontrado aqui Para executar o código mostrado neste tutorial. UM Girub O repositório também está disponível. Estrela nosso repositório se você gostou!

Para se aprofundar na maneira como a IA é usada na medicina, você não pode errar com o Você tem para medicina Curso on -line, oferecido pela Coursera. Se você deseja se concentrar na análise de imagem médica com aprendizado profundo, eu recomendo começar a partir do Curso de Udemy, baseado em Pytorch.

Começaremos com o básico da imagem da TC. Você pode pular esta seção se já estiver familiarizado com a imagem da TC.

Imagem de TC

Física das tomografias

A tomografia computadorizada (TC) usa feixes de raios-X para obter intensidades de pixels 3D do corpo humano. Um cátodo aquecido libera vigas de alta energia (elétrons), que por sua vez liberam sua energia como radiação de raios-X. Os raios X passam pelos tecidos do corpo humano e atingem um detector do outro lado. Um tecido denso (ou seja, ossos) absorverá mais radiação do que os tecidos moles (ou seja, gordura). Quando os raios X não são absorvidos do corpo (ou seja, na região do ar dentro dos pulmões) e atingem o detector, os vemos como preto, semelhante a um filme preto. No oposto, os tecidos densos são retratados como brancos.

Desta maneira, A imagem da TC é capaz de distinguir diferenças de densidade e criar uma imagem 3D do corpo.

Fonte: Christopher P. Hess, MD, Ph.D e Derk Purcell, MD, Departamento de Radiologia e Imagem Biomédica na UCSF

Aqui está um vídeo de 1 min que achei muito conciso:

https://www.youtube.com/watch?v=l9swbatrrbg

Intensidades de TC e unidades de Hounsfield

A absorção de raios-X é medida no Escala de Hounsfield. Nesta escala, fixamos a intensidade do ar para -1000 e a água para 0 intensidade. É essencial entender que Housenfield é uma escala absoluta, diferentemente da ressonância magnética, onde temos uma escala relativa de 0 a 255.

A imagem ilustra alguns dos tecidos básicos e seus valores de intensidade correspondentes. Lembre -se de que as imagens são barulhentas. Os números podem variar ligeiramente em imagens reais.

A escala de Hounsfield. Imagem por autor.

Os ossos têm alta intensidade. Normalmente, prendemos a imagem para ter um alcance máximo superior. Por exemplo, o valor máximo pode ser de 1000, por razões práticas.

O problema: as bibliotecas de visualização funcionam na escala (0,255). Não seria muito aconselhável visualizar toda a escala de Hounsfield (de -1000 a 1000+) a 256 escalas para diagnóstico médico.

Em vez disso, limitamos nossa atenção a diferentes partes desse intervalo e focamos nos tecidos subjacentes.

Visualização de dados de TC: Nível e janela

A convenção de imagem médica para prender a linha Housenfield é escolhendo uma intensidade central, chamada de nível e uma janela, conforme descrito:

Na verdade, é uma convenção bastante feia para cientistas da computação. Gostaríamos apenas do Min e do Max do intervalo:

$max = nível + janela/2$

$min = nível – janela/2$

import matplotlib.pyplot as plt
import numpy as np
def show_slice_window(slice, level, window):
   """
   Function to display an image slice
   Input is a numpy 2D array
   """
   max = level + window/2
   min = level - window/2
   slice = slice.clip(min,max)
   plt.figure()
   plt.imshow(slice.T, cmap="gray", origin="lower")
   plt.savefig('L'+str(level)+'W'+str(window))

Se você não estiver convencido, a próxima imagem o convencerá de que a mesma imagem de CT é mais rica do que um canal de imagem comum:

Janela e nivelamento no CT podem fornecer imagens bastante diferentes. Imagem por autor.

Para referência, aqui está uma lista de faixas de visualização:

Região/tecido	Janela	Nível
cérebro	80	40
pulmões	1500	-600
fígado	150	30
Tecidos moles	250	50
osso	1800	400

Hora de jogar!

Segmentação pulmonar com base em valores de intensidade

Não vamos apenas segmentar os pulmões, mas também encontraremos a área real em $mm^2$. Para fazer isso, precisamos encontrar o tamanho real das dimensões do pixel. Cada imagem pode ter um diferente (Pixdim no arquivo de cabeçalho Nifty). Vamos ver o arquivo de cabeçalho primeiro:

import nibabel as nib
ct_img = nib.load(exam_path)
print(ct_img.header)

Aqui vou mostrar apenas alguns campos importantes do cabeçalho:

<class 'nibabel.nifti1.Nifti1Header'> object, endian='<'
sizeof_hdr      : 348
dim             : (  2 512 512   1   1   1   1   1)
datatype        : int16
bitpix          : 16
pixdim          : (1.   0.78515625 0.78515625   1.   1.   1.    1.   1.  )
srow_x          : ( -0.78515625   0.           0.         206.60742   )
srow_y          : (  0.          -0.78515625   0.         405.60742   )
srow_z          : (   0.     0.     1.  -304.5)

Para o registro, SROW_X, SROW_Y, SROW_Z é o Affine Matrix da imagem. Bitpix é quantos bits usamos para representar cada intensidade de pixels.

Então, vamos definir uma função que lê essas informações do arquivo de cabeçalho. Com base no formato bacana, cada dimensão no arquivo Nifty tem uma dimensão pixel. O que precisamos é descobrir os índices de duas dimensões da imagem e suas respectivas dimensões de pix.

Etapa 1: Encontre as dimensões do pixel para calcular a área em mm^2

def find_pix_dim(ct_img):
    """
    Get the pixdim of the CT image.
    A general solution that gets the pixdim indicated from the image dimensions. From the last 2 image dimensions, we get their pixel dimension.
    Args:
        ct_img: nib image
    Returns: List of the 2 pixel dimensions
    """
    pix_dim = ct_img.header("pixdim") 
    dim = ct_img.header("dim") 
    max_indx = np.argmax(dim)
    pixdimX = pix_dim(max_indx)
    dim = np.delete(dim, max_indx)
    pix_dim = np.delete(pix_dim, max_indx)
    max_indy = np.argmax(dim)
    pixdimY = pix_dim(max_indy)
    return (pixdimX, pixdimY) 

Etapa 2: imagem binarize usando limiar de intensidade

Esperamos que os pulmões estejam no alcance da unidade de Housendfield (-1000, -300). Para esse fim, precisamos prender o intervalo de imagens para (-1000, -300) e binarize os valores para 0 e 1, então obteremos algo assim:

Imagem por autor

Etapa 3: descoberta de contorno

Vamos esclarecer o que é um contorno antes de qualquer outra coisa:

Para visão computacional, um contorno é um conjunto de pontos que descrevem uma linha ou área. Portanto, para cada contorno detectado, não obteremos uma máscara binária completa, mas um conjunto com um monte de valores X e Y.

Ok, como podemos isolar a área desejada? Hmmm .. Vamos pensar sobre isso. Nós nos preocupamos com as regiões pulmonares que são mostradas em branco. Se pudéssemos encontrar um algoritmo para identificar conjuntos fechados ou qualquer tipo de contorno na imagem que possa ajudar. Depois de algumas pesquisas online, encontrei o Método de quadrados em marcha que encontra contornos constantes de valor em uma imagem de Skimagechamado skimage.measure.find_contours().

Depois de usar esta função, visualizo os contornos detectados na imagem original do CT:

Imagem por autor

Aqui está a função!

def intensity_seg(ct_numpy, min=-1000, max=-300):
   clipped = clip_ct(ct_numpy, min, max)
   return measure.find_contours(clipped, 0.95)

Etapa 4: encontre a área pulmonar de um conjunto de possíveis contornos

Observe que eu usei uma imagem diferente para mostrar uma caixa de borda que o corpo do paciente não é um conjunto de pontos fechado. OK, não exatamente o que queremos, mas vamos ver se poderíamos resolver isso.

Como nos preocupamos apenas com os pulmões, precisamos definir algum tipo de restrições para excluir as regiões indesejadas.

Para fazer isso, extraí um polígono convexo do contorno usando Scipy. Depois de assumir 2 restrições:

Isso pode ou não incluir o contorno do corpo, resultando em mais de 3 contornos. Quando isso acontece, o corpo é facilmente descartado com o maior volume do contorno que satisfaz as suposições pré-descritas.

def find_lungs(contours):
    """
    Chooses the contours that correspond to the lungs and the body
    First, we exclude non-closed sets-contours
    Then we assume some min area and volume to exclude small contours
    Then the body is excluded as the highest volume closed set
    The remaining areas correspond to the lungs
    Args:
        contours: all the detected contours
    Returns: contours that correspond to the lung area
    """
    body_and_lung_contours = ()
    vol_contours = ()
    for contour in contours:
        hull = ConvexHull(contour)
        if hull.volume > 2000 and set_is_closed(contour):
            body_and_lung_contours.append(contour)
            vol_contours.append(hull.volume)
    if len(body_and_lung_contours) == 2:
        return body_and_lung_contours
    elif len(body_and_lung_contours) > 2:
        vol_contours, body_and_lung_contours = (list(t) for t in
                zip(*sorted(zip(vol_contours, body_and_lung_contours))))
        body_and_lung_contours.pop(-1) 
    return body_and_lung_contours 

Como um caso de borda, estou mostrando que o algoritmo não se restringe a apenas duas regiões dos pulmões. Inspecione o contorno azul abaixo:

Imagem por autor

Etapa 5: contorno para máscara binária

Em seguida, o salvamos como um arquivo bacana, para que precisemos converter o conjunto de pontos em uma máscara binária pulmonar. Para isso, usei o travesseiro Python Lib que desenha um polígono e cria uma máscara de imagem binária. Então eu fundi todas as máscaras dos contornos do pulmão já encontrados.

import numpy as np
from PIL import Image, ImageDraw
def create_mask_from_polygon(image, contours):
    """
    Creates a binary mask with the dimensions of the image and
    converts the list of polygon-contours to binary masks and merges them together
    Args:
        image: the image that the contours refer to
        contours: list of contours
    Returns:
    """
    lung_mask = np.array(Image.new('L', image.shape, 0))
    for contour in contours:
        x = contour(:, 0)
        y = contour(:, 1)
        polygon_tuple = list(zip(x, y))
        img = Image.new('L', image.shape, 0)
        ImageDraw.Draw(img).polygon(polygon_tuple, outline=0, fill=1)
        mask = np.array(img)
        lung_mask += mask
    lung_mask(lung_mask > 1) = 1  
    return lung_mask.T  

A área pulmonar desejada em $mm^2$ é simplesmente o número de elementos diferentes de zero multiplicados pelas duas dimensões de pixels da imagem correspondente.

As áreas pulmonares são salvas em um arquivo CSV junto com o nome da imagem.

Finalmente, para salvar a máscara como bacana, usei o valor de 255 para a área pulmonar em vez de 1 para poder exibir em um visualizador bacana. Além disso, salvo a imagem com a transformação afim da fatia de CT inicial para poder ser exibida de maneira significativa (alinhada sem conflitos de rotação).

def save_nifty(img_np, name, affine):
    """
    binary masks should be converted to 255 so it can be displayed in a nii viewer
    we pass the affine of the initial image to make sure it exits in the same
    image coordinate space
    Args:
        img_np: the binary mask
        name: output name
        affine: 4x4 np array
    Returns:
    """
    img_np(img_np == 1) = 255
    ni_img = nib.Nifti1Image(img_np, affine)
    nib.save(ni_img, name + '.nii.gz')

Finalmente, abri a máscara com um visualizador bacana comum para o Linux validar que tudo correu bem. Aqui estão instantâneos para a fatia número 4:

Imagem por autor

Eu usei um espectador de imagem médica gratuita chamada Aliza no Linux.

Segmentar os navios principais e calcular a proporção da área dos navios sobre os navios

Se houver um pixel com um valor de intensidade acima de -500 hu dentro da área do pulmão, consideraremos -o como um vaso.

Primeiro, fazemos multiplicação de elemento entre a imagem da CT e a máscara pulmonar para obter apenas os pulmões. Posteriormente, definimos os zeros que resultaram da multiplicação no elemento para -1000 (ar em Hu) e, finalmente, mantemos apenas as intensidades maiores que -500 como embarcações.

def create_vessel_mask(lung_mask, ct_numpy, denoise=False):
    vessels = lung_mask * ct_numpy  
    vessels(vessels == 0) = -1000
    vessels(vessels >= -500) = 1
    vessels(vessels < -500) = 0
    show_slice(vessels)
    if denoise:
        return denoise_vessels(lungs_contour, vessels)
    show_slice(vessels)
    return vessels

Uma amostra deste processo pode ser ilustrada abaixo:

A máscara de embarcação com algum barulho. Imagem por autor.

Analisar e melhorar o resultado da segmentação

Como você pode ver, temos algumas partes do contorno dos pulmões, que acredito que gostaríamos de evitar. Para esse fim, eu criei um função de denoising Isso considera a distância da máscara a todos os pontos de contorno. Se estiver abaixo de 0,1, eu defino o valor do pixel como 0 e, como resultado, os exclua dos vasos detectados.

def denoise_vessels(lung_contour, vessels):
    vessels_coords_x, vessels_coords_y = np.nonzero(vessels)  
    for contour in lung_contour:
        x_points, y_points = contour(:, 0), contour(:, 1)
        for (coord_x, coord_y) in zip(vessels_coords_x, vessels_coords_y):
            for (x, y) in zip(x_points, y_points):
                d = euclidean_dist(x - coord_x, y - coord_y)
                if d <= 0.1:
                    vessels(coord_x, coord_y) = 0
    return vessels

Abaixo, você pode ver a diferença entre a imagem denoise à direita e a máscara inicial:

Imagem por autor

Se sobrepormos a máscara na imagem CT original que obtemos:

Imagem por autor

def overlay_plot(im, mask):
    plt.figure()
    plt.imshow(im.T, 'gray', interpolation='none')
    plt.imshow(mask.T, 'jet', interpolation='none', alpha=0.5)

Agora que temos a máscara, a área da embarcação é calculada semelhante ao que eu fiz pelos pulmões, levando em consideração a dimensão de pixel de imagem individual.

def compute_area(mask, pixdim):
    """
    Computes the area (number of pixels) of a binary mask and multiplies the pixels
    with the pixel dimension of the acquired CT image
    Args:
        lung_mask: binary lung mask
        pixdim: list or tuple with two values
    Returns: the lung area in mm^2
    """
    mask(mask >= 1) = 1
    lung_pixels = np.sum(mask)
    return lung_pixels * pixdim(0) * pixdim(1)

Os índices são armazenados em um arquivo CSV no notebook.

Conclusão e leituras adicionais

Acredito que agora você tem um entendimento sólido das imagens de TC e suas particularidades. Podemos fazer muitas coisas impressionantes com informações tão ricas de imagens 3D.

Apoie -nos olhando nosso projeto em Girub. Finalmente, para tutoriais mais avançados, consulte nossa imagem médica Artigos.

Foi-me pedido um grande curso de IA prático. Estou pensando em escrever um livro sobre imagens médicas em 2021. Até então, você pode aprender com o curso Coursera Você tem para medicina.

* Divulgação: Observe que alguns dos links acima podem ser links de afiliados e, sem nenhum custo adicional, ganharemos uma comissão se você decidir fazer uma compra depois de clicar.