No Stack Trace Segmentation Fault -- reconstruction_loop Error (greyreconstruct.py)

jeffditty commented 5 years ago

Description

When running fill depressions (Python Library pysheds) on a 6 GB DEM, running into a segmentation fault error. Using IPDB, I located the segmentation fault error happens at line 205 of greyreconstruct.py. The code works fine when working on a subsetted sample of the 6 GB DEM, but when applied full scale, segmentation fault happens continuously. I am currently trying to get permissions to give you the DEM causing the error.

start = index_sorted[0]
reconstruction_loop(value_rank, prev, next, nb_strides, start,
                    image_stride)

I am running on a linux OS r5.12xlarge EC2 instance with a container.

Way to reproduce


# Place the full code we need to recreate your issue here

from pysheds.grid import Grid
grid = Grid.from_raster(ec2_local_Dem, data_name='dem')
grid.fill_depressions(data='dem', out_name='filled_dem')

# upload all necessary images to github too!
Working on getting permissions

Version information

# Paste the output of the following python commands
from __future__ import print_function

import sys; print(sys.version)
3.7.3 | packaged by conda-forge | (default, Mar 27 2019, 23:01:00) 
[GCC 7.3.0]

import platform; print(platform.platform())
Linux-4.4.0-1052-aws-x86_64-with-debian-9.7

import skimage; print("scikit-image version: {}".format(skimage.__version__))
scikit-image version: 0.15.0

import numpy; print("numpy version: {}".format(numpy.__version__))
numpy version: 1.16.3

# your output here. 

The following is about 20 lines of stepping through the greyreconstruct.py code.

> /opt/conda/lib/python3.7/site-packages/numpy/core/numeric.py(337)full()
    336     multiarray.copyto(a, fill_value, casting='unsafe')
--> 337     return a
    338 

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(171)reconstruction()
    170     images = np.full(dims, pad_value, dtype='float64')
--> 171     images[(0, *inside_slices)] = seed
    172     images[(1, *inside_slices)] = mask

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(172)reconstruction()
    171     images[(0, *inside_slices)] = seed
--> 172     images[(1, *inside_slices)] = mask
    173 

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(176)reconstruction()
    175     # a flattened array
--> 176     value_stride = np.array(images.strides[1:]) // images.dtype.itemsize
    177     image_stride = images.strides[0] // images.dtype.itemsize

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(177)reconstruction()
    176     value_stride = np.array(images.strides[1:]) // images.dtype.itemsize
--> 177     image_stride = images.strides[0] // images.dtype.itemsize
    178     selem_mgrid = np.mgrid[[slice(-o, d - o)

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(178)reconstruction()
    177     image_stride = images.strides[0] // images.dtype.itemsize
--> 178     selem_mgrid = np.mgrid[[slice(-o, d - o)
    179                             for d, o in zip(selem.shape, offset)]]

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(179)reconstruction()
    178     selem_mgrid = np.mgrid[[slice(-o, d - o)
--> 179                             for d, o in zip(selem.shape, offset)]]
    180     selem_offsets = selem_mgrid[:, selem].transpose()

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(180)reconstruction()
    179                             for d, o in zip(selem.shape, offset)]]
--> 180     selem_offsets = selem_mgrid[:, selem].transpose()
    181     nb_strides = np.array([np.sum(value_stride * selem_offset)

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(181)reconstruction()
    180     selem_offsets = selem_mgrid[:, selem].transpose()
--> 181     nb_strides = np.array([np.sum(value_stride * selem_offset)
    182                            for selem_offset in selem_offsets], np.int32)

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(182)reconstruction()
    181     nb_strides = np.array([np.sum(value_stride * selem_offset)
--> 182                            for selem_offset in selem_offsets], np.int32)
    183 

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(184)reconstruction()
    183 
--> 184     images = images.flatten()
    185 

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(187)reconstruction()
    186     # Erosion goes smallest to largest; dilation goes largest to smallest.
--> 187     index_sorted = np.argsort(images).astype(np.int32)
    188     if method == 'dilation':

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(188)reconstruction()
    187     index_sorted = np.argsort(images).astype(np.int32)
--> 188     if method == 'dilation':
    189         index_sorted = index_sorted[::-1]

ipdb> index_sorted                                                                                                                   
array([ 1117839987,  1292776875,  1292776874, ...,   744169283,
         744169305, -2001703307], dtype=int32)
ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(192)reconstruction()
    191     # Make a linked list of pixels sorted by value. -1 is the list terminator.
--> 192     prev = np.full(len(images), -1, np.int32)
    193     next = np.full(len(images), -1, np.int32)

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(193)reconstruction()
    192     prev = np.full(len(images), -1, np.int32)
--> 193     next = np.full(len(images), -1, np.int32)
    194     prev[index_sorted[1:]] = index_sorted[:-1]

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(194)reconstruction()
    193     next = np.full(len(images), -1, np.int32)
--> 194     prev[index_sorted[1:]] = index_sorted[:-1]
    195     next[index_sorted[:-1]] = index_sorted[1:]

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(195)reconstruction()
    194     prev[index_sorted[1:]] = index_sorted[:-1]
--> 195     next[index_sorted[:-1]] = index_sorted[1:]
    196 

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(198)reconstruction()
    197     # Cython inner-loop compares the rank of pixel values.
--> 198     if method == 'dilation':
    199         value_rank, value_map = rank_order(images)

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(200)reconstruction()
    199         value_rank, value_map = rank_order(images)
--> 200     elif method == 'erosion':
    201         value_rank, value_map = rank_order(-images)

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(201)reconstruction()
    200     elif method == 'erosion':
--> 201         value_rank, value_map = rank_order(-images)
    202         value_map = -value_map

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(202)reconstruction()
    201         value_rank, value_map = rank_order(-images)
--> 202         value_map = -value_map
    203 

ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(204)reconstruction()
    203 
--> 204     start = index_sorted[0]
    205     reconstruction_loop(value_rank, prev, next, nb_strides, start,

ipdb> value_map                                                                                                                      
array([ 2613.59375,  2613.3125 ,  2613.     , ...,  -884.84375,
        -885.625  , -9999.     ])
ipdb> n                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(205)reconstruction()
    204     start = index_sorted[0]
--> 205     reconstruction_loop(value_rank, prev, next, nb_strides, start,
    206                         image_stride)

ipdb> s                                                                                                                              
> /opt/conda/lib/python3.7/site-packages/skimage/morphology/greyreconstruct.py(206)reconstruction()
    205     reconstruction_loop(value_rank, prev, next, nb_strides, start,
--> 206                         image_stride)
    207 

ipdb> s                                                                                                                              
Segmentation fault (core dumped)

stefanv commented 5 years ago

Thanks for the report @jeffditty ; you should never see segfaults, so this looks like a bug. The data would be helpful, especially if you can provide a small image that triggers it.

jeffditty commented 5 years ago

Unfortunately, I tried and we cannot get around the fact that the data is proprietary. Is an image without the actual DEM useful? @stefanv

hmaarrfk commented 5 years ago

Yeah, as far as helping improve the library, if you can recreate the issue with np.rand that would be fine too ;)

grlee77 commented 2 years ago

It has been a long time, but I think #6277 is likely the same issue as this one. There is more recent discussion there and I think we will got momentum to address the issue soon.

scikit-image / scikit-image