Q:

The fastest way to "fill holes" in 3D?

I have a 3D binary image that I must fill all of the holes in the image.
"Fill holes 3D" in this case, is basically a 3D version of 
BW2 = imfill(BW,'holes') from MatLab, 
but I am not focusing on Matlab here. I am focusing on the algorithm in general.
Currently, I simply perform a "flood fill 3D" algorithm on all voxels at the border area of the 3D matrix. Anything not filled are holes so take this as a mask and remove all voxels on the real image with this mask.
Sure, the algorithm works, but somehow, I feel this is inefficient.
What is the fastest way to "fill holes" in 3D in this case?

A:

I don't have the privilege to comment, so I'm going write my comment as an answer.
I'm guessing your function, as well as imfill in Matlab, performs the "filling" iteratively.
Two basic improvements to speed these kind of iterative functions up would be to 1. increase the seed points, i.e. using other possible locations in addition to the border area, maybe base on image intensity before the binerization process, or whatever characteristics the images have, and 2. try to parallelize your function (see Dirk-Jan Kroon's Region Growing function for an example).  If this doesn't speed things up, or you just have too many seed points, then I would suggest you start implementing a function that uses your GPU for parallel processing.
Now there may be "non-iterative" functions that fill 3d holes... so I'll leave that to the experts here on dsp.se.  It would also be nice if you can show us a slice of your 3D image for a better understanding of what you are dealing with.

