3D reconstruction from images

Introduction

We needed to take slices of brain or confocal microscope images and convert them back into 3D objects. Two examples are shown.

Typically, slices are processed so that there is a scalar 3D array representing the data set.The 3D array is then converted into a surface or set of surfaces using isosurface extraction techniques.

The programs

The process of building and showing the 3D object was broken into two parts to save memory and processing time. The first program converted the images into a Matlab 3D array. The second program builds the objects, combines them and displays them.

The first example is an insect brain. This work was performed for Nichole VanderSal at Cornell. The inputs wre TIFF files. The first Matlab program to read the files folllows. The TIFF files must be sequentially numbered. After this program is run, the "slice" matrix is manually stored as a Matlab *.mat file.

clf
clear all
%process a bunch of TIFF files into a volume
%the original images were had-colored by region to one of 7 colors
%so the scalar value on the output runs from from 0 to 7
%file setup
filebase = 'e:\Audax';
startFrame = 1;
endFrame = 74;
%read frames, reduce size, show frames, and build volume
for i=startFrame:endFrame
    filename=[filebase, num2str(i,'%2d'),'.tif']
    temp=double(imresize(imread(filename), 0.5));
    slice(:,:,i) = (temp(:,:,1)==255) + 2*(temp(:,:,2)==255) + 4*(temp(:,:,3)==255);
    imagesc(slice(:,:,i));
    colormap('gray')
    drawnow
end

The display program computes an isosurface for each of the hand-traced structures, smooths the resulting isosurfaces, colors them, then displays them as a 3D object.

clf
clear all

%load the file with the preprocessed images
load 'c:\my documents\matlab\nichole\slice.mat'
%the original images were had-colored by region to one of 7 colors

%patch smoothing factor
rfactor = 0.125; 
%isosurface size adjustment
level = .8;
%useful string constants
c2 = 'facecolor';
c1 = 'edgecolor';

%build one isosurface for each of 7 different levels
%The "slice" matrix takes on one of 7 integer values,
%so each of the following converts the data to a binary
%volume variable, then computes the isosurface between the
%1 and 0 regions

p=patch(isosurface(smooth3(slice==1),level));
reducepatch(p,rfactor)
set(p,c2,[1,0,0],c1,'none');

p=patch(isosurface(smooth3(slice==2),level));
reducepatch(p,rfactor)
set(p,c2,[0,1,0],c1,'none');

p=patch(isosurface(smooth3(slice==3),level));
reducepatch(p,rfactor)
set(p,c2,[1,1,0],c1,'none');

p=patch(isosurface(smooth3(slice==4),level));
reducepatch(p,rfactor)
set(p,c2,[0,0,1],c1,'none');

p=patch(isosurface(smooth3(slice==5),level));
reducepatch(p,rfactor)
set(p,c2,[1,0,1],c1,'none');

p=patch(isosurface(smooth3(slice==6),level));
reducepatch(p,rfactor)
set(p,c2,[0,1,1],c1,'none');

p=patch(isosurface(smooth3(slice==7),level));
reducepatch(p,rfactor)
set(p,c2,0.8*[1,1,1],c1,'none');

%lighting/image control
set(gca,'projection','perspective')
box on
lighting phong
light('position',[1,1,1])
light('position',[-1,-1,-1])
light('position',[-1, 1,-1])
light('position',[ 1,-1,-1])

%set relative proportions
daspect([1,1,1])
axis on
set(gca,'color',[1,1,1]*.6)

view(-30,30)

rotate3d on

The resulting image is shown below. An animation orbits the camera around the reconstruction.

The second example is a fluorescent neuromuscular junction. This work was performed for Kelly Ryan at Cornell. The inputs were TIFF files taken from a fluroscent, confocal microscope.

The first program:

clf
clear all

%file setup
filebase = 'C:\My Documents\BigData\KellyRyan\8-23-';
startFrame = 101;
endFrame = 122;

for i=startFrame:endFrame
    filename=[filebase, num2str(i,'%3d'),'.TIF']
    slice(:,:,i-100)=imread(filename);
    %image(slice(:,:,i-100));
    %colormap('gray')
end

The scond program computed an isosurface of intensity.

clf
clear all

%stacked data
load 'c:\my documents\bigdata\kellyryan\array.mat' slice

%build the surface of constant fluro intensity
p=patch(isosurface(slice,175));
reducepatch(p, .5)
set(p,'facecolor',[1,1,1],'edgecolor','none');
set(gca,'projection','perspective')
box on
light('position',[1,1,1])
light('position',[-1,-1,-1])

%scale in the z direcion to compensate
%for slice thickness
zscale=.5
daspect([1,1,zscale])
axis on
set(gca,'color',[1,1,1]*.8)
set(gca,'xlim',[0 250], 'ylim',[0 250])

az=0
el=0
view(az,el)

The NMJ shapes are hown below.