Hi, Christian:
I'm wondering if you have any tool for estimating SNR of DTI data?
Best
--
Jiansong Xu, M.D., Ph.D.
Associate Research Scientist
Dept. of Psychiatry
Yale School of Medicine
Doctors Building, Suite 215
2 Church Street
New Haven, CT 06519
On 1/9/09 2:39 AM, "Christian Gaser" <[log in to unmask]> wrote:
> 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)
>> ----------------------------------------------------------
|