Hi Donald,

Given that you want to deal with two overlapping spheres, it is a natural
jump to want to handle any small number of overlapping spheres. The
trouble is your '-expr' will become much more complex, but the answer is
yes.

Assuming you have a dataset filled with N (possibly truncated) spheres,
having respective values of 1..N, and you want to insert sphere N+1.

You could set the new value to N+1 where the previous value is zero, and
then figure out how to fight it out for the other locations. The easy
case would be akin to just doing two spheres, where you know which sphere
you are fighting with. But if you don't know, you will have to test them
all.

In general, the goal at every voxel within the new radius is to set its
value to that of the closest center. Which means you either iterate over
all center for any newly assigned one, or you do it all at once. Let me
just present a do-it-all-at-once option.

For each voxel set the value to argmax(list of negative distances from
the centers) * step(is voxel in any shpere). So for 3 spheres with centers
at cx1, cy1, cz1, ... and squared radii rs1, rs2, rs3:

3dcalc -a zero_dset+orig -expr "argmax(-((x-$cx1)*(x-$cx1)+(y-$cy1)*(y-$cy1)+(z-$cz1)*(z-$cz1)),-((x-$cx2)*(x-$cx2)+(y-$cy2)*(y-$cy2)+(z-$cz2)*(z-$cz2)),-((x-$cx3)*(x-$cx3)+(y-$cy3)*(y-$cy3)+(z-$cz3)*(z-$cz3)))*step(step((x-$cx1)*(x-$cx1)+(y-$cy1)*(y-$cy1)+(z-$cz1)*(z-$cz1)-$rs1)+step((x-$cx2)*(x-$cx2)+(y-$cy2)*(y-$cy2)+(z-$cz2)*(z-$cz2)-$rs2)+step((x-$cx3)*(x-$cx3)+(y-$cy3)*(y-$cy3)+(z-$cz3)*(z-$cz3)-$rs3))" -prefix spheres3


Note that there is no argmin function. So argmax is taken of the negative
distances from the centers. So that part assigns a value of 1..N equal to
the index of the closest sphere. Then you multiply the results by the
step of: the sum of steps answering "is the voxel in sphere i?".

Therefore argmax() will be in 1..N, and step( sum of steps ) will be 1
if the voxel is in any sphere, and 0 if it is not. So their product will
be the index of the closest sphere if it is inside any sphere.

You can expand this to any number of spheres, as long as you don't run
out of 'expr' space, in which case you will have to iterate the steps.

Hmmmm, that's probably the most exciting 3dcalc expression I've ever
written...

And thanks to Bob for adding that argmax() function a couple of months
ago.

- rick