Dutch national flag problem – matlab

Dutch national flag problem is a classic programming problem. It groups the input array into three sub-groups in just one pass. For detailed information, please refer to http://en.wikipedia.org/wiki/Dutch_national_flag_problem

function x = category_sort(x, low, high)

len = length(x);
p = 1; q = len;
i = 1;

while i <= q
    if x(i) < low     % move the current element to the front
        tmp = x(p);
        x(p) = x(i);
        x(i) = tmp;

        p = p + 1;  
        i = i + 1;     
    elseif x(i) >= high     % move the current element to the back
        tmp = x(q);
        x(q) = x(i);
        x(i) = tmp;
        q = q - 1;
    else    % move to the next element in the array
        i = i + 1;

We can test the code by simply running the following script:

x = 3 * rand(20, 1);
x = category_sort(x, 1, 2)

Implement strtok in C

To implement the C function ‘strtok’, one doesn’t need to allocate extra memory for the input string before modifying it, since the caller function is supposed to make copy of the input string and deallocate the original input string. The function only needs to define static variables to store the starting of the string.

The caller function provides a string as the input during the first call, and ‘null’ for the following function calls. For details about how this function works, please refer to http://www.cplusplus.com/reference/cstring/strtok/.

char* sp = NULL; /* the start position of the string */

char* strtok1(char* str, const char* delimiters) {

	int i = 0;
	int len = strlen(delimiters);

	/* check in the delimiters */
	if(len == 0)
		printf("delimiters are empty\n");

	/* if the original string has nothing left */
	if(!str && !sp)
		return NULL;

	/* initialize the sp during the first call */
	if(str && !sp)
		sp = str;

	/* find the start of the substring, skip delimiters */
	char* p_start = sp;
	while(true) {
		for(i = 0; i < len; i ++) {
			if(*p_start == delimiters[i]) {
				p_start ++;

		if(i == len) {
		       sp = p_start;

	/* return NULL if nothing left */
	if(*sp == '\0') {
		sp = NULL;
		return sp;

	/* find the end of the substring, and
        replace the delimiter with null */
	while(*sp != '\0') {
		for(i = 0; i < len; i ++) {
			if(*sp == delimiters[i]) {
				*sp = '\0';

		sp ++;
		if (i < len)

	return p_start;

This following is the test code from cplusplus.com.

#include <stdio.h>
#include <string.h>

int main ()
  char str[] ="- This, a sample string.";
  char * pch;
  printf ("Splitting string \"%s\" into tokens:\n",str);
  pch = strtok1(str," ,.-");
  while (pch != NULL)
    printf ("%s\n",pch);
    pch = strtok1(NULL, " ,.-");
  return 0;

As pointed out by Ian in the following comments [9.25.2014], there is a much simpler solution to this problem. Feel free to test it and let us know if it has any bugs.

char* strtok1(char *str, const char* delim) {
	static char* _buffer;
	if(str != NULL) _buffer = str;
	if(_buffer[0] == '\0') return NULL;

	char *ret = _buffer, *b;
	const char *d;

	for(b = _buffer; *b !='\0'; b++) {
		for(d = delim; *d != '\0'; d++) {
			if(*b == *d) {
				*b = '\0';
				_buffer = b+1;

				// skip the beginning delimiters
				if(b == ret) { 
				return ret;

	return ret;

Kile can not edit tex file

Sometimes when I use Kile to open tex files, kile just simply cannot edit them. I can move the cursor freely but cannot type in anything. In this case, you just need to uncheck Tools -> Read Only Mode.

If the read-only problem is caused by ‘The file xxx.tex was opened and contained too long lines (more than 1 024 characters). Too long lines were wrapped and the document is set to read-only mode, as saving will modify its content.’, you may walk around this by setting the line lenght to 0, under Settings -> Configure Kile -> Open/Save -> Line Lenght Limit.

The wield case I met is that my kile didn’t have menu bar, although kate has that. I searched online and found out that you can re-enable the menu bar by editing the kilerc under ~/.kde/share/config/. Changing ‘menubar=Disabled’ to ‘Enabled’ would make the menu bar re-appear.

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 http://www.apteryx.fr/dicom/dicom_conformance

% 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 @ https://geuz.org/trac/gmsh/attachment/wiki/STLRemeshing/pelvis.stl.gz (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 (http://wias-berlin.de/software/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 (http://mrl.sci.utah.edu/software/febio/), a FEM analysis software for biomechanical simulations. A paper of using both TetGen and Febio can be found here.

install Fedora 16 from Ubuntu, network installation

I just switched to Fedora from Ubuntu. This time the installation is fairly easy. First, I copied vmlinuz and initrd.img from the isolinux folder of the Fedora installation CD to /boot, and then I rebooted the system to enter the grub command line. And typed:

grub> linux /vmlinuz
grub> initrd /initrd.img
grub> boot

then the Fedora installation begun. Later on you just need to choose network interface for installation. The whole installation lasted less than one hour in my case. One similar topic from my previous posts can be found at link .


You might meet the login problems, “Oh! No! Something has gone wrong. A problem has occured and the system can’t recover. Please contact a system adminstrator”, which could be solved by:

yum remove selinux-policy selinux-policy-targeted
yum update libsepol
yum --nogpgcheck install libsepol http://kojipkgs.fedoraproject.org/
      selinux-policy-3.10.0-55.fc16.noarch.rpm --enablerepo=u*g
yum --nogpgcheck install http://kojipkgs.fedoraproject.org/packages/
      selinux-policy-targeted-3.10.0-55.fc16.noarch.rpm --enablerepo=u*g


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 http://sourceforge.net/apps/mediawiki/kile/index.php?title=KileOnWindows

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.