Gli algoritmi di elaborazione del segnale definiti in tempo / spazio / frequenza continui vengono in genere implementati campionando il segnale su una griglia discreta e convertendo gli integrali in somme (e derivati in differenze). I filtri spaziali sono implementati attraverso la convoluzione con un kernel di convoluzione (cioè la somma ponderata dei vicini).
Esiste un enorme bagaglio di conoscenze sul filtraggio dei segnali del dominio del tempo campionati; i filtri nel dominio del tempo sono implementati come filtri a risposta impulsiva finita , dove il campione di output corrente viene calcolato come una somma ponderata del precedenti campioni di input N; o filtri a risposta impulsiva infinita, in cui l'uscita corrente è una somma ponderata degli ingressi precedenti e delle precedenti uscite . Formalmente, i filtri del tempo discreto sono descritti usando z-transform , che è l'analogo del tempo discreto per Trasformata di Laplace . La trasformazione bilineare mappa l'una con l'altra ( c2d e d2c in Matlab).
How would you evaluate functions at arbitrary points?
Quando hai bisogno del valore di un segnale in un punto che non giace direttamente sulla tua griglia di campionamento, interpola il suo valore da punti vicini. L'interpolazione può essere semplice come scegliere il campione più vicino, calcolare una media ponderata dei campioni più vicini, o inserire una funzione analitica arbitrariamente complicata nei dati campionati e valutare questa funzione alle coordinate necessarie. L'interpolazione su una griglia più uniforme è upsampling . Se il segnale originale (continuo) non contiene dettagli (cioè frequenze) più fini della metà della griglia di campionamento, allora la funzione continua può essere perfettamente ricostruita dalla versione campionata (il Teorema di campionamento di Nyquist-Shannon ). Per un esempio di come puoi interpolare in 2D, vedi interpolazione bilineare .
In Matlab puoi usare interp1 o interp2 per interpolare 1D o dati 2D campionati regolarmente (rispettivamente) o griddata per interpolare da dati 2D campionati irregolarmente.
Would you have a for-loop going through each voxel and computing the corresponding formula?
Sì, esattamente.
Matlab ti evita di doverlo fare tramite espliciti cicli di apertura perché è progettato per operare su matrici e vettori (cioè array multidimensionali). In Matlab questo è chiamato "vettorizzazione". Gli integrali definiti possono essere approssimati con sum , cumsum , trapz , cumtrapz , ecc.
I've read the book "Digital Image Processing" by Gonzalez and Woods but I'm still at a loss. I have also read about the Numerical Recipes book series. Would that be the correct way?
Sì, Numerical Recipes sarebbe un ottimo inizio. È molto pratico e copre la maggior parte dei metodi numerici di cui avrai bisogno. (Scoprirai che Matlab implementa già tutto ciò che ti serve, ma Ricette numeriche fornirà uno sfondo eccellente.)
I have taken an "algorithms and data structures" class, but I don't see the relation between the material presented there and implementing scientific algorithms.
Il materiale trattato nei corsi "Algoritmi e strutture dati" tende a concentrarsi su strutture come elenchi, matrici, alberi e grafici contenenti numeri interi o stringhe e operazioni come l'ordinamento e la selezione: problemi per i quali esiste in genere un singolo risultato corretto. Quando si tratta di algoritmi scientifici, questa è solo metà della storia. L'altra metà riguarda i metodi per stimare i numeri reali e le funzioni analitiche. Lo troverai in un corso su "Metodi numerici" (o "Analisi numerica"; ad esempio questo - scorri verso il basso per le diapositive): come stimare le funzioni speciali, come stimare integrali e derivati, ecc. Qui uno dei compiti principali è stimare l'accuratezza del risultato, e uno schema comune è quello di iterare una routine che migliora una stima fino a quando non è sufficientemente accurato (Potresti chiederti come Matlab sa come fare qualcosa di semplice come stimare un valore di sin(x) per qualche x .)
Come semplice esempio, ecco uno script breve che calcola una trasformazione di radon di un'immagine in Matlab. La trasformazione del radon richiede proiezioni di un'immagine su una serie di angoli di proiezione. Invece di cercare di calcolare la proiezione lungo un angolo arbitrario, ruoterò invece l'intera immagine usando imrotate , in modo che la proiezione sia sempre verticale. Quindi possiamo prendere la proiezione semplicemente usando sum , poiché sum di una matrice restituisce un vettore contenente la somma su ciascuna colonna.
Potresti scrivere il tuo imrotate se preferisci, usando interp2 .
%%# Home-made Radon Tranform
%# load a density map (image).
A = phantom;
n_pixels = size(A, 1); %# image width (assume square)
%# At what rotation angles do we want to take projections?
n_thetas = 101;
thetas = linspace(0, 180, n_thetas);
result = zeros(n_thetas, n_pixels);
%# Loop over angles
for ii=1:length(thetas)
theta = thetas(ii);
rotated_image = imrotate(A, theta, 'crop');
result(ii, :) = sum(rotated_image);
end
%# display the result
imagesc(thetas, 1:n_pixels, result.');
xlabel('projection angle [degrees]');
Quello che una volta era un integrale della densità lungo un raggio è ora una somma su una colonna di un'immagine discretamente campionata, che a sua volta è stata trovata interpolando l'immagine originale su un sistema di coordinate trasformato.