1- # -*- coding: utf-8 -*-
2-
31# Copy of the map_coordinates implementation from https://github.com/dask/dask-image/pull/237
42
53# Copyright (c) 2017-2018, dask-image Developers (see AUTHORS.rst for details)
3028# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
3129# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
3230
33- import functools
34- import math
35- from itertools import product
3631
37- from dask import delayed
3832import dask .array as da
3933import numpy as np
34+ from dask import delayed
4035from dask .base import tokenize
41- from dask .highlevelgraph import HighLevelGraph
42- import scipy
43- from scipy .ndimage import affine_transform as ndimage_affine_transform
36+ from scipy .ndimage import labeled_comprehension as ndimage_labeled_comprehension
4437from scipy .ndimage import map_coordinates as ndimage_map_coordinates
45- from scipy .ndimage import labeled_comprehension as \
46- ndimage_labeled_comprehension
47- from scipy .special import sindg , cosdg
48-
49- from ..dispatch ._dispatch_ndinterp import (
50- dispatch_affine_transform ,
51- dispatch_asarray ,
52- dispatch_spline_filter ,
53- dispatch_spline_filter1d ,
54- )
55- from ..ndfilters ._utils import _get_depth_boundary
5638
5739__all__ = [
5840 "map_coordinates" ,
5941]
6042
6143
6244def _map_single_coordinates_array_chunk (
63- input , coordinates , order = 3 , mode = ' constant' ,
64- cval = 0.0 , prefilter = False ):
45+ input , coordinates , order = 3 , mode = " constant" , cval = 0.0 , prefilter = False
46+ ):
6547 """
6648 Central helper function for implementing map_coordinates.
6749
@@ -95,23 +77,26 @@ def _map_single_coordinates_array_chunk(
9577
9678 # linearize input chunk locations
9779 coords_input_chunk_locations_linear = np .sum (
98- coords_input_chunk_locations * np .array (
99- [np .prod (input .numblocks [:dim ])
100- for dim in range (input .ndim )])[::- 1 ],
101- axis = 1 , dtype = np .int64 )
80+ coords_input_chunk_locations
81+ * np .array ([np .prod (input .numblocks [:dim ]) for dim in range (input .ndim )])[::- 1 ],
82+ axis = 1 ,
83+ dtype = np .int64 ,
84+ )
10285
10386 # determine the input chunks that have coords associated and
10487 # count how many coords map onto each input chunk
105- chunk_indices_count = np .bincount (coords_input_chunk_locations_linear ,
106- minlength = np .prod (input .numblocks ))
88+ chunk_indices_count = np .bincount (
89+ coords_input_chunk_locations_linear , minlength = np .prod (input .numblocks )
90+ )
10791 required_input_chunk_indices = np .where (chunk_indices_count > 0 )[0 ]
10892 required_input_chunks = input_chunk_locations [required_input_chunk_indices ]
10993 coord_rc_count = chunk_indices_count [required_input_chunk_indices ]
11094
11195 # inverse mapping: input chunks to coordinates
112- required_input_chunk_coords_indices = \
113- [np .where (coords_input_chunk_locations_linear == irc )[0 ]
114- for irc in required_input_chunk_indices ]
96+ required_input_chunk_coords_indices = [
97+ np .where (coords_input_chunk_locations_linear == irc )[0 ]
98+ for irc in required_input_chunk_indices
99+ ]
115100
116101 # STEP 2: for each input chunk that has been associated to at least
117102 # one coordinate, calculate the minimal slice required to map all
@@ -123,82 +108,89 @@ def _map_single_coordinates_array_chunk(
123108 # Note that this slice can be larger than the given chunk when coords
124109 # lie at chunk borders.
125110 # (probably there's a more efficient way to do this)
126- input_slices_lower = np .array ([np .clip (
127- ndimage_labeled_comprehension (
128- np .floor (coordinates [dim ] - order // 2 ),
129- coords_input_chunk_locations_linear ,
130- required_input_chunk_indices ,
131- np .min , np .int64 , 0 ),
132- 0 , input .shape [dim ] - 1 )
133- for dim in range (input .ndim )])
134-
135- input_slices_upper = np .array ([np .clip (
136- ndimage_labeled_comprehension (
137- np .ceil (coordinates [dim ] + order // 2 ) + 1 ,
138- coords_input_chunk_locations_linear ,
139- required_input_chunk_indices ,
140- np .max , np .int64 , 0 ),
141- 0 , input .shape [dim ])
142- for dim in range (input .ndim )])
143-
144- input_slices = np .array ([input_slices_lower , input_slices_upper ])\
145- .swapaxes (1 , 2 )
111+ input_slices_lower = np .array (
112+ [
113+ np .clip (
114+ ndimage_labeled_comprehension (
115+ np .floor (coordinates [dim ] - order // 2 ),
116+ coords_input_chunk_locations_linear ,
117+ required_input_chunk_indices ,
118+ np .min ,
119+ np .int64 ,
120+ 0 ,
121+ ),
122+ 0 ,
123+ input .shape [dim ] - 1 ,
124+ )
125+ for dim in range (input .ndim )
126+ ]
127+ )
128+
129+ input_slices_upper = np .array (
130+ [
131+ np .clip (
132+ ndimage_labeled_comprehension (
133+ np .ceil (coordinates [dim ] + order // 2 ) + 1 ,
134+ coords_input_chunk_locations_linear ,
135+ required_input_chunk_indices ,
136+ np .max ,
137+ np .int64 ,
138+ 0 ,
139+ ),
140+ 0 ,
141+ input .shape [dim ],
142+ )
143+ for dim in range (input .ndim )
144+ ]
145+ )
146+
147+ input_slices = np .array ([input_slices_lower , input_slices_upper ]).swapaxes (1 , 2 )
146148
147149 # STEP 3: For each previously obtained slice and its associated
148150 # coordinates, define a dask task and apply ndimage.map_coordinates
149151
150152 # prepare building dask graph
151153 # define one task per associated input chunk
152- name = "map_coordinates_chunk-%s" % tokenize (
153- input ,
154- coordinates ,
155- order ,
156- mode ,
157- cval ,
158- prefilter
159- )
154+ name = f"map_coordinates_chunk-{ tokenize (input , coordinates , order , mode , cval , prefilter )} "
160155
161156 keys = [(name , i ) for i in range (len (required_input_chunks ))]
162157
163158 # pair map_coordinates calls with input slices and mapped coordinates
164159 values = []
165160 for irc in range (len (required_input_chunks )):
166161
167- ric_slice = [slice (
168- input_slices [0 ][irc ][dim ],
169- input_slices [1 ][irc ][dim ])
170- for dim in range (input .ndim )]
162+ ric_slice = [
163+ slice (input_slices [0 ][irc ][dim ], input_slices [1 ][irc ][dim ]) for dim in range (input .ndim )
164+ ]
171165 ric_offset = input_slices [0 ][irc ]
172166
173- values .append ((
174- ndimage_map_coordinates ,
175- input [tuple (ric_slice )],
176- coordinates [:, required_input_chunk_coords_indices [irc ]]
177- - ric_offset [:, None ],
178- None ,
179- order ,
180- mode ,
181- cval ,
182- prefilter
183- ))
167+ values .append (
168+ (
169+ ndimage_map_coordinates ,
170+ input [tuple (ric_slice )],
171+ coordinates [:, required_input_chunk_coords_indices [irc ]] - ric_offset [:, None ],
172+ None ,
173+ order ,
174+ mode ,
175+ cval ,
176+ prefilter ,
177+ )
178+ )
184179
185180 # build dask graph
186- dsk = dict (zip (keys , values ))
181+ dsk = dict (zip (keys , values , strict = False ))
187182 ar = da .Array (dsk , name , tuple ([list (coord_rc_count )]), input .dtype )
188183
189184 # STEP 4: rearrange outputs of ndimage.map_coordinates
190185 # to match input order
191- orig_order = np .argsort (
192- [ic for ric_ci in required_input_chunk_coords_indices
193- for ic in ric_ci ])
186+ orig_order = np .argsort ([ic for ric_ci in required_input_chunk_coords_indices for ic in ric_ci ])
194187
195188 # compute result and reorder
196189 # (ordering first would probably unnecessarily inflate the task graph)
197190 return ar .compute ()[orig_order ]
198191
199192
200- def map_coordinates (input , coordinates , order = 3 ,
201- mode = 'constant' , cval = 0.0 , prefilter = False ):
193+ def map_coordinates (input , coordinates , order = 3 , mode = "constant" , cval = 0.0 , prefilter = False ):
202194 """
203195 Wraps ndimage.map_coordinates.
204196
@@ -250,9 +242,7 @@ def map_coordinates(input, coordinates, order=3,
250242 # make sure indices are not split across chunks, i.e. that there's
251243 # no chunking along the first dimension
252244 if len (coordinates .chunks [0 ]) > 1 :
253- coordinates = da .rechunk (
254- coordinates ,
255- (- 1 ,) + coordinates .chunks [1 :])
245+ coordinates = da .rechunk (coordinates , (- 1 ,) + coordinates .chunks [1 :])
256246
257247 # if the input array is not a dask array, convert it into one
258248 if type (input ) is not da .Array :
@@ -263,7 +253,7 @@ def map_coordinates(input, coordinates, order=3,
263253 # a dirty trick: it is split into its components and passed as a delayed
264254 # object, which reconstructs the original array when the task is
265255 # executed. Therefore two `compute` calls are required to obtain the
266- # final result, one of which is peformed by
256+ # final result, one of which is performed by
267257 # `_map_single_coordinates_array_chunk`
268258 output = da .map_blocks (
269259 _map_single_coordinates_array_chunk ,
0 commit comments