Modellizzazione matematica

Per semplicità si pone $T_{P}=1$ nell'equazione (2.5). Se si calcola la trasformata di Fourier di ambo i membri, si ottiene:

$$\hat{g_{\Delta}}(\omega_{x},\omega_{y})=4sin^{2}\frac{1}{2}\Delta\omega_{y} \hat{f}(\omega_{x},\omega_{y}) \,\,$$ (2.6)

dove $\omega_{x},\omega_{y}$ sono le frequenze associate alle variabili spaziali $x,y$.
Come già osservò Beckers  [1], l'equazione (2.6) evidenzia che l'immagine chopped and nodded, così  come l'immagine chopped, non contiene informazioni sulla trasformata di Fourier $\hat{f}(\omega_{x},\omega_{y})$ alle frequenze $\omega_{y,k}=2\pi k/ \Delta$ (con $k=0,\pm1,\pm2,...$), anche se $\hat{g_{\Delta}}(\omega_{x},\omega_{y})$ è diversa da zero in tali frequenze a causa del noise che contamina la $g_{\Delta}(x,y)$.
Per questo motivo la ricostruzione di $\hat{f}(\omega_{x},\omega_{y})$ non può essere ottenuta dividendo l'equazione (2.6) per il fattore $sin^{2}(\Delta\omega_{y}/2)$; in effetti calcolare l'immagine $f(x,y)$, partendo dall'acquisizione $g_{\Delta}(x,y)$ risulta un problema ``mal posto''.
Per ovviare a questo inconveniente si possono utilizzare metodi di regolarizzazione  [2]. In generale questi non possono far ricorso a metodi di Fourier, poiché $f$ e $g$ sono definite su domini diversi. Infatti, se la regione da prendere in considerazione corrisponde all'intervallo $[0,Y]$ nella variabile $y$, $g_{\Delta}(x,y)$ è ivi definita (per $x$ fissato) ma, come si può dedurre dall'equazione (2.5), essa riceve contributi dai valori di $f(x,y)$ nell'intervallo $[-\Delta,Y+\Delta]$, che è molto più ampio di $[0,Y]$. Un metodo che ricostruisca $f(x,y)$ in questo intervallo porterà un aumento del range. Di questo fatto si deve tenere conto quando si discretizza il problema.
Assumiamo che il piano del detector sia suddiviso in $N \times N$ pixels ciascuno con dimensione $\delta \times \delta$, etichettati con indice $j$, corrispondente alla colonna (variabile $x$) e con indice $m$ corrispondente alla riga (variabile $y$) dell'array. Inoltre supponiamo che la direzione del movimento di chopping sia parallela alle colonne e che l'ampiezza di chopping $\Delta$ sia un multiplo della distanza di campionamento $\delta$ ( $\Delta=K\delta$). Spesso consideriamo la maggior parte delle immagini delle $128\times128$ con $K$ tipicamente (ma non necessariamente) fra 30 e 50.
Se $g_{j,m}$ e $f_{j,m}$ sono rispettivamente i campioni di $g(x,y)$ e $f(x,y)$, per ogni indice $j$, i valori $g_{j,m}$ (con $m=1,...,N$) formano un vettore di lunghezza $N$ che denoteremo con $\mathbf{g}_{j}$. Esso riceve contributi dagli $N+2K$ valori di $f_{j,m}$ ($m=1,..,N+2K$) il quale forma un vettore $\mathbf{f}_{j}$ di lunghezza $N+2K$. Le componenti di $\mathbf{f}_{j}$ con $m=K+1,...K+N$ corrispondono ai punti di campionamento nella regione di interesse, che chiameremo regione d'osservazione.
Con le notazioni sopra citate, l'equazione (2.5) viene modificata, in termini discreti, nella relazione:

$$g_{j,m}=-f_{j,m}+2f_{j,m+K}-f_{j,m+2K} . \,\,$$ (2.7)

Introduciamo ora la matrice di imaging $[A]$ definita come:

$$[A]_{m,n}=-\delta_{m,n}+2\delta_{m+K,n}-\delta_{m+2K,n} \,\,$$ (2.8)

dove $m=1,2,...,N$ e $n=1,2,....N+2K$ ed inoltre vale:

$$\delta_{m,n}=\left\{ \begin{array}{l} 0 \,\, se \,\, m = n \\ 1 \,\, se \,\, m \neq n \end{array}\right. \,\,$$ (2.9)

In notazione matriciale l'equazione (2.7) diventa:

$$\mathbf{g}_{j}=[A]\mathbf{f}_{j}. \,\,$$ (2.10)

$[A]$ è una matrice rettangolare con $N$ righe e $N+2K$ colonne. Ne segue che il sistema lineare dell'equazione (2.10) è sottodeterminato e, per ogni $\mathbf{g}_{j}$, ci sono almeno $2K$ soluzioni linearmente indipendenti.
Dato che la matrice di imaging non dipende dall'indice $j$, il problema di ricostruzione di immagini si riduce al problema di risolvere l'equazione (2.10) per ogni colonna dell'immagine.
Poiché la matrice $[A]$ risulta mal condizionata, la soluzione è numericamente instabile e di ciò occorrerà tenere conto. Utilizzando il modello descritto dall'equazione (2.10), però, si commettono diversi tipi di errori:

• Il primo errore è dovuto al noise che consiste soprattutto in read-out noise e background Poisson noise (cioè rumore gaussiano bianco). Dal momento che nel caso di grande background il secondo prevale, l'immagine reale può essere rappresentata nel seguente modo:

  
  

  $$\mathbf{g}_{j}=[A]\mathbf{f}_{j}+\mathbf{w}_{j} \,\,$$ (2.11)

  
  
  
  
  dove $\mathbf{w}_{j}$ è una realizzazione di un processo gaussiano. Quando si approssima l'equazione (2.11) con la (2.10) si commette un errore che viene riprodotto e può essere amplificato dal metodo di ricostruzione a causa del mal condizionamento della matrice $[A]$. Per eliminare questo problema si fa uso di metodi che regolarizzano la soluzione, cioè che permettono di ottenere soluzioni approssimate e numericamente stabili.

• Il secondo errore può essere causato dal fatto che l'ampiezza di chopping $\Delta$ non è un multiplo intero della dimensione $\delta$ dei pixels.
  In questo caso l'espressione semplice che usiamo per la matrice $[A]$ dovrà essere sostituita da una più complicata.
  Da simulazioni numeriche si è notato che, se il valore $\Delta/\delta$ è piuttosto grande, come accade in pratica, questo problema perde importanza.
  Si può ottenere lo stesso risultato utilizzando, come valore di $K$, l'intero più vicino a $\Delta/\delta$.

• Il terzo tipo di errore che si può commettere è causato dal non allineamento del detector array con la direzione di chopping e nodding. Si tratta dell'errore più serio tra quelli esposti poiché vìola un'assunzione fondamentale, cioè che la ricostruzione dell'immagine possa essere effettuata colonna per colonna. Se si verifica tale situazione è necessario correggerla mediante una rotazione dell'immagine acquisita.