Read LiDar/LAS file in matlab.

Below is the c code for reading LiDar/LAS file into matlab. To compile the code, simply use
mex LASReader.cpp

The code reads the points coordinates only. Please see the code for details.

For details of LiDar file format, please visit this link.

De Boor’s Algorithm in C++

Here is the c++ code for the De Boor algorithm. It calculates a point C(x) on a B-Spline curve of any degree.

Point deBoor(int k,int degree, int i, double x, double* knots, Point *ctrlPoints)
{   // Please see wikipedia page for detail
	// note that the algorithm here kind of traverses in reverse order
	// comapred to that in the wikipedia page
	if( k == 0)
		return ctrlPoints[i];
	else
	{
		double alpha = (x-knots[i])/(knots[i+degree+1-k]-knots[i]);
		return (deBoor(k-1,degree, i-1, x, knots, ctrlPoints)*(1-alpha )+deBoor(k-1,degree, i, x, knots, ctrlPoints)*alpha );
	}
}

The Point class is defined as the follow:

class Point
{
public:
	Point(){x=0.;y=0.;z=0.;};
	// copy operator
	Point operator=(const Point pt) ;
	Point operator+(const Point pt) const;
	//Point operator-(const Point pt) const;
	Point operator*(double m) const;
	Point operator/(double m) const;
	double x,y,z;
};

Point Point::operator=(const Point pt)
{
	x = pt.x;
	y = pt.y;
	z = pt.z;
	return *this;
}
Point Point::operator+(const Point pt) const
{
	Point temp;
	temp.x = x + pt.x;
	temp.y = y + pt.y;
	temp.z = z + pt.z;
	return temp;
}
Point Point::operator*(double m) const
{
	Point temp;
	temp.x = x*m;
	temp.y = y*m;
	temp.z = z*m;
	return temp;
}
Point Point::operator/(double m) const
{
	Point temp;
	temp.x = x/m;
	temp.y = y/m;
	temp.z = z/m;
	return temp;
}

Include EPS graphics in pdflatex

It’s a shame that pdflatex supports jpg, png, pdf… graphic type but not EPS while latex supports EPS but not the others. From my previous post, there is an easy way to convert EPS to pdf.  But it is just annoying that every time you want to include an eps graphic you have to do the conversion manually. After googling, I found out that there is a way to let pdflatex do the conversion for you. Here is how.

Convert EPS to PDF: problem of bounding box

When you save a multiple plots/subfigures figure as an eps file in matlab, the bounding box might not be properly defined. People who use pdflatex need to convert EPS to PDF. The undefined bounding box will cause problems if you use epstopdf to convert EPS to PDF. This problem is shown in the first row in the picture below. The top-left figure is the eps file without properly defined bounding box. Because of this, epstopdf simply printed a letter size pdf file that only contains part of the figure as shown in the top-right figure. You can try to use the option –exact

epstopdf --exact bad.eps

It will probably generate an warning msg like “Warning: BoundingBox not found!” and the resulting pdf is still wrong. One way to fix this problem is to use epstool

epstool --copy --bbox bad.eps --output good.eps

The –bbox option tells epstool to generate a (correct) bounding box for the eps file. With a correctly defined bounding box, you can easily convert an eps file to a pdf file.

epstopdf --exact good.eps

The results are shown in the 2nd row of the picture below. epstopdf generated an pdf file in the same size as the eps file.

Create iPhone ringtone for FREE using iTune 8

iTune allow users to create their own ringtone from the music file downloaded via iTune store. It’s very convenient for users to create ringtones from their favorite songs  but it’s not free. This guide will show you how to convert mp3/acc files to iPhone for FREE simply using iTune 8.

Install flash player for firefox in ubuntu 8.04

1. Open firefox. Go to youtube.com and select any video. If you havn’t install flah, a message of “Hello, you either have JavaScript turned off or an ol version of Adobe’s Flash Player. Get the lastest Flash player.”. Click on the link and you will be redirected to adobe website to download it. Select .deb for Ubuntu x.xx.

2. Let the Package installer software install it. Restart firefox and it should be working.

If it’s still not working, follow the following steps.

3. Open Terminal and input
$> firefox &
to execute firefox.

4. Go to youtube.com and click on any video. It will be just like step 1 but this time you can see the error message from the Terminal. For me, the error messages say
LoadPlugin: failed to initialize shared library /usr/lib/adobe-flashplugin/libflashplayer.so [libnss3.so: cannot open shared object file: No such file or directory]
The problem is not in flash player but because of a dependency lib of flash player mssing .

5. Open Synaptic package manager (System >> Administration >> Synaptic package manager) and install the missing lib (libnss3). Restart firefox and the flash should work fine.

WebDAV client with encrypted SSL connection in ubuntu

WebDAV is a set of extensions to the HTTP protocol that allow users to edit and manage files on remote www servers. Here is the installation guide.

1. #> sudo apt-get install davfs2

2. create a folder
mkdir /media/webdav

3. Mount remote webdav server by
sudo mount -t davfs https://ip_address:port /media/webdav

The index of the largest element in a 2D array in matlab

[v,ind]=max(X);
[v1,ind1]=max(max(X));
disp(['The largest element in X is' num2str(v1) ' at (' num2str(ind(ind1)) ',' num2str(ind1) ')']);

Converting images to avi file in matlab

Here is the matlab code to convert a collection of images to a video file.

%// Will open an avi file name test.avi in local folder
aviobj = avifile('test.avi');
%// the quality of this video file
aviobj.Quality = 80;
%// compression method. See matlab manual for details.
aviobj.COMPRESSION ='None';%%color image

for i =1:numOfFrames
    %// apply image processing algorithms to the image here.
    %// image must be in size width x height x 3
    %//  in other words, color image.
    ..........
    ..........
    %// add image to the end of the avi file
    aviobj = addframe(aviobj,image);
end
%// close the file handle.
aviobj = close(aviobj);

Click here for detail instruction on addframe function.
Keywords: convert images jpb bmp to video avi compressed codec

Using rdesktop to remotely access window XP from ubuntu/linux

First, install rdesktop if you don’t already have it. rdesktop is a client software for RDP (remote desktop protocol) which is used in various microsoft OS’s.

$ sudo apt-get install rdesktop

Open a terminal and input

$ rdesktop -f -a 16 <ip/domain name>

where -f is for full screen mode and -a is the connection color depth. -a 16 means 16 bits color depth will be used for this connection. These are the two option I used. Other two useful options are -u <username> and -p <password>.

Press Ctrl-Alt-Enter to switch between full screen and window mode.

3D mesh correspondence selection tool

This is a program that allows users to manually assign corresponding triangles between two triangular meshes. Currently it only supports triangular-only meshs. To compile the code, run make PROG=hw3c.

User has to first click on a triangle from the model on the left window and select its corresponding triangle from the model on the right. The index of corresponding facets will be stored. Press ‘w’ to save the correspondences into file “cr.txt”. You will have to backup the old cr.txt yourself.

User can rotate each model individually by moving the mouse and hold the right mouse button and zoom in/out by using the mouse wheel. If the mouse wheel is not responding, use key ‘z’ and ‘x’ to zoom in and out. Press ‘Esc’ to exit the program.

Source code can be downloaded here.

Screeshot of the program

Displaying a sequence of images like video in matlab

images is an array of dim imageHeight x imageWidth x frameNumber

for i = 1:size(images,3)
    imshow(images(:,:,i));
    title(['frame# ' num2str(i)]);
    drawnow; %// this is why it works
end

Optical flow in matlab

Implementation based on Horn’s optical flow algorithm.

function [Vx,Vy] = OpticalFlow(images,alpha,iterations)
%// Calculating optical flow of a sequence of images.
%// images : 3D array that contains a sequence of images. size of images is (imageHeight, imageWidth, frame number)
%// alpha
%// iterations.
[height,width,frames]=size(images);

%//initialzation of u and v
Vx = zeros(height,width);
Vy = zeros(height,width);

for k = 1:frames-1

    % //initialization of Ex Ey and Et
    Ex = zeros(height-1,width-1,frames-1);
    Ey = zeros(height-1,width-1,frames-1);
    Et = zeros(height-1,width-1,frames-1);

    %//calculating Ex Ey and Et in frame k.
    for x = 2:width-1
        for y = 2:height-1
            Ex(y,x,k) = (images(y+1,x+1,k)-images(y+1,x,k)+images(y,x+1,k)...
                -images(y,x,k)+images(y+1,x+1,k+1)-images(y+1,x,k+1)...
                +images(y,x+1,k+1)-images(y,x,k+1))/4;

            Ey(y,x,k) = (images(y,x,k)-images(y+1,x,k)+images(y,x+1,k)...
                -images(y+1,x+1,k)+images(y,x,k+1)-images(y+1,x,k+1)...
                +images(y,x+1,k+1)-images(y+1,x+1,k+1))/4;

            Et(y,x,k) = (images(y+1,x,k+1)-images(y+1,x,k)+images(y,x,k+1)...
                -images(y,x,k)+images(y+1,x+1,k+1)-images(y+1,x+1,k)...
                +images(y,x+1,k+1)-images(y,x+1,k))/4;
        end
    end

    for nn = 1:iterations
        for x = 2:width-1
            for y = 2:height-1

                Vxbar = (Vx(y-1,x)+Vx(y,x+1)+Vx(y+1,x)+Vx(y,x-1))/6+...
                     (Vx(y-1,x-1)+Vx(y-1,x+1)+Vx(y+1,x+1)+Vx(y+1,x-1))/12;

                Vybar = (Vy(y-1,x)+Vy(y,x+1)+Vy(y+1,x)+Vy(y,x-1))/6+...
                    (Vy(y-1,x-1)+Vy(y-1,x+1)+Vy(y+1,x+1)+Vy(y+1,x-1))/12;

                %// chapter 12 of Horn's paper
                temp = (Ex(y,x,k)*Vxbar+Ey(y,x,k)*Vybar+Et(y,x,k))/(alpha^2 + Ex(y,x,k)^2 + Ey(y,x,k)^2);
                %// update u and v
                Vx(y,x) = Vxbar-Ex(y,x,k)*temp;
                Vy(y,x) = Vybar-Ey(y,x,k)*temp;
            end
        end
    end

end

Calculation of matrix inverse in C/C++

Below is the code to calculate matrix inverse of a matrix of arbitrary size (order) by using analytic solution. This method is known to be slow for very large matrix because of the recursion. However, I used this mainly for calculating inverse of 4×4 matrices and it worked just fine. You can also use CalcDeterminant and GetMinor to calculate determinant of a matrix.

MatrixInversion: the main function to calculate the inverse matrix.
CalcDeterminant: calculate the determinant of a matrix.
GetMinor: Get minor matrix.

Fast calculation of integral image in matlab

After implementing a mex-files for calculating integral image in matlab, a friend of mine told me this can be easily done with the cumsum function in matlab.

intImage = cumsum(cumsum(double(img)),2);

SOM (Self-organizing map) code in matlab

Below is the code of applying SOM on handwritten digits recongnition. It was implemented for a homework assignment in a course offered by professor Paul Gader. I used only ten handwritten digits images provided by Dr. Gader. Each image is resized to 30×30 and reshape to a 900×1 column vector. data in the code below is the training set. It’s a coolection of 900-dimensional column vectors.

Be aware of memory alignment of IplImage in OpenCV

The imageData of IplImage is aligned 4 or 8 bytes in order to enhance the processing speed. For example, let’s say you have a color image of size 98-by-98 (pixel depth 8 bits ) and imageData is aligned 4 bytes. The size of each row will not be 98×3=294 bytes but 296 bytes. The widthStep of IplImage shows the actual size of aligned image row in bytes. This can help us access imageData correctly. Here is a code that copy imageData to a user-created buffer.

//
IplImage *colorImage;
// load an color image
colorImage = cvLoadImage("LM_0273.jpg");
// create a temp buffer
unsigned char *buffer,*temp2;
buffer = new unsigned char[colorImage->width*colorImage->height*colorImage->nChannels];
temp2 = buffer;

// pointer to imageData
unsigned char *temp1 = (unsigned char*) colorImage->imageData;
// copy imagedata to buffer row by row
for(int i=0;i<colorImage->height;i++)
{
    // memory copy
    memcpy(temp2, temp1, colorImage->width*colorImage->nChannels);
    // imageData jump to next line
    temp2 = temp2 + colorImage->widthStep;
    // buffer jump to next line
    temp1 = temp1+ colorImage->width*colorImage->nChannels;
}
// ......................
//.....

Calculating principal components of data with VERY HIGH dimensionality

Suppose we have a set of data X with dimensionality D and total number of data is N. You might encounter insufficient memory problem when D is a large number or simply waste too much time calculating those redundant components when D >> N. Using singular value decomposition can prevent those problems.

From singular value decomposition, we know that a D-by-N matrix X can be decomposed as , where U is composed of the eigenvectors of and of the size D-by-D and V is composed of the eigenvectors of and of size N-by-N. Because is an identity matrix(column vectors in V are orthonormal), the first equation can be written as . From the question, the first m (m<=N) columns of matrix U are actually the principle components of X. We can actually get the first effect m columns of U without actually calculating . Here is how it works.

Histograms of two sets of data with different color in matlab

figure;
hist(data1);
hold on;
%//make data1 red
%//get the handle of the bars in a histogram
h = findobj(gca,'Type','patch');
%//color of the bar is red and the color of the border
%// of the bar is white!
set(h,'FaceColor','r','EdgeColor','w');
%//data 2 use default color!
hist(data2);

Edit: 03/30/09

As suggested by Adam in the comments, there is a better way to achieve that. Thanks Adam!!
_______________________________

The variable h contains handles to each of the histograms. Just use “set” on each element individually.

e.g.:

hist(data1);
 hold on;
 hist(data2);
 hist(data3);

h = findobj(gca,’Type’,’patch’);
display(h)

set(h(1),’FaceColor’,’r’,’EdgeColor’,’k’);
set(h(2),’FaceColor’,’g’,’EdgeColor’,’k’);
set(h(2),’FaceColor’,’b’,’EdgeColor’,’k’);

Load a sequence of images into matlab

I have implemented the following script probably every time I need to load a bunch of images into matlab. This time, I am going to save it for the future:)