plot2bmp(Z,varargin)

function plot2bmp(Z,varargin)

% PLOT2BMP   Convert 2D Matlab Arrays Into .BMP Images
%
%   x = linspace(-pi,pi,201);
%   [X,Y] = meshgrid(x);
%   Z = sin(X) + sin(Y);
%   plot2bmp(Z) 
%   % % saves a bitmap image of 2D array Z with filename "plot2bmp.bmp". Maps
%   % % the array Z to a hue scale.
%               
%   options.filename = 'my_pic';
%   plot2bmp(Z,options) 
%   % % specify a file name prefix.
%
%   options.colormap = 'bw';
%   plot2bmp(Z,options)
%   % % All possible colormaps are 'bw' (black and white) 'hot', 'two-tone', 
%   % % 'full' (full spectrum), or 'basic' (default).
%
%   options.invert=true;
%   plot2bmp(Z,options) 
%   % % Invert the colors
%
%   x = linspace(-3,3,201);
%   [X,Y] = meshgrid(x);
%   Z = (X + 1i*Y).^3;
%   plot2bmp(Z,options) 
%   % % Method for generating complex plots
%
%   plot2bmp(H,S,V) 
%   % % Specify all HSV data
%
%   plot2bmp(H,V)
%   % % Specify all HV data
%
%   options.rgbvals = true;
%   plot2bmp(R,G,B,options)
%   % % Specify all RGB data
%
%   x = linspace(-pi,pi,201);
%   [X,Y] = meshgrid(x);
%   Z = sin(X) + sin(Y);
%   options.range = [-.75,.75];
%   options.filename = 'plot2bmp1';
%   plot2bmp(Z,options)
%   options.range = [-2,2];
%   options.filename = 'plot2bmp2';
%   plot2bmp(Z,options)
%   % % Specify a range for the data
%
%   x = linspace(-pi,pi,201);
%   [X,Y] = meshgrid(x);
%   Z = sin(X) + sin(Y);
%   options.colormap = 'two-tone';
%   options.rgbperm = 5;
%   plot2bmp(Z,options)
%   % % Do 1 of 6 permutations on the RGB values


%    License: BSD
%    Copyright 2013 Bryan C. Smith
%    $Revision: 1.0 $ $Date: 27-May-2013
%    Tested in MATLAB 7.13.0.564 (R2011b)



filename = 'plot2bmp';
colormap = 'basic';
invert = false;
rgbperm = 1;
range = [];
rgbvals = false;

p = 1;
if nargin > 1
    if isreal(varargin{1})
        p = 2;
        if nargin > 2
            if isreal(varargin{2})
                p = 3;
                V = varargin{2};
                S = varargin{1};
                if nargin > 3
                    options = varargin{3};
                end
            else
                V = varargin{1};
                options = varargin{2};
            end
        else
            V = varargin{1};
            options = struct;      
        end
    else
        options = varargin{1};
    end
else
    options = struct;
end


names = fieldnames(options);
for kkk = 1:length(names)
    eval([names{kkk} ' = options.' names{kkk} ';']);
end


[m,n] = size(Z);
step = ceil(5000/m);


if p > 1
    colormap = 'full';
end
if rgbvals == true
    if ~(p==3)
        error('For RGB, must have n x m x 3 data vector');
    end
    range = [];
end

complexvals = ~isreal(Z);
if complexvals
    if p > 1
        error('No method for complex 3D array')
    end    

    V = abs(Z);    
    Z = angle(Z);
    
    vmax = max(max(V));
    vmin = min(min(V));
    fmax = pi + pi/3;
    fmin = -pi;
else
    fmax = max(max(Z(:,:,1)));
    fmin = min(min(Z(:,:,1)));
end

if length(range)==2
    if complexvals
        if (range(2)<vmax)
            V = min(V,range(2));
        end
        if (range(1)>vmin)
            V = max(V,range(1));
        end
        vmax = range(2);
        vmin = range(1);
    else
        if (range(2)<fmax)
            Z(:,:,1) = min(Z(:,:,1),range(2));
        end
        if (range(1)>fmin)
            Z(:,:,1) = max(Z(:,:,1),range(1));
        end
        fmax = range(2);
        fmin = range(1);
    end
end

if complexvals
    ss = vmax-vmin;
    if ss/abs(vmax+vmin) < 1e-14
        V = ones(size(V));
    else
        V = (V-vmin)/ss;
    end
    V = V.^.5;
end

if p==1
    % Map F to [0, 1];
    scale = (fmax-fmin);
    if scale == 0
        Z = zeros(size(Z));
    else
        Z = (Z - fmin)/scale;
    end
end

a = 1:6;
if sum(rgbperm==a)~= 1
	error('plot2bmp:IncorrectRGBPermutation',...
        '%s','Permutation of RGB colors must be an integer between 1 and 6.');
end

savfile = [filename '.bmp'];

outFile = fopen(savfile,'w','ieee-le');
C = p2bheader(m,n);
fwrite(outFile, C, 'uchar');

for count = 1:step:n
    j = count:min(n, count + step - 1);
    if complexvals
        HH = Z(:,j); 
        VV = V(:,j);
        SS = 1;
    elseif p==1
        HH = Z(:,j);     
        VV = 1;
        SS = 1;
    elseif p==2
        HH = Z(:,j); 
        VV = V(:,j);
        SS = 1;
    else
        HH = Z(:,j);
        VV = V(:,j);
        SS = S(:,j);
    end

    HH = p2bcmap(HH,colormap);

    if rgbvals == false
        [R,G,B] = getrgbvals(HH,SS,VV,colormap,rgbperm,invert);
    else
        R = HH;
        G = SS;
        B = VV;
    end    

	C = p2bdatatobmp(R,G,B,m,length(j));
	fwrite(outFile, C, 'uchar');
end

fclose(outFile);


function [R,G,B] = getrgbvals(H,S,V,colormap,rgbperm,invert)

if strcmp(colormap,'bw')
    [R,G,B] = p2bhue2rgb(1,0,H);
elseif rgbperm == 1
    [R,G,B] = p2bhue2rgb(H,S,V);
elseif rgbperm == 2
    [G,B,R] = p2bhue2rgb(H,S,V);
elseif rgbperm == 3
    [B,R,G] = p2bhue2rgb(H,S,V);
elseif rgbperm == 4
    [R,B,G] = p2bhue2rgb(H,S,V);
elseif rgbperm == 5
    [B,G,R] = p2bhue2rgb(H,S,V);
else
    [G,R,B] = p2bhue2rgb(H,S,V);
end

if invert == true;
    R = 255-R;
    G = 255-G;
    B = 255-B;
end

function F = p2bcmap(F,colormap)
% Apply a colormap
if strcmp(colormap,'basic')
    F = (5/6)*F;
elseif strcmp(colormap,'twotone')
    F = 0.667*F + 0.5*(1-F);
elseif strcmp(colormap,'hot')
    F = (1/6)*(1-F);
elseif ~( strcmp(colormap,'full')||strcmp(colormap,'bw') )
    error('Colormap Undefined')
end


function C = p2bdatatobmp(R, G, B, w, l)
pad = mod(w,4);

% size of partition data entries.
sz3 = (3*w + pad)*l;

C = zeros(sz3,1);
ii = 3*(1:w);

for j = 1:l
    jj = (3*w + pad)*(j-1);
    C(ii + jj - 2) = B(:,j);
    C(ii + jj - 1) = G(:,j);
    C(ii + jj) = R(:,j);
end

function a = p2bheader(w, l)
% Creates header info for .bmp file.

% padded size, width, length
pad = mod(w,4);

% number of data entries
sz2 = (3*w + pad)*l;

% full file size
sz = sz2 + 54;
sprintf('File size is %7.6e',sz);

a = zeros(54,1);
a(1:2,:) = [66 77]';

[hhex, length] = p2bint2hex(sz);
hex = zeros(length,1);
for i = 1:length
    hex(i) = hhex(length-i+1);
end

dist = 2 + length;
a(3:dist) = hex;

n = 8 - length;
z = zeros(n,1);

olddist = dist;
dist = dist + n + 8;
a((1+olddist):dist) = [z;  54; 0; 0; 0; 40; 0; 0; 0];

[hhex, length] = p2bint2hex(w);
hex = zeros(length,1);
for i = 1:length
    hex(i) = hhex(length-i+1);
end
olddist = dist;
dist = dist + length;
a((1+olddist):dist) = hex;

n = 4 - length;
z = zeros(n,1);
olddist = dist;
dist = dist + n;

a((1+olddist):dist) = z;

[hhex, length] = p2bint2hex(l);
hex = zeros(length,1);
for i = 1:length
    hex(i) = hhex(length-i+1);
end
olddist = dist;
dist = dist + length;

a((1+olddist):dist) = hex;

n = 4 - length;
z = zeros(n,1);
olddist = dist;
dist = dist + n + 8;
a((1+olddist):dist) = [z; 1; 0; 24; 0; 0; 0; 0; 0];

[hhex, length] = p2bint2hex(sz2);
hex = zeros(length,1);
for i = 1:length
    hex(i) = hhex(length-i+1);
end
olddist = dist;
dist = dist + length;
a((1+olddist):dist) = hex;

n = 20 - length;
z = zeros(n,1);

olddist = dist;
dist = dist + n;
a((1+olddist):dist) = z;


function [R,G,B] = p2bhue2rgb(H,varargin)

% H,S,L are matrices or vectors.

type = 1;
if nargin > 3
    if strcmp(varargin{3},'lightness')
        type = 2;
    elseif ~strcmp(varargin{3},'value')
        error('Input value unkown')
    end
end


if nargin > 2
    V = varargin{2};
else
    V = 1;
end


if nargin > 1
    S = varargin{1};
else
    S = 1;
end


d = 1/6;

R = zeros(size(H));
G = R;
B = G;

% Hue
J = (255*6)*H;
JJ = J - 510;
JJJ = J - 1020;

U = (H < d);
R = R + 255*U;
G = G + J.*U;

U = (H >= d & H < 2*d);
R = R - JJ.*U;
G = G + 255*U;

U = (H >= 2*d & H < 3*d);
G = G + 255*U;
B = JJ.*U;

U = (H >= 3*d & H <= 4*d);
G = G - JJJ.*U;
B = B + 255*U;

U = (H > 4*d & H <= 5*d);
R = R + JJJ.*U;
B = B + 255*U;

U = (H > 5*d);
R = R + 255*U;
B = B + (1530 - J).*U;

% Saturation
S = 1 - S;
R = R + (255 - R).*S;
G = G + (255 - G).*S;
B = B + (255 - B).*S;

% Value
if type == 1
    R = R.*V;
    G = G.*V;
    B = B.*V;
end

% Lightness
if type == 2
   U = (V < 0.5);
   V = V/.5;
   
   R(U) = R(U).*V(U);
   R(~U) = R(~U) + (255-R(~U)).*(V(~U)-1);
   
   G(U) = G(U).*V(U);
   G(~U) = G(~U) + (255-G(~U)).*(V(~U)-1);
   
   B(U) = B(U).*V(U);
   B(~U) = B(~U) + (255-B(~U)).*(V(~U)-1);
end

R = round(R);
G = round(G);
B = round(B);


function [hex, length] = p2bint2hex(num)

% Converts a integer to hexadecimal representation.

hex = [];
length = 1;
q = 0;
while q==0
    f = num/256^(length);
    if f < 1
        q = 1;
    else
        length = length+1;
    end
end

k = length - 1;
while k > 0
    f = floor(num/256^k);
    hex = [hex; f];
    num = num - f*256^k;
    k = k - 1;
end
hex = [hex; num];