Compute Library
 21.02
NEMagnitudePhaseKernel.cpp
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2016-2020 Arm Limited.
3  *
4  * SPDX-License-Identifier: MIT
5  *
6  * Permission is hereby granted, free of charge, to any person obtaining a copy
7  * of this software and associated documentation files (the "Software"), to
8  * deal in the Software without restriction, including without limitation the
9  * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
10  * sell copies of the Software, and to permit persons to whom the Software is
11  * furnished to do so, subject to the following conditions:
12  *
13  * The above copyright notice and this permission notice shall be included in all
14  * copies or substantial portions of the Software.
15  *
16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22  * SOFTWARE.
23  */
25 
26 #include "arm_compute/core/Error.h"
33 
34 #include <arm_neon.h>
35 #include <cstdint>
36 
37 using namespace arm_compute;
38 
39 namespace arm_compute
40 {
41 class Coordinates;
42 } // namespace arm_compute
43 
44 namespace
45 {
46 // Defines for computing atan2
47 constexpr float SCALE_FACTOR = 0.7111111111111111f;
48 constexpr float PI = 3.141592653589793f;
49 constexpr float SCALE_180 = 180.0f / PI;
50 constexpr float SCALE_360 = SCALE_180 * SCALE_FACTOR;
51 constexpr float PI_4 = 0.7853981633974483f;
52 constexpr float COEFF1 = 0.0663f;
53 constexpr float COEFF2 = 0.2447f;
54 } // namespace
55 
56 namespace
57 {
58 inline float32x4_t inv(float32x4_t x)
59 {
60  float32x4_t result = vrecpeq_f32(x);
61  result = vmulq_f32(vrecpsq_f32(x, result), result);
62  return result;
63 }
64 
65 inline float32x4_t atan2_0_360(float32x4_t gx, float32x4_t gy)
66 {
67  const float32x4_t zero = vdupq_n_f32(0.0f);
68  const float32x4_t epsilon = vdupq_n_f32(1e-9f);
69  const float32x4_t piover4 = vdupq_n_f32(PI_4);
70  const float32x4_t coeff1 = vdupq_n_f32(COEFF1);
71  const float32x4_t coeff2 = vdupq_n_f32(COEFF2);
72  const float32x4_t ninety = vdupq_n_f32(90.0f * SCALE_FACTOR);
73  const float32x4_t oneeighty = vdupq_n_f32(180.0f * SCALE_FACTOR);
74  const float32x4_t threesixty = vdupq_n_f32(360.0f * SCALE_FACTOR);
75  const float32x4_t scale = vdupq_n_f32(SCALE_360);
76 
77  float32x4_t abs_gx = vabsq_f32(gx);
78  float32x4_t abs_gy = vabsq_f32(gy);
79  float32x4_t tmin = vminq_f32(abs_gx, abs_gy);
80  float32x4_t tmax = vmaxq_f32(abs_gx, abs_gy);
81  float32x4_t z = vmulq_f32(tmin, inv(vaddq_f32(tmax, epsilon)));
82  float32x4_t absz = vabsq_f32(z);
83  float32x4_t term = vmulq_f32(z, vsubq_f32(vdupq_n_f32(1.0f), absz));
84 
85  /* Compute y = pi/4 * x - x*(abs(x)-1)*(0.2447+0.0663 * abs(x) */
86  float32x4_t result = vaddq_f32(coeff2, vmulq_f32(absz, coeff1));
87  result = vmulq_f32(result, term);
88  result = vmlaq_f32(result, piover4, z);
89 
90  /* Radians to degrees conversion with applied a scale factor in order to have the result [0, 255] */
91  result = vmulq_f32(result, scale);
92 
93  /* If z > 1, result = 90 - result */
94  result = vbslq_f32(vcgeq_f32(abs_gx, abs_gy), result, vsubq_f32(ninety, result));
95 
96  /* Choose correct quadrant */
97  result = vbslq_f32(vcltq_f32(gx, zero), vsubq_f32(oneeighty, result), result);
98  result = vbslq_f32(vcltq_f32(gy, zero), vsubq_f32(threesixty, result), result);
99 
100  return result;
101 }
102 
103 inline float32x4_t atan2_0_180(float32x4_t gx, float32x4_t gy)
104 {
105  const float32x4_t zero = vdupq_n_f32(0.0f);
106  const float32x4_t epsilon = vdupq_n_f32(1e-9f); // epsilon used to avoiding division by 0
107  const float32x4_t piover4 = vdupq_n_f32(PI_4);
108  const float32x4_t coeff1 = vdupq_n_f32(COEFF1);
109  const float32x4_t coeff2 = vdupq_n_f32(COEFF2);
110  const float32x4_t ninety = vdupq_n_f32(90.0f);
111  const float32x4_t oneeighty = vdupq_n_f32(180.0f);
112  const float32x4_t threesixty = vdupq_n_f32(360.0f);
113  const float32x4_t scale = vdupq_n_f32(SCALE_180);
114 
115  float32x4_t abs_gx = vabsq_f32(gx);
116  float32x4_t abs_gy = vabsq_f32(gy);
117  float32x4_t tmin = vminq_f32(abs_gx, abs_gy);
118  float32x4_t tmax = vmaxq_f32(abs_gx, abs_gy);
119  float32x4_t z = vmulq_f32(tmin, inv(vaddq_f32(tmax, epsilon)));
120  float32x4_t absz = vabsq_f32(z);
121 
122  /* Compute y = pi/4 * z - z*(abs(z)-1)*(0.2447+0.0663 * abs(z) */
123  float32x4_t term = vmulq_f32(z, vsubq_f32(vdupq_n_f32(1.0f), absz));
124  float32x4_t result = vaddq_f32(coeff2, vmulq_f32(absz, coeff1));
125  result = vmulq_f32(result, term);
126  result = vmlaq_f32(result, piover4, z);
127 
128  /* Radians to degrees conversion */
129  result = vmulq_f32(result, scale);
130 
131  /* If z > 1, result = 90 - result */
132  result = vbslq_f32(vcgeq_f32(abs_gx, abs_gy), result, vsubq_f32(ninety, result));
133 
134  /* Choose correct quadrant */
135  result = vbslq_f32(vcltq_f32(gx, zero), vsubq_f32(oneeighty, result), result);
136  result = vbslq_f32(vcltq_f32(gy, zero), vsubq_f32(threesixty, result), result);
137  result = vbslq_f32(vcgtq_f32(result, oneeighty), vsubq_f32(result, oneeighty), result);
138 
139  return result;
140 }
141 
142 inline float32x4_t invsqrtv(float32x4_t x)
143 {
144  float32x4_t sqrt_reciprocal = vrsqrteq_f32(x);
145 
146  sqrt_reciprocal = vmulq_f32(vrsqrtsq_f32(vmulq_f32(x, sqrt_reciprocal), sqrt_reciprocal),
147  sqrt_reciprocal);
148  sqrt_reciprocal = vmulq_f32(vrsqrtsq_f32(vmulq_f32(x, sqrt_reciprocal), sqrt_reciprocal),
149  sqrt_reciprocal);
150 
151  return sqrt_reciprocal;
152 }
153 
154 inline float32x4_t sqrtv(float32x4_t x)
155 {
156  float32x4_t res = vdupq_n_f32(0.5f);
157  return vmlaq_f32(res, x, invsqrtv(x));
158 }
159 
160 inline int16x8_t magnitude_l2(int16x8_t input1, int16x8_t input2)
161 {
162  const int32x4x2_t square_x =
163  {
164  {
165  vmull_s16(vget_low_s16(input1), vget_low_s16(input1)),
166  vmull_s16(vget_high_s16(input1), vget_high_s16(input1))
167  }
168  };
169 
170  const int32x4x2_t square_y =
171  {
172  {
173  vmull_s16(vget_low_s16(input2), vget_low_s16(input2)),
174  vmull_s16(vget_high_s16(input2), vget_high_s16(input2))
175  }
176  };
177 
178  const uint32x4x2_t sum =
179  {
180  {
181  vaddq_u32(vreinterpretq_u32_s32(square_x.val[0]), vreinterpretq_u32_s32(square_y.val[0])),
182  vaddq_u32(vreinterpretq_u32_s32(square_x.val[1]), vreinterpretq_u32_s32(square_y.val[1]))
183  }
184  };
185 
186  const float32x4x2_t res =
187  {
188  {
189  sqrtv(vcvtq_f32_u32(sum.val[0])),
190  sqrtv(vcvtq_f32_u32(sum.val[1]))
191  }
192  };
193 
194  return vcombine_s16(vqmovn_s32(vcvtq_s32_f32(res.val[0])),
195  vqmovn_s32(vcvtq_s32_f32(res.val[1])));
196 }
197 
198 inline int16x8_t magnitude_l1(int16x8_t input1, int16x8_t input2)
199 {
200  /* Saturating add */
201  return vqaddq_s16(vqabsq_s16(input1), vqabsq_s16(input2));
202 }
203 
204 inline uint8x8_t phase_signed(int16x8_t input1, int16x8_t input2)
205 {
206  const float32x4_t zeropointfive = vdupq_n_f32(0.5f);
207 
208  float32x4_t inputx_f32_high = vcvtq_f32_s32(vmovl_s16(vget_high_s16(input1)));
209  float32x4_t inputx_f32_low = vcvtq_f32_s32(vmovl_s16(vget_low_s16(input1)));
210  float32x4_t inputy_f32_high = vcvtq_f32_s32(vmovl_s16(vget_high_s16(input2)));
211  float32x4_t inputy_f32_low = vcvtq_f32_s32(vmovl_s16(vget_low_s16(input2)));
212 
213  /* Compute fast atan2 */
214  float32x4_t angle_high = atan2_0_360(inputx_f32_high, inputy_f32_high);
215  float32x4_t angle_low = atan2_0_360(inputx_f32_low, inputy_f32_low);
216 
217  angle_high = vaddq_f32(angle_high, zeropointfive);
218  angle_low = vaddq_f32(angle_low, zeropointfive);
219 
220  return vmovn_u16(vcombine_u16(vqmovun_s32(vcvtq_s32_f32(angle_low)),
221  vqmovun_s32(vcvtq_s32_f32(angle_high))));
222 }
223 
224 inline uint8x8_t phase_unsigned(int16x8_t input1, int16x8_t input2)
225 {
226  const float32x4_t zeropointfive = vdupq_n_f32(0.5f);
227 
228  float32x4_t inputx_f32_high = vcvtq_f32_s32(vmovl_s16(vget_high_s16(input1)));
229  float32x4_t inputx_f32_low = vcvtq_f32_s32(vmovl_s16(vget_low_s16(input1)));
230  float32x4_t inputy_f32_high = vcvtq_f32_s32(vmovl_s16(vget_high_s16(input2)));
231  float32x4_t inputy_f32_low = vcvtq_f32_s32(vmovl_s16(vget_low_s16(input2)));
232 
233  /* Compute fast atan2 */
234  float32x4_t angle_high = atan2_0_180(inputx_f32_high, inputy_f32_high);
235  float32x4_t angle_low = atan2_0_180(inputx_f32_low, inputy_f32_low);
236 
237  angle_high = vaddq_f32(angle_high, zeropointfive);
238  angle_low = vaddq_f32(angle_low, zeropointfive);
239 
240  return vmovn_u16(vcombine_u16(vqmovun_s32(vcvtq_s32_f32(angle_low)),
241  vqmovun_s32(vcvtq_s32_f32(angle_high))));
242 }
243 } // namespace
244 
245 template <MagnitudeType mag_type, PhaseType phase_type>
247  : _func(nullptr), _gx(nullptr), _gy(nullptr), _magnitude(nullptr), _phase(nullptr)
248 {
249 }
250 
251 template <MagnitudeType mag_type, PhaseType phase_type>
253 {
256  ARM_COMPUTE_ERROR_ON((nullptr == magnitude) && (nullptr == phase));
257 
258  const bool run_mag = magnitude != nullptr;
259  const bool run_phase = phase != nullptr;
260 
261  if(run_mag)
262  {
264  }
265 
266  if(run_phase)
267  {
269  }
270 
271  _gx = gx;
272  _gy = gy;
273  _magnitude = magnitude;
274  _phase = phase;
275 
276  if(run_mag && run_phase)
277  {
278  /* Run magnitude and phase */
280  }
281  else
282  {
283  if(run_mag)
284  {
285  /* Run magnitude */
287  }
288  else if(run_phase)
289  {
290  /* Run phase */
292  }
293  else
294  {
295  ARM_COMPUTE_ERROR("At least one output must be NOT NULL");
296  }
297  }
298 
299  constexpr unsigned int num_elems_processed_per_iteration = 16;
300 
301  // Configure kernel window
302  Window win = calculate_max_window(*gx->info(), Steps(num_elems_processed_per_iteration));
303  AccessWindowHorizontal magnitude_access(magnitude == nullptr ? nullptr : magnitude->info(), 0, num_elems_processed_per_iteration);
304  AccessWindowHorizontal phase_access(phase == nullptr ? nullptr : phase->info(), 0, num_elems_processed_per_iteration);
305 
309  magnitude_access,
310  phase_access);
311 
313  gy->info()->valid_region());
314 
315  magnitude_access.set_valid_region(win, valid_region);
316  phase_access.set_valid_region(win, valid_region);
317 
318  INEKernel::configure(win);
319 }
320 
321 template <MagnitudeType mag_type, PhaseType phase_type>
323 {
324  Iterator gx(_gx, window);
325  Iterator gy(_gy, window);
326  Iterator magnitude(_magnitude, window);
327 
328  execute_window_loop(window, [&](const Coordinates &)
329  {
330  const int16x8x2_t input1 =
331  {
332  {
333  vld1q_s16(reinterpret_cast<int16_t *>(gx.ptr())),
334  vld1q_s16(reinterpret_cast<int16_t *>(gx.ptr()) + 8)
335  }
336  };
337 
338  const int16x8x2_t input2 =
339  {
340  {
341  vld1q_s16(reinterpret_cast<int16_t *>(gy.ptr())),
342  vld1q_s16(reinterpret_cast<int16_t *>(gy.ptr()) + 8)
343  }
344  };
345 
346  /* Compute magnitude */
347  int16x8x2_t mag{ {} };
348 
349  if(MagnitudeType::L2NORM == mag_type)
350  {
351  mag.val[0] = magnitude_l2(input1.val[0], input2.val[0]);
352  mag.val[1] = magnitude_l2(input1.val[1], input2.val[1]);
353  }
354  else
355  {
356  mag.val[0] = magnitude_l1(input1.val[0], input2.val[0]);
357  mag.val[1] = magnitude_l1(input1.val[1], input2.val[1]);
358  }
359 
360  /* Store magnitude */
361  vst1q_s16(reinterpret_cast<int16_t *>(magnitude.ptr()), mag.val[0]);
362  vst1q_s16(reinterpret_cast<int16_t *>(magnitude.ptr()) + 8, mag.val[1]);
363  },
364  gx, gy, magnitude);
365 }
366 
367 template <MagnitudeType mag_type, PhaseType phase_type>
369 {
370  Iterator gx(_gx, window);
371  Iterator gy(_gy, window);
372  Iterator phase(_phase, window);
373 
374  execute_window_loop(window, [&](const Coordinates &)
375  {
376  const int16x8x2_t input1 =
377  {
378  {
379  vld1q_s16(reinterpret_cast<int16_t *>(gx.ptr())),
380  vld1q_s16(reinterpret_cast<int16_t *>(gx.ptr()) + 8)
381  }
382  };
383 
384  const int16x8x2_t input2 =
385  {
386  {
387  vld1q_s16(reinterpret_cast<int16_t *>(gy.ptr())),
388  vld1q_s16(reinterpret_cast<int16_t *>(gy.ptr()) + 8)
389  }
390  };
391 
392  /* Compute phase */
393  uint8x8x2_t vphase{ {} };
394 
395  if(PhaseType::SIGNED == phase_type)
396  {
397  vphase.val[0] = phase_signed(input1.val[0], input2.val[0]);
398  vphase.val[1] = phase_signed(input1.val[1], input2.val[1]);
399  }
400  else
401  {
402  vphase.val[0] = phase_unsigned(input1.val[0], input2.val[0]);
403  vphase.val[1] = phase_unsigned(input1.val[1], input2.val[1]);
404  }
405 
406  /* Store phase */
407  vst1q_u8(phase.ptr(), vcombine_u8(vphase.val[0], vphase.val[1]));
408  },
409  gx, gy, phase);
410 }
411 
412 template <MagnitudeType mag_type, PhaseType phase_type>
414 {
415  Iterator gx(_gx, window);
416  Iterator gy(_gy, window);
417  Iterator magnitude(_magnitude, window);
418  Iterator phase(_phase, window);
419 
420  execute_window_loop(window, [&](const Coordinates &)
421  {
422  const int16x8x2_t input1 =
423  {
424  {
425  vld1q_s16(reinterpret_cast<int16_t *>(gx.ptr())),
426  vld1q_s16(reinterpret_cast<int16_t *>(gx.ptr()) + 8)
427  }
428  };
429 
430  const int16x8x2_t input2 =
431  {
432  {
433  vld1q_s16(reinterpret_cast<int16_t *>(gy.ptr())),
434  vld1q_s16(reinterpret_cast<int16_t *>(gy.ptr()) + 8)
435  }
436  };
437 
438  /* Compute magnitude */
439  int16x8x2_t mag{ {} };
440 
441  if(MagnitudeType::L2NORM == mag_type)
442  {
443  mag.val[0] = magnitude_l2(input1.val[0], input2.val[0]);
444  mag.val[1] = magnitude_l2(input1.val[1], input2.val[1]);
445  }
446  else
447  {
448  mag.val[0] = magnitude_l1(input1.val[0], input2.val[0]);
449  mag.val[1] = magnitude_l1(input1.val[1], input2.val[1]);
450  }
451 
452  /* Store magnitude */
453  vst1q_s16(reinterpret_cast<int16_t *>(magnitude.ptr()), mag.val[0]);
454  vst1q_s16(reinterpret_cast<int16_t *>(magnitude.ptr()) + 8, mag.val[1]);
455 
456  /* Compute phase */
457  uint8x8x2_t vphase{ {} };
458 
459  if(PhaseType::SIGNED == phase_type)
460  {
461  vphase.val[0] = phase_signed(input1.val[0], input2.val[0]);
462  vphase.val[1] = phase_signed(input1.val[1], input2.val[1]);
463  }
464  else
465  {
466  vphase.val[0] = phase_unsigned(input1.val[0], input2.val[0]);
467  vphase.val[1] = phase_unsigned(input1.val[1], input2.val[1]);
468  }
469 
470  /* Store phase */
471  vst1q_u8(phase.ptr(), vcombine_u8(vphase.val[0], vphase.val[1]));
472  },
473  gx, gy, magnitude, phase);
474 }
475 
476 template <MagnitudeType mag_type, PhaseType phase_type>
478 {
479  ARM_COMPUTE_UNUSED(info);
482  ARM_COMPUTE_ERROR_ON(_func == nullptr);
483 
484  (this->*_func)(window);
485 }
486 
Window calculate_max_window(const ValidRegion &valid_region, const Steps &steps, bool skip_border, BorderSize border_size)
const Window & window() const
The maximum window the kernel can be executed on.
Definition: IKernel.cpp:28
Angle range: [0, 360].
#define ARM_COMPUTE_ERROR(msg)
Print the given message then throw an std::runtime_error.
Definition: Error.h:352
1 channel, 1 U8 per channel
DATA_TYPE sum(__global const DATA_TYPE *input)
Calculate sum of a vector.
#define ARM_COMPUTE_ERROR_ON(cond)
If the condition is true then an error message is printed and an exception thrown.
Definition: Error.h:466
SimpleTensor< uint8_t > phase(const SimpleTensor< T > &gx, const SimpleTensor< T > &gy, PhaseType phase_type)
Definition: Phase.cpp:35
uchar16 phase_signed(DATA_TYPE16 a, DATA_TYPE16 b)
Calculates signed phase between two inputs.
DATA_TYPE16 magnitude_l2(int16 a, int16 b)
Calculates L2 normalization between two inputs.
Template interface for the kernel to compute magnitude and phase.
const ValidRegion valid_region
Definition: Scale.cpp:221
Interface for Neon tensor.
Definition: ITensor.h:36
Copyright (c) 2017-2021 Arm Limited.
virtual ValidRegion valid_region() const =0
Valid region of the tensor.
bool update_window_and_padding(Window &win, Ts &&... patterns)
Update window and padding size for each of the access patterns.
Definition: WindowHelpers.h:46
#define ARM_COMPUTE_UNUSED(...)
To avoid unused variables warnings.
Definition: Error.h:152
Class to describe a number of elements in each dimension.
Definition: Steps.h:40
Coordinates of an item.
Definition: Coordinates.h:37
Implementation of a row access pattern.
uchar16 phase_unsigned(DATA_TYPE16 a, DATA_TYPE16 b)
Calculates unsigned phase between two inputs.
virtual ITensorInfo * info() const =0
Interface to be implemented by the child class to return the tensor&#39;s metadata.
constexpr uint8_t * ptr() const
Return a pointer to the current pixel.
Definition: Helpers.inl:139
ValidRegion intersect_valid_regions(const Ts &... regions)
Intersect multiple valid regions.
Definition: WindowHelpers.h:74
#define ARM_COMPUTE_ERROR_ON_UNCONFIGURED_KERNEL(k)
Definition: Validate.h:941
1 channel, 1 S16 per channel
L2 normalization type.
#define ARM_COMPUTE_ERROR_ON_DATA_TYPE_CHANNEL_NOT_IN(t, c,...)
Definition: Validate.h:790
ScaleKernelInfo info(interpolation_policy, default_border_mode, PixelValue(), sampling_policy, false)
void run(const Window &window, const ThreadInfo &info) override
Execute the kernel on the passed window.
Information about executing thread and CPU.
Definition: CPPTypes.h:235
unsigned int num_elems_processed_per_iteration
void execute_window_loop(const Window &w, L &&lambda_function, Ts &&... iterators)
Iterate through the passed window, automatically adjusting the iterators and calling the lambda_funct...
Definition: Helpers.inl:77
__kernel void magnitude_phase(__global uchar *gx_ptr, uint gx_stride_x, uint gx_step_x, uint gx_stride_y, uint gx_step_y, uint gx_offset_first_element_in_bytes, __global uchar *gy_ptr, uint gy_stride_x, uint gy_step_x, uint gy_stride_y, uint gy_step_y, uint gy_offset_first_element_in_bytes, __global uchar *magnitude_ptr, uint magnitude_stride_x, uint magnitude_step_x, uint magnitude_stride_y, uint magnitude_step_y, uint magnitude_offset_first_element_in_bytes, __global uchar *phase_ptr, uint phase_stride_x, uint phase_step_x, uint phase_stride_y, uint phase_step_y, uint phase_offset_first_element_in_bytes)
Calculate the magnitude and phase of given the gradients of an image.
Container for valid region of a window.
Definition: Types.h:188
void configure(const ITensor *gx, const ITensor *gy, ITensor *magnitude, ITensor *phase)
Initialise the kernel&#39;s input, output.
Iterator updated by execute_window_loop for each window element.
Definition: Helpers.h:46
Describe a multidimensional execution window.
Definition: Window.h:39
SimpleTensor< T > magnitude(const SimpleTensor< T > &gx, const SimpleTensor< T > &gy, MagnitudeType magnitude_type)
Definition: Magnitude.cpp:35
#define ARM_COMPUTE_ERROR_ON_INVALID_SUBWINDOW(f, s)
Definition: Validate.h:205
DATA_TYPE16 magnitude_l1(DATA_TYPE16 a, DATA_TYPE16 b)
Calculates L1 normalization between two inputs.