How to fill the holes in 3D image reconstructed from binary images of CT slides

I have around 200 CT slides from patients and need to reconstruct them and segment out the eye orbit. I read the original DICOM image, use imbinarize(I,T) to change them to binary image, and then stack them up to form the 3D orbital image (see attached images below).

However, as you can see, there are many holes inside the orbit. To cope with this, I use bwmorph() with 'thicken' and 'bridge', then use imclose() with the 'disk' parameter:

for i = 30 : 160  
    info = dicominfo(dirOutput(i).name,'UseDictionaryVR',true);
    imgTemp = dicomread(info);
    imgCropped = imcrop(imgTemp,bboxCrop);
    imgCroppedBin = imbinarize(imgCropped,0.51728);
    img_thick = bwmorph(imgCroppedBin,'thicken');
    img_bridge = bwmorph(img_thick,'bridge',Inf);
    se = strel('disk',1);
    I(:,:,i) = imclose(img_bridge,se);
end;

I(:,:,i) is the resulting 3D array. Then I use fv = isosurface(I,ISOVALUE) and patch() to reconstruct the 3D array to a 3D image (as attached above).

The resulting 3D image still has many holes unfilled. Increasing the STREL element size or decreasing the ISOVALUE will make the final model thick and coarse (e.g. with stair case inside).

May I ask how can I fill the resulting holes in the 3D model and still preserve the original smoothness of the image? Should I do this in the 3D image or do this in 2D binary image before the reconstruction?

Thanks!

Answer 1

[Matt J]

It might help to look at some of these Examples with alphaShape. Basically, if you tune the alpha radius right, I think you might be able to plug the holes.

[Ringo Cheung]

Thanks for your solution. But alphaShapes need to know the exact coordinates of the holes to be suppressed.

And it might be impossible to do so.

On the other hand, with the experiment I did before, the 3D holes should come from the disconnected region of the 2D binary images before they are stacked up :

So that's why I use bwmorph() with 'thicken' and 'bridge' to connect them. But the result ends up unsatisfactory as said above.

Any idea to improve this algorithm? Thanks!

[Matt J]

Thanks for your solution. But alphaShapes need to know the exact coordinates of the holes to be suppressed.

They shouldn't. I've done similar things with them before.

However, in lieu of bwmorph, you could also try imclose on the whole volume (not slice-by-slice) with a 3-dimensional strel (a sphere, say, and not a disk).

[Ringo Cheung]

Hi Matt,

1) So you mean I input the (x,y,z) coordinates of my 3D volume to alphaShapes(), then it will auto compute the position of the holes and then fill them?

2) I have used imclose('sphere', size) both in 2D and 3D. The result is no good. It produces a coarse inner wall of the volume with a large value of "size" (e.g. 20) and also, it is very computationally intensive, requiring 10 to 20 min. to complete the hole-filling process. I used to think that gpuarray can help. However, it only supports flat STREL and so a "sphere" shape cannot take advantage of the GPU. Thanks!

[Sean de Wolski]

1) yes. Use find and ind2sub to build up your row/column/page indices and then call alphaShape on that.

[Matt J]

As for computation, I wonder if you have to be processing the whole volume. If you can compute approximately where the eye socket is within your volume, maybe extract a cropped sub-volume including only that area.

[Ringo Cheung]

Hi Sean,

My code to construct alphaShape :

[r,c,v] = find(I); 
x1 = c(:); 
y1 = r(:); 
z1 = v(:); 
P = [x1 y1 z1]; 
shp = alphaShape(P); 
figure; 
plot(shp);

But the output of plot has nothing. Just a 3D axis displayed.

[Ringo Cheung]

Code amended :

ind = find(I);
[x, y, z] = ind2sub(size(I), ind);
P = [x y z];
shp = alphaShape(P);
figure;
plot(shp);

Result plotted as attached.

Thanks!

[Ringo Cheung]

So, seems that the figure output was successful using alphaShape(). But any chance of exporting it to .STL file?

I have tried to study surf3stl(), patch2stl() and stlwrite(). But all 3 require the input to have faces and vertices coming from isosurface(), NOT plot().

Grateful if anyone of you can have any idea? Thanks!

[Matt J]

I think

   [bf,P] = boundaryFacets(shp)

will get you the necessary faces and vertices.

[Ringo Cheung]

Hi Matt,

The function is OK, and the output parameters can be passed to stlwrite() without problem.

To Matt and Sean,

Thanks for you two's great help, you really save my day!

[Sean de Wolski]

Looks good Ringo!

[Ringo Cheung]

Thanks, although there are still things to warp up like smoothing and cleaning up the model etc.

[Ringo Cheung]

One point :

Should I adjust the alpha radius or the HoleThreshold values so that I can fill the holes easier?

Thanks!

[Matt J]

I think ideally the radius should be just large enough that a sphere of that radius couldn't fit through any of the holes.

[Sean de Wolski]

You can use the alphaSpectrum to identify the largest radius below what you'd expect that to be.

[Ringo Cheung]

Dear both,

Got it. Noted with thanks!