Interpolating Between Two Planes in 3d space

I'm developing a tool that lets you circle/enclose things on a 3d "volume." I want to save time by labelling "slices" 1 and 3, and "filling in" slice 2 from that information.

Two easy solutions are to:

1. slice2 = slice1 AND slice3 (gets the overlap between the two)
2. slice2 = slice2 OR  slice3 (true for any pixel true in either image)

These are ok and fast, but I would prefer to do something more intelligent by making the shape some sort of average/ interpolation between the two. You can imagine it as trying to find the cliff face that connects a plane at sea level and some plateau in the air.

Example: Fill in slices 2-4 from this 3d matrix. (Created using montage)

Feel free to come up with totally new ideas. I'm going to put my thoughts so far below.

Some stuff I've thought about that might help you, the answerer, but that I haven't been able to use successfully.

• You can do a bwperim on each image.
  
• You can try to "average" the images (or a weighted average).
  

The best I've got so far:

Add the images. Gives you overlap and two perimeters:
-an inner perimeter (inside of which will definitely be 1)
-and an outer perimeter (inside of which is questionable).
You can also mask the area that is >0 AND <2, which is a mask of this questionable area.
Run a bwdist on the two-perimeter image and mask that:

![masked bwdist image] 2

Not sure how to go from here, though. A line that took the "maximum" contour along that area would work, but I'm not sure how to do that robustly.

Any ideas on fixing my idea or any other ideas are welcome!

Thanks.

I figured this out today, after reading: "Efficient Semiautomatic Segmentation of 3D Objects in Medical Images" by Schenk, et. al.

here's the function I wrote:

function out = interp_shape(top,bottom,num)


if nargin<2;
    error('not enough args');
end
if nargin<3;
    num = 1;
end
if ~num>0 && round(num)== num; 
    error('number of slices to be interpolated must be integer >0');
end

top = signed_bwdist(top); % see local function below
bottom = signed_bwdist(bottom);

r = size(top,1);
c = size(top,2);
t = num+2;

[x y z] = ndgrid(1:r,1:c,[1 t]); % existing data
[xi yi zi] = ndgrid(1:r,1:c,1:t); % including new slice

out = interpn(x,y,z,cat(3,bottom,top),xi,yi,zi);
out = out(:,:,2:end-1)>=0;

function im = signed_bwdist(im)
im = -bwdist(bwperim(im)).*~im + bwdist(bwperim(im)).*im;