ACloudViewer  3.9.4
A Modern Library for 3D Data Processing
ContinuousConvBackpropFilter.h
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // - CloudViewer: www.cloudViewer.org -
3 // ----------------------------------------------------------------------------
4 // Copyright (c) 2018-2024 www.cloudViewer.org
5 // SPDX-License-Identifier: MIT
6 // ----------------------------------------------------------------------------
7 
8 #pragma once
9 
10 #include <tbb/parallel_for.h>
11 
12 #include <mutex>
13 
14 #include "ml/impl/continuous_conv/CoordinateTransformation.h"
15 
16 namespace cloudViewer {
17 namespace ml {
18 namespace impl {
19 
20 // Implementation of CConvBackropFilterCPU
21 template <class TFeat,
22  class TOut,
23  class TReal,
24  class TIndex,
25  InterpolationMode INTERPOLATION,
26  CoordinateMapping MAPPING,
27  bool ALIGN_CORNERS,
28  bool INDIVIDUAL_EXTENT,
29  bool ISOTROPIC_EXTENT,
30  bool POINT_IMPORTANCE>
31 void _CConvBackropFilterCPU(TOut* filter_backprop,
32  const std::vector<int>& filter_dims,
33  size_t num_out,
34  const TReal* out_positions,
35  size_t num_inp,
36  const TReal* inp_positions,
37  const TFeat* inp_features,
38  const TFeat* inp_importance,
39  size_t neighbors_index_size,
40  const TIndex* neighbors_index,
41  const TFeat* neighbors_importance,
42  const int64_t* neighbors_row_splits,
43  const TReal* extents,
44  const TReal* offsets,
45  const TFeat* out_features_gradient,
46  bool normalize) {
47  const bool NEIGHBORS_IMPORTANCE = neighbors_importance;
48  // const int VECSIZE = 32;
49 #define VECSIZE 32
50  typedef Eigen::Array<TReal, VECSIZE, 1> Vec_t;
51  typedef InterpolationVec<TReal, VECSIZE, INTERPOLATION> InterpolationVec_t;
52  InterpolationVec_t interpolation;
53 
54  const int in_channels = filter_dims[filter_dims.size() - 2];
55  const int out_channels = filter_dims[filter_dims.size() - 1];
56 
57  int spatial_filter_size = 1;
58  for (int i = 0; i < 3; ++i) spatial_filter_size *= filter_dims[i];
59  const int total_filter_size =
60  spatial_filter_size * in_channels * out_channels;
61  Eigen::Array<int, 3, 1> filter_size_xyz(filter_dims[2], filter_dims[1],
62  filter_dims[0]);
63 
64  memset(filter_backprop, 0, sizeof(TOut) * total_filter_size);
65  std::mutex filter_backprop_mutex;
66 
67  typedef Eigen::Array<TFeat, VECSIZE, Eigen::Dynamic> Matrix;
68  typedef Eigen::Array<TReal, VECSIZE, 3> Matrix3C;
69 
70  tbb::parallel_for(
71  tbb::blocked_range<size_t>(0, num_out, 32),
72  [&](const tbb::blocked_range<size_t>& r) {
73  int range_length = r.end() - r.begin();
74 
75  Eigen::Matrix<TFeat, Eigen::Dynamic, Eigen::Dynamic> B(
76  in_channels * spatial_filter_size, range_length);
77  B.setZero();
78  Eigen::Matrix<TFeat, Eigen::Dynamic, Eigen::Dynamic> C(
79  out_channels, range_length);
80 
81  Matrix infeat(VECSIZE, in_channels);
82 
83  Eigen::Array<TReal, 3, 1> offsets_(offsets[0], offsets[1],
84  offsets[2]);
85 
86  Matrix3C inv_extents;
87  if (INDIVIDUAL_EXTENT == false) {
88  if (ISOTROPIC_EXTENT) {
89  inv_extents = 1 / extents[0];
90  } else {
91  inv_extents.col(0) = 1 / extents[0];
92  inv_extents.col(1) = 1 / extents[1];
93  inv_extents.col(2) = 1 / extents[2];
94  }
95  }
96 
97  for (size_t out_idx = r.begin(); out_idx != r.end();
98  ++out_idx) {
99  const int out_col = out_idx - r.begin();
100  const size_t neighbor_start = neighbors_row_splits[out_idx];
101  const size_t neighbor_end =
102  neighbors_row_splits[out_idx + 1];
103  TOut normalizer(0);
104 
105  if (INDIVIDUAL_EXTENT) {
106  if (ISOTROPIC_EXTENT) {
107  inv_extents = 1 / extents[out_idx];
108  } else {
109  inv_extents.col(0) = 1 / extents[3 * out_idx + 0];
110  inv_extents.col(1) = 1 / extents[3 * out_idx + 1];
111  inv_extents.col(2) = 1 / extents[3 * out_idx + 2];
112  }
113  }
114 
115  typename InterpolationVec_t::Weight_t interp_weights;
116  typename InterpolationVec_t::Idx_t interp_indices;
117 
118  int vec_valid_count = 0;
119  Vec_t x, y, z;
120 
121  // set to zero to avoid problems with vectors with less than
122  // VECSIZE valid entries
123  x.setZero();
124  y.setZero();
125  z.setZero();
126  for (size_t n = neighbor_start; n < neighbor_end; ++n) {
127  const size_t inp_idx = neighbors_index[n];
128  const int i = vec_valid_count;
129  x(i) = inp_positions[inp_idx * 3 + 0] -
130  out_positions[out_idx * 3 + 0];
131  y(i) = inp_positions[inp_idx * 3 + 1] -
132  out_positions[out_idx * 3 + 1];
133  z(i) = inp_positions[inp_idx * 3 + 2] -
134  out_positions[out_idx * 3 + 2];
135 
136  const TFeat n_importance =
137  (NEIGHBORS_IMPORTANCE ? neighbors_importance[n]
138  : TFeat(1));
139  normalizer += TOut(n_importance);
140 
141  for (int ic = 0; ic < in_channels; ++ic)
142  infeat(i, ic) =
143  inp_features[inp_idx * in_channels + ic];
144 
145  TFeat importance = TFeat(1);
146  if (POINT_IMPORTANCE)
147  importance = inp_importance[inp_idx];
148  if (NEIGHBORS_IMPORTANCE) importance *= n_importance;
149 
150  if (POINT_IMPORTANCE || NEIGHBORS_IMPORTANCE) {
151  for (int ic = 0; ic < in_channels; ++ic)
152  infeat(i, ic) *= importance;
153  }
154 
155  ++vec_valid_count;
156  if (vec_valid_count == VECSIZE) {
157  ComputeFilterCoordinates<ALIGN_CORNERS, MAPPING,
158  TReal, VECSIZE>(
159  x, y, z, filter_size_xyz, inv_extents,
160  offsets_);
161  interpolation.Interpolate(
162  interp_weights, interp_indices, x, y, z,
163  filter_size_xyz, in_channels);
164  for (int k = 0; k < VECSIZE; ++k)
165  for (int j = 0; j < InterpolationVec_t::Size();
166  ++j) {
167  for (int ic = 0; ic < in_channels; ++ic)
168  B(interp_indices(j, k) + ic, out_col) +=
169  TFeat(interp_weights(j, k)) *
170  infeat(k, ic);
171  }
172  vec_valid_count = 0;
173  }
174  }
175  if (vec_valid_count) {
176  ComputeFilterCoordinates<ALIGN_CORNERS, MAPPING, TReal,
177  VECSIZE>(
178  x, y, z, filter_size_xyz, inv_extents,
179  offsets_);
180  interpolation.Interpolate(interp_weights,
181  interp_indices, x, y, z,
182  filter_size_xyz, in_channels);
183  for (int k = 0; k < vec_valid_count; ++k)
184  for (int j = 0; j < InterpolationVec_t::Size();
185  ++j) {
186  for (int ic = 0; ic < in_channels; ++ic)
187  B(interp_indices(j, k) + ic, out_col) +=
188  TFeat(interp_weights(j, k)) *
189  infeat(k, ic);
190  }
191  }
192 
193  C.col(out_col) = Eigen::Map<
194  const Eigen::Array<TFeat, Eigen::Dynamic, 1>>(
195  out_features_gradient + out_idx * out_channels,
196  out_channels, 1);
197 
198  if (normalize && normalizer != TOut(0))
199  C.col(out_col) /= TFeat(normalizer);
200 
201  } // out_idx
202 
203  Eigen::Matrix<TOut, Eigen::Dynamic, Eigen::Dynamic> A(
204  out_channels, spatial_filter_size * in_channels);
205 
206  A = (C * B.transpose()).template cast<TOut>();
207 
208  {
209  std::lock_guard<std::mutex> lock(filter_backprop_mutex);
210  int linear_i = 0;
211  for (int j = 0; j < spatial_filter_size * in_channels; ++j)
212  for (int i = 0; i < out_channels; ++i, ++linear_i) {
213  filter_backprop[linear_i] += A(i, j);
214  }
215  }
216  });
217 #undef VECSIZE
218 }
219 
292 template <class TFeat, class TOut, class TReal, class TIndex>
293 void CConvBackpropFilterCPU(TOut* filter_backprop,
294  const std::vector<int>& filter_dims,
295  size_t num_out,
296  const TReal* out_positions,
297  size_t num_inp,
298  const TReal* inp_positions,
299  const TFeat* inp_features,
300  const TFeat* inp_importance,
301  size_t neighbors_index_size,
302  const TIndex* neighbors_index,
303  const TFeat* neighbors_importance,
304  const int64_t* neighbors_row_splits,
305  const TReal* extents,
306  const TReal* offsets,
307  const TFeat* out_features_gradient,
308  InterpolationMode interpolation,
309  CoordinateMapping coordinate_mapping,
310  bool align_corners,
311  bool individual_extent,
312  bool isotropic_extent,
313  bool normalize) {
314  bool has_importance = inp_importance;
315 
316 #define FN_PARAMETERS \
317  filter_backprop, filter_dims, num_out, out_positions, num_inp, \
318  inp_positions, inp_features, inp_importance, neighbors_index_size, \
319  neighbors_index, neighbors_importance, neighbors_row_splits, \
320  extents, offsets, out_features_gradient, normalize
321 
322 #define CALL_TEMPLATE(INTERPOLATION, MAPPING, ALIGN_CORNERS, \
323  INDIVIDUAL_EXTENT, ISOTROPIC_EXTENT, HAS_IMPORTANCE) \
324  if (INTERPOLATION == interpolation && MAPPING == coordinate_mapping && \
325  ALIGN_CORNERS == align_corners && \
326  INDIVIDUAL_EXTENT == individual_extent && \
327  ISOTROPIC_EXTENT == isotropic_extent && \
328  HAS_IMPORTANCE == has_importance) \
329  _CConvBackropFilterCPU<TFeat, TOut, TReal, TIndex, INTERPOLATION, \
330  MAPPING, ALIGN_CORNERS, INDIVIDUAL_EXTENT, \
331  ISOTROPIC_EXTENT, HAS_IMPORTANCE>( \
332  FN_PARAMETERS);
333 
334 #define CALL_TEMPLATE2(INTERPOLATION, MAPPING) \
335  CALL_TEMPLATE(INTERPOLATION, MAPPING, true, true, true, true) \
336  CALL_TEMPLATE(INTERPOLATION, MAPPING, true, true, true, false) \
337  CALL_TEMPLATE(INTERPOLATION, MAPPING, true, true, false, true) \
338  CALL_TEMPLATE(INTERPOLATION, MAPPING, true, true, false, false) \
339  CALL_TEMPLATE(INTERPOLATION, MAPPING, true, false, true, true) \
340  CALL_TEMPLATE(INTERPOLATION, MAPPING, true, false, true, false) \
341  CALL_TEMPLATE(INTERPOLATION, MAPPING, true, false, false, true) \
342  CALL_TEMPLATE(INTERPOLATION, MAPPING, true, false, false, false) \
343  CALL_TEMPLATE(INTERPOLATION, MAPPING, false, true, true, true) \
344  CALL_TEMPLATE(INTERPOLATION, MAPPING, false, true, true, false) \
345  CALL_TEMPLATE(INTERPOLATION, MAPPING, false, true, false, true) \
346  CALL_TEMPLATE(INTERPOLATION, MAPPING, false, true, false, false) \
347  CALL_TEMPLATE(INTERPOLATION, MAPPING, false, false, true, true) \
348  CALL_TEMPLATE(INTERPOLATION, MAPPING, false, false, true, false) \
349  CALL_TEMPLATE(INTERPOLATION, MAPPING, false, false, false, true) \
350  CALL_TEMPLATE(INTERPOLATION, MAPPING, false, false, false, false)
351 
352 #define CALL_TEMPLATE3(INTERPOLATION) \
353  CALL_TEMPLATE2(INTERPOLATION, CoordinateMapping::BALL_TO_CUBE_RADIAL) \
354  CALL_TEMPLATE2(INTERPOLATION, \
355  CoordinateMapping::BALL_TO_CUBE_VOLUME_PRESERVING) \
356  CALL_TEMPLATE2(INTERPOLATION, CoordinateMapping::IDENTITY)
357 
358 #define CALL_TEMPLATE4 \
359  CALL_TEMPLATE3(InterpolationMode::LINEAR) \
360  CALL_TEMPLATE3(InterpolationMode::LINEAR_BORDER) \
361  CALL_TEMPLATE3(InterpolationMode::NEAREST_NEIGHBOR)
362 
363  CALL_TEMPLATE4
364 
365 #undef CALL_TEMPLATE
366 #undef CALL_TEMPLATE2
367 #undef CALL_TEMPLATE3
368 #undef CALL_TEMPLATE4
369 
370 #undef FN_PARAMETERS
371 }
372 
373 } // namespace impl
374 } // namespace ml
375 } // namespace cloudViewer
CALL_TEMPLATE4
#define CALL_TEMPLATE4
VECSIZE
#define VECSIZE
normalize
__host__ __device__ float2 normalize(float2 v)
Definition: cutil_math.h:1179
cloudViewer::ml::impl::_CConvBackropFilterCPU
void _CConvBackropFilterCPU(TOut *filter_backprop, const std::vector< int > &filter_dims, size_t num_out, const TReal *out_positions, size_t num_inp, const TReal *inp_positions, const TFeat *inp_features, const TFeat *inp_importance, size_t neighbors_index_size, const TIndex *neighbors_index, const TFeat *neighbors_importance, const int64_t *neighbors_row_splits, const TReal *extents, const TReal *offsets, const TFeat *out_features_gradient, bool normalize)
Definition: ContinuousConvBackpropFilter.h:31
cloudViewer::ml::impl::CConvBackpropFilterCPU
void CConvBackpropFilterCPU(TOut *filter_backprop, const std::vector< int > &filter_dims, size_t num_out, const TReal *out_positions, size_t num_inp, const TReal *inp_positions, const TFeat *inp_features, const TFeat *inp_importance, size_t neighbors_index_size, const TIndex *neighbors_index, const TFeat *neighbors_importance, const int64_t *neighbors_row_splits, const TReal *extents, const TReal *offsets, const TFeat *out_features_gradient, InterpolationMode interpolation, CoordinateMapping coordinate_mapping, bool align_corners, bool individual_extent, bool isotropic_extent, bool normalize)
Definition: ContinuousConvBackpropFilter.h:293
cloudViewer::ml::impl::ComputeFilterCoordinates
void ComputeFilterCoordinates(Eigen::Array< T, VECSIZE, 1 > &x, Eigen::Array< T, VECSIZE, 1 > &y, Eigen::Array< T, VECSIZE, 1 > &z, const Eigen::Array< int, 3, 1 > &filter_size, const Eigen::Array< T, VECSIZE, 3 > &inv_extents, const Eigen::Array< T, 3, 1 > &offset)
Definition: CoordinateTransformation.h:107
cloudViewer
Generic file read and write utility for python interface.
Definition: AutoSegmentationTools.h:16
cloudViewer::ml::impl::InterpolationVec
Class for computing interpolation weights.
Definition: CoordinateTransformation.h:185