Laboratori Matlab - Informatica medica

provate

% ogni volta a scrivere char(ans) sotto ogni fread, vedrete che troverete

% delle informazioni o delle lettere che per voi non vogliono dire niente

% ma sono indicazioni precise dentro il file originale

g= fread(f,8);

char(g); %troverete SIEMENS

% ES 3

function posizione = trovaelemento (f, group, element)

%codifico il TAG in binario da esadecimale

%tag di esempio 0008, 0070) trovo 8 e 0 come se leggessi da destra il mio

%numero

V = [ hex2dec(group(3:4)) hex2dec(group(1:2)) hex2dec(element(3:4))

hex2dec(element(1:2))];

%mi posiziono all'inizio, mi sposto di 0 byte rispetto all'origine del file

%che per convenzione è rappresentato dal -1

fseek(f, 0, -1);

A = fread(f);

%è un vettore posizione e lo sto inizializzando a 0 e ogni volta lo aggiorno

posizione=[];

for i = 1:length(A)-3

%verifica che la dimensione di A sia uguale al vettore V

if isequal (A(i:i+3)',V)

%se è vera, concatena al vettore posizione il vettore corrente di i

posizione = [posizione i];

end

end

% ES 4

posizione = trovaelemento (f,'0008', '0020');

fseek(f, 536, -1);

fread(f, 4);

fread(f,2);

char(ans)';

fread(f,2);

char(ans)';

fread(f,8);

char(ans)'; %esce la data "20120516", 16/05/2012

posizione = trovaelemento (f, '7FE0', '0010') % questa è l'immagine in se, ve l'ha

messa

% qui più che altro per trarre in inganno. L'informazione che vi da è che

% troverete l'immagine al byte 68031, è anche un suggerimento per la domanda n. 6,

% che comunque vedremo prima come l'ha vista lei, ricordatevi di tornare

% qui più tardi

% ES 5

posizione = trovaelemento (f, '0020', '000D');

% Uso il primo numero che esce da "posizione" nel comando sottostante

fseek(f, 2168, -1);

fread(f,4);

dec2hex(ans); % questa è la stringa che corrisponde a 0020,000D, non ho capito se

richiedeva lei

% o se chiedeva questa

fread(f,2);

char(ans)'; %da cui esce "UI"

fread(f,2);

char(ans)';

fread(f, 54); %è la dimensione del contenuto "scritto" di quel tag

char(ans)';

% ES 6

I=dicomread('image.dcm');

imshow(I, 'DisplayRange' , []); %ti fa vedere l'immagine

dicomdisp('image.dcm'); %QUESTO COMANDO TI FA VEDERE LA STRUTTURA DEL FILE DI TESTO

% SU CUI HAI LAVORATO PRATICAMENTE A VUOTO FINO AD ORA, TE LO METTO QUI NON

% PERCHè SONO STRONZA, MA ADESSO PUOI RIVEDERE TUTTO QUELLO FATTO AVENDO A

% MENTE QUESTO SCHEMA, A ME è SERVITO MOLTO IMPARARE PRIMA I COMANDI E LA

% SEQUENZA, E POI CAPIRE SU COSA STESSI EFFETTIVAMENTE LAVORANDO PER POI

% RIVEDERE IL TUTTO

% Comunque, ti va direttamente alla fine e nell'ultima riga puoi vedere "Pixel-

Data" che è l'effettiva immagine e ti

% dice che occupa 524288 bytes, quindiper sapere quanti bytes occupa il file

% di testo senza immagine non faccio altro che togliere questo numero alla

% lunghezza totale del file

length(A) - 524288; % viene fuori 68042 occupati dal file, che, come visto

nell'esercizio 4

% corrisponde al byte al quale inizia l'immagine trovato tramite il comando

% "posizione=trovaelemento(..)"

% ES 7

load('ecg.mat');

figure;

plot(ECGtm, ECG);

hold on;

plot(ECGtm(ann), ones(1,10), 'ro');

plot(ECGtm(ann), ECG(ann), 'ro');

figure;

subplot(1,2,1), plot(ECGtm(1:1000), ECG(1:1000));

subplot(1,2,2), plot(ECGtm(1001:2000), ECG(1001:2000));

LABORATORIO 3

% ES 1

%estrarre dal file ICD10.tt i codici e le realative descrizioni per la

%categoria G000-G09: "Infiammatory disases of the central nervous sistem"

f=fopen('ICD10.txt');

lista=['G00'; 'G01'; 'G02'; 'G03'; 'G04'; 'G05'; 'G06'; 'G07'; 'G08'; 'G09'];

n=0;

%ottengo una matrice di 10 rigehe per 3 colonne????

%inizializzo a 0 una variabile che conterà i termini che effettivamente

%andiamo a estrarre in quanto facenti parte della categoria voluta

%utilizzo il while per leggere l'intero documento, quella funzione valuta

%se siamo al termine della lettura del file

while feof(f)==0

line=fgetl(f);

if isempty(strmatch(line(2:4),lista)) == 0

%prendo i primi caratteri delle nostra riga appena letta e la

%confronto con la lista, la funzione strmatch mi consente

%questo e mi resituisce la posizione dell'occorrenza della

%nostra stringa, se il risultato è un vettore non vuoto la

%funzione isempty diversa da 0, non è vuoto hai trovato

%qualcosa allora ho trovato un'istanza dei termini che cercavo

%e vado a contarli

if isempty(findstr('-', line))

n=n+1;

[code, descrizione] = strtok(line);

%sto utilizzando la funzione prende la stringa in ingresso,

%line nel nostro caso, salta eventuali spazi all'inizio, nel

%primo argomento mette tutti i caratteri prima del carattere di

%tab, nel secondo argomento vengono memorizzati i caratteri restanti

s(n) = struct('code', code, ...

'descrizione', descrizione);

end

end

end

fclose(f)

%ho delle righe di intestazione che voglio eliminare, contengono il

%trattino qindi aggiungo un if

% ES 2

%creare un array strutture contenente tutti i termini che contengono la

%parola 'deficiency'

f=fopen('ICD10.txt');

n=0;

while feof(f)==0

line=fgetl(f);

if isempty(findstr(' deficiency', ...

line))==0;

n=n+1;

%estrae e copia il code e descrizione

[code, descrizione]= strtok(line);

d(n) = struct('code', code, 'descrizione', descrizione);

end

end

% ES 3 – ECG

load('ecg.mat')

%trovare e aggiungere al grafico i relativi marcatori del picco S che segue

%il picco R nel battito

D=diff(ECG); %segnale differenziale

S=[]; %vettore che conterrà le posizioni dei picchi S

for i=1:length(ann)

%ann contiene i picchi R annotati manualmente da un esperto

d=D(ann(i) + 3:end);

positivi=find(d>0);

s=positivi(1)-1; %prendo il primo dei positivi che troviamo

%parto dal picco R, il segnale scende, valuto la differenza tra un

%campionee il precedente, questa sarà sempre negativa fin quando non

%trovo il minimo, mi interessa il primo minimo e questo sarà la

%posizione del mio s

s=s+ann(i)+3;

S= [S s];

end figure

plot (ECGtm, ECG)

hold on

plot(ECGtm(S), ECG(S), 'mx')

LABORATORIO 4

% ROC

%PUNTO 1

%il soggetto è considerato affetto se il marcatore è >= 1,3, rispetto allo

%standard del disease calcolare sensitività, specificità e accuratezza del

%test

load('ex1.mat');

POS=find(disease);

NEG=find(disease==0);

n=length(disease);

TP= length(find(concentration(POS)>1.3));

TN=length(find(concentration(NEG)<=1.3));

SENSITIVITY = TP/length(POS);

SPECIFICITY = TN/length(NEG);

ACCURACY = (TN + TN)/n;

% PUNTO 2 - Vado a fare la stessa cosa, ma andando a variare la soglia nel

% ciclo for, questo per valutare il marcatore (vedere teoria info med +

% segnali)

soglie =unique(concentration);

ROCx=[];

ROCy=[];

for i=1:length(soglie)

TP= length(find(concentration(POS)>soglie(i)));

TN=length(find(concentration(NEG)<=soglie(i)));

SENSITIVITY = TP/length(POS);

SPECIFICITY = TN/length(NEG);

ROCx=[ROCx, 1-SPECIFICITY];

ROCy=[ROCy, SENSITIVITY];

end

figure

plot(ROCx, ROCy)

% stiamo così valutando il marcatore in base a TUTTI i valori soglia

xlabel('1-specificity')

ylabel('sensitivity')

% Kaplan-Meier

Appunto personale a posteriori, questo è il laboratorio più difficile, questa prima parte in verde è presa dagli audio che

avevo del laboratorio, è sicuramente poco chiara, ma ho preferito comunque metterla a disposizione, ma se vi fa

confusione provata a partire direttamente dal testo.

% censored fa riferimento ai soggetti che sono all'interno della casistica

% all'inizio del periodo di osservazione, ma poi non si hanno informazioni

% relative a questi soggetti, questi partecipano al compito della curva di

% Kaplan-Meier ma vengono esclusi nel momento in cui non si hanno più

% informazioni rispetto al loro follow-up

% tempi in settimane per l'evento critico

% man mano che ho eventi il numero dei soggetti a rischio diminuisce, anche

% a causa dei censored

% la curva è una produttoria delle probabilità che riguardano la sopravvivenza

% a un certo tempo

% ciascuno avrà il proprio tempo 0, che può essere una certa diagnosi, tutti

% sono riferiti allo stesso istante iniziale rispetto al quale parte il

% calcolo del tempo, quindi avrò tempi cronologici diversi per i quali ho lo

% stesso tipo di informazione a partire da una certa diagnosi > evoluzione

% partiamo da 21 soggetti

N=21;

EVENT=[6 6 6 7 10 13 16 22 23];

censored = [6 9 10 11 17 19 20 25 25 32 32 34];

% tempi che ho nella casistica + 0 iniziale

t1=[0 unique(EVENT)];

% numero di eventi per ogni tempo

eventi1=[];

for i = 1:length(t1), eventi1=[eventi1 length(find(EVENT==t1(i)))]; end

%aggiorna il vettore che è in t1 della lunghezza del vettore trovato

% ora consideriamo i soggetti censored e li andiamo a contare intervallo per

% intervallo, inizializziamo a vuoto il vettore c1

c1=[];

for i = 1:length(t1) - 1, c1=[c1 length(find(t1(i)<= censored) && (censored <

t1(i+1)))]; end

% quanto la curva scende è solo in corrispondenza solo dei tempi critici

% memorizzati in event, i censored sono all'interno del gradino e vengono

% raffigurati lungo esso. Gradino che corrisponde ad un intervallo tra due

% eventi critici. Inoltre andiamo ad aggiungere a c1

c1= [c1 length(find(censored>=23))]; %23 tempo max di event

% parto con un 100% di sogetti, questo può diminuire in corrispondenza degli

% eventi critici

nt=0;

% soggetti a rischio nel j-esimo intervallo di tempo

nt=nt + N*ones(1, length(t1));

% qui vado a vedere anche i soggetti censored, diminuisco la quantità sia

% degli esclusi sia dell'evento critico

for j=2:length(t1), nt(j)=nt(j-1)-eventi1(j-1)-c1(j-1); end

% probabilità di eventi nel j-esimo intervallo

f=eventi1./nt;

robabilità di sopravvivenza

s=1-f;

% qualcosa per costruire kaplan meier

preKM = s(1);

KM=1;

% ciclo for per le sopravvivenze nei singoli intervalli

% produttoria delle frequenze calcolate fino a quel tempo

% devo moltiplicare le probabilità da 0 fino a quel tempo

for k=1:7, preKM= preKM*s(k