Category Archives: Techniques

write DICOM images from a 3D array

To write a DICOM image sequence from a 3D data array using MATLAB, you will need to specify ‘PixelSpacing’ and the ‘SliceThickness’, otherwise the 3D reconstruction software would have no idea to determine the length of the volume data. One way is to copy a sample metadata structure from another volume, update it with the new  ‘PixelSpacing’ and ‘SliceThickness’,  and pass it to dicomwrite. dicomwrite can write the correct ‘PixelSpacing’ and ‘SliceThickness’ into the dicom file, however, the reconstruction software (for example, slicer) would still treat ‘SliceThickness’ as 1 even it displays the correct ‘SliceThickness’ value. Another way is to directly specify these matadata values as the following. You just need to find the corresponding SOP Class UID for your dicom data from

% T defines the 'PixelSpacing' and the 'SliceThickness'

%% save dicom images
for i = 1 : n_slices

    status = dicomwrite(uint16(data(:,:,i)), ...
    file_name, ...
    'SliceThickness', T(3), ...
    'PixelSpacing', [T(1); T(2)], ...
    'SliceLocation', location_base - (i-1)*T(3), ...
    'ImageOrientationPatient', [1 0 0 0 1 0], ...
    'CreateMode', 'Copy', ...
    'SOPClassUID', '1.2.840.10008.');


Show the input data for the accumulating function of Matlab accumarray

When using MATLAB accumarray, if the subscripts in 'subs' are not sorted, the accumulating function ‘fun‘ should not depend on the order of the values in its input data. Otherwise you would have weird results. To show the input data for the accumulating function, one can define the function as ‘@(x){x}’. This would save the input data for you. For example,

>> val = 101:105;
>> subs = [1; 2; 4; 2; 4]
>> A = accumarray(subs, val, [], @(x){x})

A =

[       101]
[2x1 double]
[2x1 double]

Converting STL surface mesh to volume mesh using GMSH

Just three simple steps:
> File -> Merge -> “filename.stl”
> Geometry -> Elementary entities -> Add -> New -> Volume, click on the surface to build the volume
> Mesh -> 3D

Sample STL file can be downloaded @ (user name and password are ‘gmsh’ ‘gmsh’.)

If some errors happened, it is the problem of your input stl surface mesh. One common error is “self intersecting surface mesh, computing intersections (this could tak a while)”. If your surface is generated by some segmentation software and too complicated to fix. You may want to convert your surface mesh to a solid volume and then rebuild the surface mesh through isosurface reconstruction.

Or you can use meshlab to find those problematic faces and edges. Use the tools from Filters/Selection/Select Self Intersecting Faces (Select non Manifold Edges) to select these regions, delete them, and then use hole filling tool to generate the final surface mesh. This sometimes can solve some weird volume mesh problems, for example, volume mesh outside the surface mesh.


[10/27/2014] Instead of using GMSH, I actually used TetGen ( for several of my projects for volume mesh generation. It is a great tool for its simple syntax and lots of easy to understand examples. It was also integrated into Febio (, a FEM analysis software for biomechanical simulations. A paper of using both TetGen and Febio can be found here.

Kile on Windows

Installing Kile on Windows now is very simple, thanks to the effort of Windows KDE team. You just need to select Kile when you install Windows KDE packages. Detailed manual could be found at

Please remember to install a Latex distribution, i.e., miktex or texlive, before installing kile. I personally like texlive better, simply because dviout integrated in texlive loads dvi files much faster than YAP.


dviout tmpps1.pbm problem

Sometimes dviout maybe not responding and you have the message “gswin32.exe is making tmpps1.pbm” from the program status bar.  This is an old problem with dviout windows. As long as the package ‘hyperref’ is used in your tex file, you would have this problem.

tsearchn and pointLocation

To use pointLocation() to find the enclosing simplex (triangle/tetrahedron) for each query point location in Matlab, one has to supply this function with a ‘DelaunayTri’ structure as an input. This is a newly introduced class for performing delaunay triangulation. While in my case I need to reuse my saved delaunay triangulation results and assigned it to the Triangulation component of the ‘DelaunayTri’ structure. However I got this error,

??? Error using ==> subsasgn
Cannot assign values to the Triangulation.

I think this is because Triangulation is a private member of the class DelaunayTri, and every time one revises the values of X, this Triangulation  is updated automatically. If you are still allowed, you can use the old function ‘tsearchn’ to find the enclosing simplex of the delaunay triangulation. It is also very easy to write a simple function to compute the barycentric coordinate for each query point location. In case ‘tsearchn’ is not supported by the future Matlab, you can also write your own version of ‘tsearchn’. The brute force way is to check the barycentric coordinate of that point for each triangle. If one of the coordinates is negative, then the point is outside of that triangle.

Two video clips side by side using imagemagick and ffmpeg

Just copy the following code into a script file and run it at your bash. Please remmeber to “d2u your_file_name” otherwise you would have the error “syntax error near unexpected token $’\r’ “.


n=900 #image number

for(( i=0; i<n; i++))
montage -label ‘%f’ “xrCor$(printf “%04d” $i).png”  “xrSag$(printf “%04d” $i).png” -tile x1 -geometry +4+4 “montage$(printf “%04d” $i).bmp”

ffmpeg -b 7000k -r 10 -i montage%04d.bmp -b 7000k -r 10 -vcodec mpeg4 -vtag DIVX montage.avi