Questo algoritmo di image processing nasce dalla necessità di avere un volume 3D a partire dalle slices di una TC. Online si trovano diversi algoritmi ma la maggior parte richiedono immagini DICOM. In questo caso, il dataset di riferimento, è composto da 970 immagini in png. Inoltre, nonostante l’immagine sia effettivamente in scala di grigi, il file contiene tutti e tre i canali RGB, vedremo infatti che verrà convertito. Non trovando ne tool ne algoritmi disponibili ho deciso di crearmelo da solo.
Per creare l’algoritmo ho utilizzato Matlab e i tool di elaborazione di immagini. Le immagini sono prese da un dataset pubblicato dall’ESA a seguito di uno studio sugli effetti della gravità sulla massa ossea. Dovrò sfruttare alcuni dei dati estratti per un progetto di omogenizzazione dell’osso trabecolare.
In generale, l’algoritmo prevede un pre-processing delle singole immagini volto a filtrare il rumore e rendere più smooth le impronte delle trabecole. Successivamente si procede a binarizzare effettivamente l’immagine e poi si va a rendere l’immagine 3D impilando le diverse fette. Poi si va ad elaborare l’immagine 3D. Inoltre, ho aggiunto anche un esempio di meshing sfruttando Mathematica e relative funzionalità. Questo permette, almeno in linea teorica, di implementare la mesh anche in un’analisi strutturale FEM. Chiaramente il limite è la potenza di calcolo e l’elevatissimo numero di gradi di libertà che si genera, ma di questo ne parlo meglio nel progetto precedentemente citato.
Image Processing
Le immagini presenti nel dataset purtroppo non sono di ottima qualità. Hanno una buona risoluzione ma presentano tanto rumore. Inoltre, sono immagini raster, e questo rende complesso passare ad un volume 3D. In questo caso, per ottenere il volume 3D, sarebbe meglio avere un set di immagini binarie che rappresentano il vuoto (pixel==0) e le trabecole ossee(pixel==1). L’immagine ha una risoluzione di 2048×2048 ottenuta da una CT con uno Scanco μCT Scanner dove il singolo voxel ha una dimensione di 37 μm.
Prima di partire con il pre-processing possiamo importare tutte le immagini e visualizzarne la sequenza:
Per farlo andiamo a considerare le 10 cartelle e carichiamo tutti i file all’interno di una matrice tridimensionale. Quindi esportiamo la sequenza come video.
selpath = uigetdir; %select image processing folder
folderpath=strcat(selpath,'/file_raw/slice*');
folderlist=dir(folderpath);
for k=1:length(folderlist)
temp_list=dir(strcat(selpath,'/file_raw/',folderlist(k).name,'/ESA*.png'));
if k==1
file_list=append(strcat(selpath,'/file_raw/',folderlist(k).name,'/'),{temp_list.name});
else
file_list=[file_list, append(strcat(selpath,'/file_raw/',folderlist(k).name,'/'),{temp_list.name})];
end
end
xlsfiles=sort(file_list); %sort by name ascending
file1=imread(xlsfiles{1}); %first file for store measures
imm=[zeros(numel(file1(1,:)),numel(file1(:,1)))];
for i=1:length(xlsfiles) %pre allocate volumetrix matrix
imm(:,:,i)=0;
end
imload=imm;
parfor i=1:length(xlsfiles) %load images
imload(:,:,i)=mat2gray(imread(xlsfiles{i}));
end
% export video with tic-toc timing (for benchmarking)
video=VideoWriter(strcat(selpath,'/flow.mp4'),'MPEG-4'); %default 30 fps
open(video); %open the file for writing
tic
for i=1:length(imload(1,1,:))
writeVideo(video,imload(:,:,i));
end
close(video);
tocAllora possiamo caricare tutti i dati all’interno della matrice imload(i,j,k), dove sul piano i,j abbiamo la singola immagine, ESA29-99-L3-Byte.0***.png.
Preprocessing
Sulla base dell’obiettivo dell’algoritmo possiamo iniziare aumentando il contrasto. Considerando i tre indici di contrasto noti in letteratura, è evidente che soltanto CI2 può esserci utile. Questo perchè il primo e il terzo soffrono dei pixel in saturazione che sono evidenti nel set, quindi saranno entrambi stabili tra le immagini.
\operatorname{CI1}={\max ({I})-\min( I) }\\\:\\ \operatorname{CI2}={\sigma}(I)\\\:\\\operatorname{CI3}={\max ({I})-\min( I)\over\max ({I})+\min( I) }Il secondo indice ci permette invece di scegliere una, o più immagini per fare alcune considerazioni.
Quindi c’è un ciclo che va a trovare dei valori ottimi di aumento di contrasto. Il metodo scelto per aumentare il contrasto è lo stretching dell’istogramma a cui passiamo le due soglie, minima e massima, calcolandole dai quantili e il grado di non linearità.
g(x,y)=\begin{cases}
0\quad \text{se }I(x,y)< th_1\\
1\quad \text{se }I(x,y)> th_2\\
\left(I(x,y)-th_1\over th_2-th_1\right)^\gamma \quad \text{se }th_1\leq (x,y)\leq th_2
\end{cases}Nel fare questa operazione non possiamo tenere conto solo della deviazione standard ma dobbiamo tenere conto anche, e soprattutto, del risultato desiderato. Questo porta alla selezione dei parametri di partenza per l’ottimizzazione.
Chiaramente è stata ottimizzata anche la dimensione della finestra di filtraggio.
%%
clear CII2 row col val val_g gamma m M MAT val_g val_g_i
gamma=1.2:0.2:2;
m=1:2:40;
M=97:1:99;
CII2=zeros(m(end),M(end),10*gamma(end));
for i=1:length(m) %lower threeshold
for j=1:length(M) %upper threeshold
for k=1:length(gamma) %linearity coefficient
ima_c=imadjust(ima,[m(i)/100 M(j)/100],[0 1],gamma(k));
CII2(m(i),M(j),k)=std(ima_c(:))/std(ima(:));
end
end
end
%search best index
for k=1:length(gamma)
MAT=CII2(:,:,k);
[row(k),col(k)]=find(MAT==max(max(MAT)));
end
for k=1:length(gamma)
val_g(k)=CII2(row(k),col(k),k);
end
val_g_i=find(val_g==(max(val_g)));
val1=row(val_g_i)
val2=col(val_g_i)
val3=gamma(val_g_i)
%test contrast
contrastim=imadjust(ima,[row(val_g_i)/100 col(val_g_i)/100],[0 1],gamma(val_g_i));
figure;
imshowpair(imload(:,:,best_contrast_index),contrastim,'montage')
mex=strcat('m=',string(val1),', M=',string(val2),', gamma=',string(val3));
title(mex)
%% contrast set
imcontr=imm;
for k=1:length(imload(1,1,:))
imcontr(:,:,k)=imadjust(imload(:,:,k),[row(val_g_i)/100 col(val_g_i)/100],[0 1],gamma(val_g_i));
end
%% filtering test
L=6;
mask1=fspecial('average',L);
mask2=fspecial('gaussian',L,(L-1)/6);
ima_f1=imfilter(ima,mask1);
figure;
imshowpair(ima,ima_f1,'montage')
ima_f2=imfilter(ima,mask2);
figure;
imshowpair(ima,ima_f2,'montage')
%% filtering
imfilt=imm;
for k=1:length(imload(1,1,:))
imfilt(:,:,k)=imfilter(imcontr(:,:,k),mask1);
endQuesto aumento di contrasto va nella direzione opposta alla seconda operazione, ovvero il denoising. Volendo rimuovere del rumore possiamo provare con un filtro gaussiano o un filtro di media. Generata la maschera vediamo come otteniamo un risultato migliore con il filtro di averaging. Non serve esagerare, anzi, esagerando andremo a rimuovere l’aumento di contrasto. Nonostante possa sembrare un controsenso è bene sottolineare che l’insieme di queste due procedure permette di avere un rendering 3D più smooth e meno spigoloso. Inoltre, permette di avere più margine sulla soglia di binarizzazione.
Binarizzazione
Allora è possibile procedere a binarizzare l’immagine considerando una soglia tale da garantire sia una buona qualità delle trabecole sia una sufficiente rimozione del background. Una buona idea per ottimizzare questa soglia è considerare i dati in letteratura sulle massa trabecolare, corrispondente al 20% del volume (vedi progetto di omogenizzazione).
figure()
imshow(imfilt(:,:,best_contrast_index));
taglio=drawfreehand(gca);
BW=createMask(taglio);
figure;
imshow(BW);
%...%
imbin=imm;
parfor i=1:length(imload(1,1,:))
imbin(:,:,i)=imbinarize(imfilt(:,:,i),0.46);
imbin(:,:,i)=imbin(:,:,i).*BW;
imwrite(imbin(:,:,i),strcat(selpath,'\binarized\',string(i),'.png'))
endNello stesso ciclo usiamo anche una maschera disegnata manualmente per filtrare gli elementi di bordo e esportiamo le singole immagini, per usarle poi in Mathematica.
3D rendering
A questo punto, avendo ottimizzato tutti i passaggi precedenti, è possibile procedere al rendering 3D sfruttando il toolbox VolumeViewer. Il tool è abbastanza intuito ed è dotato anche di app per settare i diversi parametri tramite GUI.
figure()
volshow(imbin,config)
Pore fraction
Possiamo andare anche a verificare l’effettiva frazione di pori rispetto le trabecole. Per farlo possiamo sfruttare il tool Image Region Analyzer che ci permettere di individuare le proprietà di singole aree nell’immagine binarizzata. Prendiamo quindi la slice no. 650 e ripuliamola con una maschera BW. Leggiamone le proprietà con il tool e andiamo a prendere la somma delle singole aree identificate, questo ci darà il numero totale di pixel di osso trabecolare (area).
slice=imread('binarized/600.png');
sliceb=imbinarize(slice); %effective logic image
figure()
imshow(sliceb);
taglio=drawfreehand(gca);
BW=createMask(taglio);
sliceb=sliceb.*BW; %clean image
area=sum(propsTable.Area);%total area calculated by tool exported properties
AA=imclose(sliceb,strel('disk',20));
AA1=imfill(AA,8,'holes');
figure;imshowpair(AA,AA1);
total=numel(AA1(find(AA1==1)));
fraction=area/total
Quindi possiamo usare dei comandi per riempire l’area interna le trabecole e questo ci darà il numero totale di pixel (total). Quindi il loro rapporto del 27.98% si ritrova conforme ai dati in letteratura.
Mesh
Lo stesso risultato è possibile ottenerlo anche in Mathematica. Inoltre, un algoritmo più semplice anche se meno preciso, potrebbe essere quello di passare direttamente alla binarizzazione, saltando il preprocessing. Il risultato è sicuramente diverso e meno preciso ma può essere una buona alternativa per velocizzare i tempi.
files = FileNames["*.png",
"C:\\Users\\bigba\\OneDrive\\Documenti\\GitHub\\bone\
homogenization\\image_processing\\file_raw\\slices600-699\\"];
DDD = Table[Import[files[[n]]], {n, 1, 90}];
(* riprendi dalla binarizzazione in matlab*)
Do[DDD[[n]] = Binarize[ColorConvert[DDD[[n]], "Grayscale"], 0.8], {n, 1, 90}];
(*
Imm=Binarize[ColorConvert[Imm,"Grayscale"],0.8]
Imm=Binarize[ColorConvert[Imm2,"Grayscale"],0.8];
*)
A = Image3D[DDD]
res = Timing[MESH = ImageMesh[A]];
Export[FileNameJoin[{NotebookDirectory[], "mesh.png"}], MESH]
NotebookDirectory[]
FileNameJoin[{NotebookDirectory[], "mesh.png"}]
res[[1]]/60(*minute to generate mesh*)
Timing[Export[FileNameJoin[{NotebookDirectory[], "mesh.png"}], MESH,ImageResolution -> 300]]
Timing[Export[FileNameJoin[{NotebookDirectory[], "mesh.svg"}], MESH]]
Timing[Export[FileNameJoin[{NotebookDirectory[], "mesh.obj"}], MESH]]
Benchmarking
Questi algoritmi non sono molto leggeri vista la quantità di dati analizzati (2048x2048x1000 livelli) e della diverse operazioni. Ad esempio, la creazione del video ha richiesto 150,25 sec. Per l’aumento di contrasto 225,40 sec e per il filtraggio 220,11 sec. Le parti computazionalmente più pesanti sono state il caricamento con 714,27 sec e le riassegnazioni. È stato necessario anche cancellare più volte alcune variabili per liberare memoria.
In Mathematica sono stati necessari 10.07 secondi per realizzare la mesh su 200 slices (500÷700) ma oltre 20 minuti per generare il render grafico e l’export. Questo a causa della scheda video che risultava continuamente saturata.
Questi dati sono relativi ad una Workstation con i7-108050H, 32 GB RAM DDR4 2933 MHz, SSD NVMe, NVIDIA Quadro P620.
Scopri di più sull’osso trabecolare
Il codice completo è disponibile alla repository del progetto:
