Sorry, to be clear I ways proposing two different ways of doing to same thing. The first is to use `numpy.where`. I like this method better because it leaves vol unchanged. This method will broadcast the mask against the volume, but probably not the way you'd like, it applies numpys broadcasting rules so if mas is (3, 2) and vol is (4, 3, 2), newvol will be like `array([mask * subvol for subvol in vol])`. Because probably want to pad mask after the last dimension we might want to do something like:

```
# Check that vol has more dims than mask
mask = mask.reshape(mask.shape + (vol.ndim - mask.ndim) * (1,))
newvol = npwhere(mask, vol, 0)
return newvol
```

If one wants to modify vol inplace we can simply use `*= mask`, but again we should decide if we want to add new dimensions at the begging of shape (this is the behavior of numpy's broadcasting) or if we want to add them at the end. If we add them at the end the code would look like this:

```
# Check that vol has more dims than mask
# Check/force mask to be of type bool
mask = mask.reshape(mask.shape + (vol.ndim - mask.ndim) * (1,))
vol *= mask
```

I hope that's more clear.