basically you want to "normalize" the terrain Z value to the ground surface.
You create a ground model using your ground points - for example a grid with cells - aka raster - holding an interpolated value of the ground Z value, and then subtracting all the points to that value. Just build a grid over your XY bounding box; below an example of a 90 X 100 grid:
int nRowCells = 100;
int nColCells = 90;
vector< vector <float> > grid;
for(int i=0; i < nRowCells; ++i)
{
std::vector<coord> row(nColCells, .0f);
grid.push_back( row );
}
Then you assign each of your ground point to a cell using its coordinate and cell resolution.
By the way, in LasTools this is possible to the "lasheight" module, if you do some Googleing you will read about how it works.