%  Copyright (c) 2010, Fred Sigworth
%  All rights reserved.
%
%  Redistribution and use in source and binary forms, with or without
%  modification, are permitted provided that the following conditions are
%  met:
%
%      * Redistributions of source code must retain the above copyright
%        notice, this list of conditions and the following disclaimer.
%      * Redistributions in binary form must reproduce the above copyright
%        notice, this list of conditions and the following disclaimer in
%        the documentation and/or other materials provided with the distribution
%
%  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
%  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
%  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
%  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
%  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
%  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
%  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
%  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
%  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
%  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
%  POSSIBILITY OF SUCH DAMAGE.
%
% <http://www.mathworks.com/matlabcentral/fileexchange/27021-imagic--mrc-and-dm3-file-i-o>

function WriteMRC(map,rez,filename)
% function WriteMRC(map,rez,filename)
% Write out a 2D image or a 3D volume as an MRC map file, for example for viewing in
% Chimera.  'map' is the 3D array, rez is the voxel size in angstroms.
%
% fixed for 2D output as well.  fs 28 Aug 07
% 
% % Test program: create an ellipsoid and store it as a map file.
%   [x y z]=ndgrid(-32:31);
%   m=(x.^2+(y*1.5).^2+z.^2)>20^2;
%   WriteMRC(m,1,'maptest.mrc');
  
f=WriteMRCHeader(map,rez,filename);
count2=fwrite(f,map,'single');
fclose(f);

return;


function handle=WriteMRCHeader(map,rez,filename,nim)
% function handle=WriteMRCHeader(sizes,rez,filename,nim)
% Write out the header of an MRC map file, and leave the file open,
% returning the file handle, so that data can be written sequentially into
% the file and then the file closed.  Data are written in little-endian style,
% as float32 values.
% If you want to write more images than are contained in map, give the
% number in the optional nim argument.
%
% Example: write out 10,000 images.
%     images=randn(64,64,1000);
%     f=WriteMRCHeader(images,2.8,'test.mrc',10000);
%     fwrite(f,images,'float32');
%     for i=2:10
%         images=randn(64,64,1000);
%         fwrite(f,images,'float32');
%     end;
%     fclose(f);

% Files are always written in little-ended format.
% Figure out if we have a little-ended machine.
q=typecast(int32(1),'uint8');
machineLE=(q(1)==1);  % true for little-endian machine

hdr=int32(zeros(256,1));

sizes=size(map);

if numel(sizes)<3
    sizes(3)=1;
end;
if nargin >3
    sizes(3)=nim;
end;

% Get statistics

map=reshape(map,numel(map),1);  % convert it into a 1D vector
theMean=mean(map);
theSD=std(map);
theMax=max(map);
theMin=min(map);


hdr(1:3)=sizes; % number of columns, rows, sections
hdr(4)=2;  % mode: real, float values
hdr(8:10)=hdr(1:3);  % number of intervals along x,y,z
hdr(11:13)=typecast(single(single(hdr(1:3))*rez),'int32');  % Cell dimensions
hdr(14:16)=typecast(single([90 90 90]),'int32');   % Angles
hdr(17:19)=(1:3)';  % Axis assignments
hdr(20:22)=typecast(single([theMin theMax theMean]'),'int32');
hdr(23)=0;  % Space group 0 (default)
if machineLE
    hdr(53)=typecast(uint8('MAP '),'int32');
    hdr(54)=typecast(uint8([68 65 0 0]),'int32');  % LE machine stamp.
else
    hdr(53)=typecast(uint8(' PAM'),'int32');  % LE machine stamp, for writing with BE machine.
    hdr(54)=typecast(uint8([0 0 65 68]),'int32');
end

hdr(55)=typecast(single(theSD),'int32');

handle=fopen(filename,'w','ieee-le');
count1=fwrite(handle,hdr,'int32');
