Repository: scipy/scipy
Pull request: BUG: segfault in ndimage.fourier_ellipsoid with length-1 dims
#### Reference issue
Closes gh-12763

#### What does this implement/fix?
Replaces a segmentation fault with a `ValueError`.

The three python functions `ndimage.fourier_ellipsoid`, `fourier_gaussian` and `fourier_uniform` call the C function `NI_FourierFilter` in `ndimage/src/ni_fourier.c`. If the size of any of the dimensions of the supplied input array is less than 2 a segmentation fault occurs in the C code. This PR checks the input array dimension sizes and raises a `ValueError` instead.

The new code is in the C function rather than repeated 3 times in the python functions.

#### Additional information
Test added and docstrings updated.

Problem statement:
ndimage.fourier_ellipsoid segmentation fault
On Scipy 1.5.2. Tested on mac and linux:

```
Python 3.6.8 |Anaconda, Inc.| (default, Dec 30 2018, 01:22:34) 
[GCC 7.3.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> import numpy as np
>>> x = np.ones((1,128))
>>> from scipy import ndimage
>>> ndimage.fourier_ellipsoid(x,3)
Segmentation fault (core dumped)
```

It does work fine if x is defined to be a 2x128 instead of 1x128.
I also tried to reduce the size to a 1 and keep it as a 1x128, still segfaulted in this case.

Hints from issue discussion:
This problem originates in the C file `ni_fourier.c`:

https://github.com/scipy/scipy/blob/2850365e2e447b8e4c1d5ede5193afa5bc682804/scipy/ndimage/src/ni_fourier.c#L227-L229

For each dimension in the `input` array that has a size greater than 1 it allocates a buffer to store calculations. If the size of the dimension is 1 (or 0) it does not do so and the corresponding `params[kk]` is left as a null pointer. A segmentation fault occurs when this null pointer is subsequently deferenced.

Because of the explicit check that the dimension size is `> 1` it looks like the code is not expecting such values. Hence the fix is to check for dimension sizes less than 2 and fail gracefully rather than cause a segmentation fault. I am happy to submit a PR to do this.

Commit message: BUG: avoid segfault in ndimage.fourier_ellipsoid with length-1 dims.

Closes #12763.

Patch to review:
diff --git a/scipy/ndimage/src/ni_fourier.c b/scipy/ndimage/src/ni_fourier.c
index 700878305c45..a018e3641e67 100644
--- a/scipy/ndimage/src/ni_fourier.c
+++ b/scipy/ndimage/src/ni_fourier.c
@@ -225,12 +225,10 @@ int NI_FourierFilter(PyArrayObject *input, PyArrayObject* parameter_array,
         params[kk] = NULL;
     }
     for (kk = 0; kk < PyArray_NDIM(input); kk++) {
-        if (PyArray_DIM(input, kk) > 1) {
-            params[kk] = malloc(PyArray_DIM(input, kk) * sizeof(double));
-            if (!params[kk]) {
-                PyErr_NoMemory();
-                goto exit;
-            }
+        params[kk] = malloc(PyArray_DIM(input, kk) * sizeof(double));
+        if (!params[kk]) {
+            PyErr_NoMemory();
+            goto exit;
         }
     }
 
diff --git a/scipy/ndimage/tests/test_fourier.py b/scipy/ndimage/tests/test_fourier.py
index f7a7e212a456..a4a57352cdaf 100644
--- a/scipy/ndimage/tests/test_fourier.py
+++ b/scipy/ndimage/tests/test_fourier.py
@@ -9,7 +9,7 @@
 
 class TestNdimageFourier:
 
-    @pytest.mark.parametrize('shape', [(32, 16), (31, 15)])
+    @pytest.mark.parametrize('shape', [(32, 16), (31, 15), (1, 10)])
     @pytest.mark.parametrize('dtype, dec',
                              [(numpy.float32, 6), (numpy.float64, 14)])
     def test_fourier_gaussian_real01(self, shape, dtype, dec):
@@ -35,7 +35,7 @@ def test_fourier_gaussian_complex01(self, shape, dtype, dec):
         a = fft.ifft(a, shape[0], 0)
         assert_almost_equal(ndimage.sum(a.real), 1.0, decimal=dec)
 
-    @pytest.mark.parametrize('shape', [(32, 16), (31, 15)])
+    @pytest.mark.parametrize('shape', [(32, 16), (31, 15), (1, 10)])
     @pytest.mark.parametrize('dtype, dec',
                              [(numpy.float32, 6), (numpy.float64, 14)])
     def test_fourier_uniform_real01(self, shape, dtype, dec):
@@ -93,7 +93,7 @@ def test_fourier_shift_complex01(self, shape, dtype, dec):
         assert_array_almost_equal(a.imag, numpy.zeros(shape),
                                   decimal=dec)
 
-    @pytest.mark.parametrize('shape', [(32, 16), (31, 15)])
+    @pytest.mark.parametrize('shape', [(32, 16), (31, 15), (1, 10)])
     @pytest.mark.parametrize('dtype, dec',
                              [(numpy.float32, 5), (numpy.float64, 14)])
     def test_fourier_ellipsoid_real01(self, shape, dtype, dec):


Code review comments:
diff_hunk: @@ -225,12 +225,10 @@ int NI_FourierFilter(PyArrayObject *input, PyArrayObject* parameter_array,
         params[kk] = NULL;
     }
     for (kk = 0; kk < PyArray_NDIM(input); kk++) {
-        if (PyArray_DIM(input, kk) > 1) {
-            params[kk] = malloc(PyArray_DIM(input, kk) * sizeof(double));
path: scipy/ndimage/src/ni_fourier.c
text: @user1:
I think this might create new issues with length-0 axes. Specifically, `_NI_UNIFORM` on line 270 has:
```c
if (params[hh]) {
	params[hh][0] = 1.0;
    // etc...
}
```
which is invalid if `params[hh]` is a non-null pointer to a 0-sized array.


@author:
Yes, you're right. I didn't notice as there are no tests for length-0 dimensions, perhaps not surprisingly as there were no tests for length-1 dimensions.

So, should length-0 dimensions be supported? In general this seems like a bad idea but I am new around here.

@user1:
I would say it's better to support empty dimensions unless there's a good reason not to. In this case, the result should just be an empty array with the same shape as the input. So long as the parameter initialization doesn't crash it should be fine.

Also note places where the code already checks for `> 0`, so again it's making an effort to support it.