function V = spm_write_vol(V,Y)
% Write an image volume to disk, setting scales and offsets as appropriate
% FORMAT V = spm_write_vol(V,Y)
% V (input)  - a structure containing image volume information (see spm_vol)
% Y          - a one, two or three dimensional matrix containing the image voxels
% V (output) - data structure after modification for writing.
%
% Note that if there is no 'pinfo' field, then SPM will figure out the
% max and min values from the data and use these to automatically determine
% scalefactors.  If 'pinfo' exists, then the scalefactor in this is used.
%_______________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% John Ashburner
% $Id$

use_offset = false;

if ndims(Y)>3, error('Can only handle a maximum of 3 dimensions.'), end
dim = [size(Y) 1 1 1];
if ~all(dim(1:3) == V.dim(1:3)),
    error('Incompatible dimensions.');
end

if ~isfield(V,'pinfo'),
    V.pinfo = [1;0;0];
    rescal  = true;
elseif ~all(isfinite(V.pinfo(1:2))) || V.pinfo(1) == 0,
    V.pinfo(1:2) = [1;0];
    rescal  = true;
else
    rescal = false;
end

if rescal,
    % Set scalefactors and offsets
    %-----------------------------------------------------------------------
    dt           = V.dt(1);
    s            = find(dt == [2 4 8 256 512 768]);
    if isempty(s),
        V.pinfo(1:2) = [1;0];
    else
        dmnmx        = [0 -2^15 -2^31 -2^7 0 0 ; 2^8-1 2^15-1 2^31-1 2^7-1 2^16-1 2^32-1];
        dmnmx        = dmnmx(:,s);
        mxs          = zeros(dim(3),1)+NaN;
        mns          = zeros(dim(3),1)+NaN;

        for p=1:dim(3),
            tmp    = double(Y(:,:,p));
            tmp    = tmp(isfinite(tmp));
            if ~isempty(tmp),
                mxs(p) = max(tmp);
                mns(p) = min(tmp);
            end
        end

        mx = max(mxs(isfinite(mxs)));
        mn = min(mns(isfinite(mns)));
        if isempty(mx), mx = 0; end;
        if isempty(mn), mn = 0; end;
        if mx~=mn,
            if use_offset,
                V.pinfo(1,1) = (mx-mn)/(dmnmx(2)-dmnmx(1));
                V.pinfo(2,1) = (dmnmx(2)*mn-dmnmx(1)*mx)/(dmnmx(2)-dmnmx(1));
            else
                if dmnmx(1)<0,
                    V.pinfo(1) = max(mx/dmnmx(2),mn/dmnmx(1));
                else
                    V.pinfo(1) = mx/dmnmx(2);
                end
                V.pinfo(2) = 0;
            end
        else
            V.pinfo(1,1) = mx/dmnmx(2);
            V.pinfo(2,1) = 0;
        end
    end
end

%-Create and write image
%-----------------------------------------------------------------------
V = spm_create_vol(V);
V = spm_write_plane(V,Y,':');

