create 3D image from coordinates and intensity values

I am trying to create a 3D image of size 1000x1000x1000 voxels with all the voxels being zero and then assign a random value in the 2000 to 2001 range instead of 0 to some specific elements in the array and finally store it as a binary file.

The array named "coord" is the Nx3 matrix coordinates (x,y,z) of the points that I need them to be assigned the random value in the 3D array.))

I should mention that all the x,y,z values of the coordinate matrix are floating point numbers with: 0<=x<=1000 0<=y<=1000 0<=z<=1000

My aim is to export the 3D matrix in a binary format (other than MATLAB's default binary format) so that I can use it with other programs. Here is what I've been up to so far:

load coord;
a=coord(:,1);
b=coord(:,2);
c=coord(:,3);
d=rand(1000,1)*2000;
dd = 0:2:1000;
[xq,yq,zq] = meshgrid(dd,dd,dd);
vq = griddata3(a,b,c,d,xq,yq,zq,'nearest');
h=figure;
plot3(a,b,c,'ro')
%=========================================%
fid=fopen('data.bin','w');
fwrite(fid,vq,'single');
fclose(fid);

In the above code a, b and c are the x, y and z coordinates of each point and d is the corresponding intensity values for the desired range. While it is possible to create a 3D mesh (using meshgrid) and then interpolate the intensity values for mesh points (using griddata3), the final result (vq) would not be the actual points (ai,bi,ci) and corresponding intensities , but rather an interpolated set of intensities assigned to the mesh points which is pretty useful for visualization purposes (for instance if you like to fit a 3D surface which fits through actual data).

I am simply trying to find a way to store the actual data-points and their intensities into a file and export it. Any help is highly appreciated.

[Walter Roberson]

If you want the actual point locations and intensities to be stored, you need to define the order, and define what should happen if multiple points end up in the same cuboid.

With the information you have (not) given,

sortrows([a,b,c,d])

would seem to be a plausible output.

[Mehdi]

Hi Walter. I am not quite sure if I understand your point. Could you explain further on "multiple point on the same cuboid" ?

Each point has a unique (xi,yi,zi) coordinate in the 3D space and I guess this is sufficient info to locate it. There are no duplicate points in the data either.

I am not quite sure about the ordering... Could you give an example so that I can understand it .

Answer 1

[Walter Roberson]

When you use griddata3() or anything similar, you are interpolating positions into a 3D grid, and the coarseness of that grid is determined by your mesh, not by how close the points are. As a result, if two points happen to be sufficiently close together, they end up gridded into the same cuboid even though they are distinct coordinates.

You are not necessarily required to use a grid, but even if not then your other program might expect some particular sequence of data, such as that the data be sorted first on x, then on y within x, and then on z within x and y. If there are restrictions such as that, then you need to build them into the way you order the data.

If the other program does not care what spacial order the points are in, then you need to know what binary data format the other program expects.

Here is one plausible output routine:

load coord;

d = rand(size(coords,1),1)*2000;

outdata = [coords(:,1:3), d(:)] .';  %transposed!

fid = fopen('data.bin','w');
fwrite(fid, outdata, 'single');
fclose(fid);

This writes the data in binary, first x, first y, first z, first d, second x, second y, second z, second d, and so on. The array length is self-defining by the number of entries you read.

Another plausible routine is

   load coord;

d = rand(size(coords,1),1)*2000;

outdata = [coords(:,1:3), d(:)];  %not transposed!

fid = fopen('data.bin','w');
fwrite(fid, outdata, 'single');
fclose(fid);

This writes the data in binary, all x, all y, all z, all d. It is not obvious at the time of reading what each value corresponds to. On the other hand, it is easy to deal with in MATLAB:

resize( fread(inputfid, 'single=>double'), [], 4)

A problem with both of these two files is that they are not voxelized. As soon as your voxelize, you need to deal with the problem of multiple values falling into the same voxel.

[Mehdi]

@Walter Thanks a lot. I am more clear on this now. However, with regards to griddata3(), I reckon it does not interpolate the grid-point coordinates but their assigned values based on those of the actual data using one of the interpolation methods (linear, nearest neighbor etc).

Could you hint on how I can voxelate the created data? I guess I may end up rounding up the coordinate x,y,z as in the process of voxelating the binary file, this should happen (?)

[Walter Roberson]

I do not recall how griddata3() itself works. It is recommended that griddata3() be replaced with methods in the TriScatteredInterp classes, http://www.mathworks.com/help/techdoc/ref/triscatteredinterp.html

Roughly speaking, the approach there is to do a Delaunay triangulation of the coordinates, and use that to create an interpolating function, which you then apply your voxel grid in order to determine the voxel values. With such a process, if multiple input points happen to influence the same point being interpolated (e.g., because the points were closer together than the gridding distance) then the multiple points would have an influence.

Note that for both griddata3() and the TriScatteredInterp approaches, if a particular voxel does not happen to contain a point, then its value is interpolated from the nearby actual points that do exist. I'm not sure that is what you want.

Now, if you could be sure that you did not have multiple input points that end up in the same voxel, I could suggest an output format:

Create your voxel grid not of intensities, but rather of the row index of the contained point.

I'd have to think a bit about the most efficient way of doing that. I don't think it should be too bad.

[Mehdi]

Thanks Walter. As you mentioned there is quite a few issues to be taken care of before getting the voxelation right.

I'd appreciate it if you could suggest a method that voxelates the binary file. I'm not familiar with this part.

[Walter Roberson]

Please expand on what you mean by "voxelates the binary file" ?

[Mehdi]

Referring to your first method above, I end up with a binary image which needs to be voxelized. Could you explain on how this can be achieved? is there particular functions in MATLAB to do this.

[Walter Roberson]

I'm not sure which possibility of mine you are referring to as "first" ? There is more than one meaning of "binary" and "image" involved in this discussion.

Answer 2

[Image Analyst]

I'm wondering if you're getting confused by thinking that what you get with griddata3 (now superseded by TriScatteredInterp) is a 3D volumetric image, when it's not. They give you a surface, which is kind of like a 2.5 D data set. You have an x, and a y, and for each of those two dimensions, there is a value z. The value z depends on the x and y coordinate and is not an independent dimension. It's actually a 2D dataset but being displayed as a surface. Just because the surface looks 3D does not mean that the data being plotted is 3D. A true volumetric data set, like you'd get from MRI or CT, has full independent x, y, and z coordinates and a value. I don't think that's what you're after is you're trying to use griddata3.

If you did want a true 3D dataset, you'd do this:

load coord;
width = 1000;
x=coord(:,1);
y=coord(:,2);
z=coord(:,3);
numberOfCoords = length(x);
array3D=rand(width ,width ,width );
for xk= 1 : numberOfCoords
   thisX = round(x(xk));
   for yk= 1 : numberOfCoords 
     thisY = round(y(yk));
     for zk= 1 : numberOfCoords 
          thisZ = round(z(zk));
          array3D(thisX, thisY, thisZ) = 2000 + rand(1);
      end
    end
end

But again, I don't think you want that, or maybe you do - I'm not sure. Or maybe you want it both ways?

[Mehdi]

@Image Analyst. Thanks a lot for your helpful info. I am actually after creating a 3D (CT SCAN type) image from a set of 3D coordinates and intensity values. I have already written something with a similar approach to your code above to create the image. However, If I am right, you have also rounded off the coordinate values to make them integer (hence losing the floating point precision).

I guess it is not possible to incorporate the coordinates in floating mode in a tomogram. Please correct me if you think that's not the case. The reason I think so is that in a voxelized image, usually each voxel will have an integer set of x,y,z, coordinate values.

[Walter Roberson]

griddata3() _does_ fit values to three independent coordinates. griddata() is the routine that you are thinking of.

TriScatteredInterp is perhaps more commonly used over a planar grid, but it is fully capable of working in 3 space.

[Walter Roberson]

CT Scans always involve regular coordinate intervals (though not necessarily the same spacing in each direction.) The coordinates do not need to be integers, but they do need to be equal intervals in any one direction.

Often for something like a CT Scan or an MRI, one would use a DICOM file, with the DICOM metadata set to indicate the X-Y-Z coordinate origin, the pixel size (possibly non-square), and the distance between slices.

Writing out something like that is quite different than trying to keep the full original coordinate locations of each point.