function drc_vbm_smooth_native(P,T,Q)
% Smooth native segmentations like spm_preproc_write does when creating
% warped versions of them. 
% Smoothed ims are saved as float32 (like warped ims) to avoid data-loss.
%
% P  - native seg image to be smoothed
% T  - template image (e.g. pick tpm/grey.nii)
% Q  - optional new filename for smoothed image (otherwise overwrite orig)
%
% (modified from spm_smooth and spm_preproc_write)

dtype = spm_type('float32');
if ischar(P); P = spm_vol(P); end
if nargin < 3; Q = P.fname; end
if ischar(T); T = spm_vol(T); end
if isstruct(P),
	ivx = sqrt(sum(P.mat(1:3,1:3).^2));
else
	ivx = [1 1 1];
end;

% determine amount of smoothing (from spm_preproc_write)
ovx = abs(det(T.mat(1:3,1:3)))^(1/3);
fwhm = max(ovx./ivx - 1, 0.1);

% setup output volume
if ischar(Q) && isstruct(P),
	[pth,nam,ext,num] = spm_fileparts(Q);
	q         = fullfile(pth,[nam,ext]);
	Q         = P;
	Q.fname   = q;
	if ~isempty(num),
		Q.n       = str2double(num);
	end;
	if ~isfield(Q,'descrip'), Q.descrip = sprintf('SPM compatible'); end;
	% Q.descrip = sprintf('%s - conv(%g,%g,%g)',Q.descrip, fwhm);
	if dtype~=0 && (finite(spm_type(P.dt(1),'maxval')) || ~finite(spm_type(Q.dt(1),'maxval'))),
		Q.dt(1) = dtype;
		if ~finite(spm_type(dtype,'maxval')),
			Q.pinfo = [1 0 0]';
		else
			maxv = spm_type(P.dt(1),'maxval')*P.pinfo(1) + P.pinfo(2);
			sf   = maxv/spm_type(Q.dt(1),'maxval');
			Q.pinfo = [sf 0 0]';
		end;
	end;
end;

img = spm_read_vols(P); % can now safely overwrite P

% compute parameters for spm_conv_vol (from spm_preproc_write)
lim = ceil(2*fwhm);
x  = -lim(1):lim(1); x = spm_smoothkern(fwhm(1),x); x  = x/sum(x);
y  = -lim(2):lim(2); y = spm_smoothkern(fwhm(2),y); y  = y/sum(y);
z  = -lim(3):lim(3); z = spm_smoothkern(fwhm(3),z); z  = z/sum(z);
i  = (length(x) - 1)/2;
j  = (length(y) - 1)/2;
k  = (length(z) - 1)/2;

% smooth
% (Uncomment the two data-type conversion lines below to exactly match
% spm5's modulated-warped output images.)
% img = uint8(img*255); % (as spm_preproc_write -- apply on uint8 data!)
spm_conv_vol(img,img,x,y,z,-[i,j,k]); % (works in-place)
% img = single(img)/255; % (convert back to float32)

spm_write_vol(Q,img);