X-Ray Guided: a coronary stenting simulator on one fluoro pedal

Raffaele
/
X-Ray Guided: a coronary stenting simulator on one fluoro pedal
about 19 hours ago
3
28
9
2
9984
game.m
function g = game()
%game  X-RAY GUIDED - a coronary stenting simulator on one fluoro pedal
%
%   Stati (s.mode):
%     0 titolo   1 fermo (cieca)   2 spinta alla cieca   3 scopia in navigazione
%     4 referto  5 complicanza / fine turno
%     6 angiografia / roadmap   7 posizionamento stent   8 gonfiaggio
%     9 scheda paziente
%
%   Turno in sala: pazienti in sequenza, ognuno con 2 caratteristiche casuali
%   (s.f), velocita crescente (s.v), fine turno alla terza complicanza.
%
%   NB: il limite di caratteri e' contato su tree2str(mtree(...)), che aggiunge
%   spazi e indentazione: i commenti sono gratis, annidamenti e righe no.
%   Conteggio locale: python3 _dev/t2s.py game.m
g.init = @initState;
g.update = @updateState;
g.render = @renderFrame;
end

function s = initState(ax)
% stile da monitor di angiografo: fondo nero, testi grigio chiaro, accenti gialli
G = [.8 .8 .8];
D = [.3 .3 .3];
% helper di creazione (tutti gli oggetti nascono qui, una volta sola)
I = @(x,y,c) image(ax,x,y,c,'CDataMapping','scaled');
T = @(x,y,z,t) text(ax,x,y,t,'Color',G,'FontSize',z);
R = @(p,c) rectangle(ax,'Position',p,'EdgeColor',c);
% geometria immagine RX 128x128 (origine al centro, Y verso l'alto)
[X,Y] = meshgrid(-63.5:63.5);
s.X = X; s.Y = Y;
s.M = 1-(X.^2+Y.^2)/3e4;                % flat panel quadrato, leggera vignettatura
s.B = .3*exp(-((X+34)/13).^2);          % colonna vertebrale
% monitor, simbolo raggi a pixel, testi in sovrimpressione
h.im = I([32 508],[92 568],0);
h.ref = I([36 176],[424 564],0);          % immagine di riferimento (secondo monitor)
R([28 88 484 484],D);
[x,y] = meshgrid(-7:7); r = hypot(x,y);
t = r<1.6 | r>2.6 & r<7.2 & mod(atan2(y,x),2*pi/3)<pi/3;
h.rad = I([454 500],[514 560],cat(3,t,.85*t,0*t));
h.lih = T(190,550,14,'');
h.par = T(40,112,10,'');
set([h.lih h.par],'BackgroundColor','k','FontWeight','bold');   % leggibili sopra l'immagine chiara
% schema del percorso
line(ax,[548 548],[130 560],'Color',D,'LineWidth',4);
h.pos = line(ax,548,130,'Marker','s','MarkerSize',11,'Color','w');
n = {'FEMORAL','ILIAC','AORTA','ARCH','OSTIUM'};
for k = 1:5
    T(562,22.5+107.5*k,9,n{k});
end
% barre: resistenza / pressione palloncino, dose, danno
n = {'FEEL','DOSE','HARM'};
for k = 1:3
    x = 615+50*k;
    R([x 130 30 430],D);
    h.bar(k) = R([x 130 30 1],'none');
    set(T(x+15,112,9,n{k}),'HorizontalAlignment','center');
end
h.tg = R([660 385 40 27],'y');          % finestra 95-105 % della nominale
h.ped = T(30,50,13,'');
h.ped.Color = [1 .8 .2];                % riga del pedale in giallo
h.msg = T(270,330,16,'');
set(h.msg,'HorizontalAlignment','center','FontWeight','bold');
h.hud = T(540,585,10,'');
% assi impostati dopo aver creato le immagini, cosi image() non li puo cambiare
set(ax,'YDir','normal','Color','k','CLim',[0 1.05]);
colormap(ax,gray(256));
set(ax.Children,'HitTest','off','PickableParts','none');
s.h = h;
s.best = 0; s.life = 0;
s = newCase(s);
s.mode = 0; s.n = 0;
end

function s = newCase(s)
% Nuovo paziente (e nuovo turno se le complicanze sono finite)
if s.life<1
    % statistiche del turno: stent riusciti, DAP e scopia totali
    [s.n,s.tot,s.sk,s.ns,s.td,s.tf] = deal(0);
    s.cl = [];                       % elenco delle complicanze
    s.life = 3;
end
s.n = s.n+1;
% 2 caratteristiche su 6: 1 obeso  2 iliache tortuose  3 tachicardico
%                         4 vaso piccolo  5 insuff. renale  6 calcifico
f = ismember(1:6,randperm(6,2));
s.f = f;
s.v = .88+.12*s.n;                   % fattore di velocita crescente
s.nz = .09+.07*f(1);                 % rumore quantico
s.vw = 6-2*f(4);                     % semi-larghezza della coronaria
% ostacoli casuali + ostio (ultimo, passabile solo in scopia)
n = 4+2*f(2);
s.op = [.08+.8*((1:n)-.5+.6*(rand(1,n)-.5))/n, 1];     % posizione
s.os = [min(.95,.2+.5*rand(1,n)+.25*f(2)), 1];         % forza 0..1
s.oa = (20+22*rand(1,n+1)).*sign(randn(1,n+1));        % direzione del lume (gradi)
[s.p,s.m,s.t,s.ft,s.dap,s.harm,s.feel,s.th,s.P,s.cov,s.sc,s.stall,s.b0,s.p0,s.why,s.cine] = deal(0);
s.k = 1; s.fade = 1; s.q = -1; s.sp = -50;
s.ph = 6*rand;
s.ca = 6+12*rand;                    % ampiezza delle curve della coronaria
s.cw = 16+12*rand;                   % passo delle curve
s.d = zeros(128);
s.ref = s.d;
s.ys = 30*rand-5;                    % centro della stenosi (pixel)
s.ls = 10+8*rand;                    % lunghezza della stenosi
s.pn = 10+round(4*rand);             % pressione nominale (atm)
s.mode = 9;
end

function s = updateState(s, dt, pr)
m = s.mode;
if any(m==[0 4 5 9])
    if pr && m==9
        s.mode = 1;
    elseif pr
        s = newCase(s);
    end
elseif s.dap>=60 || s.harm>=100
    s = fin(s,1+5*(s.dap>=60));         % dissezione, o limite di dose (uguale per tutti)
elseif m>5
    s = cor(s,dt,pr);
else
    s = nav(s,dt,pr);
end
end

function s = nav(s, dt, pr)
% NAVIGAZIONE femorale -> ostio. m: 1 fermo, 2 spinta alla cieca, 3 scopia
m = s.mode; k = s.k; z = s.stall;       % z: fermo su un ostacolo a inizio frame
a = s.os(k); b = s.oa(k);
X = s.X; Y = s.Y;
s.t = s.t+dt;
% avanzamento (cieca .07, in scopia .05) e distanza dal prossimo ostacolo
s.p = min(s.p+(m>1)*~z*(.11-.02*m)*s.v*dt, s.op(k));
q = s.op(k)-s.p;
s.stall = q<=0;
% resistenza "in mano": fondo, falsi picchi, rampa prima dell'ostacolo
s.feel = (m>1)*(.2+.25*max(0,sin(23*s.p+s.ph))^8+.4*max(0,1-q/.03)+z*(.1+.3*a));
s.harm = s.harm+(m==2)*z*18*a*dt;       % forzare alla cieca fa danno
s.fade = max(.25+.75*(m>2), s.fade-.5*dt);             % LIH che si spegne
s.dap = s.dap+(m>2)*(1+.5*s.f(1))*dt;
s.ft = s.ft+(m>2)*dt;
% punta: ruota ampia se ferma sull'ostacolo, piccola oscillazione altrimenti
w = (1.6+.8*a)*min(s.v,1.5)*s.t;        % velocita della punta, con tetto
s.th = (.15+.85*z)*(40*sin(w)+12*sin(2.3*w+s.ph));
s.m = s.m+(m>2)*z*s.v*dt/(3+6*a);       % micromovimenti (piu rapidi coi pazienti)
if m>2
    % frame RX centrato sulla punta. Asse del vaso: dritto nel tratto attuale,
    % piegato di b sopra il prossimo ostacolo e di b0 sotto l'ultimo superato;
    % la guida segue l'asse nel tratto gia percorso.
    c = cosd(s.th); e = sind(s.th);
    u = X*e+Y*c;
    xc = X-tand(b)*max(0,Y-1000*q)+tand(s.b0)*min(0,Y+1000*(s.p-s.p0));
    s = shot(s, s.B.*(1+.3*cos((Y+900*s.p)/5)) ...
        + .12*exp(-(xc/8).^2) ...
        + .8*(Y<=0 & abs(xc)<1) ...
        + 1.2*(u>=0 & u<=10+10*s.m & abs(X*c-Y*e)<1.3));
end
if z && (m>2 && (pr && abs(s.th-b)<14 || s.m>=1) || m==2 && ~pr && rand<.5*(1-a)*dt)
    % ostacolo superato: tap allineato, micromovimenti, o forzatura riuscita
    s.k = k+1; s.m = 0; s.stall = 0;
    s.b0 = b; s.p0 = s.p;               % la curva appena presa resta alle spalle
    s.mode = 2+4*(k==numel(s.op));      % dopo l'ostio si passa alla coronaria
elseif pr
    s.harm = s.harm+12*(m>2)*z;         % tap fuori allineamento
    s.mode = 2+(m==2)+(m>2)*z;          % 1->2, 2->3, 3->2 (o resta 3 se sull'ostacolo)
end
end

function s = cor(s, dt, pr)
% CORONARIA. m: 6 angiografia, 7 posizionamento stent, 8 gonfiaggio
m = s.mode; f = s.f; q = s.q;
j = m==6 && q>=0;                       % iniezione in corso (cine)
s.t = s.t+dt;
s.ft = s.ft+dt;
s.fade = 1;
g = s.P/s.pn;
s.feel = g/1.6;
yb = s.Y-(3+2*f(3))*sin((7.5+5*f(3))*s.t);             % battito cardiaco
xc = s.X-s.ca*sin(yb/s.cw+s.ph);                           % distanza dall'asse del vaso
x = @(y) max(0,(q-(y+64)/90)/1.6);                     % bolo che risale il vaso
b = @(x) (x.*exp(1-x)).^2;
k = abs(yb-s.sp);
L = s.ls/2+3;                           % semi-lunghezza dello stent (su misura)
u = 14*s.v*q;                           % corsa dello stent
% cine: durante l'iniezione e durante ogni puff (1.5 s dopo il fine corsa)
s.cine = j || m==7 && u>95 && mod(u,95)<21*s.v;
s.dap = s.dap+(1+4*s.cine)*(1+.5*f(1))*dt;
% contrasto: bolo del cine, oppure puff automatico a ogni fine corsa dello stent
C = j*(1.3-.7*f(5))*b(x(yb))+(m==7 && u>95)*(.6-.3*f(5))*b(mod(u,95)/7/s.v);
wv = s.vw*(1-.63*(1-min(g,1))*exp(-((yb-s.ys)/s.ls*2).^4));   % lume con stenosi
% senza contrasto il vaso non si vede: restano guida, marker e palloncino
s = shot(s, s.B.*(1+.3*cos(s.Y/5)) ...
    + C.*exp(-(xc./wv).^2) ...
    + .8*(abs(xc)<1) ...
    + (m>6)*(1.5*(abs(k-L)<1.5 & abs(xc)<2.5) + .3*(k<L & abs(xc)<1.5+(s.vw-1.5)*g)));
q = q+(m>6 || j)*dt;
s.q = q;
if m==6 && ~j
    s.q = pr-1;                         % tap: parte l'iniezione
elseif m==6 && (pr || q>6)
    s.ref = s.d;                        % il frame del tap diventa il riferimento
    s.mode = 7; s.q = 0;
elseif m==7
    s.sp = 45-abs(mod(u,190)-95);       % stent avanti e indietro
    s.mode = 7+pr; s.q = q*~pr;
elseif m>7
    s.P = s.P+(2+6*g+(1+3*f(6))*sin((5+4*f(6))*q))*(.8+.2*s.v)*dt; % gonfiaggio che accelera
    g = s.P/s.pn;
    if pr || g>=1.6
        s.cov = max(0,min(1,(L+s.ls/2-abs(s.sp-s.ys))/s.ls));
        s.harm = max(s.harm+40*(g>1.2), 100*(g>=1.6)); % sovradistensione / rottura
        s.sc = max(0,round(1000*s.cov*max(0,1-2*abs(g-1))-4*s.dap-5*s.harm-2*s.t));
        % causa: 2 stent fuori lesione, 3 sotto-espanso, 4 sovradisteso, 5 rottura
        s = fin(s, max([2*(s.cov<.8), 3*(g<.8), 4*(g>1.2), 5*(g>=1.6)]));
    end
end
end

function s = fin(s, c)
% Fine paziente: c = 0 riuscito (punteggio con moltiplicatore di serie),
% altrimenti codice della complicanza
ok = ~c;
s.why = c;
s.cl = [s.cl c(~ok)];
s.ns = s.ns+ok;
s.td = s.td+s.dap;
s.tf = s.tf+s.ft;
s.sk = ok*(s.sk+1);
s.tot = s.tot+ok*round(s.sc*(.75+s.sk/4));
s.life = s.life-~ok;
s.mode = 5-ok;
s.best = max(s.best,s.tot);
end

function s = shot(s, A)
% Frame RX dalla mappa di attenuazione A: rumore quantico + filtro ricorsivo
I = exp(-A);
s.d = .5*s.d+.5*s.M.*(I+s.nz*sqrt(I).*randn(128));
end

function renderFrame(~, s)
h = s.h; m = s.mode; z = s.stall;
g = s.P/s.pn;
h.im.CData = s.d*s.fade*(m<4 || m>5);
h.pos.YData = 130+430*s.p;
% barre; la prima e rossa oltre soglia, in gonfiaggio gialla/verde/rossa
v = min([s.feel, s.dap/60, s.harm/100],1);
% grigio normale, giallo sotto pressione, rosso oltre soglia, verde in finestra
c = {[.8 .8 .8],[1 .8 .2],[1 .3 .25],[.3 1 .5]};
c{1} = c{(m>7)*(2+2*(g>=.95)-(g>1.05))+(m<8)*(1+2*(v(1)>.6))};
for k = 1:3
    set(h.bar(k),'Position',[615+50*k 130 30 max(430*v(k),1)],'FaceColor',c{k});
end
h.tg.Visible = m==8;
set(h.ref,'CData',s.ref,'Visible',m==7 || m==8);
% riga del pedale
p = {'start','push wire, blind','stop + fluoro on', ...
     'fluoro off, push blind','when aligned / WAIT: creep', ...
     'next patient','next patient','inject contrast', ...
     'at peak, save reference','deploy stent on lesion', ...
     'stop balloon in yellow box','patient in'};
h.ped.String = ['TAP: ' p{m+1+(m>3 || m==3 && z)+(m>6 || m==6 && s.q>=0)}];
% overlay da angiografo
% intestazione del monitor: nome della fase + stato
e = 'WIRE NAVIGATION';
o = [e newline 'X-RAY OFF - LIH'];
q = {'',o,o,e,'','','ANGIOGRAPHY','STENT POSITIONING',sprintf('BALLOON INFLATION\n%.1f / %d atm',s.P,s.pn),''};
h.lih.String = q{m+1};
x = m==3 || m>5 && m<9;                 % raggi accesi
h.rad.Visible = x;
j = s.cine;
kv = 66+8*s.f(1)+round(9*s.p+2*sin(s.t))+6*j;
a = sprintf('\n%d kV  %.0f mA  %g fps',kv,kv/9*(1+60*j),7.5+7.5*j);
h.par.String = sprintf('DAP %.1f Gy.cm2  FT %d s%s',s.dap,round(s.ft),a(1:numel(a)*x));
% messaggi a centro monitor
T = {'obese: noise + dose','tortuous: more resistance','tachycardic: fast heart', ...
     'small vessel: faint','renal: weak contrast','calcified: jumpy balloon'};
y = {'VESSEL DISSECTION','LESION MISSED','UNDER-EXPANSION', ...
     'OVERSTRETCH','CORONARY RUPTURE','SKIN DOSE INJURY'};
% complicanza durante il turno, oppure sunto di fine turno
w = {sprintf('COMPLICATION\n%s\nstrike %d of 3\n\nSHIFT %d',y{max(s.why,1)},3-s.life,s.tot), ...
     sprintf(['SHIFT OVER\n%d patients, %d stents placed\n\ncomplications:\n%s\n' ...
     'total DAP %.0f Gy.cm2  fluoro %.1f min\n\nSHIFT %d  BEST %d'], ...
     s.n,s.ns,sprintf('%s\n',y{s.cl}),s.td,s.tf/60,s.tot,s.best)};
t = {sprintf('X-RAY GUIDED\na coronary stenting simulator\n\nThe tap is your X-ray pedal.\nThe line at the bottom tells you\nwhat it does at each step.\nWaiting is safer but costs dose.\nFull HARM or DOSE = complication.\n3 complications end your shift.\n\nSPACE or click = pedal'), ...
     '',repmat('RESISTANCE',1,1*z),'', ...
     sprintf('STENT DEPLOYED\n\ncoverage %.0f %%  expansion %.0f %%\nDAP %.0f Gy.cm2  trauma %.0f %%\n\nSCORE +%d x%.2f', ...
         100*s.cov,100*g,s.dap,s.harm,s.sc,.75+s.sk/4), ...
     w{1+(s.life<1)}, ...
     '','','',sprintf('PATIENT %d\n\n%s\n%s',s.n,T{s.f})};
h.msg.String = t{m+1};
h.hud.String = sprintf('PT %d  SCORE %d  STRIKES %d/3',s.n,s.tot,3-s.life);
end
Play Game
Preview

Play Game opens this in MATLAB Drive. Click Run this app there to play.