Adaptive 3D grid for Surface Reconstruction

Question

I have implemented the algorithm described in http://physbam.stanford.edu/~fedkiw/papers/stanford2001-03.pdf to perform Surface Reconstruction from 3D unorganized data sets (Point Clouds). The main characteristic of this method is to represent a surface implicitly (aka volumetric representation) on a 3D grid and solve a partial differential equation to make the surface fit the data. The code I wrote (C++) is working fine and gives me nice results. But it's quite naive since I compute the solution iterating through each point of the 3D grid (so for a grid of 128 points per dimension I have more than 2 millions of point).

Since, I'm a complete newbie in this field, I'd be very glad if someone could give some reference for implementing adaptive grids (like octrees maybe) or whatever strategy to make more efficient the code!!!

Thanks.