There are a few differences. The main ones I see in the code are that:
* the smoothing of the histogram no longer assumes that the histogram consists
of a load of stick functions, and instead assumes that it is continuously
defined function.
* hot spots in the images are ignored when computing the image maximum (for
working out the binning). Here in the FIL, there is an artifact at the edge
of some of our images, whereby a few voxels contain extremely high values.
Now only the smallest 99.99% of voxels are used for figuring out the binning.
* The Powell's method optimisation should also be a little more stable in
SPM5.
Best regards,
-John
john@ash:/local/spm> diff -db spm2/spm_coreg.m spm5/spm_coreg.m
75c75,79
< % @(#)spm_coreg.m 2.4 John Ashburner 03/03/12
---
> % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
>
> % John Ashburner
> % $Id: spm_coreg.m 1007 2007-11-21 12:37:33Z john $
>
86c90
< else,
---
> else
90c94
< if ~isfield(flags,fnms{i}), flags =
setfield(flags,fnms{i},getfield(def_flags,fnms{i})); end;
---
> if ~isfield(flags,fnms{i}), flags.(fnms{i}) =
def_flags.(fnms{i}); end;
96,97c100,101
< VG = spm_vol(spm_get(1,'IMAGE','Select reference image'));
< else,
---
> VG = spm_vol(spm_select(1,'image','Select reference image'));
> else
102,103c106,107
< VF = spm_vol(spm_get(1,'IMAGE','Select moved image'));
< else,
---
> VF = spm_vol(spm_select(Inf,'image','Select moved image(s)'));
> else
105c109
< if ischar(VF), VF = spm_vol(VF); end;
---
> if ischar(VF) || iscellstr(VF), VF = spm_vol(strvcat(VF)); end;
108,116d111
<
< % Voxel sizes (mm)
< vxg = sqrt(sum(VG.mat(1:3,1:3).^2));
< vxf = sqrt(sum(VF.mat(1:3,1:3).^2));
<
< % Smoothnesses
< fwhmg = sqrt(max([1 1 1]*flags.sep(end)^2 - vxg.^2, [0 0 0]))./vxg;
< fwhmf = sqrt(max([1 1 1]*flags.sep(end)^2 - vxf.^2, [0 0 0]))./vxf;
<
118a114,115
> vxg = sqrt(sum(VG.mat(1:3,1:3).^2));
> fwhmg = sqrt(max([1 1 1]*flags.sep(end)^2 - vxg.^2, [0 0 0]))./vxg;
121,125d117
< if ~isfield(VF, 'uint8'),
< VF.uint8 = loaduint8(VF);
< VF = smooth_uint8(VF,fwhmf); % Note side effects
< end;
<
130,136c122,139
< x = flags.params(:);
< for samp=flags.sep(:)',
< [x,fval] = spm_powell(x(:),
xi,sc,mfilename,VG,VF,samp,flags.cost_fun,flags.fwhm);
< x = x(:)';
< end;
< if flags.graphics,
< display_results(VG,VF,x,flags);
---
>
> for k=1:numel(VF),
> VFk = VF(k);
> if ~isfield(VFk, 'uint8'),
> VFk.uint8 = loaduint8(VFk);
> vxf = sqrt(sum(VFk.mat(1:3,1:3).^2));
> fwhmf = sqrt(max([1 1 1]*flags.sep(end)^2 - vxf.^2, [0 0
0]))./vxf;
> VFk = smooth_uint8(VFk,fwhmf); % Note side effects
> end;
>
> xk = flags.params(:);
> for samp=flags.sep(:)',
> xk = spm_powell(xk(:),
xi,sc,mfilename,VG,VFk,samp,flags.cost_fun,flags.fwhm);
> x(k,:) = xk(:)';
> end;
> if flags.graphics,
> display_results(VG(1),VFk(1),xk(:)',flags);
> end;
155,160c158,160
< sm = max(fwhm/sqrt(8*log(2)),0.001); % FWHM -> Gaussian param
< t = max(round(3*sm(1)),0); krn1 = exp(-([-t:t].^2)/sm(1)^2) ; krn1 =
krn1/sum(krn1) ; H = conv2(H,krn1);
< t = max(round(3*sm(2)),0); krn2 = exp(-([-t:t].^2)/sm(2)^2)'; krn2 =
krn2/sum(krn2)'; H = conv2(H,krn2);
< %d = 32;
< %H = sum(reshape(H,[256/d d 256]),1);
< %H = reshape(sum(reshape(H,[d 256/d d]),2),[d d]);
---
> lim = ceil(2*fwhm);
> krn1 = smoothing_kernel(fwhm(1),-lim(1):lim(1)) ; krn1 = krn1/sum(krn1); H =
conv2(H,krn1);
> krn2 = smoothing_kernel(fwhm(2),-lim(2):lim(2))'; krn2 = krn2/sum(krn2); H =
conv2(H,krn2);
212c212
< if size(V.pinfo,2)==1 & V.pinfo(1) == 2,
---
> if size(V.pinfo,2)==1 && V.pinfo(1) == 2,
215c215
< else,
---
> else
226a227,247
>
> % Another pass to find a maximum that allows a few hot-spots in the data.
> spm_progress_bar('Init',V.dim(3),...
> ['2nd pass max/min of ' spm_str_manip(V.fname,'t')],...
> 'Planes complete');
> nh = 2048;
> h = zeros(nh,1);
> for p=1:V.dim(3),
> img = spm_slice_vol(V,spm_matrix([0 0 p]),V.dim(1:2),1);
> img = img(isfinite(img));
> img = round((img+((mx-mn)/(nh-1)-mn))*((nh-1)/(mx-mn)));
> if spm_matlab_version_chk('7.0')>=0,
> h = h + accumarray(img,1,[nh 1]);
> else
> h = h + full(sparse(img,1,1,nh,1));
> end
> spm_progress_bar('Set',p);
> end;
> tmp = [find(cumsum(h)/sum(h)>0.9999); nh];
> mx = (mn*nh-mx+tmp(1)*(mx-mn))/(nh-1);
>
230a252
> %udat = zeros(V.dim,'uint8'); Needs MATLAB 7 onwards
232c254,255
< udat(V.dim(1),V.dim(2),V.dim(3))=0;
---
> udat(V.dim(1),V.dim(2),V.dim(3)) = 0;
>
238,239c261,262
< udat(:,:,p) = uint8(round((img-mn)*(255/(mx-mn))));
< else,
---
> udat(:,:,p) =
uint8(max(min(round((img-mn)*(255/(mx-mn))),255),0));
> else
242c265
< udat(:,:,p) = uint8(round((img+r-mn)*(255/(mx-mn))));
---
> udat(:,:,p) =
uint8(max(min(round((img+r-mn)*(255/(mx-mn))),255),0));
250c273
< if ~spm_type(V.dim(4),'intt'),
---
> if ~spm_type(V.dt(1),'intt'),
252c275
< else,
---
> else
255c278
< else,
---
> else
264,274c287,290
< s = fwhm/sqrt(8*log(2));
< x = round(6*s(1)); x = [-x:x];
< y = round(6*s(2)); y = [-y:y];
< z = round(6*s(3)); z = [-z:z];
< x = exp(-x.^2/(2*s(1).^2+eps));
< y = exp(-y.^2/(2*s(2).^2+eps));
< z = exp(-z.^2/(2*s(3).^2+eps));
< x = x/sum(x);
< y = y/sum(y);
< z = z/sum(z);
<
---
> lim = ceil(2*fwhm);
> x = -lim(1):lim(1); x = smoothing_kernel(fwhm(1),x); x = x/sum(x);
> y = -lim(2):lim(2); y = smoothing_kernel(fwhm(2),y); y = y/sum(y);
> z = -lim(3):lim(3); z = smoothing_kernel(fwhm(3),z); z = z/sum(z);
282a299,332
> function krn = smoothing_kernel(fwhm,x)
>
> % Variance from FWHM
> s = (fwhm/sqrt(8*log(2)))^2+eps;
>
> % The simple way to do it. Not good for small FWHM
> % krn = (1/sqrt(2*pi*s))*exp(-(x.^2)/(2*s));
>
> % For smoothing images, one should really convolve a Gaussian
> % with a sinc function. For smoothing histograms, the
> % kernel should be a Gaussian convolved with the histogram
> % basis function used. This function returns a Gaussian
> % convolved with a triangular (1st degree B-spline) basis
> % function.
>
> % Gaussian convolved with 0th degree B-spline
> % int(exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= -0.5..0.5)
> % w1 = 1/sqrt(2*s);
> % krn = 0.5*(erf(w1*(x+0.5))-erf(w1*(x-0.5)));
>
> % Gaussian convolved with 1st degree B-spline
> % int((1-t)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= 0..1)
> % +int((t+1)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t=-1..0)
> w1 = 0.5*sqrt(2/s);
> w2 = -0.5/s;
> w3 = sqrt(s/2/pi);
> krn = 0.5*(erf(w1*(x+1)).*(x+1) + erf(w1*(x-1)).*(x-1) - 2*erf(w1*x ).*
x)...
> +w3*(exp(w2*(x+1).^2) + exp(w2*(x-1).^2) - 2*exp(w2*x.^2));
>
> krn(krn<0) = 0;
> return;
> %_______________________________________________________________________
>
> %_______________________________________________________________________
289c339
< txt = 'Information Theoretic Coregistration';
---
> %txt = 'Information Theoretic Coregistration';
336,337c386,387
< h1 = spm_orthviews('Image',VG.fname,[0.01 0.01 .48 .49]);
< h2 = spm_orthviews('Image',VF.fname,[.51 0.01 .48 .49]);
---
> spm_orthviews('Image',VG,[0.01 0.01 .48 .49]);
> h2 = spm_orthviews('Image',VF,[.51 0.01 .48 .49]);
340c390
< spm_orthviews('Space',h1);
---
> spm_orthviews('Space');
john@ash:/local/spm> diff -db spm2/spm_powell.m spm5/spm_powell.m
17,18c17,18
< % Method is based on Powell's optimisation method from Numerical Recipes
< % in C (Press, Flannery, Teukolsky & Vetterling).
---
> % Method is based on Powell's optimisation method described in
> % Numerical Recipes (Press, Flannery, Teukolsky & Vetterling).
20c20,24
< % @(#)spm_powell.m 2.3 John Ashburner 01/09/28
---
> % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
>
> % John Ashburner
> % $Id: spm_powell.m 691 2006-11-22 17:44:19Z john $
>
24,25c28,30
< for iter=1:128,
< fprintf('iteration %d...\n', iter);
---
> for iter=1:512,
> if numel(p)>1, fprintf('iteration %d...\n', iter); end;
> ibig = numel(p);
31c36
< [p,junk,f] = linmin(p,xi(:,i),func,f,tolsc,varargin{:});
---
> [p,junk,f] = min1d(p,xi(:,i),func,f,tolsc,varargin{:});
37c42
< if sqrt(sum(((p(:)-pp(:))./tolsc(:)).^2))<1, return; end;
---
> if numel(p)==1 || sqrt(sum(((p(:)-pp(:))./tolsc(:)).^2))<1, return;
end;
40c45
< [p,xi(:,ibig),f] = linmin(p,p-pp,func,f,tolsc,varargin{:});
---
> [p,xi(:,ibig),f] = min1d(p,p-pp,func,f,tolsc,varargin{:});
48c53
< function [p,pi,f] = linmin(p,pi,func,f,tolsc,varargin)
---
> function [p,pi,f] = min1d(p,pi,func,f,tolsc,varargin)
51c56
< global lnm % used in linmineval
---
> global lnm % used in funeval
55,56c60,61
< linmin_plot('Init', 'Line Minimisation','Function','Parameter Value');
< linmin_plot('Set', 0, f);
---
> min1d_plot('Init', 'Line Minimisation','Function','Parameter Value');
> min1d_plot('Set', 0, f);
60c65
< [f,pmin] = brents(t,tol);
---
> [f,pmin] = search(t,tol);
66c71
< linmin_plot('Clear');
---
> min1d_plot('Clear');
72c77
< function f = linmineval(p)
---
> function f = funeval(p)
75c80
< global lnm % defined in linmin
---
> global lnm % defined in min1d
78c83
< linmin_plot('Set',p,f);
---
> min1d_plot('Set',p,f);
84c89
< % Bracket the minimum (t(2)) between t(1) and t(2)
---
> % Bracket the minimum (t(2)) between t(1) and t(3)
90c95
< t(2).f = linmineval(t(2).p);
---
> t(2).f = funeval(t(2).p);
92c97
< % if not better then swap
---
> % if t(2) not better than t(1) then swap
94c99
< tmp = t(1);
---
> t(3) = t(1);
96c101
< t(2) = tmp;
---
> t(2) = t(3);
100c105
< t(3).f = linmineval(t(3).p);
---
> t(3).f = funeval(t(3).p);
109a115,116
> if pol(3)>0,
> % minimum is when gradient of polynomial is zero
110a118,122
>
> % A very conservative constraint on the displacement
> if d > (1+gold)*(t(3).p-t(2).p),
> d = (1+gold)*(t(3).p-t(2).p);
> end;
111a124,133
> else,
> % sign of pol(3) (the 2nd deriv) is not +ve
> % so extend out by golden ratio instead
> u.p = t(3).p+gold*(t(3).p-t(2).p);
> end;
>
> % FUNCTION EVALUATION
> u.f = funeval(u.p);
>
> if (t(2).p < u.p) == (u.p < t(3).p),
113,114d134
< ulim = t(2).p+100*(t(3).p-t(2).p);
< if (t(2).p-u.p)*(u.p-t(3).p) > 0.0,
116d135
< u.f = linmineval(u.p);
127,151d145
< % try golden search instead
< u.p = t(3).p+gold*(t(3).p-t(2).p);
< u.f = linmineval(u.p);
<
< elseif (t(3).p-u.p)*(u.p-ulim) > 0.0
< % u is between t(3) and ulim
< u.f = linmineval(u.p);
< if u.f < t(3).f,
< % still no minimum as function is still decreasing
< % t(1) = t(2);
< t(2) = t(3);
< t(3) = u;
< u.p = t(3).p+gold*(t(3).p-t(2).p);
< u.f = linmineval(u.p);
< end;
<
< elseif (u.p-ulim)*(ulim-t(3).p) >= 0.0,
< % gone too far - constrain it
< u.p = ulim;
< u.f = linmineval(u.p);
<
< else,
< % try golden search instead
< u.p = t(3).p+gold*(t(3).p-t(2).p);
< u.f = linmineval(u.p);
163c157
< function [f,p] = brents(t, tol)
---
> function [f,p] = search(t, tol)
166,167c160
< % 1 - golden ratio
< Cgold = 1 - (sqrt(5)-1)/2;
---
> gold1 = 1-(sqrt(5)-1)/2;
173,181d165
< % t(1) and t(3) bracket the minimum
< if t(1).p>t(3).p,
< brk(1) = t(3).p;
< brk(2) = t(1).p;
< else,
< brk(1) = t(1).p;
< brk(2) = t(3).p;
< end;
<
183,190c167,169
< tmp = t(1);
< t(1) = t(2);
< t(2) = tmp;
< if t(2).f>t(3).f,
< tmp = t(2);
< t(2) = t(3);
< t(3) = tmp;
< end;
---
> [junk,ind] = sort(cat(1,t.f));
> t = t(ind);
> brk = [min(cat(1,t.p)) max(cat(1,t.p))];
214,215c193
< % and (not sure if it is necessary) that the solution is a minimum
< % rather than a maximum
---
> % and that the solution is a minimum rather than a maximum
217c195
< if abs(d) >= abs(ppd)/2 | u.p <= brk(1)+eps2 | u.p >= brk(2)-eps2 |
pol(3)<=0,
---
> if abs(d) > abs(ppd)/2 | u.p < brk(1)+eps2 | u.p > brk(2)-eps2 |
pol(3)<=0,
220c198
< d = Cgold*(brk(1)-t(1).p);
---
> d = gold1*(brk(1)-t(1).p);
222c200
< d = Cgold*(brk(2)-t(1).p);
---
> d = gold1*(brk(2)-t(1).p);
228c206
< u.f = linmineval(u.p);
---
> u.f = funeval(u.p);
239c217
< if u.f <= t(2).f | t(1).p==t(2).p,
---
> if u.f <= t(2).f,
242c220
< elseif u.f <= t(3).f | t(1).p==t(3).p | t(2).p==t(3).p,
---
> elseif u.f <= t(3).f,
247d224
< warning('Too many iterations in Brents');
252c229
< function linmin_plot(action,arg1,arg2,arg3,arg4)
---
> function min1d_plot(action,arg1,arg2,arg3,arg4)
254c231
< global linminplot
---
> persistent min1dplot
257c234
< linmin_plot('Init');
---
> min1d_plot('Init');
274,275c251,252
< linminplot =
struct('pointer',get(fg,'Pointer'),'name',get(fg,'Name'),'ax',[]);
< linmin_plot('Clear');
---
> min1dplot =
struct('pointer',get(fg,'Pointer'),'name',get(fg,'Name'),'ax',[]);
> min1d_plot('Clear');
278c255
< linminplot.ax = axes('Position', [0.15 0.1 0.8
0.75],...
---
> min1dplot.ax = axes('Position', [0.15 0.1 0.8
0.75],...
280c257
< lab = get(linminplot.ax,'Xlabel');
---
> lab = get(min1dplot.ax,'Xlabel');
282c259
< lab = get(linminplot.ax,'Ylabel');
---
> lab = get(min1dplot.ax,'Ylabel');
284c261
< lab = get(linminplot.ax,'Title');
---
> lab = get(min1dplot.ax,'Title');
287c264
<
'LineWidth',2,'Tag','LinMinPlot','Parent',linminplot.ax,...
---
>
'LineWidth',2,'Tag','LinMinPlot','Parent',min1dplot.ax,...
309,312c286,289
< if isstruct(linminplot),
< if ishandle(linminplot.ax), delete(linminplot.ax);
end;
< set(fg,'Pointer',linminplot.pointer);
< set(fg,'Name',linminplot.name);
---
> if isstruct(min1dplot),
> if ishandle(min1dplot.ax), delete(min1dplot.ax); end;
> set(fg,'Pointer',min1dplot.pointer);
> set(fg,'Name',min1dplot.name);
On Thursday 26 June 2008 00:36, Veni, Gopalkrishna wrote:
> Hello SPM group,
>
>
>
> I wonder if there are any differences in processing between
> SPM5 and SPM2 for co-registration. If so, can you please let me know?
>
>
>
> Thanks in advance,
>
> Gopal
>
>
>
> Gopalkrishna Veni | Imaging Programmer
>
> Nevada Cancer Institute | Imaging
>
> One Breakthrough Way, Las Vegas, NV 89135
> T: 702.822.5221 | F: 702.944.2376| C: 702.217.6257
> [log in to unmask] <mailto:[log in to unmask]> |
> www.nevadacancerinstitute.org <http://www.nevadacancerinstitute.org>
>
>
>
>
>
> --------------------------------------------------------------------------
> CONFIDENTIALITY NOTICE: This e-mail message, including any
> attachments, is for the sole use of the intended
> recipient(s) and may contain confidential, proprietary,
> and/or privileged information protected by law. If you are
> not the intended recipient, you may not use, copy, or
> distribute this e-mail message or its attachments. If you
> believe you have received this e-mail message in error,
> please contact the sender by reply e-mail and destroy all
> copies of the original message
|
