Thread Subject:

Efficient way for Creating 3D Matrix? (Cubes in a sphere)

Hi All:

I want to create a sphere filled by voxels (unit cube). I'v written the below small code.
e = 20;
r = 1;
t = linspace(-r,+r,e);
for i = 1:e
    for j = 1:e
        for k= 1:e
            if (t(i))^2 + (t(j))^2 + (t(k))^2 <= r^2
                lp = length(pore_voxel);
                sphere_voxel(1+lp,:) = [t(i) t(j) t(k)];
            end
        end
    end
end
Where t(i), t(j) and t(k) are the center coordinate of the cubes and r is the radius of sphere.

Questions:
1) Is there any other way to get rid of these "for loops" ?
2) Is there any way to create a 3D matrix ? Because, The "parfor" can't be used in the above code as long as "lp" and "sphere_voxel" matrices are changing in each loop.
3) What is your recommendation to speed up this code ? (it is about 2 or 3 million cubes inside the sphere).
It should be noted that I simplified the code and "knnsearch" is used within the code which drastically decreases the code speed. That's why I'm curious to find an efficient way.

Sorry in above code "sphere_voxel(l+lp,:) = [t(i) t(j) t(k)];" should be replaced by sphere_voxel(e+lp,:) = [t(i) t(j) t(k)];

I dont know why I did many mistakes to write this little code. Here hopefully is the correct one. Pl help me with that?

sphere_voxel = [];
e = 20;
r = 1;
t = linspace(-r,+r,e);
for i = 1:e
    for j = 1:e
        for k= 1:e
            if (t(i))^2 + (t(j))^2 + (t(k))^2 <= r^2
                lp = length(sphere_voxel);
                sphere_voxel(e+lp,:) = [t(i) t(j) t(k)];
            end
        end
    end
end

 

"mahdi roozbhani" <m.m.roozbahani@gmail.com> wrote in message <jnhs4e$lg8$1@newscl01ah.mathworks.com>...
> I dont know why I did many mistakes to write this little code. Here hopefully is the correct one. Pl help me with that?
>
> sphere_voxel = [];
> e = 20;
> r = 1;
> t = linspace(-r,+r,e);
> for i = 1:e
> for j = 1:e
> for k= 1:e
> if (t(i))^2 + (t(j))^2 + (t(k))^2 <= r^2
> lp = length(sphere_voxel);
> sphere_voxel(e+lp,:) = [t(i) t(j) t(k)];
> end
> end
> end
> end
- - - - - - - - - - - -
 [TK,TJ,TI] = ndgrid(linspace(-r,+r,e));
 S = TI.^2+TJ.^2+TK.^2<=r^2;
 sphere_voxel = [TI(S),TJ(S),TK(S)];

Roger Stafford