Print

---

Hi Jürgen,

to check several T1-sequences on our scanner I have written a small tool to estimate SNR 
and CNR for _ spatially normalized_ images. I use the tissue priors (probability > 85%) to create a mask of GM/WM/CSF/Background, which is maybe not perfect, but works very 
well to roughly compare a bunch of images in terms of SNR/CNR. Of course the best 
would be to use individual segmentations for masking...

Regards,

Christian


____________________________________________________________________________

Christian Gaser, Ph.D.
Assistant Professor of Computational Neuroscience
Department of Psychiatry
Friedrich-Schiller-University of Jena
Jahnstrasse 3, D-07743 Jena, Germany
Tel: ++49-3641-934752	Fax:   ++49-3641-934755
e-mail: [log in to unmask]
http://dbm.neuro.uni-jena.de

-----------------------------------------------------------------------------------------------------
function CG_SNR_CNR
%CG_SNR_CNR	to check SNR and CNR in T1 images
%
% Images have to be in the same orientation with same voxel size
% and dimension (e.g. normalized images)

TH = 0.85;

if strcmp(spm('ver'),'SPM5') | strcmp(spm('ver'),'SPM8')
	P = spm_select(Inf,'IMAGE','Select T1-images');
	priors = str2mat(...
		fullfile(spm('Dir'),'tpm','grey.nii'),...
		fullfile(spm('Dir'),'tpm','white.nii'),...
		fullfile(spm('Dir'),'tpm','csf.nii'));
else
	P = spm_get(Inf,'IMAGE','Select T1-images');
	priors = str2mat(...
		fullfile(spm('Dir'),'apriori','gray.mnc'),...
		fullfile(spm('Dir'),'apriori','white.mnc'),...
		fullfile(spm('Dir'),'apriori','csf.mnc'));
end
Vp = spm_vol(priors);

V = spm_vol(P);
n = size(P,1);

hold = 1;

if length(V)>1 & any(any(diff(cat(1,V.dim),1,1),1))
	error('images don''t all have same dimensions'), end
if any(any(any(diff(cat(3,V.mat),1,3),3)))
	error('images don''t all have same orientation & voxel size'), end

M = V(1).mat;

gm  = zeros(V(1).dim(1:3));
wm  = zeros(V(1).dim(1:3));
csf = zeros(V(1).dim(1:3));
bkg = zeros(V(1).dim(1:3));

%-Start progress plot
%-----------------------------------------------------------------------
for z=1:V(1).dim(3)
	Mi  = spm_matrix([0 0 z]);

	% Load slice z from all apriori images
	M1  = M\Vp(1).mat\Mi;
	gm(:,:,z)  = spm_slice_vol(Vp(1),M1,V(1).dim(1:2),[1 0]);
	wm(:,:,z) = spm_slice_vol(Vp(2),M1,V(1).dim(1:2),[1 0]);
	csf(:,:,z)   = spm_slice_vol(Vp(3),M1,V(1).dim(1:2),[1 0]);
end

bkg = 1 - (gm + wm + csf);

% threshold prior maps
gm  = gm > TH;
wm  = wm > TH;
csf = csf > TH;
bkg = bkg > TH;

snr_cnr = zeros(n,4);
spm_progress_bar('Init',n,'Volumes completed');
% load all images.
for i=1:n
    vol = spm_read_vols(V(i));

    vol_gm  = vol(find(gm));
    avg_gm  = mean(vol_gm);
    vol_wm  = vol(find(wm));
    avg_wm  = mean(vol_wm);
    vol_csf = vol(find(csf));
    avg_csf = mean(vol_csf);

    % check order of means
    [tmp, ind] = sort([avg_csf avg_gm avg_wm]);
    if sum(abs(ind - (1:3)))
        fprintf(['Warning: Intensity values are unusual for T1-images. Please check: ' 
V(i).fname '\n']);
    end

	bkg_lt_csf  = vol(find(bkg & (vol < avg_csf)));
	avg_noise = mean(bkg_lt_csf);
	std_noise = std(bkg_lt_csf) + eps;

    % check for zero means
    if ~isempty(find([avg_gm avg_wm avg_csf avg_noise] == 0))
        fprintf(['Warning: Mean intensity is zero. Please check: ' V(i).fname '\n']);
    end

    snr_cnr(i,:) = [avg_gm/std_noise,avg_wm/std_noise,abs(avg_wm - 
avg_gm)/std_noise,abs(avg_gm - avg_csf)/std_noise];

    spm_progress_bar('Set',i);
end

spm_progress_bar('Clear');

sum_snr_cnr = zeros(n,1);
for i=1:n
    sum_snr_cnr(i) = sum(snr_cnr(i,:)./mean(snr_cnr));
end
[new_sum_snr_cnr, ind] = sort(sum_snr_cnr);

fprintf('SNR GM\tSNR WM\tCNR GM/WM\tCNR GM/CSF\tFilename\n');

for i=1:n
    j = ind(i);
    [pth, name] = fileparts(V(j).fname);
    
fprintf('%6.2f\t%6.2f\t%10.2f\t%10.2f\t%s\n',snr_cnr(j,1),snr_cnr(j,2),snr_cnr(j,3),snr_c
nr(j,4),name);
end

if n > 2
	show = spm_input('Show files (beginning from bad)?','+1','yes|no',[1 0],2);
	if show
		number = min([length(ind) 15]);
		number = spm_input('How many files ?','+1','e',number);

		list = str2mat(V(fliplr(ind)).fname);
		list2 = list(1:number,:);
		spm_check_registration(list2)
	end
end

return

On Thu, 8 Jan 2009 14:51:24 +0100, Jürgen Hänggi <[log in to unmask]> 
wrote:

>Dear SPM experts
>
>I would like to compute the SNR of two T1-w. MRI scans of the same subject
>derived from two different scanners. In the past, I used three VOIs. One VOI
>was placed on the head of the caudate nucleus (grey matter mean signal), one
>on the splenium of the corpus callosum (white matter mean signal), and the
>last VOI was placed outside the brain (standard deviation of noise). This
>seems to be the standard procedure to compute SNR.
>
>But, I'am interested in the cortical sheed (neocortex) only. Therefore my
>questions:
>
>1. Do the measures of SNR (CNR) in subcortical structures (caudate nucleus
>and corpus callosum) is also representative for the cortical sheed?
>
>2. There is a function in SPM (img_snr.m) that computes the SNR voxel-wise
>across the whole image. With this function, the standard deviation used to
>divide the mean signal is computed across the whole image (brain and
>background) and not only across the background.
>
>Does this kind of computing the SNR is also meaningful?
>
>3. Are there better procedures to compute the SNR across the whole brain
>such as using a mask for the brain to get the mean signal of tissue and one
>mask for the background to get the SD of the noise?
>
>Thanks in advance
>Kindly regards
>Jürgen 
>
>
>
>----------------------------------------------------------
>Juergen Haenggi
>Ph.D. (Dr. des.)
>Division Neuropsychology
>Institute of Psychology
>University of Zurich
>Binzmuehlestrasse 14, PO Box 25
>8050 Zurich, Switzerland
>0041 44 635 73 97 (phone office)
>0041 76 445 86 84 (phone mobile)
>0041 44 635 74 09 (fax office)
>BIN 4.D.04 (office room number)
>j.haenggi[at]psychologie.uzh.ch (email)
>http://www.psychologie.uzh.ch/neuropsy/ (website)
>http://www.juergenhaenggi.ch (private website)
>----------------------------------------------------------