One easy solution would be to use a dilation morphological operator. Dilate your image by one voxel, then substract the original image from your dilated image. The result is the boundary region. You can then multiply by 0.5 and sum with your original image to have the boundary region in gray:

ana = analyze75read('test.hdr');

dil = imdilate(ana, ones(3,3,3)); 

boundary = dil-ana;

output = 0.5*boundary + ana; %may need to convert to float, since it'll round up if you're using ints 

Matlab also has a tool for 2D images called bwmorph. If a slice-by-slice approach is acceptable, you can run:

ana = analyze75read('test.hdr');

slice = ana(:,:,30);

bwmorph(slice, 'remove')

And then just repeat per slice.

Are you sure you can't do that with bwboundaries???? if you call bwboundaries with a binary image it will return an Nx1 cell array where N represents the amount of objects  in the segmented image...and for each N you'll have the x and y coordinates for each pixel in the object boundary.....so if for example you do something like this:

B = bwboundaries(BW);

for k = 1:length(B)

       boundary = B{k};

       plot(boundary(:,2), boundary(:,1), 'w', 'LineWidth', 2)

       ////////////you are ploting the pixel in the boundary of object B{k}

      ////////////boundary(:,2) = X axis and boundary(:,1) = Y axis

end