%  cubic_manifold_general
%
%   Hardcoded scale similar coordination of a 3x3x3 set.
%   Outputs an nth order scale similar Space Volume.
%
%   Calls:
%       write_Space_volume
%
%   Mark Dow, 2/5/04
%   modified by Mark Dow, 12/30/04


function cubic_manifold_general()

filestem = 'Tetra_5';
nBorderWidth = 4;
nOrder = 5;             % volume size is O( 3^nOrder ), 15 MB for nOrder = 5

nSize = 3^nOrder

vol = int8(0);
vol( 1:nSize + 2*nBorderWidth, 1:nSize + 2*nBorderWidth, 1:nSize + 2*nBorderWidth ) = 0;

for i = 1 + nBorderWidth : nSize + nBorderWidth
    
    ii = i - nBorderWidth;
    
    for j = 1 + nBorderWidth : nSize + nBorderWidth
        
        jj = j - nBorderWidth;
        
        for k = 1 + nBorderWidth : nSize + nBorderWidth
            
            kk = k - nBorderWidth;
            
            bFill = 1;
            bEnd = 0;
            
            t = 0;
            iii = ii-1;
            jjj = jj-1;
            kkk = kk-1;
            iit = mod( iii, 3 ) - 1;
            jjt = mod( jjj, 3 ) - 1;
            kkt = mod( kkk, 3 ) - 1;
            
            while bFill == 1 & bEnd == 0
                
                t = t + 1;
                
                % Faces and center.
  %              if    (  abs(iit) + abs(jjt) + abs(kkt) ==  1  ) | ( iit==  0 & jjt ==  0 & kkt ==  0  )
  
                % Tetrahedral. 
                if     (  iit ==  -1 & jjt ==  -1 & kkt ==  -1  )...
                    |  (  iit ==   0 & jjt ==   0 & kkt ==  -1  )...
                    |  (  iit ==   1 & jjt ==   1 & kkt ==  -1  )...
                    |  (  iit ==   0 & jjt ==  -1 & kkt ==   0  )...
                    |  (  iit ==   1 & jjt ==   0 & kkt ==   0  )...
                    |  (  iit ==   1 & jjt ==  -1 & kkt ==   1  )...
                    |  (  iit ==   0 & jjt ==   1 & kkt ==   0  )...
                    |  (  iit ==  -1 & jjt ==   1 & kkt ==   1  )...
                    |  (  iit ==  -1 & jjt ==   0 & kkt ==   0  )...
                    |  (  iit ==   0 & jjt ==   0 & kkt ==   1  )
                    
                % Serpinski. 
%                 if     (  iit ==  -1 & jjt ==  -1 & kkt ==  -1  )...
%                     |  (  iit ==   0 & jjt ==  -1 & kkt ==  -1  )...
%                     |  (  iit ==   1 & jjt ==  -1 & kkt ==  -1  )...
%                     |  (  iit ==  -1 & jjt ==   0 & kkt ==  -1  )...
%                     |  (  iit ==   1 & jjt ==   0 & kkt ==  -1  )...
%                     |  (  iit ==  -1 & jjt ==   1 & kkt ==  -1  )...
%                     |  (  iit ==   0 & jjt ==   1 & kkt ==  -1  )...
%                     |  (  iit ==   1 & jjt ==   1 & kkt ==  -1  )...
%                     |  (  iit ==  -1 & jjt ==  -1 & kkt ==   0  )...
%                     |  (  iit ==   1 & jjt ==  -1 & kkt ==   0  )...
%                     |  (  iit ==  -1 & jjt ==   1 & kkt ==   0  )...
%                     |  (  iit ==   1 & jjt ==   1 & kkt ==   0  )...
%                     |  (  iit ==  -1 & jjt ==  -1 & kkt ==   1  )...
%                     |  (  iit ==   0 & jjt ==  -1 & kkt ==   1  )...
%                     |  (  iit ==   1 & jjt ==  -1 & kkt ==   1  )...
%                     |  (  iit ==  -1 & jjt ==   0 & kkt ==   1  )...
%                     |  (  iit ==   1 & jjt ==   0 & kkt ==   1  )...
%                     |  (  iit ==  -1 & jjt ==   1 & kkt ==   1  )...
%                     |  (  iit ==   0 & jjt ==   1 & kkt ==   1  )...
%                     |  (  iit ==   1 & jjt ==   1 & kkt ==   1  )...
                   
                     iii = round( floor(iii/3) );
                     jjj = round( floor(jjj/3) );
                     kkk = round( floor(kkk/3) );
                     iit = mod( iii , 3 ) - 1;
                     jjt = mod( jjj, 3 ) - 1;
                     kkt = mod( kkk, 3 ) - 1;                    
               else  
                    
                  % Don't fill, step out of while loop.
                    bFill = 0;       
                end
                
                if     t >= nOrder
                    bEnd = 1;
                end
                
            end
            
            % If a voxel is never excluded, mark it.
            if bFill == 1
                
                vol( i, j, k ) = 192;
            end
            
        end
    end
end

write_Space_volume( vol, filestem )