GCC Code Coverage Report
Directory: ./ Exec Total Coverage
File: src/limn/qn.c Lines: 0 168 0.0 %
Date: 2017-05-26 Branches: 0 257 0.0 %

Line Branch Exec Source
1
/*
2
  Teem: Tools to process and visualize scientific data and images             .
3
  Copyright (C) 2013, 2012, 2011, 2010, 2009  University of Chicago
4
  Copyright (C) 2008, 2007, 2006, 2005  Gordon Kindlmann
5
  Copyright (C) 2004, 2003, 2002, 2001, 2000, 1999, 1998  University of Utah
6
7
  This library is free software; you can redistribute it and/or
8
  modify it under the terms of the GNU Lesser General Public License
9
  (LGPL) as published by the Free Software Foundation; either
10
  version 2.1 of the License, or (at your option) any later version.
11
  The terms of redistributing and/or modifying this software also
12
  include exceptions to the LGPL that facilitate static linking.
13
14
  This library is distributed in the hope that it will be useful,
15
  but WITHOUT ANY WARRANTY; without even the implied warranty of
16
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17
  Lesser General Public License for more details.
18
19
  You should have received a copy of the GNU Lesser General Public License
20
  along with this library; if not, write to Free Software Foundation, Inc.,
21
  51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22
*/
23
24
25
#include "limn.h"
26
27
/*
28
** "16simple": 16 bit representation based on
29
** <http://www.gamedev.net/reference/articles/article1252.asp>
30
**
31
** info: [z sign] [y sign] [x sign] [y component] [x component]
32
** bits:    1        1        1          7             6
33
**
34
** despite appearances, this is isotropic in X and Y
35
*/
36
37
void
38
_limnQN16simple_QNtoV_f(float *vec, unsigned int qn) {
39
  int xi, yi;
40
  float x, y, z, n;
41
42
  xi = qn & 0x3F;
43
  yi = (qn >> 6) & 0x7F;
44
  if (xi + yi >= 127) {
45
    xi = 127 - xi;
46
    yi = 127 - yi;
47
  }
48
  x = xi/126.0f;
49
  y = yi/126.0f;
50
  z = 1.0f - x - y;
51
  x = (qn & 0x2000) ? -x : x;
52
  y = (qn & 0x4000) ? -y : y;
53
  z = (qn & 0x8000) ? -z : z;
54
  n = AIR_CAST(float, 1.0/sqrt(x*x + y*y + z*z));
55
  vec[0] = x*n;
56
  vec[1] = y*n;
57
  vec[2] = z*n;
58
}
59
60
unsigned int
61
_limnQN16simple_VtoQN_f(const float *vec) {
62
  float x, y, z, L, w;
63
  unsigned int sgn = 0, xi, yi;
64
65
  x = vec[0];
66
  y = vec[1];
67
  z = vec[2];
68
  if (x < 0) {
69
    sgn |= 0x2000;
70
    x = -x;
71
  }
72
  if (y < 0) {
73
    sgn |= 0x4000;
74
    y = -y;
75
  }
76
  if (z < 0) {
77
    sgn |= 0x8000;
78
    z = -z;
79
  }
80
  L = x + y + z;
81
  if (!L) {
82
    return 0;
83
  }
84
  w = 126.0f/L;
85
  xi = AIR_CAST(unsigned int, x*w);
86
  yi = AIR_CAST(unsigned int, y*w);
87
  if (xi >= 64) {
88
    xi = 127 - xi;
89
    yi = 127 - yi;
90
  }
91
  return sgn | (yi << 6) | xi;
92
}
93
94
/* ----------------------------------------------------------------  */
95
96
void
97
_limnQN16border1_QNtoV_f(float *vec, unsigned int qn) {
98
  unsigned int ui, vi;
99
  float u, v, x, y, z, n;
100
101
  ui = qn & 0xFF;
102
  vi = qn >> 8;
103
  u = AIR_CAST(float, AIR_AFFINE(0.5, ui, 254.5, -0.5, 0.5));
104
  v = AIR_CAST(float, AIR_AFFINE(0.5, vi, 254.5, -0.5, 0.5));
105
  x =  u + v;
106
  y =  u - v;
107
  z = 1 - AIR_ABS(x) - AIR_ABS(y);
108
  z *= AIR_CAST(int, ((ui&1) ^ (vi&1)) << 1) - 1;
109
  n = AIR_CAST(float, 1.0/sqrt(x*x + y*y + z*z));
110
  vec[0] = x*n;
111
  vec[1] = y*n;
112
  vec[2] = z*n;
113
}
114
115
unsigned int
116
_limnQN16border1_VtoQN_f(const float *vec) {
117
  float L, u, v, x, y, z;
118
  unsigned int ui, vi, zi;
119
  char me[]="limnQNVto16PB1";
120
121
  x = vec[0];
122
  y = vec[1];
123
  z = vec[2];
124
  L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z);
125
  if (!L) {
126
    return 0;
127
  }
128
  x /= L;
129
  y /= L;
130
  u = x + y;
131
  v = x - y;
132
  ui = airIndex(-1, u, 1, 254); ui++;
133
  vi = airIndex(-1, v, 1, 254); vi++;
134
  zi = (ui ^ vi) & 0x01;
135
  if (!zi && z > 1.0/128.0) {
136
    ui -= (((ui >> 7) & 0x01) << 1) - 1;
137
  }
138
  else if (zi && z < -1.0/128.0) {
139
    vi -= (((vi >> 7) & 0x01) << 1) - 1;
140
  }
141
  zi = (ui ^ vi) & 0x01;
142
  if (!zi && z > 1.0/127.0) {
143
    printf("%s: panic01\n", me);
144
  }
145
  else if (zi && z < -1.0/127.0) {
146
    printf("%s: panic02\n", me);
147
  }
148
  return (vi << 8) | ui;
149
}
150
151
/* ----------------------------------------------------------------  */
152
153
  /*
154
    x =  [  1.0   1.0 ] u
155
    y    [  1.0  -1.0 ] v
156
  */
157
158
  /*
159
    u =  [  0.5   0.5 ] x
160
    v    [  0.5  -0.5 ] y
161
  */
162
163
  /* xor of low bits == 0 --> z<0 ; xor == 1 --> z>0 */
164
165
/* May 11 2003 GK (visually) verified that except at equator seam,
166
   the error due to quantization is unbiased */
167
168
#define _EVEN_CHECK_QtoV(HNB, TT, vec, qn) \
169
  ui = (qn) & ((1<<HNB)-1);                \
170
  vi = ((qn) >> HNB) & ((1<<HNB)-1);              \
171
  u = AIR_AFFINE(0, ui, ((1<<HNB)-1), -0.5, 0.5); \
172
  v = AIR_AFFINE(0, vi, ((1<<HNB)-1), -0.5, 0.5); \
173
  x =  u + v;                                     \
174
  y =  u - v;                                     \
175
  z = 1 - AIR_ABS(x) - AIR_ABS(y);                \
176
  z *= AIR_CAST(int, ((ui & 1) ^ (vi & 1)) << 1) - 1; \
177
  n = 1.0/sqrt(x*x + y*y + z*z); \
178
  (vec)[0] = AIR_CAST(TT, x*n); \
179
  (vec)[1] = AIR_CAST(TT, y*n); \
180
  (vec)[2] = AIR_CAST(TT, z*n)
181
182
#define _EVEN_CHECK_VtoQ(HNB, vec) \
183
  x = (vec)[0]; \
184
  y = (vec)[1]; \
185
  z = (vec)[2]; \
186
  L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z); \
187
  if (!L) { \
188
    return 0; \
189
  } \
190
  x /= L; \
191
  y /= L; \
192
  if (z > 0) { \
193
    xi = airIndex(-1.0, x, 1.0, ((1<<HNB)-1)); \
194
    yi = airIndex(-1.0 - 1.0/((1<<HNB)-1), y, 1.0 + 1.0/((1<<HNB)-1), \
195
                  (1<<HNB)); \
196
    ui = xi + yi - ((1<<(HNB-1))-1); \
197
    vi = xi - yi + (1<<(HNB-1)); \
198
  } \
199
  else { \
200
    xi = airIndex(-1.0 - 1.0/((1<<HNB)-1), x, 1.0 + 1.0/((1<<HNB)-1), \
201
                  (1<<HNB)); \
202
    yi = airIndex(-1, y, 1, ((1<<HNB)-1)); \
203
    ui = xi + yi - ((1<<(HNB-1))-1); \
204
    vi = xi - yi + ((1<<(HNB-1))-1); \
205
  } \
206
  return (vi << HNB) | ui
207
208
#define _EVEN_OCTA_QtoV(HNB, TT, vec, qn) \
209
  xi = (qn) & ((1<<HNB)-1); \
210
  yi = ((qn) >> HNB) & ((1<<HNB)-1); \
211
  x = AIR_AFFINE(-0.5, xi, ((1<<HNB)-1)+0.5, -1, 1); \
212
  y = AIR_AFFINE(-0.5, yi, ((1<<HNB)-1)+0.5, -1, 1); \
213
  z = 1 - AIR_ABS(x) - AIR_ABS(y); \
214
  if (z < 0) { \
215
    if (x > 0) { \
216
      x += z; \
217
    } else { \
218
      x -= z; \
219
    } \
220
    if (y > 0) { \
221
      y += z; \
222
    } else { \
223
      y -= z; \
224
    } \
225
  } \
226
  n = 1.0/sqrt(x*x + y*y + z*z); \
227
  (vec)[0] = AIR_CAST(TT, x*n); \
228
  (vec)[1] = AIR_CAST(TT, y*n); \
229
  (vec)[2] = AIR_CAST(TT, z*n)
230
231
#define _EVEN_OCTA_VtoQ(HNB, vec) \
232
  x = (vec)[0]; \
233
  y = (vec)[1]; \
234
  z = (vec)[2]; \
235
  L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z); \
236
  if (!L) { \
237
    return 0; \
238
  } \
239
  x /= L; \
240
  y /= L; \
241
  z /= L; \
242
  if (z < 0) { \
243
    if (x > 0) { \
244
      x -= z; \
245
    } else { \
246
      x += z; \
247
    } \
248
    if (y > 0) { \
249
      y -= z; \
250
    } else { \
251
      y += z; \
252
    } \
253
  } \
254
  xi = airIndex(-1, x, 1, (1<<HNB)); \
255
  yi = airIndex(-1, y, 1, (1<<HNB)); \
256
  return (yi << HNB) | xi
257
258
/* _16_CHECK_QtoV and _16_CHECK_VtoQ are not actually used */
259
#define _16_CHECK_QtoV(vec, qn) \
260
  ui = (qn) & 0xFF; \
261
  vi = ((qn) >> 8) & 0xFF; \
262
  u = AIR_AFFINE(0, ui, 255, -0.5, 0.5); \
263
  v = AIR_AFFINE(0, vi, 255, -0.5, 0.5); \
264
  x =  u + v; \
265
  y =  u - v; \
266
  z = 1 - AIR_ABS(x) - AIR_ABS(y); \
267
  /* would this be better served by a branch? */ \
268
  z *= (((ui ^ vi) & 0x01) << 1) - 1; \
269
  n = 1.0/sqrt(x*x + y*y + z*z); \
270
  (vec)[0] = x*n; \
271
  (vec)[1] = y*n; \
272
  (vec)[2] = z*n
273
274
#define _16_CHECK_VtoQ(vec) \
275
  x = (vec)[0]; \
276
  y = (vec)[1]; \
277
  z = (vec)[2]; \
278
  L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z); \
279
  if (!L) { \
280
    return 0; \
281
  } \
282
  x /= L; \
283
  y /= L; \
284
  if (z > 0) { \
285
    xi = airIndex(-1.0, x, 1.0, 255); \
286
    yi = airIndex(-1.0 - 1.0/255, y, 1.0 + 1.0/255, 256); \
287
    ui = xi + yi - 127; \
288
    vi = xi - yi + 128; \
289
  } \
290
  else { \
291
    xi = airIndex(-1.0 - 1.0/255, x, 1.0 + 1.0/255, 256); \
292
    yi = airIndex(-1, y, 1, 255); \
293
    ui = xi + yi - 127; \
294
    vi = xi - yi + 127; \
295
  }
296
297
void
298
_limnQN16checker_QNtoV_f(float *vec, unsigned int qn) {
299
  unsigned int ui, vi;
300
  double u, v, x, y, z, n;
301
302
  _EVEN_CHECK_QtoV(8, float, vec, qn);
303
}
304
305
void
306
_limnQN16checker_QNtoV_d(double *vec, unsigned int qn) {
307
  unsigned int ui, vi;
308
  double u, v, x, y, z, n;
309
310
  _EVEN_CHECK_QtoV(8, double, vec, qn);
311
}
312
313
unsigned int
314
_limnQN16checker_VtoQN_f(const float *vec) {
315
  double L, x, y, z;
316
  unsigned int xi, yi, ui, vi;
317
318
  _EVEN_CHECK_VtoQ(8, vec);
319
}
320
321
unsigned int
322
_limnQN16checker_VtoQN_d(const double *vec) {
323
  double L, x, y, z;
324
  unsigned int xi, yi, ui, vi;
325
326
  _EVEN_CHECK_VtoQ(8, vec);
327
}
328
329
/* ----------------------------------------------------------------  */
330
331
void
332
_limnQN16octa_QNtoV_f(float *vec, unsigned int qn) {
333
  unsigned int xi, yi;
334
  double x, y, z, n;
335
336
  _EVEN_OCTA_QtoV(8, float, vec, qn);
337
}
338
339
void
340
_limnQN16octa_QNtoV_d(double *vec, unsigned int qn) {
341
  unsigned int xi, yi;
342
  double x, y, z, n;
343
344
  _EVEN_OCTA_QtoV(8, double, vec, qn);
345
}
346
347
unsigned int
348
_limnQN16octa_VtoQN_f(const float *vec) {
349
  double L, x, y, z;
350
  unsigned int xi, yi;
351
352
  _EVEN_OCTA_VtoQ(8, vec);
353
}
354
355
unsigned int
356
_limnQN16octa_VtoQN_d(const double *vec) {
357
  double L, x, y, z;
358
  unsigned int xi, yi;
359
360
  _EVEN_OCTA_VtoQ(8, vec);
361
}
362
363
/* ----------------------------------------------------------------  */
364
365
/* 15 bit --> HNB == 7 */
366
367
#define _ODD_OCTA_QtoV(HNB, TT, vec, qn) \
368
  ui = qn & ((1<<HNB)-1); \
369
  vi = (qn >> HNB) & ((1<<HNB)-1); \
370
  zi = (qn >> (2*HNB)) & 0x01; \
371
  u = AIR_CAST(TT, AIR_AFFINE(-0.5, ui, ((1<<HNB)-1)+0.5, -0.5, 0.5)); \
372
  v = AIR_CAST(TT, AIR_AFFINE(-0.5, vi, ((1<<HNB)-1)+0.5, -0.5, 0.5)); \
373
  x =  u + v; \
374
  y =  u - v; \
375
  z = 1 - AIR_ABS(x) - AIR_ABS(y); \
376
  z *= AIR_CAST(int, zi << 1) - 1;    \
377
  n = AIR_CAST(TT, 1.0/sqrt(x*x + y*y + z*z)); \
378
  vec[0] = AIR_CAST(TT, x*n); \
379
  vec[1] = AIR_CAST(TT, y*n); \
380
  vec[2] = AIR_CAST(TT, z*n)
381
382
#define _ODD_OCTA_VtoQ(HNB, vec) \
383
  x = vec[0]; \
384
  y = vec[1]; \
385
  z = vec[2]; \
386
  L = AIR_ABS(x) + AIR_ABS(y) + AIR_ABS(z); \
387
  if (!L) { \
388
    return 0; \
389
  } \
390
  x /= L; \
391
  y /= L; \
392
  u = x + y; \
393
  v = x - y; \
394
  ui = airIndex(-1, u, 1, (1<<HNB)); \
395
  vi = airIndex(-1, v, 1, (1<<HNB)); \
396
  zi = (z > 0); \
397
  return (zi << (2*HNB)) | (vi << HNB) | ui
398
399
400
void
401
_limnQN15octa_QNtoV_f(float *vec, unsigned int qn) {
402
  unsigned int ui, vi, zi;
403
  float u, v, x, y, z, n;
404
405
  _ODD_OCTA_QtoV(7, float, vec, qn);
406
}
407
408
void
409
_limnQN15octa_QNtoV_d(double *vec, unsigned int qn) {
410
  unsigned int ui, vi, zi;
411
  double u, v, x, y, z, n;
412
413
  _ODD_OCTA_QtoV(7, double, vec, qn);
414
}
415
416
unsigned int
417
_limnQN15octa_VtoQN_f(const float *vec) {
418
  float L, u, v, x, y, z;
419
  unsigned int ui, vi, zi;
420
421
  _ODD_OCTA_VtoQ(7, vec);
422
}
423
424
unsigned int
425
_limnQN15octa_VtoQN_d(const double *vec) {
426
  double L, u, v, x, y, z;
427
  unsigned int ui, vi, zi;
428
429
  _ODD_OCTA_VtoQ(7, vec);
430
}
431
432
/* ----------------------------------------------------------- */
433
434
void
435
_limnQN14checker_QNtoV_f(float *vec, unsigned int qn) {
436
  unsigned int ui, vi;
437
  double u, v, x, y, z, n;
438
439
  _EVEN_CHECK_QtoV(7, float, vec, qn);
440
}
441
442
void
443
_limnQN14checker_QNtoV_d(double *vec, unsigned int qn) {
444
  unsigned int ui, vi;
445
  double u, v, x, y, z, n;
446
447
  _EVEN_CHECK_QtoV(7, double, vec, qn);
448
}
449
450
unsigned int
451
_limnQN14checker_VtoQN_f(const float *vec) {
452
  double L, x, y, z;
453
  unsigned int xi, yi, ui, vi;
454
455
  _EVEN_CHECK_VtoQ(7, vec);
456
}
457
458
unsigned int
459
_limnQN14checker_VtoQN_d(const double *vec) {
460
  double L, x, y, z;
461
  unsigned int xi, yi, ui, vi;
462
463
  _EVEN_CHECK_VtoQ(7, vec);
464
}
465
466
/* ----------------------------------------------------------------  */
467
468
void
469
_limnQN14octa_QNtoV_f(float *vec, unsigned int qn) {
470
  unsigned int xi, yi;
471
  double x, y, z, n;
472
473
  _EVEN_OCTA_QtoV(7, float, vec, qn);
474
}
475
476
void
477
_limnQN14octa_QNtoV_d(double *vec, unsigned int qn) {
478
  unsigned int xi, yi;
479
  double x, y, z, n;
480
481
  _EVEN_OCTA_QtoV(7, double, vec, qn);
482
}
483
484
unsigned int
485
_limnQN14octa_VtoQN_f(const float *vec) {
486
  double L, x, y, z;
487
  unsigned int xi, yi;
488
489
  _EVEN_OCTA_VtoQ(7, vec);
490
}
491
492
unsigned int
493
_limnQN14octa_VtoQN_d(const double *vec) {
494
  double L, x, y, z;
495
  unsigned int xi, yi;
496
497
  _EVEN_OCTA_VtoQ(7, vec);
498
}
499
500
/* ----------------------------------------------------------- */
501
502
void
503
_limnQN13octa_QNtoV_f(float *vec, unsigned int qn) {
504
  unsigned int ui, vi, zi;
505
  float u, v, x, y, z, n;
506
507
  _ODD_OCTA_QtoV(6, float, vec, qn);
508
}
509
510
void
511
_limnQN13octa_QNtoV_d(double *vec, unsigned int qn) {
512
  unsigned int ui, vi, zi;
513
  double u, v, x, y, z, n;
514
515
  _ODD_OCTA_QtoV(6, double, vec, qn);
516
}
517
518
unsigned int
519
_limnQN13octa_VtoQN_f(const float *vec) {
520
  float L, u, v, x, y, z;
521
  unsigned int ui, vi, zi;
522
523
  _ODD_OCTA_VtoQ(6, vec);
524
}
525
526
unsigned int
527
_limnQN13octa_VtoQN_d(const double *vec) {
528
  double L, u, v, x, y, z;
529
  unsigned int ui, vi, zi;
530
531
  _ODD_OCTA_VtoQ(6, vec);
532
}
533
534
/* ----------------------------------------------------------- */
535
536
void
537
_limnQN12checker_QNtoV_f(float *vec, unsigned int qn) {
538
  unsigned int ui, vi;
539
  double u, v, x, y, z, n;
540
541
  _EVEN_CHECK_QtoV(6, float, vec, qn);
542
}
543
544
void
545
_limnQN12checker_QNtoV_d(double *vec, unsigned int qn) {
546
  unsigned int ui, vi;
547
  double u, v, x, y, z, n;
548
549
  _EVEN_CHECK_QtoV(6, double, vec, qn);
550
}
551
552
unsigned int
553
_limnQN12checker_VtoQN_f(const float *vec) {
554
  double L, x, y, z;
555
  unsigned int xi, yi, ui, vi;
556
557
  _EVEN_CHECK_VtoQ(6, vec);
558
}
559
560
unsigned int
561
_limnQN12checker_VtoQN_d(const double *vec) {
562
  double L, x, y, z;
563
  unsigned int xi, yi, ui, vi;
564
565
  _EVEN_CHECK_VtoQ(6, vec);
566
}
567
568
/* ----------------------------------------------------------------  */
569
570
void
571
_limnQN12octa_QNtoV_f(float *vec, unsigned int qn) {
572
  unsigned int xi, yi;
573
  double x, y, z, n;
574
575
  _EVEN_OCTA_QtoV(6, float, vec, qn);
576
}
577
578
void
579
_limnQN12octa_QNtoV_d(double *vec, unsigned int qn) {
580
  unsigned int xi, yi;
581
  double x, y, z, n;
582
583
  _EVEN_OCTA_QtoV(6, double, vec, qn);
584
}
585
586
unsigned int
587
_limnQN12octa_VtoQN_f(const float *vec) {
588
  double L, x, y, z;
589
  unsigned int xi, yi;
590
591
  _EVEN_OCTA_VtoQ(6, vec);
592
}
593
594
unsigned int
595
_limnQN12octa_VtoQN_d(const double *vec) {
596
  double L, x, y, z;
597
  unsigned int xi, yi;
598
599
  _EVEN_OCTA_VtoQ(6, vec);
600
}
601
602
/* ----------------------------------------------------------- */
603
604
void
605
_limnQN11octa_QNtoV_f(float *vec, unsigned int qn) {
606
  unsigned int ui, vi, zi;
607
  float u, v, x, y, z, n;
608
609
  _ODD_OCTA_QtoV(5, float, vec, qn);
610
}
611
612
void
613
_limnQN11octa_QNtoV_d(double *vec, unsigned int qn) {
614
  unsigned int ui, vi, zi;
615
  double u, v, x, y, z, n;
616
617
  _ODD_OCTA_QtoV(5, double, vec, qn);
618
}
619
620
unsigned int
621
_limnQN11octa_VtoQN_f(const float *vec) {
622
  float L, u, v, x, y, z;
623
  unsigned int ui, vi, zi;
624
625
  _ODD_OCTA_VtoQ(5, vec);
626
}
627
628
unsigned int
629
_limnQN11octa_VtoQN_d(const double *vec) {
630
  double L, u, v, x, y, z;
631
  unsigned int ui, vi, zi;
632
633
  _ODD_OCTA_VtoQ(5, vec);
634
}
635
636
/* ----------------------------------------------------------- */
637
638
void
639
_limnQN10checker_QNtoV_f(float *vec, unsigned int qn) {
640
  unsigned int ui, vi;
641
  double u, v, x, y, z, n;
642
643
  _EVEN_CHECK_QtoV(5, float, vec, qn);
644
}
645
646
void
647
_limnQN10checker_QNtoV_d(double *vec, unsigned int qn) {
648
  unsigned int ui, vi;
649
  double u, v, x, y, z, n;
650
651
  _EVEN_CHECK_QtoV(5, double, vec, qn);
652
}
653
654
unsigned int
655
_limnQN10checker_VtoQN_f(const float *vec) {
656
  double L, x, y, z;
657
  unsigned int xi, yi, ui, vi;
658
659
  _EVEN_CHECK_VtoQ(5, vec);
660
}
661
662
unsigned int
663
_limnQN10checker_VtoQN_d(const double *vec) {
664
  double L, x, y, z;
665
  unsigned int xi, yi, ui, vi;
666
667
  _EVEN_CHECK_VtoQ(5, vec);
668
}
669
670
/* ----------------------------------------------------------------  */
671
672
void
673
_limnQN10octa_QNtoV_f(float *vec, unsigned int qn) {
674
  unsigned int xi, yi;
675
  double x, y, z, n;
676
677
  _EVEN_OCTA_QtoV(5, float, vec, qn);
678
}
679
680
void
681
_limnQN10octa_QNtoV_d(double *vec, unsigned int qn) {
682
  unsigned int xi, yi;
683
  double x, y, z, n;
684
685
  _EVEN_OCTA_QtoV(5, double, vec, qn);
686
}
687
688
unsigned int
689
_limnQN10octa_VtoQN_f(const float *vec) {
690
  double L, x, y, z;
691
  unsigned int xi, yi;
692
693
  _EVEN_OCTA_VtoQ(5, vec);
694
}
695
696
unsigned int
697
_limnQN10octa_VtoQN_d(const double *vec) {
698
  double L, x, y, z;
699
  unsigned int xi, yi;
700
701
  _EVEN_OCTA_VtoQ(5, vec);
702
}
703
704
/* ----------------------------------------------------------- */
705
void
706
_limnQN9octa_QNtoV_f(float *vec, unsigned int qn) {
707
  unsigned int ui, vi, zi;
708
  float u, v, x, y, z, n;
709
710
  _ODD_OCTA_QtoV(4, float, vec, qn);
711
}
712
713
void
714
_limnQN9octa_QNtoV_d(double *vec, unsigned int qn) {
715
  unsigned int ui, vi, zi;
716
  double u, v, x, y, z, n;
717
718
  _ODD_OCTA_QtoV(4, double, vec, qn);
719
}
720
721
unsigned int
722
_limnQN9octa_VtoQN_f(const float *vec) {
723
  float L, u, v, x, y, z;
724
  unsigned int ui, vi, zi;
725
726
  _ODD_OCTA_VtoQ(4, vec);
727
}
728
729
unsigned int
730
_limnQN9octa_VtoQN_d(const double *vec) {
731
  double L, u, v, x, y, z;
732
  unsigned int ui, vi, zi;
733
734
  _ODD_OCTA_VtoQ(4, vec);
735
}
736
737
/* ----------------------------------------------------------- */
738
739
void
740
_limnQN8checker_QNtoV_f(float *vec, unsigned int qn) {
741
  unsigned int ui, vi;
742
  double u, v, x, y, z, n;
743
744
  _EVEN_CHECK_QtoV(4, float, vec, qn);
745
}
746
747
void
748
_limnQN8checker_QNtoV_d(double *vec, unsigned int qn) {
749
  unsigned int ui, vi;
750
  double u, v, x, y, z, n;
751
752
  _EVEN_CHECK_QtoV(4, double, vec, qn);
753
}
754
755
unsigned int
756
_limnQN8checker_VtoQN_f(const float *vec) {
757
  double L, x, y, z;
758
  unsigned int xi, yi, ui, vi;
759
760
  _EVEN_CHECK_VtoQ(4, vec);
761
}
762
763
unsigned int
764
_limnQN8checker_VtoQN_d(const double *vec) {
765
  double L, x, y, z;
766
  unsigned int xi, yi, ui, vi;
767
768
  _EVEN_CHECK_VtoQ(4, vec);
769
}
770
771
/* ----------------------------------------------------------------  */
772
773
void
774
_limnQN8octa_QNtoV_f(float *vec, unsigned int qn) {
775
  unsigned int xi, yi;
776
  double x, y, z, n;
777
778
  _EVEN_OCTA_QtoV(4, float, vec, qn);
779
}
780
781
void
782
_limnQN8octa_QNtoV_d(double *vec, unsigned int qn) {
783
  unsigned int xi, yi;
784
  double x, y, z, n;
785
786
  _EVEN_OCTA_QtoV(4, double, vec, qn);
787
}
788
789
unsigned int
790
_limnQN8octa_VtoQN_f(const float *vec) {
791
  double L, x, y, z;
792
  unsigned int xi, yi;
793
794
  _EVEN_OCTA_VtoQ(4, vec);
795
}
796
797
unsigned int
798
_limnQN8octa_VtoQN_d(const double *vec) {
799
  double L, x, y, z;
800
  unsigned int xi, yi;
801
802
  _EVEN_OCTA_VtoQ(4, vec);
803
}
804
805
/* ----------------------------------------------------------- */
806
807
void (*
808
limnQNtoV_f[LIMN_QN_MAX+1])(float *, unsigned int) = {
809
  NULL,
810
  _limnQN16simple_QNtoV_f,
811
  _limnQN16border1_QNtoV_f,
812
  _limnQN16checker_QNtoV_f,
813
  _limnQN16octa_QNtoV_f,
814
  _limnQN15octa_QNtoV_f,
815
  _limnQN14checker_QNtoV_f,
816
  _limnQN14octa_QNtoV_f,
817
  _limnQN13octa_QNtoV_f,
818
  _limnQN12checker_QNtoV_f,
819
  _limnQN12octa_QNtoV_f,
820
  _limnQN11octa_QNtoV_f,
821
  _limnQN10checker_QNtoV_f,
822
  _limnQN10octa_QNtoV_f,
823
  _limnQN9octa_QNtoV_f,
824
  _limnQN8checker_QNtoV_f,
825
  _limnQN8octa_QNtoV_f
826
};
827
828
void (*
829
limnQNtoV_d[LIMN_QN_MAX+1])(double *, unsigned int) = {
830
  NULL,
831
  NULL,
832
  NULL,
833
  _limnQN16checker_QNtoV_d,
834
  _limnQN16octa_QNtoV_d,
835
  _limnQN15octa_QNtoV_d,
836
  _limnQN14checker_QNtoV_d,
837
  _limnQN14octa_QNtoV_d,
838
  _limnQN13octa_QNtoV_d,
839
  _limnQN12checker_QNtoV_d,
840
  _limnQN12octa_QNtoV_d,
841
  _limnQN11octa_QNtoV_d,
842
  _limnQN10checker_QNtoV_d,
843
  _limnQN10octa_QNtoV_d,
844
  _limnQN9octa_QNtoV_d,
845
  _limnQN8checker_QNtoV_d,
846
  _limnQN8octa_QNtoV_d
847
};
848
849
unsigned int (*
850
limnVtoQN_f[LIMN_QN_MAX+1])(const float *vec) = {
851
  NULL,
852
  _limnQN16simple_VtoQN_f,
853
  _limnQN16border1_VtoQN_f,
854
  _limnQN16checker_VtoQN_f,
855
  _limnQN16octa_VtoQN_f,
856
  _limnQN15octa_VtoQN_f,
857
  _limnQN14checker_VtoQN_f,
858
  _limnQN14octa_VtoQN_f,
859
  _limnQN13octa_VtoQN_f,
860
  _limnQN12checker_VtoQN_f,
861
  _limnQN12octa_VtoQN_f,
862
  _limnQN11octa_VtoQN_f,
863
  _limnQN10checker_VtoQN_f,
864
  _limnQN10octa_VtoQN_f,
865
  _limnQN9octa_VtoQN_f,
866
  _limnQN8checker_VtoQN_f,
867
  _limnQN8octa_VtoQN_f
868
};
869
870
unsigned int (*
871
limnVtoQN_d[LIMN_QN_MAX+1])(const double *vec) = {
872
  NULL,
873
  NULL,
874
  NULL,
875
  _limnQN16checker_VtoQN_d,
876
  _limnQN16octa_VtoQN_d,
877
  _limnQN15octa_VtoQN_d,
878
  _limnQN14checker_VtoQN_d,
879
  _limnQN14octa_VtoQN_d,
880
  _limnQN13octa_VtoQN_d,
881
  _limnQN12checker_VtoQN_d,
882
  _limnQN12octa_VtoQN_d,
883
  _limnQN11octa_VtoQN_d,
884
  _limnQN10checker_VtoQN_d,
885
  _limnQN10octa_VtoQN_d,
886
  _limnQN9octa_VtoQN_d,
887
  _limnQN8checker_VtoQN_d,
888
  _limnQN8octa_VtoQN_d
889
};
890
891
unsigned int
892
limnQNBins[LIMN_QN_MAX+1] = {
893
  0,
894
  (1 << 16),
895
  (1 << 16),
896
  (1 << 16),
897
  (1 << 16),
898
  (1 << 15),
899
  (1 << 14),
900
  (1 << 14),
901
  (1 << 13),
902
  (1 << 12),
903
  (1 << 12),
904
  (1 << 11),
905
  (1 << 10),
906
  (1 << 10),
907
  (1 << 9),
908
  (1 << 8),
909
  (1 << 8)
910
};
911
912
/*
913
** can use via test/tqn:
914
  foreach Q ( 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 )
915
    echo $Q
916
    test/tqn -s 400 -q $Q \
917
     | unu tile -a 2 0 1 -s 2 3 \
918
     | unu tile -a 2 0 1 -s 2 1 | unu quantize -b 8 -o q${Q}.png
919
  end
920
*/
921
int
922
limnQNDemo(Nrrd *nqn, unsigned int reso, int qni) {
923
  static const char me[]="limnQNDemo";
924
  unsigned int ui, vi, oi;
925
  double *qdata, ll, uu, vv, ww, vecd[3], unqd[3];
926
  float vecf[3], unqf[3];
927
928
  if (!nqn) {
929
    biffAddf(LIMN, "%s: got NULL pointer", me);
930
    return 1;
931
  }
932
  if (nrrdMaybeAlloc_va(nqn, nrrdTypeDouble, 4,
933
                        AIR_CAST(size_t, reso),
934
                        AIR_CAST(size_t, reso),
935
                        AIR_CAST(size_t, 6),
936
                        AIR_CAST(size_t, 2))){
937
    biffMovef(LIMN, NRRD, "%s: couldn't alloc output", me);
938
    return 1;
939
  }
940
  if (!(limnQNUnknown < qni && qni < limnQNLast)) {
941
    biffAddf(LIMN, "%s: qni %d not in valid range [%d,%d]", me,
942
             qni, limnQNUnknown+1, limnQNLast-1);
943
    return 1;
944
  }
945
  qdata = AIR_CAST(double *, nqn->data);
946
  for (oi=0; oi<6; oi++) {
947
    for (vi=0; vi<reso; vi++) {
948
      vv = AIR_AFFINE(0, vi, reso-1, 1, -1);
949
      for (ui=0; ui<reso; ui++) {
950
        uu = AIR_AFFINE(0, ui, reso-1, -1, 1);
951
        ll = uu*uu + vv*vv;
952
        if (ll <= 1) {
953
          ww = sqrt(1 - ll);
954
          if (oi % 2) {
955
            ww *= -1;
956
          }
957
        } else {
958
          continue;
959
        }
960
        switch (oi) {
961
        case 0:
962
        case 1:
963
          ELL_3V_SET(vecd, uu, vv, ww);
964
          ELL_3V_SET_TT(vecf, float, uu, vv, ww);
965
          break;
966
        case 2:
967
        case 3:
968
          ELL_3V_SET(vecd, uu, ww, vv);
969
          ELL_3V_SET_TT(vecf, float, uu, ww, vv);
970
          break;
971
        case 4:
972
        case 5:
973
          ELL_3V_SET(vecd, ww, uu, vv);
974
          ELL_3V_SET_TT(vecf, float, ww, uu, vv);
975
          break;
976
        }
977
        if (limnVtoQN_d[qni] && limnQNtoV_d[qni]) {
978
          (limnQNtoV_d[qni])(unqd, (limnVtoQN_d[qni])(vecd));
979
          qdata[ui + reso*(vi + reso*(oi + 6))] = ell_3v_angle_d(unqd, vecd);
980
        }
981
        if (limnVtoQN_f[qni] && limnQNtoV_f[qni]) {
982
          (limnQNtoV_f[qni])(unqf, (limnVtoQN_f[qni])(vecf));
983
          qdata[ui + reso*(vi + reso*(oi + 0))] = ell_3v_angle_f(unqf, vecf);
984
        }
985
      }
986
    }
987
  }
988
  return 0;
989
}