ndimage.shift results in incorrect values for inplace for float arrays

[JoOkuma]

Description

When using in-place operations for float arrays, the image.shift returns incorrect results.
I noticed this bug when getting incorrect results and not an empty result as my example, but that's how I reproduced in a minimal example. Because it wasn't happening every time.

To Reproduce

import cupy  as cp
import napari
from cupyx.scipy import ndimage as ndi

def main() -> None:

    # Define the volume size
    volume_shape = (512, 512, 512)
    size = 500
    center = 256

    # integer works
    arr = cp.zeros(volume_shape, dtype=cp.float32)
    arr[tuple(slice(center - size // 2, center + size // 2) for _ in range(3))] = 1

    # Apply a shift using scipy.ndimage.shift
    shift_vector = (5, -68, -130)  # Shift in (z, y, x) order

    print("Before", arr.min(), arr.max())

    # not in-place works
    # arr = ndi.shift(
    #     arr,
    #     shift_vector,
    #     order=1,
    # )
    ndi.shift(
        arr,
        shift_vector,
        order=1,
        output=arr,
    )

    print("After", arr.min(), arr.max())

    napari.view_image(arr.get())
    napari.run()


if __name__ == "__main__":
    main()

Installation

Conda-Forge (conda install ...)

Environment

OS                           : Linux-5.15.0-92-generic-x86_64-with-glibc2.31
Python Version               : 3.10.12
CuPy Version                 : 12.2.0
CuPy Platform                : NVIDIA CUDA
NumPy Version                : 1.24.4
SciPy Version                : 1.11.2
Cython Build Version         : 0.29.36
Cython Runtime Version       : None
CUDA Root                    : /home/jordao/miniconda3/envs/imgproc
nvcc PATH                    : None
CUDA Build Version           : 11080
CUDA Driver Version          : 12000
CUDA Runtime Version         : 11080
cuBLAS Version               : (available)
cuFFT Version                : 10900
cuRAND Version               : 10300
cuSOLVER Version             : (11, 4, 1)
cuSPARSE Version             : (available)
NVRTC Version                : (11, 8)
Thrust Version               : 101501
CUB Build Version            : 101501
Jitify Build Version         : 11d783c
cuDNN Build Version          : 8800
cuDNN Version                : 8800
NCCL Build Version           : 21903
NCCL Runtime Version         : 21903
cuTENSOR Version             : 10700
cuSPARSELt Build Version     : None
Device 0 Name                : NVIDIA GeForce RTX 3090
Device 0 Compute Capability  : 86
Device 0 PCI Bus ID          : 0000:03:00.0
Device 1 Name                : NVIDIA GeForce RTX 3090
Device 1 Compute Capability  : 86
Device 1 PCI Bus ID          : 0000:04:00.0

Additional Information

No response

[kmaehashi]

Thanks for reporting @JoOkuma! The shift operation (and many ndimage features) rely on raw CUDA kernel generated on-the-fly. I haven't fully checked the code but it looks it does not support cases s.t. input and output memory overlap.

Supporting in-place will need almost full rewrite of the kernel so I'm thinking of adding a check and making a copy automatically if there is an overlap.

cupy/cupyx/scipy/ndimage/_interpolation.py

Lines 640 to 655 in 8fb6d04

	output = _util._get_output(output, input)
	if input.dtype.kind in 'iu':
	input = input.astype(cupy.float32)
	filtered, nprepad = _filter_input(input, prefilter, mode, cval, order)
	integer_output = output.dtype.kind in 'iu'
	_util._check_cval(mode, cval, integer_output)
	large_int = _prod(input.shape) > 1 << 31
	kern = _interp_kernels._get_shift_kernel(
	input.ndim, large_int, input.shape, mode, cval=cval, order=order,
	integer_output=integer_output, nprepad=nprepad)
	shift = cupy.asarray(shift, dtype=cupy.float64, order='C')
	if shift.ndim != 1:
	raise ValueError('shift must be 1d')
	if shift.size != filtered.ndim:
	raise ValueError('len(shift) must equal input.ndim')
	kern(filtered, shift, output)

[JoOkuma]

Making a copy and updating the docs to say it's not in place would be excellent!
Thanks for the quick reply @kmaehashi