SplineInterpolation2D.h
Go to the documentation of this file.
1#ifndef SHARK_CORE_IMAGES_CPU_SPLINE_INTERPOLATION_2D_H
2#define SHARK_CORE_IMAGES_CPU_SPLINE_INTERPOLATION_2D_H
3
4#include <shark/LinAlg/Base.h>
5#include <shark/Core/Shape.h>
6#include <shark/Core/OpenMP.h>
7namespace shark{
8namespace image{
9template<class T>
10void splineInterpolation2D(
11 blas::dense_matrix_adaptor<T const, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> images,
12 Shape const& shape,
13 blas::dense_matrix_adaptor<T const, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> points, std::size_t pointsPerImage,
14 blas::dense_matrix_adaptor<T, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> values
15){
16 std::size_t height = shape[0];
17 std::size_t width = shape[1];
18 std::size_t numChannels = shape[2];
19 std::size_t stride = (pointsPerImage == points.size1())? 0: pointsPerImage;
20
21 values.clear();
22
23 SHARK_PARALLEL_FOR(int im = 0; im < (int)images.size1(); ++im){
24 auto image = to_matrix(row(images,im), width * height, numChannels);
25 auto v = to_matrix(row(values,im), pointsPerImage, numChannels);
26 for(std::size_t p = 0; p != pointsPerImage; ++p){
27 std::size_t pointIdx = p + im * stride;
28 using std::min; using std::max;
29 T basex = std::floor(points(pointIdx,1) * width);
30 T t=points(pointIdx, 1) * width - basex;
31 T t2=t*t;
32 T t3=t2*t;
33
34 T x[4]={-t3+3*t2-3*t+1, 3*t3-6*t2+4, -3*t3+3*t2+3*t+1, t3};
35 std::size_t px[4]={
36 (std::size_t)min<T>(max<T>(basex-1,0),width-1),
37 (std::size_t)max<T>(min<T>(basex,width-1),0),
38 (std::size_t)min<T>(max<T>(basex+1,0),width-1),
39 (std::size_t) max<T>(min<T>(basex+2,width-1),0)
40 };
41
42 T basey = std::floor(points(pointIdx,0) * height);
43 t=points(pointIdx, 0) * height - basey;
44 t2=t*t;
45 t3=t2*t;
46
47 T y[4]={-t3+3*t2-3*t+1, 3*t3-6*t2+4, -3*t3+3*t2+3*t+1, t3};
48 std::size_t py[4]={
49 (std::size_t)min<T>(max<T>(basey-1,0),height-1),
50 (std::size_t)max<T>(min<T>(basey,height-1),0),
51 (std::size_t)min<T>(max<T>(basey+1,0),height-1),
52 (std::size_t) max<T>(min<T>(basey+2,height-1),0)
53 };
54 for(std::size_t k = 0;k < 4; k++){
55 for(std::size_t l = 0;l < 4; l++){
56 std::size_t index = px[l] + width * py[k];
57 noalias(row(v,p)) += (x[l]*y[k]/36) * row(image,index);
58 }
59 }
60 }
61 }
62}
63
64template<class T>
65void splineInterpolation2D(
66 blas::dense_matrix_adaptor<T const, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> images,
67 Shape const& shape,
68 blas::dense_matrix_adaptor<T const, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> points, std::size_t pointsPerImage,
69 blas::dense_matrix_adaptor<T, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> values,
70 blas::dense_matrix_adaptor<T, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> valuesdx,
71 blas::dense_matrix_adaptor<T, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> valuesdy
72){
73 std::size_t height = shape[0];
74 std::size_t width = shape[1];
75 std::size_t numChannels = shape[2];
76 std::size_t stride = (pointsPerImage == points.size1())? 0: pointsPerImage;
77
78 values.clear();
79 valuesdx.clear();
80 valuesdy.clear();
81 SHARK_PARALLEL_FOR(int im = 0; im < (int)images.size1(); ++im){
82 auto image = to_matrix(row(images,im), width * height, numChannels);
83 auto v = to_matrix(row(values,im), pointsPerImage, numChannels);
84 auto vdx = to_matrix(row(valuesdx,im), pointsPerImage, numChannels);
85 auto vdy = to_matrix(row(valuesdy,im), pointsPerImage, numChannels);
86 for(std::size_t p = 0; p != points.size1(); ++p){
87 std::size_t pointIdx = p + im * stride;
88 using std::min; using std::max;
89 T basex = std::floor(points(pointIdx,1) * width);
90 T t=points(pointIdx, 1) * width - basex;
91 T t2=t*t;
92 T t3=t2*t;
93
94 T x[4]={-t3+3*t2-3*t+1, 3*t3-6*t2+4, -3*t3+3*t2+3*t+1, t3};
95 T dxdt[4]={-3*t2+6*t-3, 9*t2-12*t, -9*t2+6*t+3, 3*t2};
96 std::size_t px[4]={
97 (std::size_t)min<T>(max<T>(basex-1,0),width-1),
98 (std::size_t)max<T>(min<T>(basex,width-1),0),
99 (std::size_t)min<T>(max<T>(basex+1,0),width-1),
100 (std::size_t) max<T>(min<T>(basex+2,width-1),0)
101 };
102
103 T basey = std::floor(points(pointIdx,0) * height);
104 t=points(pointIdx, 0) * height - basey;
105 t2=t*t;
106 t3=t2*t;
107
108 T y[4]={-t3+3*t2-3*t+1, 3*t3-6*t2+4, -3*t3+3*t2+3*t+1, t3};
109 T dydt[4]={-3*t2+6*t-3, 9*t2-12*t, -9*t2+6*t+3, 3*t2};
110 std::size_t py[4]={
111 (std::size_t)min<T>(max<T>(basey-1,0),height-1),
112 (std::size_t)max<T>(min<T>(basey,height-1),0),
113 (std::size_t)min<T>(max<T>(basey+1,0),height-1),
114 (std::size_t) max<T>(min<T>(basey+2,height-1),0)
115 };
116
117 for(std::size_t k = 0;k < 4; k++){
118 for(std::size_t l = 0;l < 4; l++){
119 std::size_t index = px[l] + width * py[k];
120 T vdI =x[l]*y[k]/36;
121 T dfdx1 = dxdt[l] * y[k]/36;
122 T dfdx2 = x[l] * dydt[k]/36;
123
124 noalias(row(v,p)) += vdI * row(image,index);
125 noalias(row(vdx,p)) += dfdx1 * row(image,index);
126 noalias(row(vdy,p)) += dfdx2 * row(image,index);
127 }
128 }
129 }
130 }
131}
132template<class T>
133void splineInterpolation2DDerivative(
134 blas::dense_matrix_adaptor<T const, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> imageDerivatives,
135 Shape const& shape,
136 blas::dense_matrix_adaptor<T const, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> points, std::size_t pointsPerImage,
137 blas::dense_matrix_adaptor<T, blas::row_major, blas::continuous_dense_tag, blas::cpu_tag> results
138){
139 std::size_t height = shape[0];
140 std::size_t width = shape[1];
141 std::size_t numChannels = shape[2];
142 std::size_t stride = (pointsPerImage == points.size1())? 0: pointsPerImage;
143
144 results.clear();
145
146 SHARK_PARALLEL_FOR(int im = 0; im < (int)imageDerivatives.size1(); ++im){
147 auto imageDer = to_matrix(row(imageDerivatives,im), pointsPerImage, numChannels);
148 auto result = to_matrix(row(results,im), width * height, numChannels);
149 for(std::size_t p = 0; p != points.size1(); ++p){
150 std::size_t pointIdx = p + im * stride;
151 using std::min; using std::max;
152 T basex = std::floor(points(pointIdx,1) * width);
153 T t=points(pointIdx, 1) * width - basex;
154 T t2=t*t;
155 T t3=t2*t;
156
157 T x[4]={-t3+3*t2-3*t+1, 3*t3-6*t2+4, -3*t3+3*t2+3*t+1, t3};
158 std::size_t px[4]={
159 (std::size_t)min<T>(max<T>(basex-1,0),width-1),
160 (std::size_t)max<T>(min<T>(basex,width-1),0),
161 (std::size_t)min<T>(max<T>(basex+1,0),width-1),
162 (std::size_t) max<T>(min<T>(basex+2,width-1),0)
163 };
164
165 T basey = std::floor(points(pointIdx,0) * height);
166 t=points(pointIdx, 0) * height - basey;
167 t2=t*t;
168 t3=t2*t;
169
170 T y[4]={-t3+3*t2-3*t+1, 3*t3-6*t2+4, -3*t3+3*t2+3*t+1, t3};
171 std::size_t py[4]={
172 (std::size_t)min<T>(max<T>(basey-1,0),height-1),
173 (std::size_t)max<T>(min<T>(basey,height-1),0),
174 (std::size_t)min<T>(max<T>(basey+1,0),height-1),
175 (std::size_t) max<T>(min<T>(basey+2,height-1),0)
176 };
177 for(std::size_t k = 0;k < 4; k++){
178 for(std::size_t l = 0;l < 4; l++){
179 std::size_t index = px[l] + width * py[k];
180 T vdI = x[l]*y[k]/36;
181 noalias(row(result,index)) += vdI * row(imageDer,p);
182 }
183 }
184 }
185 }
186}
187
188}}
189
190#endif