nmervis(sequence,n,varargin)

function nmervis(sequence,n,varargin)

% Set default display options
viewstyle = '3D';
transp = 'off';

% Read and set any varargin options
propertyArgIn = varargin;
while length(propertyArgIn) >= 2,
    prop = propertyArgIn{1};
    val = propertyArgIn{2};
    propertyArgIn = propertyArgIn(3:end);
    switch lower(prop)
        case 'view'
            viewstyle = val;
        case 'scalealpha'
            transp = val;
        otherwise
        error('nmervis:InvalidProperty', ...
            [prop,' is not a valid nmervis option'])
    end
end

% Get nmercount
counts = nmercount(sequence,n);

% Display Sierpinski gasket of depth n
figure('Color','w');
tetrahandles = biosierpinski(n);

% Set default look of tetrahedra
set(tetrahandles,'EdgeColor','none');
set(tetrahandles,'FaceColor','flat');
set(tetrahandles,'FaceAlpha',0);

% Find greatest count
greatestcount = 0;
for i=1:n
    if counts{i,2}>greatestcount
        greatestcount = counts{i,2};
    end
end

% Set color and transparency of tetrahedra according
% to sequence counts
for i=1:size(counts)
    tetra = findobj('Tag',counts{i,1});
    set(tetra,'FaceVertexCData',counts{i,2},...
        'CDataMapping','direct');
    if strcmp(transp,'on')
        set(tetra,'FaceAlpha',(counts{i,2}/greatestcount)^2);
    else
        set(tetra,'FaceAlpha',0.1);
    end
end

% Add labels
text(0,0,0,'A');
text(0,1,1,'C');
text(1,0,1,'G');
text(1,1,0,'T');

% Set viewpoint
if strcmp(viewstyle,'flat')
    view(30,90);
elseif strcmp(viewstyle,'3D')
    view(30,45);
end

% Add colorbar
map=jet(greatestcount);
colormap(map);
colorbar;

% Make axes invisible
set(gca,'Visible','off');

%-----------------------------------------------------------
% Display a Sierpinski gasket, setting the Tag of each tetrahedron
% to the appropriate sequence.

% biosierpinski calls itself and increasedepth recursively to
% create a gasket of arbitrary depth.
function tetrahandles = biosierpinski(depth)

if depth ~= 0
    tetrahandles = biosierpinski(depth-1);
    newtetrahandles = increasedepth(tetrahandles);
    tetrahandles = newtetrahandles;
else
    vert = [0 0 0;0 1 1;1 0 1;1 1 0];
    faces = [1 2 3;1 2 4;1 3 4;2 3 4];
    tcolor = [0 0 0;0 0 0;0 0 0;0 0 0];
    tetrahandles = patch('Faces',faces,...
        'Vertices',vert,...
        'FaceVertexCData',tcolor,...
        'Tag','');
end

%-----------------------------------------------------------
function newtetrahandles = increasedepth(tetrahandles)

for h = 1:size(tetrahandles,2)
    vert = get(tetrahandles(h),'Vertices');
    tag = get(tetrahandles(h),'Tag');
    for i = 1:4
        for j= 1:4
            if i==j
                newvert(j,:)=vert(j,:);
            else
                newvert(j,:)=(vert(j,:)+vert(i,:))/2;
            end
        end
        newtag = [tag, int2nt(i)];
        faces = [1 2 3;1 2 4;1 3 4;2 3 4];
        tcolor = [0 0 0;0 0 0;0 0 0;0 0 0];
        newtetrahandles((h*4)+(i-4)) = patch('Faces',faces,...
            'Vertices',newvert,...
            'FaceVertexCData',tcolor,...
            'Tag',newtag);
    end
end
delete(tetrahandles);