how to Drawing 3d Voxel data

I want to draw voxels in 3d space (I have them stored in a 3d matrix). i used patch to do this (or voxel function from file exchange, draws every face of each voxel). here is my code

        function  PlotVoxelinput()
        clear all;
        clc;
      %PlotVoxelInput Summary of this function goes here
      %   Detailed explanation goes here

    voxel_init_loc=[0,0,0];
    voxel_loc=voxel_init_loc;
    voxel_size=[1,1,1];
    file='Other inputs/nucleon.raw';
    %file='Lev input/Model 128/bone_model_3061.raw';
    fileID=fopen(file);
    [A, count] = fread(fileID);
    m_hight=nthroot(count,3);
    m_width=nthroot(count,3);
    m_depth=nthroot(count,3);
    %set figure parameters
    xlim([0 m_hight]);
    ylim([0 m_width]);
    zlim([0 m_depth]);
    view(45,45);
    daspect([1 1 1])
    %get 3d data
    VoxelMat= multibandread(file,[m_hight,m_width,m_depth],'uint8',0,'bip','ieee-le');
    for i=1:m_hight
        for j=1:m_width
            for k=1:m_depth
                if VoxelMat(i,j,k)==1
                voxel(voxel_loc,voxel_size,'b',1)
                end
                voxel_loc=[voxel_loc(1),voxel_loc(2),voxel_loc(3)+1];
            end
            voxel_loc=[voxel_loc(1),voxel_loc(2)+1,voxel_init_loc(3)];
        end
        voxel_loc=[voxel_loc(1)+1,voxel_init_loc(2),voxel_init_loc(3)];
    end

    fclose(fileID);
    end

here is voxel function from file exchange

     function voxel(i,d,c,alpha);

      %VOXEL function to draw a 3-D voxel in a 3-D plot
      %
      %Usage
      %   voxel(start,size,color,alpha);
      %
      %   will draw a voxel at 'start' of size 'size' of color 'color' and
      %   transparency alpha (1 for opaque, 0 for transparent)
      %   Default size is 1
      %   Default color is blue
      %   Default alpha value is 1
      %
      %   start is a three element vector [x,y,z]
      %   size the a three element vector [dx,dy,dz]
      %   color is a character string to specify color 
      %       (type 'help plot' to see list of valid colors)
      %   
      %
      %   voxel([2 3 4],[1 2 3],'r',0.7);
      %   axis([0 10 0 10 0 10]);
      %

      %   Suresh Joel Apr 15,2003
      %           Updated Feb 25, 2004

      switch(nargin),
      case 0
          disp('Too few arguements for voxel');
          return;
      case 1
          l=1;    %default length of side of voxel is 1
          c='b';  %default color of voxel is blue
      case 2,
          c='b';
      case 3,
          alpha=1;
      case 4,
          %do nothing
      otherwise
          disp('Too many arguements for voxel');
      end;

      x=[i(1)+[0 0 0 0 d(1) d(1) d(1) d(1)]; ...
              i(2)+[0 0 d(2) d(2) 0 0 d(2) d(2)]; ...
              i(3)+[0 d(3) 0 d(3) 0 d(3) 0 d(3)]]';
      for n=1:3,
          if n==3,
              x=sortrows(x,[n,1]);
          else
              x=sortrows(x,[n n+1]);
          end;
          temp=x(3,:);
          x(3,:)=x(4,:);
          x(4,:)=temp;
          h=patch(x(1:4,1),x(1:4,2),x(1:4,3),c);
          set(h,'FaceAlpha',alpha);
          temp=x(7,:);
          x(7,:)=x(8,:);
          x(8,:)=temp;
          h=patch(x(5:8,1),x(5:8,2),x(5:8,3),c);
          set(h,'FaceAlpha',alpha);
      end;

The problem is that matlab crashes if i try to plot a 128X128X128 voxel file (only 2Mb). when i say crashes i mean its simply stuck. nothing happens and it seems most of my computer resources are trying to plot this I am pretty sure its not my pc thats the problem (i7, 8Gb RAM, bench results: 0.2806 0.2195 0.1979 0.2687 0.6053 0.6086). anyone know how can i plot these voxels?

[Walter Roberson]

How much of the matrix is occupied? e.g., what is sum(VoxelMat(:)) ?

If everything was occupied you would be trying to draw two million patches, which would drag down MATLAB graphics for sure.

[Itzik Ben Shabat]

well its not sparse, pretty dense actually. so is there a way to overcome this? a different function maybe?

Answer 1

[Walter Roberson]

See the Faces and Vertices properties of patch() at http://www.mathworks.com/help/matlab/ref/patch_props.html

That is, instead of creating one patch per voxel, create one (or more) patches that has face and vertex information appropriate for a number of voxels. The voxels do not need to be adjacent to each other. You do not need to construct a left-face if the left voxel is set, not create a bottom face if the bottom voxel is set, not create a front face if the voxel in front is set.

It might be a bit much to create a single patch with all of the voxel information. One patch per slice might be better.

Are you sure you want to create all of this at once?? You are using alpha = 1 for all patches, so it is not possible to see through any face to any detail below, so usually for graphics purpose it is only worth creating the voxels that would be on the "outside" from the line of sight. You might want to set all of the "outside" at once so that you can rotate without recalculating. If your thought was to allow zooming in to the object to see interior detail, then instead of calculating everything at once, please consider using the zoom PreAction or PostAction callback to determine what the new coordinates will be, and calculate what will be visible from the new view and create the patches for that. It is perfectly acceptable for zoom (or pan) callbacks to change the graphics object that are visible, adding or removing, so as to provide only the relevant information.

[Itzik Ben Shabat]

1. thanks! your ideas it helps a lot

2. if i understand what you mean about creating less patches. lets say i have a full 128X128X128 voxel 3d binary matrix. this means that currently I am drawing 12582912 patches (6 per voxel- top,bottom,front,back,right,left, most of them overlapping). what did you mean by "one patch per slice" ? to how many patches does that decrease? as i see it i need at least 3 for all interior ones and a bit more for exterior ones (that sums up to half the number above)

3. regarding alpha and the line of sight. basically I wish to manipulate the data and then redraw it - I am creating different levels of resolution, havnt got to the programming of this yet). however, for graphic purposes i guess its fine if i just draw what is in the line of sight. is there a MATLAB function for this?

[Walter Roberson]

The voxel() code above creates two patches object per filled voxel; more than one face is drawn at a time.

One patch per slice would be drawing (for example) all the faces and vertices for one particular Z level, draw VoxelMat(:,:,1) as a single patch, draw VoxelMat(:,:,2) as another patch, and so on. If the time for each is practical, this would take you down to 128 patches (or possible 128 + 3 with the 3 being to close off the exposed ends.)

I am not sure if there is a MATLAB routine for calculating what is in the line of sight. Possibly not built-in. I think it likely that someone has built such a routine.

[Itzik Ben Shabat]

i think you are mistaken regarding voxel(). it plots 6 faces for each filled pixel (2 at a time - top&bottom, front&back, right&left). i dont understand what you meant regarding plotting only 1 patch per voxel (what if the voxel is in the corner of the volume ? i have to plot a minimum of 3 patches).

regardless your answer has pointed me in the right direction. i have already written a code to identify which voxels are external (no big improvement), now im writing a code to detect which faces are external and i will only plot them. if this still fails to work, i guess i cant avoid moving to c and opengl which should be faster. thanks for all the help.

[Walter Roberson]

You are right, voxel() does use 6 patches per voxel, even worse than the 2 I had noticed before.

When you patch(), even if you use the simple call which lists only vertices in order, you can create more than one visual area. You can insert NaN to indicate a discontinuous jump. So for example, you could have the x coordinates

[x000 x001 x001 x000 x000 nan x225 x226 x226 x225 x225]

to define two disconnected squares.

When you define a patch() object in terms of vertices and faces, you can work even more directly, defining a list of vertices and then using their index numbers to define the faces, any reasonable number of faces at a time. This should translate pretty directly into opengl calls in the graphics driver.

[Itzik Ben Shabat]

got it. i improved the code like you said (attached below for anyone who will ever search for this), calling patch() only once for each voxel(but still drawing 6 faces). this improved the speed tremendously! now im working on drawing only the external faces. Thanks for all the help! :)

function  [ExternalIndexes]=PlotVoxelinput()
clear all;
clc;
%PlotVoxelInput opens a RAW voxel file and plots it. 
%   

voxel_init_loc=[0,0,0];
voxel_loc=voxel_init_loc;
voxel_size=[1,1,1];

%file='Other inputs/nucleon.raw';
file='Lev input/Model 128/bone_model_3061.raw';
fileID=fopen(file);
[A, count] = fread(fileID);
m_hight=nthroot(count,3);
m_width=nthroot(count,3);
m_depth=nthroot(count,3);
%set figure parameters
xlim([0 m_hight]);
ylim([0 m_width]);
zlim([0 m_depth]);
view(45,45);
daspect([1 1 1])
%get 3d data
VoxelMat= multibandread(file,[m_hight,m_width,m_depth],'uint8',0,'bip','ieee-be');
%detect the external voxels
[ExternalIndexes,ExtenalNum]=FindExternalVoxels(VoxelMat);

%plot only external voxels
for i=1:ExtenalNum
    FV.vertices=[ExternalIndexes(i,1)+[0 voxel_size(1) voxel_size(1) 0 0 voxel_size(1) voxel_size(1) 0]; ...
       ExternalIndexes(i,2)+[0 0 0 0 voxel_size(2) voxel_size(2) voxel_size(2) voxel_size(2)]; ...
        ExternalIndexes(i,3)+[0 0 voxel_size(3) voxel_size(3) 0 0 voxel_size(3) voxel_size(3)]]';
    FV.faces=[1 2 3 4; 2 6 7 3 ; 6 5 8 7; 5 1 4 8; 4 3 7 8 ; 1 2 6 5];
   patch(FV,'FaceColor','r');
end
fclose(fileID);
end

[Mehmet Mutlu]

This single improvement is indeed much faster. Thanks a lot. Drawing only external ones should accelerate further. If someone wants to replace voxel.m from filexchange:

function voxel_improved(i,d,c,alpha)
%VOXEL function to draw a 3-D voxel in a 3-D plot
%Usage
%   voxel_improved(start,size,color,alpha);
%
%   will draw a voxel at 'start' of size 'size' of color 'color' and
%   transparency alpha (1 for opaque, 0 for transparent)
%   Default size is 1
%   Default color is blue
%   Default alpha value is 1
%
%   start is a three element vector [x,y,z]
%   size the a three element vector [dx,dy,dz]
%   color is a character string to specify color 
%       (type 'help plot' to see list of valid colors)
%   voxel_improved([2 3 4],[1 2 3],'r',0.7);
%   Based on Suresh Joel Apr 15,2003 and Itzik's answer

switch(nargin),
case 0
    disp('Too few arguements for voxel');
    return;
case 1
    l=1;    %default length of side of voxel is 1
    c='b';  %default color of voxel is blue
case 2,
    c='b';
case 3,
    alpha=1;
case 4,
    %do nothing
otherwise
    disp('Too many arguements for voxel');
end;

FV.vertices=[i(1)+[0 d(1) d(1) 0 0 d(1) d(1) 0]; ...
        i(2)+[0 0 0 0 d(2) d(2) d(2) d(2)]; ...
        i(3)+[0 0 d(3) d(3) 0 0 d(3) d(3)]]';
FV.faces=[1 2 3 4; 2 6 7 3 ; 6 5 8 7; 5 1 4 8; 4 3 7 8 ; 1 2 6 5];
patch(FV,'FaceColor',c,'FaceAlpha',alpha);

end

[Badreddine Alane]

I'm using the voxel_improved to skip the Matalb stuck while trying to draw the voxel form,but it's not working, i got this every time :''Too many output arguments''.

[Walter Roberson]

It looks to me as if the code is not initializing alpha if you pass in only one or two arguments, so be sure to pass in alpha.

[Badreddine Alane]

i don't think it's about initializing alpha, i have intialize it with different values but im always getting the same error

[Itzik Ben Shabat]

If I understand correctly, this error usually occure when you assign a functions output to more variables than it actually outputs. So, in this case, voxel_impoved does not return any value so make sure you did not assign it to anything.

Regardless, for plotting voxels I have this file exchange entry that may be usefull:

https://www.mathworks.com/matlabcentral/fileexchange/50802-voxelplotter

[Walter Roberson]

Please show the complete error message, showing which line the error is occuring on, and the code for that line.

Answer 2

[Badreddine Alane]

what im doing exactly is that i I have a matrix(i extracted it from 3d point cloud) with 3 columns (x, y and z coordinates) and 41472 lines. So, each line represents one point in a volume .i want to plot a voxel for this point cloud from this matrix, but i think it's not what i should give as an input for the voxel_impoved.(i have tried to draw a voxel for each point,and also tried it for the fulll matrix ) .it's first time for me working on the voxel form, hope someone will help .

[Itzik Ben Shabat]

OK. This is deffinitely thewrong input to the function.

You first have to do a voxelization process. This means that you need to convert the point cloud into a voxel grid. a voxel grid is a 3D matrix that subdivides the 3D space into occupied\non occupied small volumes called voxels. Usually the subdivisions are equal ( i.e. the matrix is of size n x n x n, here n is the number of voxels in each axis).

So, you basically need to check if a voxel contains a point (or a certain number of points) and generate this matrix before visualizing.

[Walter Roberson]

You can use discretize() on the X, Y, and Z coordinate separately. Then you can accumarray() to generate the volume. (If you were working in 2D you could sparse() instead.)

This presumes that each coordinate triple can be approximated by the cuboid it falls into, which should typically be acceptable for a point cloud.

[M.S. Khan]

i have used it but it tell that memory is out range

Answer 3

[Badreddine Alane]

so if i understand your answers, my first task is that i have to voxleize the point cloud ?; and if it's can you please guide me how can i voxleize a point cloud

[Walter Roberson]

I described the steps at https://www.mathworks.com/matlabcentral/answers/55633-how-to-drawing-3d-voxel-data#comment_704380

xbin = discretize(X, x_edges_go_here);
ybin = discretize(Y, y_edges_go_here);
zbin = discretize(Z, z_edges_go_here);
voxelized = accumarray([xbin(:), ybin(:), zbin(:)], 1);

[Walter Roberson]

x_edges_go_here might be something like

dx = 0.01;
minx = min(x);
maxx = max(x);
first_x = floor(minx / dx) * dx;
last_x = ceil(maxx / dx) * dx;
x_edges_go_here = first_x : dx : last_x;

[M.S. Khan]

Hi Walter. its good but tell me why you took dx =0.01. how if i took dx=0.0001

please guide me. i also want to use voxel to measure volume.

regards,

M.shafiq

Answer 4

[Badreddine Alane]

thank you so much for your help