Krita Source Code Documentation
Loading...
Searching...
No Matches
multi_bspline_eval_std_c.h
Go to the documentation of this file.
1
2// einspline: a library for creating and evaluating B-splines //
3// Copyright (C) 2007 Kenneth P. Esler, Jr. //
4// //
5// This program is free software; you can redistribute it and/or modify //
6// it under the terms of the GNU General Public License as published by //
7// the Free Software Foundation; either version 2 of the License, or //
8// (at your option) any later version. //
9// //
10// This program is distributed in the hope that it will be useful, //
11// but WITHOUT ANY WARRANTY; without even the implied warranty of //
12// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
13// GNU General Public License for more details. //
14// //
15// You should have received a copy of the GNU General Public License //
16// along with this program; if not, write to the Free Software //
17// Foundation, Inc., 51 Franklin Street, Fifth Floor, //
18// Boston, MA 02110-1301 USA //
20
21#ifndef MULTI_BSPLINE_EVAL_STD_C_H
22#define MULTI_BSPLINE_EVAL_STD_C_H
23
24#include <math.h>
25#include <stdio.h>
26#include "multi_bspline_structs.h"
27
28extern const float* restrict Af;
29extern const float* restrict dAf;
30extern const float* restrict d2Af;
31extern const float* restrict d3Af;
32
33/************************************************************/
34/* 1D float-precision, complex evaulation functions */
35/************************************************************/
36inline void
37eval_multi_UBspline_1d_c (multi_UBspline_1d_c *spline,
38 double x,
39 complex_float* restrict vals)
40{
41 x -= spline->x_grid.start;
42 float ux = x*spline->x_grid.delta_inv;
43 float ipartx, tx;
44 tx = modff (ux, &ipartx); int ix = (int) ipartx;
45
46 float tpx[4], a[4];
47 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
48 complex_float* restrict coefs = spline->coefs;
49
50 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
51 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
52 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
53 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
54
55 int xs = spline->x_stride;
56
57 for (int n=0; n<spline->num_splines; n++)
58 vals[n] = 0.0;
59
60 for (int i=0; i<4; i++) {
61 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs);
62 for (int n=0; n<spline->num_splines; n++)
63 vals[n] += a[i] * coefs[n];
64 }
65}
66
67
68
69inline void
70eval_multi_UBspline_1d_c_vg (multi_UBspline_1d_c *spline,
71 double x,
72 complex_float* restrict vals,
73 complex_float* restrict grads)
74{
75 x -= spline->x_grid.start;
76 float ux = x*spline->x_grid.delta_inv;
77 float ipartx, tx;
78 tx = modff (ux, &ipartx); int ix = (int) ipartx;
79
80 float tpx[4], a[4], da[4];
81 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
82 complex_float* restrict coefs = spline->coefs;
83
84 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
85 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
86 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
87 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
88 da[0] = (dAf[ 0]*tpx[0] + dAf[ 1]*tpx[1] + dAf[ 2]*tpx[2] + dAf[ 3]*tpx[3]);
89 da[1] = (dAf[ 4]*tpx[0] + dAf[ 5]*tpx[1] + dAf[ 6]*tpx[2] + dAf[ 7]*tpx[3]);
90 da[2] = (dAf[ 8]*tpx[0] + dAf[ 9]*tpx[1] + dAf[10]*tpx[2] + dAf[11]*tpx[3]);
91 da[3] = (dAf[12]*tpx[0] + dAf[13]*tpx[1] + dAf[14]*tpx[2] + dAf[15]*tpx[3]);
92
93 int xs = spline->x_stride;
94
95 for (int n=0; n<spline->num_splines; n++) {
96 vals[n] = 0.0;
97 grads[n] = 0.0;
98 }
99
100 for (int i=0; i<4; i++) {
101 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs);
102 for (int n=0; n<spline->num_splines; n++) {
103 vals[n] += a[i] * coefs[n];
104 grads[n] += da[i] * coefs[n];
105 }
106 }
107
108 float dxInv = spline->x_grid.delta_inv;
109 for (int n=0; n<spline->num_splines; n++)
110 grads[n] *= dxInv;
111}
112
113
114inline void
115eval_multi_UBspline_1d_c_vgl (multi_UBspline_1d_c *spline,
116 double x,
117 complex_float* restrict vals,
118 complex_float* restrict grads,
119 complex_float* restrict lapl)
120{
121 x -= spline->x_grid.start;
122 float ux = x*spline->x_grid.delta_inv;
123 float ipartx, tx;
124 tx = modff (ux, &ipartx); int ix = (int) ipartx;
125
126 float tpx[4], a[4], da[4], d2a[4];
127 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
128 complex_float* restrict coefs = spline->coefs;
129
130 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
131 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
132 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
133 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
134 da[0] = (dAf[ 0]*tpx[0] + dAf[ 1]*tpx[1] + dAf[ 2]*tpx[2] + dAf[ 3]*tpx[3]);
135 da[1] = (dAf[ 4]*tpx[0] + dAf[ 5]*tpx[1] + dAf[ 6]*tpx[2] + dAf[ 7]*tpx[3]);
136 da[2] = (dAf[ 8]*tpx[0] + dAf[ 9]*tpx[1] + dAf[10]*tpx[2] + dAf[11]*tpx[3]);
137 da[3] = (dAf[12]*tpx[0] + dAf[13]*tpx[1] + dAf[14]*tpx[2] + dAf[15]*tpx[3]);
138 d2a[0] = (d2Af[ 0]*tpx[0] + d2Af[ 1]*tpx[1] + d2Af[ 2]*tpx[2] + d2Af[ 3]*tpx[3]);
139 d2a[1] = (d2Af[ 4]*tpx[0] + d2Af[ 5]*tpx[1] + d2Af[ 6]*tpx[2] + d2Af[ 7]*tpx[3]);
140 d2a[2] = (d2Af[ 8]*tpx[0] + d2Af[ 9]*tpx[1] + d2Af[10]*tpx[2] + d2Af[11]*tpx[3]);
141 d2a[3] = (d2Af[12]*tpx[0] + d2Af[13]*tpx[1] + d2Af[14]*tpx[2] + d2Af[15]*tpx[3]);
142
143 int xs = spline->x_stride;
144
145 for (int n=0; n<spline->num_splines; n++) {
146 vals[n] = 0.0;
147 grads[n] = 0.0;
148 lapl[n] = 0.0;
149 }
150
151 for (int i=0; i<4; i++) {
152 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs);
153 for (int n=0; n<spline->num_splines; n++) {
154 vals[n] += a[i] * coefs[n];
155 grads[n] += da[i] * coefs[n];
156 lapl[n] += d2a[i] * coefs[n];
157 }
158 }
159
160 float dxInv = spline->x_grid.delta_inv;
161 for (int n=0; n<spline->num_splines; n++) {
162 grads[n] *= dxInv;
163 lapl [n] *= dxInv*dxInv;
164 }
165}
166
167
168inline void
169eval_multi_UBspline_1d_c_vgh (multi_UBspline_1d_c *spline,
170 double x,
171 complex_float* restrict vals,
172 complex_float* restrict grads,
173 complex_float* restrict hess)
174{
175 eval_multi_UBspline_1d_c_vgl (spline, x, vals, grads, hess);
176}
177
178
179/************************************************************/
180/* 2D float-precision, complex evaulation functions */
181/************************************************************/
182inline void
183eval_multi_UBspline_2d_c (multi_UBspline_2d_c *spline,
184 double x, double y,
185 complex_float* restrict vals)
186{
187 x -= spline->x_grid.start;
188 y -= spline->y_grid.start;
189 float ux = x*spline->x_grid.delta_inv;
190 float uy = y*spline->y_grid.delta_inv;
191 float ipartx, iparty, tx, ty;
192 tx = modff (ux, &ipartx); int ix = (int) ipartx;
193 ty = modff (uy, &iparty); int iy = (int) iparty;
194
195 float tpx[4], tpy[4], a[4], b[4];
196 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
197 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
198 complex_float* restrict coefs = spline->coefs;
199
200 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
201 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
202 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
203 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
204
205 b[0] = (Af[ 0]*tpy[0] + Af[ 1]*tpy[1] + Af[ 2]*tpy[2] + Af[ 3]*tpy[3]);
206 b[1] = (Af[ 4]*tpy[0] + Af[ 5]*tpy[1] + Af[ 6]*tpy[2] + Af[ 7]*tpy[3]);
207 b[2] = (Af[ 8]*tpy[0] + Af[ 9]*tpy[1] + Af[10]*tpy[2] + Af[11]*tpy[3]);
208 b[3] = (Af[12]*tpy[0] + Af[13]*tpy[1] + Af[14]*tpy[2] + Af[15]*tpy[3]);
209
210 int xs = spline->x_stride;
211 int ys = spline->y_stride;
212
213 for (int n=0; n<spline->num_splines; n++)
214 vals[n] = 0.0;
215
216 for (int i=0; i<4; i++)
217 for (int j=0; j<4; j++) {
218 float prefactor = a[i]*b[j];
219 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys);
220 for (int n=0; n<spline->num_splines; n++)
221 vals[n] += prefactor*coefs[n];
222 }
223}
224
225
226inline void
227eval_multi_UBspline_2d_c_vg (multi_UBspline_2d_c *spline,
228 double x, double y,
229 complex_float* restrict vals,
230 complex_float* restrict grads)
231{
232 x -= spline->x_grid.start;
233 y -= spline->y_grid.start;
234 float ux = x*spline->x_grid.delta_inv;
235 float uy = y*spline->y_grid.delta_inv;
236 float ipartx, iparty, tx, ty;
237 tx = modff (ux, &ipartx); int ix = (int) ipartx;
238 ty = modff (uy, &iparty); int iy = (int) iparty;
239
240 float tpx[4], tpy[4], a[4], b[4], da[4], db[4];
241 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
242 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
243 complex_float* restrict coefs = spline->coefs;
244
245 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
246 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
247 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
248 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
249 da[0] = (dAf[ 0]*tpx[0] + dAf[ 1]*tpx[1] + dAf[ 2]*tpx[2] + dAf[ 3]*tpx[3]);
250 da[1] = (dAf[ 4]*tpx[0] + dAf[ 5]*tpx[1] + dAf[ 6]*tpx[2] + dAf[ 7]*tpx[3]);
251 da[2] = (dAf[ 8]*tpx[0] + dAf[ 9]*tpx[1] + dAf[10]*tpx[2] + dAf[11]*tpx[3]);
252 da[3] = (dAf[12]*tpx[0] + dAf[13]*tpx[1] + dAf[14]*tpx[2] + dAf[15]*tpx[3]);
253
254 b[0] = (Af[ 0]*tpy[0] + Af[ 1]*tpy[1] + Af[ 2]*tpy[2] + Af[ 3]*tpy[3]);
255 b[1] = (Af[ 4]*tpy[0] + Af[ 5]*tpy[1] + Af[ 6]*tpy[2] + Af[ 7]*tpy[3]);
256 b[2] = (Af[ 8]*tpy[0] + Af[ 9]*tpy[1] + Af[10]*tpy[2] + Af[11]*tpy[3]);
257 b[3] = (Af[12]*tpy[0] + Af[13]*tpy[1] + Af[14]*tpy[2] + Af[15]*tpy[3]);
258 db[0] = (dAf[ 0]*tpy[0] + dAf[ 1]*tpy[1] + dAf[ 2]*tpy[2] + dAf[ 3]*tpy[3]);
259 db[1] = (dAf[ 4]*tpy[0] + dAf[ 5]*tpy[1] + dAf[ 6]*tpy[2] + dAf[ 7]*tpy[3]);
260 db[2] = (dAf[ 8]*tpy[0] + dAf[ 9]*tpy[1] + dAf[10]*tpy[2] + dAf[11]*tpy[3]);
261 db[3] = (dAf[12]*tpy[0] + dAf[13]*tpy[1] + dAf[14]*tpy[2] + dAf[15]*tpy[3]);
262
263 int xs = spline->x_stride;
264 int ys = spline->y_stride;
265
266 for (int n=0; n<spline->num_splines; n++) {
267 vals[n] = 0.0;
268 grads[2*n+0] = grads[2*n+1] = grads[2*n+2] = 0.0;
269 }
270
271 for (int i=0; i<4; i++)
272 for (int j=0; j<4; j++) {
273 float ab = a[i]*b[j];
274 float dab[2];
275 dab[0] = da[i]* b[j];
276 dab[1] = a[i]*db[j];
277
278 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys);
279 for (int n=0; n<spline->num_splines; n++) {
280 vals [n] += ab *coefs[n];
281 grads[2*n+0] += dab[0]*coefs[n];
282 grads[2*n+1] += dab[1]*coefs[n];
283 }
284 }
285
286 float dxInv = spline->x_grid.delta_inv;
287 float dyInv = spline->y_grid.delta_inv;
288 for (int n=0; n<spline->num_splines; n++) {
289 grads[2*n+0] *= dxInv;
290 grads[2*n+1] *= dyInv;
291 }
292}
293
294inline void
295eval_multi_UBspline_2d_c_vgl (multi_UBspline_2d_c *spline,
296 double x, double y,
297 complex_float* restrict vals,
298 complex_float* restrict grads,
299 complex_float* restrict lapl)
300{
301 x -= spline->x_grid.start;
302 y -= spline->y_grid.start;
303 float ux = x*spline->x_grid.delta_inv;
304 float uy = y*spline->y_grid.delta_inv;
305 float ipartx, iparty, tx, ty;
306 tx = modff (ux, &ipartx); int ix = (int) ipartx;
307 ty = modff (uy, &iparty); int iy = (int) iparty;
308
309 float tpx[4], tpy[4], a[4], b[4], da[4], db[4], d2a[4], d2b[4];
310 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
311 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
312 complex_float* restrict coefs = spline->coefs;
313
314 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
315 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
316 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
317 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
318 da[0] = (dAf[ 0]*tpx[0] + dAf[ 1]*tpx[1] + dAf[ 2]*tpx[2] + dAf[ 3]*tpx[3]);
319 da[1] = (dAf[ 4]*tpx[0] + dAf[ 5]*tpx[1] + dAf[ 6]*tpx[2] + dAf[ 7]*tpx[3]);
320 da[2] = (dAf[ 8]*tpx[0] + dAf[ 9]*tpx[1] + dAf[10]*tpx[2] + dAf[11]*tpx[3]);
321 da[3] = (dAf[12]*tpx[0] + dAf[13]*tpx[1] + dAf[14]*tpx[2] + dAf[15]*tpx[3]);
322 d2a[0] = (d2Af[ 0]*tpx[0] + d2Af[ 1]*tpx[1] + d2Af[ 2]*tpx[2] + d2Af[ 3]*tpx[3]);
323 d2a[1] = (d2Af[ 4]*tpx[0] + d2Af[ 5]*tpx[1] + d2Af[ 6]*tpx[2] + d2Af[ 7]*tpx[3]);
324 d2a[2] = (d2Af[ 8]*tpx[0] + d2Af[ 9]*tpx[1] + d2Af[10]*tpx[2] + d2Af[11]*tpx[3]);
325 d2a[3] = (d2Af[12]*tpx[0] + d2Af[13]*tpx[1] + d2Af[14]*tpx[2] + d2Af[15]*tpx[3]);
326
327 b[0] = (Af[ 0]*tpy[0] + Af[ 1]*tpy[1] + Af[ 2]*tpy[2] + Af[ 3]*tpy[3]);
328 b[1] = (Af[ 4]*tpy[0] + Af[ 5]*tpy[1] + Af[ 6]*tpy[2] + Af[ 7]*tpy[3]);
329 b[2] = (Af[ 8]*tpy[0] + Af[ 9]*tpy[1] + Af[10]*tpy[2] + Af[11]*tpy[3]);
330 b[3] = (Af[12]*tpy[0] + Af[13]*tpy[1] + Af[14]*tpy[2] + Af[15]*tpy[3]);
331 db[0] = (dAf[ 0]*tpy[0] + dAf[ 1]*tpy[1] + dAf[ 2]*tpy[2] + dAf[ 3]*tpy[3]);
332 db[1] = (dAf[ 4]*tpy[0] + dAf[ 5]*tpy[1] + dAf[ 6]*tpy[2] + dAf[ 7]*tpy[3]);
333 db[2] = (dAf[ 8]*tpy[0] + dAf[ 9]*tpy[1] + dAf[10]*tpy[2] + dAf[11]*tpy[3]);
334 db[3] = (dAf[12]*tpy[0] + dAf[13]*tpy[1] + dAf[14]*tpy[2] + dAf[15]*tpy[3]);
335 d2b[0] = (d2Af[ 0]*tpy[0] + d2Af[ 1]*tpy[1] + d2Af[ 2]*tpy[2] + d2Af[ 3]*tpy[3]);
336 d2b[1] = (d2Af[ 4]*tpy[0] + d2Af[ 5]*tpy[1] + d2Af[ 6]*tpy[2] + d2Af[ 7]*tpy[3]);
337 d2b[2] = (d2Af[ 8]*tpy[0] + d2Af[ 9]*tpy[1] + d2Af[10]*tpy[2] + d2Af[11]*tpy[3]);
338 d2b[3] = (d2Af[12]*tpy[0] + d2Af[13]*tpy[1] + d2Af[14]*tpy[2] + d2Af[15]*tpy[3]);
339
340 int xs = spline->x_stride;
341 int ys = spline->y_stride;
342
343 complex_float lapl2[2*spline->num_splines];
344 for (int n=0; n<spline->num_splines; n++) {
345 vals[n] = 0.0;
346 grads[2*n+0] = grads[2*n+1] = 0.0;
347 lapl2[2*n+0] = lapl2[2*n+1] = 0.0;
348 }
349
350 for (int i=0; i<4; i++)
351 for (int j=0; j<4; j++) {
352 float ab = a[i]*b[j];
353 float dab[2], d2ab[2];
354 dab[0] = da[i]* b[j];
355 dab[1] = a[i]*db[j];
356 d2ab[0] = d2a[i]* b[j];
357 d2ab[1] = a[i]*d2b[j];
358
359 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys);
360 for (int n=0; n<spline->num_splines; n++) {
361 vals[n] += ab *coefs[n];
362 grads[2*n+0] += dab[0]*coefs[n];
363 grads[2*n+1] += dab[1]*coefs[n];
364 lapl2[2*n+0] += d2ab[0]*coefs[n];
365 lapl2[2*n+1] += d2ab[1]*coefs[n];
366 }
367 }
368
369 float dxInv = spline->x_grid.delta_inv;
370 float dyInv = spline->y_grid.delta_inv;
371 for (int n=0; n<spline->num_splines; n++) {
372 grads[2*n+0] *= dxInv;
373 grads[2*n+1] *= dyInv;
374 lapl2[2*n+0] *= dxInv*dxInv;
375 lapl2[2*n+1] *= dyInv*dyInv;
376 lapl[n] = lapl2[2*n+0] + lapl2[2*n+1];
377 }
378}
379
380inline void
381eval_multi_UBspline_2d_c_vgh (multi_UBspline_2d_c *spline,
382 double x, double y,
383 complex_float* restrict vals,
384 complex_float* restrict grads,
385 complex_float* restrict hess)
386{
387 x -= spline->x_grid.start;
388 y -= spline->y_grid.start;
389 float ux = x*spline->x_grid.delta_inv;
390 float uy = y*spline->y_grid.delta_inv;
391 float ipartx, iparty, tx, ty;
392 tx = modff (ux, &ipartx); int ix = (int) ipartx;
393 ty = modff (uy, &iparty); int iy = (int) iparty;
394
395 float tpx[4], tpy[4], a[4], b[4], da[4], db[4], d2a[4], d2b[4];
396 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
397 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
398 complex_float* restrict coefs = spline->coefs;
399
400 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
401 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
402 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
403 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
404 da[0] = (dAf[ 0]*tpx[0] + dAf[ 1]*tpx[1] + dAf[ 2]*tpx[2] + dAf[ 3]*tpx[3]);
405 da[1] = (dAf[ 4]*tpx[0] + dAf[ 5]*tpx[1] + dAf[ 6]*tpx[2] + dAf[ 7]*tpx[3]);
406 da[2] = (dAf[ 8]*tpx[0] + dAf[ 9]*tpx[1] + dAf[10]*tpx[2] + dAf[11]*tpx[3]);
407 da[3] = (dAf[12]*tpx[0] + dAf[13]*tpx[1] + dAf[14]*tpx[2] + dAf[15]*tpx[3]);
408 d2a[0] = (d2Af[ 0]*tpx[0] + d2Af[ 1]*tpx[1] + d2Af[ 2]*tpx[2] + d2Af[ 3]*tpx[3]);
409 d2a[1] = (d2Af[ 4]*tpx[0] + d2Af[ 5]*tpx[1] + d2Af[ 6]*tpx[2] + d2Af[ 7]*tpx[3]);
410 d2a[2] = (d2Af[ 8]*tpx[0] + d2Af[ 9]*tpx[1] + d2Af[10]*tpx[2] + d2Af[11]*tpx[3]);
411 d2a[3] = (d2Af[12]*tpx[0] + d2Af[13]*tpx[1] + d2Af[14]*tpx[2] + d2Af[15]*tpx[3]);
412
413 b[0] = (Af[ 0]*tpy[0] + Af[ 1]*tpy[1] + Af[ 2]*tpy[2] + Af[ 3]*tpy[3]);
414 b[1] = (Af[ 4]*tpy[0] + Af[ 5]*tpy[1] + Af[ 6]*tpy[2] + Af[ 7]*tpy[3]);
415 b[2] = (Af[ 8]*tpy[0] + Af[ 9]*tpy[1] + Af[10]*tpy[2] + Af[11]*tpy[3]);
416 b[3] = (Af[12]*tpy[0] + Af[13]*tpy[1] + Af[14]*tpy[2] + Af[15]*tpy[3]);
417 db[0] = (dAf[ 0]*tpy[0] + dAf[ 1]*tpy[1] + dAf[ 2]*tpy[2] + dAf[ 3]*tpy[3]);
418 db[1] = (dAf[ 4]*tpy[0] + dAf[ 5]*tpy[1] + dAf[ 6]*tpy[2] + dAf[ 7]*tpy[3]);
419 db[2] = (dAf[ 8]*tpy[0] + dAf[ 9]*tpy[1] + dAf[10]*tpy[2] + dAf[11]*tpy[3]);
420 db[3] = (dAf[12]*tpy[0] + dAf[13]*tpy[1] + dAf[14]*tpy[2] + dAf[15]*tpy[3]);
421 d2b[0] = (d2Af[ 0]*tpy[0] + d2Af[ 1]*tpy[1] + d2Af[ 2]*tpy[2] + d2Af[ 3]*tpy[3]);
422 d2b[1] = (d2Af[ 4]*tpy[0] + d2Af[ 5]*tpy[1] + d2Af[ 6]*tpy[2] + d2Af[ 7]*tpy[3]);
423 d2b[2] = (d2Af[ 8]*tpy[0] + d2Af[ 9]*tpy[1] + d2Af[10]*tpy[2] + d2Af[11]*tpy[3]);
424 d2b[3] = (d2Af[12]*tpy[0] + d2Af[13]*tpy[1] + d2Af[14]*tpy[2] + d2Af[15]*tpy[3]);
425
426 int xs = spline->x_stride;
427 int ys = spline->y_stride;
428
429 for (int n=0; n<spline->num_splines; n++) {
430 vals[n] = 0.0;
431 grads[2*n+0] = grads[2*n+1] = 0.0;
432 for (int i=0; i<4; i++)
433 hess[4*n+i] = 0.0;
434 }
435
436 for (int i=0; i<4; i++)
437 for (int j=0; j<4; j++){
438 float ab = a[i]*b[j];
439 float dab[2], d2ab[3];
440 dab[0] = da[i]* b[j];
441 dab[1] = a[i]*db[j];
442 d2ab[0] = d2a[i] * b[j];
443 d2ab[1] = da[i] * db[j];
444 d2ab[2] = a[i] * d2b[j];
445
446 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys);
447 for (int n=0; n<spline->num_splines; n++) {
448 vals[n] += ab *coefs[n];
449 grads[2*n+0] += dab[0]*coefs[n];
450 grads[2*n+1] += dab[1]*coefs[n];
451 hess [4*n+0] += d2ab[0]*coefs[n];
452 hess [4*n+1] += d2ab[1]*coefs[n];
453 hess [4*n+3] += d2ab[2]*coefs[n];
454 }
455 }
456
457 float dxInv = spline->x_grid.delta_inv;
458 float dyInv = spline->y_grid.delta_inv;
459 for (int n=0; n<spline->num_splines; n++) {
460 grads[2*n+0] *= dxInv;
461 grads[2*n+1] *= dyInv;
462 hess[4*n+0] *= dxInv*dxInv;
463 hess[4*n+1] *= dxInv*dyInv;
464 hess[4*n+3] *= dyInv*dyInv;
465 // Copy hessian elements into lower half of 3x3 matrix
466 hess[4*n+2] = hess[4*n+1];
467 }
468}
469
470
471/************************************************************/
472/* 3D float-precision, complex evaulation functions */
473/************************************************************/
474inline void
475eval_multi_UBspline_3d_c (multi_UBspline_3d_c *spline,
476 double x, double y, double z,
477 complex_float* restrict vals)
478{
479 x -= spline->x_grid.start;
480 y -= spline->y_grid.start;
481 z -= spline->z_grid.start;
482 float ux = x*spline->x_grid.delta_inv;
483 float uy = y*spline->y_grid.delta_inv;
484 float uz = z*spline->z_grid.delta_inv;
485 float ipartx, iparty, ipartz, tx, ty, tz;
486 tx = modff (ux, &ipartx); int ix = (int) ipartx;
487 ty = modff (uy, &iparty); int iy = (int) iparty;
488 tz = modff (uz, &ipartz); int iz = (int) ipartz;
489
490 float tpx[4], tpy[4], tpz[4], a[4], b[4], c[4];
491 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
492 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
493 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
494 complex_float* restrict coefs = spline->coefs;
495
496 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
497 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
498 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
499 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
500
501 b[0] = (Af[ 0]*tpy[0] + Af[ 1]*tpy[1] + Af[ 2]*tpy[2] + Af[ 3]*tpy[3]);
502 b[1] = (Af[ 4]*tpy[0] + Af[ 5]*tpy[1] + Af[ 6]*tpy[2] + Af[ 7]*tpy[3]);
503 b[2] = (Af[ 8]*tpy[0] + Af[ 9]*tpy[1] + Af[10]*tpy[2] + Af[11]*tpy[3]);
504 b[3] = (Af[12]*tpy[0] + Af[13]*tpy[1] + Af[14]*tpy[2] + Af[15]*tpy[3]);
505
506 c[0] = (Af[ 0]*tpz[0] + Af[ 1]*tpz[1] + Af[ 2]*tpz[2] + Af[ 3]*tpz[3]);
507 c[1] = (Af[ 4]*tpz[0] + Af[ 5]*tpz[1] + Af[ 6]*tpz[2] + Af[ 7]*tpz[3]);
508 c[2] = (Af[ 8]*tpz[0] + Af[ 9]*tpz[1] + Af[10]*tpz[2] + Af[11]*tpz[3]);
509 c[3] = (Af[12]*tpz[0] + Af[13]*tpz[1] + Af[14]*tpz[2] + Af[15]*tpz[3]);
510
511 int xs = spline->x_stride;
512 int ys = spline->y_stride;
513 int zs = spline->z_stride;
514
515 for (int n=0; n<spline->num_splines; n++)
516 vals[n] = 0.0;
517
518 for (int i=0; i<4; i++)
519 for (int j=0; j<4; j++)
520 for (int k=0; k<4; k++) {
521 float prefactor = a[i]*b[j]*c[k];
522 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
523 for (int n=0; n<spline->num_splines; n++)
524 vals[n] += prefactor*coefs[n];
525 }
526}
527
528
529inline void
530eval_multi_UBspline_3d_c_vg (multi_UBspline_3d_c *spline,
531 double x, double y, double z,
532 complex_float* restrict vals,
533 complex_float* restrict grads)
534{
535 x -= spline->x_grid.start;
536 y -= spline->y_grid.start;
537 z -= spline->z_grid.start;
538 float ux = x*spline->x_grid.delta_inv;
539 float uy = y*spline->y_grid.delta_inv;
540 float uz = z*spline->z_grid.delta_inv;
541 float ipartx, iparty, ipartz, tx, ty, tz;
542 tx = modff (ux, &ipartx); int ix = (int) ipartx;
543 ty = modff (uy, &iparty); int iy = (int) iparty;
544 tz = modff (uz, &ipartz); int iz = (int) ipartz;
545
546 float tpx[4], tpy[4], tpz[4], a[4], b[4], c[4],
547 da[4], db[4], dc[4];
548 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
549 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
550 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
551 complex_float* restrict coefs = spline->coefs;
552
553 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
554 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
555 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
556 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
557 da[0] = (dAf[ 0]*tpx[0] + dAf[ 1]*tpx[1] + dAf[ 2]*tpx[2] + dAf[ 3]*tpx[3]);
558 da[1] = (dAf[ 4]*tpx[0] + dAf[ 5]*tpx[1] + dAf[ 6]*tpx[2] + dAf[ 7]*tpx[3]);
559 da[2] = (dAf[ 8]*tpx[0] + dAf[ 9]*tpx[1] + dAf[10]*tpx[2] + dAf[11]*tpx[3]);
560 da[3] = (dAf[12]*tpx[0] + dAf[13]*tpx[1] + dAf[14]*tpx[2] + dAf[15]*tpx[3]);
561
562 b[0] = (Af[ 0]*tpy[0] + Af[ 1]*tpy[1] + Af[ 2]*tpy[2] + Af[ 3]*tpy[3]);
563 b[1] = (Af[ 4]*tpy[0] + Af[ 5]*tpy[1] + Af[ 6]*tpy[2] + Af[ 7]*tpy[3]);
564 b[2] = (Af[ 8]*tpy[0] + Af[ 9]*tpy[1] + Af[10]*tpy[2] + Af[11]*tpy[3]);
565 b[3] = (Af[12]*tpy[0] + Af[13]*tpy[1] + Af[14]*tpy[2] + Af[15]*tpy[3]);
566 db[0] = (dAf[ 0]*tpy[0] + dAf[ 1]*tpy[1] + dAf[ 2]*tpy[2] + dAf[ 3]*tpy[3]);
567 db[1] = (dAf[ 4]*tpy[0] + dAf[ 5]*tpy[1] + dAf[ 6]*tpy[2] + dAf[ 7]*tpy[3]);
568 db[2] = (dAf[ 8]*tpy[0] + dAf[ 9]*tpy[1] + dAf[10]*tpy[2] + dAf[11]*tpy[3]);
569 db[3] = (dAf[12]*tpy[0] + dAf[13]*tpy[1] + dAf[14]*tpy[2] + dAf[15]*tpy[3]);
570
571 c[0] = (Af[ 0]*tpz[0] + Af[ 1]*tpz[1] + Af[ 2]*tpz[2] + Af[ 3]*tpz[3]);
572 c[1] = (Af[ 4]*tpz[0] + Af[ 5]*tpz[1] + Af[ 6]*tpz[2] + Af[ 7]*tpz[3]);
573 c[2] = (Af[ 8]*tpz[0] + Af[ 9]*tpz[1] + Af[10]*tpz[2] + Af[11]*tpz[3]);
574 c[3] = (Af[12]*tpz[0] + Af[13]*tpz[1] + Af[14]*tpz[2] + Af[15]*tpz[3]);
575 dc[0] = (dAf[ 0]*tpz[0] + dAf[ 1]*tpz[1] + dAf[ 2]*tpz[2] + dAf[ 3]*tpz[3]);
576 dc[1] = (dAf[ 4]*tpz[0] + dAf[ 5]*tpz[1] + dAf[ 6]*tpz[2] + dAf[ 7]*tpz[3]);
577 dc[2] = (dAf[ 8]*tpz[0] + dAf[ 9]*tpz[1] + dAf[10]*tpz[2] + dAf[11]*tpz[3]);
578 dc[3] = (dAf[12]*tpz[0] + dAf[13]*tpz[1] + dAf[14]*tpz[2] + dAf[15]*tpz[3]);
579
580 int xs = spline->x_stride;
581 int ys = spline->y_stride;
582 int zs = spline->z_stride;
583
584 for (int n=0; n<spline->num_splines; n++) {
585 vals[n] = 0.0;
586 grads[3*n+0] = grads[3*n+1] = grads[3*n+2] = 0.0;
587 }
588
589 for (int i=0; i<4; i++)
590 for (int j=0; j<4; j++)
591 for (int k=0; k<4; k++) {
592 float abc = a[i]*b[j]*c[k];
593 float dabc[3];
594 dabc[0] = da[i]* b[j]* c[k];
595 dabc[1] = a[i]*db[j]* c[k];
596 dabc[2] = a[i]* b[j]*dc[k];
597
598 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
599 for (int n=0; n<spline->num_splines; n++) {
600 vals[n] += abc *coefs[n];
601 grads[3*n+0] += dabc[0]*coefs[n];
602 grads[3*n+1] += dabc[1]*coefs[n];
603 grads[3*n+2] += dabc[2]*coefs[n];
604 }
605 }
606
607 float dxInv = spline->x_grid.delta_inv;
608 float dyInv = spline->y_grid.delta_inv;
609 float dzInv = spline->z_grid.delta_inv;
610 for (int n=0; n<spline->num_splines; n++) {
611 grads[3*n+0] *= dxInv;
612 grads[3*n+1] *= dyInv;
613 grads[3*n+2] *= dzInv;
614 }
615}
616
617
618
619inline void
620eval_multi_UBspline_3d_c_vgl (multi_UBspline_3d_c *spline,
621 double x, double y, double z,
622 complex_float* restrict vals,
623 complex_float* restrict grads,
624 complex_float* restrict lapl)
625{
626 x -= spline->x_grid.start;
627 y -= spline->y_grid.start;
628 z -= spline->z_grid.start;
629 float ux = x*spline->x_grid.delta_inv;
630 float uy = y*spline->y_grid.delta_inv;
631 float uz = z*spline->z_grid.delta_inv;
632 float ipartx, iparty, ipartz, tx, ty, tz;
633 tx = modff (ux, &ipartx); int ix = (int) ipartx;
634 ty = modff (uy, &iparty); int iy = (int) iparty;
635 tz = modff (uz, &ipartz); int iz = (int) ipartz;
636
637 float tpx[4], tpy[4], tpz[4], a[4], b[4], c[4],
638 da[4], db[4], dc[4], d2a[4], d2b[4], d2c[4];
639 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
640 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
641 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
642 complex_float* restrict coefs = spline->coefs;
643
644 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
645 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
646 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
647 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
648 da[0] = (dAf[ 0]*tpx[0] + dAf[ 1]*tpx[1] + dAf[ 2]*tpx[2] + dAf[ 3]*tpx[3]);
649 da[1] = (dAf[ 4]*tpx[0] + dAf[ 5]*tpx[1] + dAf[ 6]*tpx[2] + dAf[ 7]*tpx[3]);
650 da[2] = (dAf[ 8]*tpx[0] + dAf[ 9]*tpx[1] + dAf[10]*tpx[2] + dAf[11]*tpx[3]);
651 da[3] = (dAf[12]*tpx[0] + dAf[13]*tpx[1] + dAf[14]*tpx[2] + dAf[15]*tpx[3]);
652 d2a[0] = (d2Af[ 0]*tpx[0] + d2Af[ 1]*tpx[1] + d2Af[ 2]*tpx[2] + d2Af[ 3]*tpx[3]);
653 d2a[1] = (d2Af[ 4]*tpx[0] + d2Af[ 5]*tpx[1] + d2Af[ 6]*tpx[2] + d2Af[ 7]*tpx[3]);
654 d2a[2] = (d2Af[ 8]*tpx[0] + d2Af[ 9]*tpx[1] + d2Af[10]*tpx[2] + d2Af[11]*tpx[3]);
655 d2a[3] = (d2Af[12]*tpx[0] + d2Af[13]*tpx[1] + d2Af[14]*tpx[2] + d2Af[15]*tpx[3]);
656
657 b[0] = (Af[ 0]*tpy[0] + Af[ 1]*tpy[1] + Af[ 2]*tpy[2] + Af[ 3]*tpy[3]);
658 b[1] = (Af[ 4]*tpy[0] + Af[ 5]*tpy[1] + Af[ 6]*tpy[2] + Af[ 7]*tpy[3]);
659 b[2] = (Af[ 8]*tpy[0] + Af[ 9]*tpy[1] + Af[10]*tpy[2] + Af[11]*tpy[3]);
660 b[3] = (Af[12]*tpy[0] + Af[13]*tpy[1] + Af[14]*tpy[2] + Af[15]*tpy[3]);
661 db[0] = (dAf[ 0]*tpy[0] + dAf[ 1]*tpy[1] + dAf[ 2]*tpy[2] + dAf[ 3]*tpy[3]);
662 db[1] = (dAf[ 4]*tpy[0] + dAf[ 5]*tpy[1] + dAf[ 6]*tpy[2] + dAf[ 7]*tpy[3]);
663 db[2] = (dAf[ 8]*tpy[0] + dAf[ 9]*tpy[1] + dAf[10]*tpy[2] + dAf[11]*tpy[3]);
664 db[3] = (dAf[12]*tpy[0] + dAf[13]*tpy[1] + dAf[14]*tpy[2] + dAf[15]*tpy[3]);
665 d2b[0] = (d2Af[ 0]*tpy[0] + d2Af[ 1]*tpy[1] + d2Af[ 2]*tpy[2] + d2Af[ 3]*tpy[3]);
666 d2b[1] = (d2Af[ 4]*tpy[0] + d2Af[ 5]*tpy[1] + d2Af[ 6]*tpy[2] + d2Af[ 7]*tpy[3]);
667 d2b[2] = (d2Af[ 8]*tpy[0] + d2Af[ 9]*tpy[1] + d2Af[10]*tpy[2] + d2Af[11]*tpy[3]);
668 d2b[3] = (d2Af[12]*tpy[0] + d2Af[13]*tpy[1] + d2Af[14]*tpy[2] + d2Af[15]*tpy[3]);
669
670 c[0] = (Af[ 0]*tpz[0] + Af[ 1]*tpz[1] + Af[ 2]*tpz[2] + Af[ 3]*tpz[3]);
671 c[1] = (Af[ 4]*tpz[0] + Af[ 5]*tpz[1] + Af[ 6]*tpz[2] + Af[ 7]*tpz[3]);
672 c[2] = (Af[ 8]*tpz[0] + Af[ 9]*tpz[1] + Af[10]*tpz[2] + Af[11]*tpz[3]);
673 c[3] = (Af[12]*tpz[0] + Af[13]*tpz[1] + Af[14]*tpz[2] + Af[15]*tpz[3]);
674 dc[0] = (dAf[ 0]*tpz[0] + dAf[ 1]*tpz[1] + dAf[ 2]*tpz[2] + dAf[ 3]*tpz[3]);
675 dc[1] = (dAf[ 4]*tpz[0] + dAf[ 5]*tpz[1] + dAf[ 6]*tpz[2] + dAf[ 7]*tpz[3]);
676 dc[2] = (dAf[ 8]*tpz[0] + dAf[ 9]*tpz[1] + dAf[10]*tpz[2] + dAf[11]*tpz[3]);
677 dc[3] = (dAf[12]*tpz[0] + dAf[13]*tpz[1] + dAf[14]*tpz[2] + dAf[15]*tpz[3]);
678 d2c[0] = (d2Af[ 0]*tpz[0] + d2Af[ 1]*tpz[1] + d2Af[ 2]*tpz[2] + d2Af[ 3]*tpz[3]);
679 d2c[1] = (d2Af[ 4]*tpz[0] + d2Af[ 5]*tpz[1] + d2Af[ 6]*tpz[2] + d2Af[ 7]*tpz[3]);
680 d2c[2] = (d2Af[ 8]*tpz[0] + d2Af[ 9]*tpz[1] + d2Af[10]*tpz[2] + d2Af[11]*tpz[3]);
681 d2c[3] = (d2Af[12]*tpz[0] + d2Af[13]*tpz[1] + d2Af[14]*tpz[2] + d2Af[15]*tpz[3]);
682
683 int xs = spline->x_stride;
684 int ys = spline->y_stride;
685 int zs = spline->z_stride;
686
687 complex_float lapl3[3*spline->num_splines];
688 for (int n=0; n<spline->num_splines; n++) {
689 vals[n] = 0.0;
690 grads[3*n+0] = grads[3*n+1] = grads[3*n+2] = 0.0;
691 lapl3[3*n+0] = lapl3[3*n+1] = lapl3[3*n+2] = 0.0;
692 }
693
694 for (int i=0; i<4; i++)
695 for (int j=0; j<4; j++)
696 for (int k=0; k<4; k++) {
697 float abc = a[i]*b[j]*c[k];
698 float dabc[3], d2abc[3];
699 dabc[0] = da[i]* b[j]* c[k];
700 dabc[1] = a[i]*db[j]* c[k];
701 dabc[2] = a[i]* b[j]*dc[k];
702 d2abc[0] = d2a[i]* b[j]* c[k];
703 d2abc[1] = a[i]*d2b[j]* c[k];
704 d2abc[2] = a[i]* b[j]*d2c[k];
705
706 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
707 for (int n=0; n<spline->num_splines; n++) {
708 vals[n] += abc *coefs[n];
709 grads[3*n+0] += dabc[0]*coefs[n];
710 grads[3*n+1] += dabc[1]*coefs[n];
711 grads[3*n+2] += dabc[2]*coefs[n];
712 lapl3[3*n+0] += d2abc[0]*coefs[n];
713 lapl3[3*n+1] += d2abc[1]*coefs[n];
714 lapl3[3*n+2] += d2abc[2]*coefs[n];
715 }
716 }
717
718 float dxInv = spline->x_grid.delta_inv;
719 float dyInv = spline->y_grid.delta_inv;
720 float dzInv = spline->z_grid.delta_inv;
721 for (int n=0; n<spline->num_splines; n++) {
722 grads[3*n+0] *= dxInv;
723 grads[3*n+1] *= dyInv;
724 grads[3*n+2] *= dzInv;
725 lapl3[3*n+0] *= dxInv*dxInv;
726 lapl3[3*n+1] *= dyInv*dyInv;
727 lapl3[3*n+2] *= dzInv*dzInv;
728 lapl[n] = lapl3[3*n+0] + lapl3[3*n+1] + lapl3[3*n+2];
729 }
730}
731
732
733
734inline void
735eval_multi_UBspline_3d_c_vgh (multi_UBspline_3d_c *spline,
736 double x, double y, double z,
737 complex_float* restrict vals,
738 complex_float* restrict grads,
739 complex_float* restrict hess)
740{
741 x -= spline->x_grid.start;
742 y -= spline->y_grid.start;
743 z -= spline->z_grid.start;
744 float ux = x*spline->x_grid.delta_inv;
745 float uy = y*spline->y_grid.delta_inv;
746 float uz = z*spline->z_grid.delta_inv;
747 float ipartx, iparty, ipartz, tx, ty, tz;
748 tx = modff (ux, &ipartx); int ix = (int) ipartx;
749 ty = modff (uy, &iparty); int iy = (int) iparty;
750 tz = modff (uz, &ipartz); int iz = (int) ipartz;
751
752 float tpx[4], tpy[4], tpz[4], a[4], b[4], c[4],
753 da[4], db[4], dc[4], d2a[4], d2b[4], d2c[4];
754 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
755 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
756 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
757 complex_float* restrict coefs = spline->coefs;
758
759 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
760 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
761 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
762 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
763 da[0] = (dAf[ 1]*tpx[1] + dAf[ 2]*tpx[2] + dAf[ 3]*tpx[3]);
764 da[1] = (dAf[ 5]*tpx[1] + dAf[ 6]*tpx[2] + dAf[ 7]*tpx[3]);
765 da[2] = (dAf[ 9]*tpx[1] + dAf[10]*tpx[2] + dAf[11]*tpx[3]);
766 da[3] = (dAf[13]*tpx[1] + dAf[14]*tpx[2] + dAf[15]*tpx[3]);
767 d2a[0] = (d2Af[ 2]*tpx[2] + d2Af[ 3]*tpx[3]);
768 d2a[1] = (d2Af[ 6]*tpx[2] + d2Af[ 7]*tpx[3]);
769 d2a[2] = (d2Af[10]*tpx[2] + d2Af[11]*tpx[3]);
770 d2a[3] = (d2Af[14]*tpx[2] + d2Af[15]*tpx[3]);
771
772 b[0] = (Af[ 0]*tpy[0] + Af[ 1]*tpy[1] + Af[ 2]*tpy[2] + Af[ 3]*tpy[3]);
773 b[1] = (Af[ 4]*tpy[0] + Af[ 5]*tpy[1] + Af[ 6]*tpy[2] + Af[ 7]*tpy[3]);
774 b[2] = (Af[ 8]*tpy[0] + Af[ 9]*tpy[1] + Af[10]*tpy[2] + Af[11]*tpy[3]);
775 b[3] = (Af[12]*tpy[0] + Af[13]*tpy[1] + Af[14]*tpy[2] + Af[15]*tpy[3]);
776 db[0] = (dAf[ 1]*tpy[1] + dAf[ 2]*tpy[2] + dAf[ 3]*tpy[3]);
777 db[1] = (dAf[ 5]*tpy[1] + dAf[ 6]*tpy[2] + dAf[ 7]*tpy[3]);
778 db[2] = (dAf[ 9]*tpy[1] + dAf[10]*tpy[2] + dAf[11]*tpy[3]);
779 db[3] = (dAf[13]*tpy[1] + dAf[14]*tpy[2] + dAf[15]*tpy[3]);
780 d2b[0] = (d2Af[ 2]*tpy[2] + d2Af[ 3]*tpy[3]);
781 d2b[1] = (d2Af[ 6]*tpy[2] + d2Af[ 7]*tpy[3]);
782 d2b[2] = (d2Af[10]*tpy[2] + d2Af[11]*tpy[3]);
783 d2b[3] = (d2Af[14]*tpy[2] + d2Af[15]*tpy[3]);
784
785 c[0] = (Af[ 0]*tpz[0] + Af[ 1]*tpz[1] + Af[ 2]*tpz[2] + Af[ 3]*tpz[3]);
786 c[1] = (Af[ 4]*tpz[0] + Af[ 5]*tpz[1] + Af[ 6]*tpz[2] + Af[ 7]*tpz[3]);
787 c[2] = (Af[ 8]*tpz[0] + Af[ 9]*tpz[1] + Af[10]*tpz[2] + Af[11]*tpz[3]);
788 c[3] = (Af[12]*tpz[0] + Af[13]*tpz[1] + Af[14]*tpz[2] + Af[15]*tpz[3]);
789 dc[0] = (dAf[ 1]*tpz[1] + dAf[ 2]*tpz[2] + dAf[ 3]*tpz[3]);
790 dc[1] = (dAf[ 5]*tpz[1] + dAf[ 6]*tpz[2] + dAf[ 7]*tpz[3]);
791 dc[2] = (dAf[ 9]*tpz[1] + dAf[10]*tpz[2] + dAf[11]*tpz[3]);
792 dc[3] = (dAf[13]*tpz[1] + dAf[14]*tpz[2] + dAf[15]*tpz[3]);
793 d2c[0] = (d2Af[ 2]*tpz[2] + d2Af[ 3]*tpz[3]);
794 d2c[1] = (d2Af[ 6]*tpz[2] + d2Af[ 7]*tpz[3]);
795 d2c[2] = (d2Af[10]*tpz[2] + d2Af[11]*tpz[3]);
796 d2c[3] = (d2Af[14]*tpz[2] + d2Af[15]*tpz[3]);
797
798 int xs = spline->x_stride;
799 int ys = spline->y_stride;
800 int zs = spline->z_stride;
801
802 for (int n=0; n<spline->num_splines; n++) {
803 vals[n] = 0.0;
804 grads[3*n+0] = grads[3*n+1] = grads[3*n+2] = 0.0;
805 for (int i=0; i<9; i++)
806 hess[9*n+i] = 0.0;
807 }
808
809 for (int i=0; i<4; i++)
810 for (int j=0; j<4; j++)
811 for (int k=0; k<4; k++) {
812 float abc = a[i]*b[j]*c[k];
813 float dabc[3], d2abc[6];
814 dabc[0] = da[i]* b[j]* c[k];
815 dabc[1] = a[i]*db[j]* c[k];
816 dabc[2] = a[i]* b[j]*dc[k];
817 d2abc[0] = d2a[i]* b[j]* c[k];
818 d2abc[1] = da[i]* db[j]* c[k];
819 d2abc[2] = da[i]* b[j]* dc[k];
820 d2abc[3] = a[i]*d2b[j]* c[k];
821 d2abc[4] = a[i]* db[j]* dc[k];
822 d2abc[5] = a[i]* b[j]*d2c[k];
823
824 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
825 for (int n=0; n<spline->num_splines; n++) {
826 vals[n] += abc *coefs[n];
827 grads[3*n+0] += dabc[0]*coefs[n];
828 grads[3*n+1] += dabc[1]*coefs[n];
829 grads[3*n+2] += dabc[2]*coefs[n];
830 hess [9*n+0] += d2abc[0]*coefs[n];
831 hess [9*n+1] += d2abc[1]*coefs[n];
832 hess [9*n+2] += d2abc[2]*coefs[n];
833 hess [9*n+4] += d2abc[3]*coefs[n];
834 hess [9*n+5] += d2abc[4]*coefs[n];
835 hess [9*n+8] += d2abc[5]*coefs[n];
836 }
837 }
838
839 float dxInv = spline->x_grid.delta_inv;
840 float dyInv = spline->y_grid.delta_inv;
841 float dzInv = spline->z_grid.delta_inv;
842 for (int n=0; n<spline->num_splines; n++) {
843 grads[3*n+0] *= dxInv;
844 grads[3*n+1] *= dyInv;
845 grads[3*n+2] *= dzInv;
846 hess[9*n+0] *= dxInv*dxInv;
847 hess[9*n+4] *= dyInv*dyInv;
848 hess[9*n+8] *= dzInv*dzInv;
849 hess[9*n+1] *= dxInv*dyInv;
850 hess[9*n+2] *= dxInv*dzInv;
851 hess[9*n+5] *= dyInv*dzInv;
852 // Copy hessian elements into lower half of 3x3 matrix
853 hess[9*n+3] = hess[9*n+1];
854 hess[9*n+6] = hess[9*n+2];
855 hess[9*n+7] = hess[9*n+5];
856 }
857}
858
859
860inline void
861eval_multi_UBspline_3d_s_vghgh (multi_UBspline_3d_c *spline,
862 double x, double y, double z,
863 complex_float* restrict vals,
864 complex_float* restrict grads,
865 complex_float* restrict hess,
866 complex_float* restrict gradhess)
867{
868 x -= spline->x_grid.start;
869 y -= spline->y_grid.start;
870 z -= spline->z_grid.start;
871 float ux = x*spline->x_grid.delta_inv;
872 float uy = y*spline->y_grid.delta_inv;
873 float uz = z*spline->z_grid.delta_inv;
874 float ipartx, iparty, ipartz, tx, ty, tz;
875 tx = modff (ux, &ipartx); int ix = (int) ipartx;
876 ty = modff (uy, &iparty); int iy = (int) iparty;
877 tz = modff (uz, &ipartz); int iz = (int) ipartz;
878
879 float tpx[4], tpy[4], tpz[4], a[4], b[4], c[4],
880 da[4], db[4], dc[4], d2a[4], d2b[4], d2c[4],
881 d3a[4], d3b[4], d3c[4];
882 tpx[0] = tx*tx*tx; tpx[1] = tx*tx; tpx[2] = tx; tpx[3] = 1.0;
883 tpy[0] = ty*ty*ty; tpy[1] = ty*ty; tpy[2] = ty; tpy[3] = 1.0;
884 tpz[0] = tz*tz*tz; tpz[1] = tz*tz; tpz[2] = tz; tpz[3] = 1.0;
885 complex_float* restrict coefs = spline->coefs;
886
887 a[0] = (Af[ 0]*tpx[0] + Af[ 1]*tpx[1] + Af[ 2]*tpx[2] + Af[ 3]*tpx[3]);
888 a[1] = (Af[ 4]*tpx[0] + Af[ 5]*tpx[1] + Af[ 6]*tpx[2] + Af[ 7]*tpx[3]);
889 a[2] = (Af[ 8]*tpx[0] + Af[ 9]*tpx[1] + Af[10]*tpx[2] + Af[11]*tpx[3]);
890 a[3] = (Af[12]*tpx[0] + Af[13]*tpx[1] + Af[14]*tpx[2] + Af[15]*tpx[3]);
891 da[0] = (dAf[ 0]*tpx[0] + dAf[ 1]*tpx[1] + dAf[ 2]*tpx[2] + dAf[ 3]*tpx[3]);
892 da[1] = (dAf[ 4]*tpx[0] + dAf[ 5]*tpx[1] + dAf[ 6]*tpx[2] + dAf[ 7]*tpx[3]);
893 da[2] = (dAf[ 8]*tpx[0] + dAf[ 9]*tpx[1] + dAf[10]*tpx[2] + dAf[11]*tpx[3]);
894 da[3] = (dAf[12]*tpx[0] + dAf[13]*tpx[1] + dAf[14]*tpx[2] + dAf[15]*tpx[3]);
895 d2a[0] = (d2Af[ 0]*tpx[0] + d2Af[ 1]*tpx[1] + d2Af[ 2]*tpx[2] + d2Af[ 3]*tpx[3]);
896 d2a[1] = (d2Af[ 4]*tpx[0] + d2Af[ 5]*tpx[1] + d2Af[ 6]*tpx[2] + d2Af[ 7]*tpx[3]);
897 d2a[2] = (d2Af[ 8]*tpx[0] + d2Af[ 9]*tpx[1] + d2Af[10]*tpx[2] + d2Af[11]*tpx[3]);
898 d2a[3] = (d2Af[12]*tpx[0] + d2Af[13]*tpx[1] + d2Af[14]*tpx[2] + d2Af[15]*tpx[3]);
899 d3a[0] = (/*d2Af[ 0]*tpx[0] + d2Af[ 1]*tpx[1] + d2Af[ 2]*tpx[2] +*/ d3Af[ 3]*tpx[3]);
900 d3a[1] = (/*d2Af[ 4]*tpx[0] + d2Af[ 5]*tpx[1] + d2Af[ 6]*tpx[2] +*/ d3Af[ 7]*tpx[3]);
901 d3a[2] = (/*d2Af[ 8]*tpx[0] + d2Af[ 9]*tpx[1] + d2Af[10]*tpx[2] +*/ d3Af[11]*tpx[3]);
902 d3a[3] = (/*d2Af[12]*tpx[0] + d2Af[13]*tpx[1] + d2Af[14]*tpx[2] +*/ d3Af[15]*tpx[3]);
903
904 b[0] = (Af[ 0]*tpy[0] + Af[ 1]*tpy[1] + Af[ 2]*tpy[2] + Af[ 3]*tpy[3]);
905 b[1] = (Af[ 4]*tpy[0] + Af[ 5]*tpy[1] + Af[ 6]*tpy[2] + Af[ 7]*tpy[3]);
906 b[2] = (Af[ 8]*tpy[0] + Af[ 9]*tpy[1] + Af[10]*tpy[2] + Af[11]*tpy[3]);
907 b[3] = (Af[12]*tpy[0] + Af[13]*tpy[1] + Af[14]*tpy[2] + Af[15]*tpy[3]);
908 db[0] = (dAf[ 0]*tpy[0] + dAf[ 1]*tpy[1] + dAf[ 2]*tpy[2] + dAf[ 3]*tpy[3]);
909 db[1] = (dAf[ 4]*tpy[0] + dAf[ 5]*tpy[1] + dAf[ 6]*tpy[2] + dAf[ 7]*tpy[3]);
910 db[2] = (dAf[ 8]*tpy[0] + dAf[ 9]*tpy[1] + dAf[10]*tpy[2] + dAf[11]*tpy[3]);
911 db[3] = (dAf[12]*tpy[0] + dAf[13]*tpy[1] + dAf[14]*tpy[2] + dAf[15]*tpy[3]);
912 d2b[0] = (d2Af[ 0]*tpy[0] + d2Af[ 1]*tpy[1] + d2Af[ 2]*tpy[2] + d2Af[ 3]*tpy[3]);
913 d2b[1] = (d2Af[ 4]*tpy[0] + d2Af[ 5]*tpy[1] + d2Af[ 6]*tpy[2] + d2Af[ 7]*tpy[3]);
914 d2b[2] = (d2Af[ 8]*tpy[0] + d2Af[ 9]*tpy[1] + d2Af[10]*tpy[2] + d2Af[11]*tpy[3]);
915 d2b[3] = (d2Af[12]*tpy[0] + d2Af[13]*tpy[1] + d2Af[14]*tpy[2] + d2Af[15]*tpy[3]);
916 d3b[0] = (/*d2Af[ 0]*tpx[0] + d2Af[ 1]*tpx[1] + d2Af[ 2]*tpx[2] +*/ d3Af[ 3]*tpy[3]);
917 d3b[1] = (/*d2Af[ 4]*tpx[0] + d2Af[ 5]*tpx[1] + d2Af[ 6]*tpx[2] +*/ d3Af[ 7]*tpy[3]);
918 d3b[2] = (/*d2Af[ 8]*tpx[0] + d2Af[ 9]*tpx[1] + d2Af[10]*tpx[2] +*/ d3Af[11]*tpy[3]);
919 d3b[3] = (/*d2Af[12]*tpx[0] + d2Af[13]*tpx[1] + d2Af[14]*tpx[2] +*/ d3Af[15]*tpy[3]);
920
921 c[0] = (Af[ 0]*tpz[0] + Af[ 1]*tpz[1] + Af[ 2]*tpz[2] + Af[ 3]*tpz[3]);
922 c[1] = (Af[ 4]*tpz[0] + Af[ 5]*tpz[1] + Af[ 6]*tpz[2] + Af[ 7]*tpz[3]);
923 c[2] = (Af[ 8]*tpz[0] + Af[ 9]*tpz[1] + Af[10]*tpz[2] + Af[11]*tpz[3]);
924 c[3] = (Af[12]*tpz[0] + Af[13]*tpz[1] + Af[14]*tpz[2] + Af[15]*tpz[3]);
925 dc[0] = (dAf[ 0]*tpz[0] + dAf[ 1]*tpz[1] + dAf[ 2]*tpz[2] + dAf[ 3]*tpz[3]);
926 dc[1] = (dAf[ 4]*tpz[0] + dAf[ 5]*tpz[1] + dAf[ 6]*tpz[2] + dAf[ 7]*tpz[3]);
927 dc[2] = (dAf[ 8]*tpz[0] + dAf[ 9]*tpz[1] + dAf[10]*tpz[2] + dAf[11]*tpz[3]);
928 dc[3] = (dAf[12]*tpz[0] + dAf[13]*tpz[1] + dAf[14]*tpz[2] + dAf[15]*tpz[3]);
929 d2c[0] = (d2Af[ 0]*tpz[0] + d2Af[ 1]*tpz[1] + d2Af[ 2]*tpz[2] + d2Af[ 3]*tpz[3]);
930 d2c[1] = (d2Af[ 4]*tpz[0] + d2Af[ 5]*tpz[1] + d2Af[ 6]*tpz[2] + d2Af[ 7]*tpz[3]);
931 d2c[2] = (d2Af[ 8]*tpz[0] + d2Af[ 9]*tpz[1] + d2Af[10]*tpz[2] + d2Af[11]*tpz[3]);
932 d2c[3] = (d2Af[12]*tpz[0] + d2Af[13]*tpz[1] + d2Af[14]*tpz[2] + d2Af[15]*tpz[3]);
933 d3c[0] = (/*d2Af[ 0]*tpx[0] + d2Af[ 1]*tpx[1] + d2Af[ 2]*tpx[2] +*/ d3Af[ 3]*tpz[3]);
934 d3c[1] = (/*d2Af[ 4]*tpx[0] + d2Af[ 5]*tpx[1] + d2Af[ 6]*tpx[2] +*/ d3Af[ 7]*tpz[3]);
935 d3c[2] = (/*d2Af[ 8]*tpx[0] + d2Af[ 9]*tpx[1] + d2Af[10]*tpx[2] +*/ d3Af[11]*tpz[3]);
936 d3c[3] = (/*d2Af[12]*tpx[0] + d2Af[13]*tpx[1] + d2Af[14]*tpx[2] +*/ d3Af[15]*tpz[3]);
937
938 int xs = spline->x_stride;
939 int ys = spline->y_stride;
940 int zs = spline->z_stride;
941
942 for (int n=0; n<spline->num_splines; n++) {
943 vals[n] = 0.0;
944 grads[3*n+0] = grads[3*n+1] = grads[3*n+2] = 0.0;
945 for (int i=0; i<9; i++)
946 hess[9*n+i] = 0.0;
947 for (int i=0; i<27; i++)
948 gradhess[27*n+i] = 0.0;
949 }
950
951 for (int i=0; i<4; i++)
952 for (int j=0; j<4; j++)
953 for (int k=0; k<4; k++) {
954 float abc = a[i]*b[j]*c[k];
955 float dabc[3], d2abc[6], d3abc[10];
956 dabc[0] = da[i]* b[j]* c[k];
957 dabc[1] = a[i]*db[j]* c[k];
958 dabc[2] = a[i]* b[j]*dc[k];
959 d2abc[0] = d2a[i]* b[j]* c[k];
960 d2abc[1] = da[i]* db[j]* c[k];
961 d2abc[2] = da[i]* b[j]* dc[k];
962 d2abc[3] = a[i]*d2b[j]* c[k];
963 d2abc[4] = a[i]* db[j]* dc[k];
964 d2abc[5] = a[i]* b[j]*d2c[k];
965
966 d3abc[0] = d3a[i]* b[j]* c[k];
967 d3abc[1] = d2a[i]* db[j]* c[k];
968 d3abc[2] = d2a[i]* b[j]* dc[k];
969 d3abc[3] = da[i]*d2b[j]* c[k];
970 d3abc[4] = da[i]* db[j]* dc[k];
971 d3abc[5] = da[i]* b[j]*d2c[k];
972 d3abc[6] = a[i]*d3b[j]* c[k];
973 d3abc[7] = a[i]*d2b[j]* dc[k];
974 d3abc[8] = a[i]* db[j]*d2c[k];
975 d3abc[9] = a[i]* b[j]*d3c[k];
976
977 complex_float* restrict coefs = spline->coefs + ((ix+i)*xs + (iy+j)*ys + (iz+k)*zs);
978 for (int n=0; n<spline->num_splines; n++) {
979 vals[n] += abc *coefs[n];
980 grads[3*n+0] += dabc[0]*coefs[n];
981 grads[3*n+1] += dabc[1]*coefs[n];
982 grads[3*n+2] += dabc[2]*coefs[n];
983 hess [9*n+0] += d2abc[0]*coefs[n];
984 hess [9*n+1] += d2abc[1]*coefs[n];
985 hess [9*n+2] += d2abc[2]*coefs[n];
986 hess [9*n+4] += d2abc[3]*coefs[n];
987 hess [9*n+5] += d2abc[4]*coefs[n];
988 hess [9*n+8] += d2abc[5]*coefs[n];
989
990 gradhess [27*n+0 ] += d3abc[0]*coefs[n];
991 gradhess [27*n+1 ] += d3abc[1]*coefs[n];
992 gradhess [27*n+2 ] += d3abc[2]*coefs[n];
993 gradhess [27*n+4 ] += d3abc[3]*coefs[n];
994 gradhess [27*n+5 ] += d3abc[4]*coefs[n];
995 gradhess [27*n+8 ] += d3abc[5]*coefs[n];
996 gradhess [27*n+13] += d3abc[6]*coefs[n];
997 gradhess [27*n+14] += d3abc[7]*coefs[n];
998 gradhess [27*n+17] += d3abc[8]*coefs[n];
999 gradhess [27*n+26] += d3abc[9]*coefs[n];
1000 }
1001 }
1002
1003 float dxInv = spline->x_grid.delta_inv;
1004 float dyInv = spline->y_grid.delta_inv;
1005 float dzInv = spline->z_grid.delta_inv;
1006 for (int n=0; n<spline->num_splines; n++) {
1007 grads[3*n+0] *= dxInv;
1008 grads[3*n+1] *= dyInv;
1009 grads[3*n+2] *= dzInv;
1010 hess [9*n+0] *= dxInv*dxInv;
1011 hess [9*n+4] *= dyInv*dyInv;
1012 hess [9*n+8] *= dzInv*dzInv;
1013 hess [9*n+1] *= dxInv*dyInv;
1014 hess [9*n+2] *= dxInv*dzInv;
1015 hess [9*n+5] *= dyInv*dzInv;
1016 // Copy hessian elements into lower half of 3x3 matrix
1017 hess [9*n+3] = hess[9*n+1];
1018 hess [9*n+6] = hess[9*n+2];
1019 hess [9*n+7] = hess[9*n+5];
1020
1021 gradhess [27*n+0 ] *= dxInv*dxInv*dxInv;
1022 gradhess [27*n+1 ] *= dxInv*dxInv*dyInv;
1023 gradhess [27*n+2 ] *= dxInv*dxInv*dzInv;
1024 gradhess [27*n+4 ] *= dxInv*dyInv*dyInv;
1025 gradhess [27*n+5 ] *= dxInv*dyInv*dzInv;
1026 gradhess [27*n+8 ] *= dxInv*dzInv*dzInv;
1027 gradhess [27*n+13] *= dyInv*dyInv*dyInv;
1028 gradhess [27*n+14] *= dyInv*dyInv*dzInv;
1029 gradhess [27*n+17] *= dyInv*dzInv*dzInv;
1030 gradhess [27*n+26] *= dzInv*dzInv*dzInv;
1031
1032 // Copy gradhess elements into rest of tensor
1033 gradhess [27*n+9 ] = gradhess [27*n+3 ] = gradhess [27*n+1 ];
1034 gradhess [27*n+18 ] = gradhess [27*n+6 ] = gradhess [27*n+2 ];
1035 gradhess [27*n+22 ] = gradhess [27*n+16 ] = gradhess [27*n+14];
1036 gradhess [27*n+12 ] = gradhess [27*n+10 ] = gradhess [27*n+4 ];
1037 gradhess [27*n+24 ] = gradhess [27*n+20 ] = gradhess [27*n+8 ];
1038 gradhess [27*n+25 ] = gradhess [27*n+23 ] = gradhess [27*n+17];
1039 gradhess [27*n+21 ] = gradhess [27*n+19 ] = gradhess [27*n+15] = gradhess [27*n+11 ] = gradhess [27*n+7 ] = gradhess [27*n+5];
1040
1041 }
1042}
1043
1044#endif
complex_float
complex float complex_float
Definition bspline_base.h:33
restrict
#define restrict
Definition local_definitions.h:4
eval_multi_UBspline_2d_c_vg
void eval_multi_UBspline_2d_c_vg(multi_UBspline_2d_c *spline, double x, double y, complex_float *restrict vals, complex_float *restrict grads)
Definition multi_bspline_eval_std_c.h:227
eval_multi_UBspline_1d_c_vg
void eval_multi_UBspline_1d_c_vg(multi_UBspline_1d_c *spline, double x, complex_float *restrict vals, complex_float *restrict grads)
Definition multi_bspline_eval_std_c.h:70
eval_multi_UBspline_2d_c_vgl
void eval_multi_UBspline_2d_c_vgl(multi_UBspline_2d_c *spline, double x, double y, complex_float *restrict vals, complex_float *restrict grads, complex_float *restrict lapl)
Definition multi_bspline_eval_std_c.h:295
eval_multi_UBspline_3d_c_vgl
void eval_multi_UBspline_3d_c_vgl(multi_UBspline_3d_c *spline, double x, double y, double z, complex_float *restrict vals, complex_float *restrict grads, complex_float *restrict lapl)
Definition multi_bspline_eval_std_c.h:620
d3Af
const float *restrict d3Af
Definition bspline_data.cpp:171
eval_multi_UBspline_2d_c_vgh
void eval_multi_UBspline_2d_c_vgh(multi_UBspline_2d_c *spline, double x, double y, complex_float *restrict vals, complex_float *restrict grads, complex_float *restrict hess)
Definition multi_bspline_eval_std_c.h:381
eval_multi_UBspline_3d_c
void eval_multi_UBspline_3d_c(multi_UBspline_3d_c *spline, double x, double y, double z, complex_float *restrict vals)
Definition multi_bspline_eval_std_c.h:475
eval_multi_UBspline_1d_c_vgl
void eval_multi_UBspline_1d_c_vgl(multi_UBspline_1d_c *spline, double x, complex_float *restrict vals, complex_float *restrict grads, complex_float *restrict lapl)
Definition multi_bspline_eval_std_c.h:115
eval_multi_UBspline_3d_c_vg
void eval_multi_UBspline_3d_c_vg(multi_UBspline_3d_c *spline, double x, double y, double z, complex_float *restrict vals, complex_float *restrict grads)
Definition multi_bspline_eval_std_c.h:530
eval_multi_UBspline_3d_s_vghgh
void eval_multi_UBspline_3d_s_vghgh(multi_UBspline_3d_c *spline, double x, double y, double z, complex_float *restrict vals, complex_float *restrict grads, complex_float *restrict hess, complex_float *restrict gradhess)
Definition multi_bspline_eval_std_c.h:861
eval_multi_UBspline_3d_c_vgh
void eval_multi_UBspline_3d_c_vgh(multi_UBspline_3d_c *spline, double x, double y, double z, complex_float *restrict vals, complex_float *restrict grads, complex_float *restrict hess)
Definition multi_bspline_eval_std_c.h:735
eval_multi_UBspline_1d_c
void eval_multi_UBspline_1d_c(multi_UBspline_1d_c *spline, double x, complex_float *restrict vals)
Definition multi_bspline_eval_std_c.h:37
dAf
const float *restrict dAf
Definition bspline_data.cpp:157
eval_multi_UBspline_1d_c_vgh
void eval_multi_UBspline_1d_c_vgh(multi_UBspline_1d_c *spline, double x, complex_float *restrict vals, complex_float *restrict grads, complex_float *restrict hess)
Definition multi_bspline_eval_std_c.h:169
d2Af
const float *restrict d2Af
Definition bspline_data.cpp:164
eval_multi_UBspline_2d_c
void eval_multi_UBspline_2d_c(multi_UBspline_2d_c *spline, double x, double y, complex_float *restrict vals)
Definition multi_bspline_eval_std_c.h:183
Af
const float *restrict Af
Definition bspline_data.cpp:150
multi_bspline_structs.h
Ugrid::start
double start
Definition bspline_base.h:85
Ugrid::delta_inv
double delta_inv
Definition bspline_base.h:89
multi_UBspline_1d_c
Definition multi_bspline_structs.h:104
multi_UBspline_1d_c::coefs
complex_float *restrict coefs
Definition multi_bspline_structs.h:107
multi_UBspline_1d_c::num_splines
int num_splines
Definition multi_bspline_structs.h:111
multi_UBspline_1d_c::x_stride
intptr_t x_stride
Definition multi_bspline_structs.h:108
multi_UBspline_1d_c::x_grid
Ugrid x_grid
Definition multi_bspline_structs.h:109
multi_UBspline_2d_c
Definition multi_bspline_structs.h:115
multi_UBspline_2d_c::coefs
complex_float *restrict coefs
Definition multi_bspline_structs.h:118
multi_UBspline_2d_c::y_stride
intptr_t y_stride
Definition multi_bspline_structs.h:119
multi_UBspline_2d_c::num_splines
int num_splines
Definition multi_bspline_structs.h:122
multi_UBspline_2d_c::y_grid
Ugrid y_grid
Definition multi_bspline_structs.h:120
multi_UBspline_2d_c::x_grid
Ugrid x_grid
Definition multi_bspline_structs.h:120
multi_UBspline_2d_c::x_stride
intptr_t x_stride
Definition multi_bspline_structs.h:119
multi_UBspline_3d_c
Definition multi_bspline_structs.h:128
multi_UBspline_3d_c::z_grid
Ugrid z_grid
Definition multi_bspline_structs.h:133
multi_UBspline_3d_c::coefs
complex_float *restrict coefs
Definition multi_bspline_structs.h:131
multi_UBspline_3d_c::x_grid
Ugrid x_grid
Definition multi_bspline_structs.h:133
multi_UBspline_3d_c::x_stride
intptr_t x_stride
Definition multi_bspline_structs.h:132
multi_UBspline_3d_c::y_stride
intptr_t y_stride
Definition multi_bspline_structs.h:132
multi_UBspline_3d_c::z_stride
intptr_t z_stride
Definition multi_bspline_structs.h:132
multi_UBspline_3d_c::y_grid
Ugrid y_grid
Definition multi_bspline_structs.h:133
multi_UBspline_3d_c::num_splines
int num_splines
Definition multi_bspline_structs.h:135