ENH: indexing an image should return an image, not a numpy array

[ncullen93]

When we index an image, it gets converted to numpy. I actually don't believe this is really in the spirit of ants, as operations should generally return ants images.

img = ants.image_read(ants.get_data('r16'))

# numpy array, but maybe this should be a cropped image?
img2 = img[:20, :20] 

Should we change the functionality to be a cropped image instead? This would actually make cropping images much easier and natural than it is now.

The functionality comes from this method in the AntsImage class:

    def __getitem__(self, idx):
        if self._array is None:
            self._array = self.numpy()

        if isinstance(idx, AntsImage):
            if not image_physical_space_consistency(self, idx):
                raise ValueError('images do not occupy same physical space')
            return self._array.__getitem__(idx.numpy().astype('bool'))
        else:
            return self._array.__getitem__(idx)

[ntustison]

Yeah, I'd opt for returning an ANTs image. But I'm fine either way.

[ncullen93]

Thanks, I will work on changing this.

[cookpa]

If the image is cropped, you'll potentially need to update its origin.

What about slicing? If you do

im2 = im[:,:,64]

currently you would get a 2D numpy array. The same functionality is available with ants.slice_image, but the user has to specify additional options. So maybe it would be best to return a 3D image with a dimension of 1 along the slice axis if it's done by array indexing?

[ncullen93]

Yeah I was thinking to basically replicate slice_image with origin and everything else updated the same way. Indexing is just a more natural way to do this operation. You don't have to call some function called np.slice_array to slice a numpy array, for instance.

I agree slicing with one value is tricky... I would personally expect a 2D image to be returned there. But I'm not sure how the origin and such should be updated in this case.

Leaving this for reference on how to extract a region: https://examples.itk.org/src/core/common/imageregionoverlap/documentation

[cookpa]

I think following slice_image makes sense. It projects the origin into the plane, ie if you slice along z, the origin (X,Y,Z) becomes (X,Y). Subsetting the orientation axes is tricky. But when we transition from 3D to 2D we lose information anyway, the 2D image does not retain knowledge of where it is placed in 3D space.

If a crop and a slicing are applied simultaneously,

cropped_vol = timeseries[:,:,min:max,0] cropped_axial_slice = im[xmin:xmax,ymin:ymax,64] cropped_sagittal_slice = im[64,ymin:ymax,zmin:zmax]

I think the correct thing to do is slice then crop.