get / set points corresponding to a surface within a 3D volume

[Tasos Papastylianou]

Let V be a volume (e.g. a 3D MRI dataset) with size(v) == [M N P]
Let X and Y form a 2D grid, where size(X) == size(Y) == [M N]
Let Z be such that size(Z) = [M N] and Z defines a surface over the grid X,Y, which is within the bounds of V.

Is there an 'indexing' operation IND with respect to the points on the surface Z such that I could do the conceptual equivalent of this?
>> U = V(IND);
and even
>> V(IND) = V(IND) .^ 2;

A 'hackish' way of achieving the first is as follows:

Given, say:
V = repmat(permute([1:100],[1 3 2]),[100 100 1]) .* repmat(magic(100),[1 1 100]) .* repmat(magic(100)',[1 1 100]);
[X Y] = ndgrid(1:100,1:100);
Z = X/2 + Y/2;

We can do: 
h = slice(V,X,Y,Z);
U = get(h,'CData');

But I do not consider this an 'indexing' operation. Nor does it allow me to 'set' that volume to something, 
only to 'get' the values. (right?)

I'd be grateful for any insight on this :)

Tas.