Filing the entire volume of a cube with small cube

I have built a hollow cube in MATLAB, I want to completely fill its volume with small cubes. Then I want to find a way to access these cubes and make paths through them, i.e if cube x is currently accessed, there should be a way to know what is its right, left, top, bottom, front, and back closest neighbors (closest neighbors= the cubes directly beside the current cube). I think we have 6 neighbors, because we have 6 different faces of the cube.

By knowing the nearest cube at every direction, a path through the cubes can be defined as a series of steps (eg, right, left, left, top, right, front). I think to be able to access every small cube and move to near ones we need to represent the small cubes in a matrix (maybe 3D), where if a small cube has a neighbor cube x to its right, then in the matrix, x will appear in the column next to the current column of the small cube. Also, if there is a direct neighbor at another depth layer (same x, y coordinates but different z coordinate, eg, the front and back neighbors) it should be indicated. Is there an easier way to identify the neighbors?

I've obtained a code by rayryeng (Building a hollow cube and filing it with small cubes in MATLAB) to fill many small cubes randomly inside the big one and build a 3D matrix, where every slice of the matrix (depth) represents one of the small cubes and the rows and columns of every slice (8 rows and 3 columns) represent the x y z coordinates of the vertices of every small cube. Please take a look at the question link that I provided to see the code.

I want to make two modifications to the code,

1- fill the big cube with small cubes in an organized way not randomly.

2- adjust the 3D matrix to represent how the small cubes are neighbors to each other.

I've tried to adjust the code in the linked question to fill the cube in an organized way. This is my trial (I added the if else in the for loop),

   clf;
figure(1);
format compact 
h(1) = axes('Position',[0.2 0.2 0.6 0.6]);

%These are the different 8 vertices of the cube, each is defined by its 3 x
%y z coordinates:
vert = [1 1 -1; 
        -1 1 -1; 
        -1 1 1; 
        1 1 1; 
        -1 -1 1;
        1 -1 1; 
        1 -1 -1;
        -1 -1 -1];
%These are the 6 faces of the cube, each is defined by connecting 4 of the
%available vertices:
fac = [1 2 3 4; 
       4 3 5 6; 
       6 7 8 5; 
       1 2 8 7; 
       6 7 1 4; 
       2 3 5 8];

% I defined a new cube whose length is 1 and centers at the origin.
vert2 = vert * .05;  
fac2 = fac;


patch('Faces',fac,'Vertices',vert,'Facecolor', 'w');  % patch function for the first big cube. 
axis([-1, 1, -1, 1, -1, 1]);
axis equal;

hold on;

patch('Faces', fac2, 'Vertices', vert2, 'FaceColor', 'r', 'EdgeColor', 'none');
material metal;
alpha('color');
alphamap('rampdown');
%view(3);


hold on;
rng(123); %// Set seed for reproducibility
num_squares = 1000; %// Set total number of squares

%// New - to store the coordinates
coords = [];

%// For remembering the colours
colors = [];
%// For each square...

for idx = 1 : num_squares


    %// Take the base cube and add an offset to each coordinate
    %// Each coordinate will range from [-1,1]
    if (idx==1)
    vert_new = bsxfun(@plus, vert2, [.01 .01 .01]);

    else 
    vert_new = bsxfun (@plus, vert_new,[.01 .01 .01] );
    end
    %// New - For the coordinates matrix
    coords = cat(3, coords, vert_new);

    %// Generate a random colour for each cube
    color = rand(1,3);

    %// New - Save the colour
    colors = cat(1, colors, color);

    %// Draw the cube
    patch('Faces', fac, 'Vertices', vert_new, 'FaceColor', color,'EdgeColor', 'none');
end




%// Post processing
material metal;
alpha('color');
alphamap('rampdown');
view(3);

But I'm getting a result like the one in the following picture,

Can anyone please tell me how to solve this problem and build the 3D matrix (or any other easier way to represent the neighbors of every cube )?

EDIT: To elaborate more on the small cube neighbor issue. Consider having a cube C anywhere inside the big cube. Let the cube position be (.5,.5,.5) which we may think of as the center of the cube. Then, this cube will have 6 neighbor cubes (positioned directly beside it) with 2 neighbors per axis. So, for the cube (.5,.5,.5) we have,

x axis neighbor to the right of C (.5 +offset, .5, .5)

x axis neighbor to the left of C (.5 -offset, .5, .5)

y axis neighbor to the top of C (.5 , .5+offset, .5)

y axis neighbor to the bottom of C (.5 , .5-offset, .5)

z axis neighbor one depth after C (.5 , .5, .5+offset)

z axis neighbor one depth before C (.5 , .5, .5-offset)

where offset may be seen as the distance between the center of a cube and any of its faces. This is just an explanation to clarify the idea, it doesn't need to be implemented the same way. I hope this is clear, I would appreciate any help in constructing this neighborhood matrix.

Thanks.

Same principle, we build a big cube, then a small one in the corner, then we repeat the small cube building with a small offset until we are full. The main difference with the old code is that this time the step change of each coordinate set is controlled (function of x,y,z coordinates of the small cube) instead of being random.

%%
clf; figure(1); format compact 
h(1) = axes('Position',[0.2 0.2 0.6 0.6]);

%These are the different 8 vertices of the cube, each is defined by its 3 x y z coordinates:
vert = [ 1  1 -1; -1  1 -1; -1  1  1; 1  1  1; -1 -1  1; 1 -1  1; 1 -1 -1; -1 -1 -1];

%These are the 6 faces of the cube, each is defined by connecting 4 of the available vertices:
fac = [1 2 3 4; 4 3 5 6; 6 7 8 5; 1 2 8 7; 6 7 1 4; 2 3 5 8];

%// How many small cube do we want
MainCubeSide = 2 ;              %// dimension of the side of the main cube
nCubeOnSide = 5 ;               %// number of small cube in one "row/column" of the main cube
nCubesTotal = nCubeOnSide^3  ;  %// total number of small cube

% define the Main container cube
MainCube.Vertices = vert *(2/MainCubeSide) ; %// because the cube as defined above has already a side=2
MainCube.Faces = fac ;
MainCube.FaceColor = 'w' ;

hMainCube = patch(MainCube);  %// patch function for the first big cube. 
axis([-1, 1, -1, 1, -1, 1]);
axis equal;
hold on;
material metal;
alpha('color');
alphamap('rampdown');
view(138,24)
%view(3);


%% // generate all the coordinates of each cube first
dstep = MainCubeSide / nCubeOnSide ;                                                 %// step size for small cube vertices
vElem = bsxfun(@plus, vert / nCubeOnSide , -( MainCubeSide/2 - dstep/2)*[1 1 1] )  ; %// elementary cube vertices

%%
hold on;
coords = zeros( size(vElem,1),size(vElem,2), nCubesTotal ) ;  %// To store the coordinates
colors = zeros( nCubesTotal , 3 ) ;                           %// To store the colours
hcube  = zeros( nCubesTotal , 1 ) ;                           %// To store the handles of the patch objects

iNeighbour = zeros( nCubesTotal , 6 ) ;   %// To save the index of the neighbours
idc = permute( reshape(1:nCubesTotal,nCubeOnSide,nCubeOnSide,nCubeOnSide) , [3 2 1] ) ;

%// For each cube ...
iCube = 0 ;
for iline=1:nCubeOnSide         %// Lines
    for icol=1:nCubeOnSide      %// Columns
        for ih=1:nCubeOnSide    %// Slice (height)
            iCube = iCube + 1 ;

            %// Take the base corner coordinates and add an offset to each coordinate
            coords(:,:,iCube) = bsxfun(@plus, vElem , dstep*[(iline-1) (icol-1) (ih-1)]);

            %// Save the colour
            colors(iCube,:) = rand(1,3) ; 

            %// Draw the cube
            hcube(iCube) = patch('Faces', fac, 'Vertices', coords(:,:,iCube), 'FaceColor', colors(iCube,:) ) ;

            drawnow     %// just for intermediate display, you can comment these 2 lines
            pause(0.05) %// just for intermediate display, you can comment these 2 lines

            %// save adjacent cubes indices
            ixAdj = [iline-1 iline+1 icol-1 icol+1 ih-1 ih+1] ;  %// indices of adjacent cubes
            idxFalse = (ixAdj<1) | (ixAdj>nCubeOnSide) ;  %// detect cube which would be "out" of the main cube
            ixAdj(idxFalse) = 1 ;                                %// just to not get an "indexing" error at this stage
            iNeighbour(iCube,:) = [idc(ixAdj(1),icol,ih)    idc(ixAdj(2),icol,ih) ...
                                   idc(iline,ixAdj(3),ih)   idc(iline,ixAdj(4),ih) ...
                                   idc(iline,icol,ixAdj(5)) idc(iline,icol,ixAdj(6)) ] ;
            iNeighbour(iCube,idxFalse) = NaN ;
        end
    end
end

This code saves the handle of each cube in the variable hcube so you can do assignment of property in bulk on all the cubes if you want. For example delete(hcube) will delete all the small cubes in one go, or set(hcube,'Facealpha',0.5) will make all cubes half transparent.

You can also set/change properties on only a subset of them hcube(idx_subset) = .... This is where knowing the adjacent cubes by their index could be useful, but your adjacency question is not completely defined yet.

Edit: I have added the neighbour tracking in the main loop. It is probably not the most efficient way to do it but it does keep an index of all the neighbour fro each elementary cube. The iNeighbour variable (size: nCubesx6) hold the handle index of each neighbour (6 possible neighbours). When a neighbour didn't exist, I opted to place a NaN instead. To retrieve the index of the neighbour directly without the NaNs, I defined a helper anonymous function:

getNeighbourIndex = @(idx) iNeighbour(idx,~isnan(iNeighbour(idx,:))) ;

Which now help you track all the neighbour of a given cube. For example:

set(hcube,'Visible','off')  %// turn off all small cubes
CubeOfInterest = 111 ;      %// select one cube
%// display the main cube of interest, and it's neighbours in transparency
set(hcube(CubeOfInterest),'Visible','on','FaceColor','r','FaceAlpha',1) 
set(hcube(getNeighbourIndex(CubeOfInterest)),'Visible','on','FaceColor','g','FaceAlpha',.05)