1 /* 2 Teem: Tools to process and visualize scientific data and images . 3 Copyright (C) 2011, 2010, 2009, 2008 Thomas Schultz 4 Copyright (C) 2010, 2009, 2008 Gordon Kindlmann 5 6 This library is free software; you can redistribute it and/or 7 modify it under the terms of the GNU Lesser General Public License 8 (LGPL) as published by the Free Software Foundation; either 9 version 2.1 of the License, or (at your option) any later version. 10 The terms of redistributing and/or modifying this software also 11 include exceptions to the LGPL that facilitate static linking. 12 13 This library is distributed in the hope that it will be useful, 14 but WITHOUT ANY WARRANTY; without even the implied warranty of 15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 16 Lesser General Public License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with this library; if not, write to Free Software Foundation, Inc., 20 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 21 */ 22 23 /* Implementation of spherical harmonics and their relation to tensors */ 24 25 #include "tijk.h" 26 27 #include "convertQuietPush.h" 28 29 #define TIJK_TABLE_TYPE 0 /* create double version */ 30 #include "shtables.h" 31 #undef TIJK_TABLE_TYPE /* explicitly undef to avoid compiler warnings */ 32 #define TIJK_TABLE_TYPE 1 /* create float version */ 33 #include "shtables.h" 34 #undef TIJK_TABLE_TYPE 35 36 const unsigned int tijk_max_esh_order=8; 37 /* for order 8, the maximum number of coefficients is 45 */ 38 #define _TIJK_MAX_ESH_LEN 45 39 /* number of coefficients for order i/2 */ 40 const unsigned int tijk_esh_len[5]={1,6,15,28,45}; 41 42 #define TIJK_EVAL_ESH_BASIS(TYPE, SUF) \ 43 unsigned int \ 44 tijk_eval_esh_basis_##SUF(TYPE *res, unsigned int order, \ 45 TYPE theta, TYPE phi) { \ 46 TYPE stheta=sin(theta); \ 47 TYPE ctheta=cos(theta); \ 48 TYPE stheta2=stheta*stheta, stheta4=stheta2*stheta2; \ 49 TYPE ctheta2=ctheta*ctheta, ctheta4=ctheta2*ctheta2; \ 50 res[0]=1.0/sqrt(4.0*AIR_PI); \ 51 if (order<2) \ 52 return 0; \ 53 res[5]=res[1]=sqrt(15.0/(16.0*AIR_PI))*stheta2; \ 54 res[5]*=sin(2*phi); res[1]*=cos(2*phi); \ 55 res[4]=res[2]=-sqrt(15.0/(4.0*AIR_PI))*ctheta*stheta; \ 56 res[4]*=sin(phi); res[2]*=cos(phi); \ 57 res[3]=sqrt(5.0/(16.0*AIR_PI))*(3.0*ctheta2-1.0); \ 58 if (order<4) \ 59 return 2; \ 60 res[14]=res[6]=sqrt(315.0/(256.0*AIR_PI))*stheta4; \ 61 res[14]*=sin(4*phi); res[6]*=cos(4*phi); \ 62 res[13]=res[7]=-sqrt(315.0/(32.0*AIR_PI))*ctheta*stheta*stheta2; \ 63 res[13]*=sin(3*phi); res[7]*=cos(3*phi); \ 64 res[12]=res[8]=sqrt(45.0/(64.0*AIR_PI))*(7*ctheta2-1)*stheta2; \ 65 res[12]*=sin(2*phi); res[8]*=cos(2*phi); \ 66 res[11]=res[9]=-sqrt(45.0/(32.0*AIR_PI))*(7*ctheta2*ctheta-3*ctheta)*stheta; \ 67 res[11]*=sin(phi); res[9]*=cos(phi); \ 68 res[10]=sqrt(9.0/(256.0*AIR_PI))*(35.0*ctheta4-30.0*ctheta2+3.0); \ 69 if (order<6) \ 70 return 4; \ 71 res[27]=res[15]=1.0/64.0*sqrt(6006.0/AIR_PI)*stheta4*stheta2; \ 72 res[27]*=sin(6*phi); res[15]*=cos(6*phi); \ 73 res[26]=res[16]=-3.0/32.0*sqrt(2002.0/AIR_PI)*stheta4*stheta*ctheta; \ 74 res[26]*=sin(5*phi); res[16]*=cos(5*phi); \ 75 res[25]=res[17]=3.0/32.0*sqrt(91.0/AIR_PI)*stheta4*(11*ctheta2-1.0); \ 76 res[25]*=sin(4*phi); res[17]*=cos(4*phi); \ 77 res[24]=res[18]=-1.0/32.0*sqrt(2730.0/AIR_PI)*stheta2*stheta*(11*ctheta2*ctheta-3*ctheta); \ 78 res[24]*=sin(3*phi); res[18]*=cos(3*phi); \ 79 res[23]=res[19]=1.0/64.0*sqrt(2730.0/AIR_PI)*stheta2*(33*ctheta4-18*ctheta2+1.0); \ 80 res[23]*=sin(2*phi); res[19]*=cos(2*phi); \ 81 res[22]=res[20]=-1.0/16.0*sqrt(273.0/AIR_PI)*stheta*(33*ctheta4*ctheta-30.0*ctheta2*ctheta+5*ctheta); \ 82 res[22]*=sin(phi); res[20]*=cos(phi); \ 83 res[21]=1.0/32.0*sqrt(13.0/AIR_PI)*(231*ctheta4*ctheta2-315*ctheta4+105*ctheta2-5.0); \ 84 if (order<8) \ 85 return 6; \ 86 res[44]=res[28]=3.0/256.0*sqrt(12155.0/AIR_PI)*stheta4*stheta4; \ 87 res[44]*=sin(8*phi); res[28]*=cos(8*phi); \ 88 res[43]=res[29]=-3.0/64.0*sqrt(12155.0/AIR_PI)*stheta4*stheta2*stheta*ctheta; \ 89 res[43]*=sin(7*phi); res[29]*=cos(7*phi); \ 90 res[42]=res[30]=1.0/128.0*sqrt(2*7293.0/AIR_PI)*stheta4*stheta2*(15*ctheta2-1); \ 91 res[42]*=sin(6*phi); res[30]*=cos(6*phi); \ 92 res[41]=res[31]=-3.0/64.0*sqrt(17017.0/AIR_PI)*stheta4*stheta*ctheta*(5*ctheta2-1); \ 93 res[41]*=sin(5*phi); res[31]*=cos(5*phi); \ 94 res[40]=res[32]=3.0/128.0*sqrt(1309.0/AIR_PI)*stheta4*(65*ctheta4-26*ctheta2+1); \ 95 res[40]*=sin(4*phi); res[32]*=cos(4*phi); \ 96 res[39]=res[33]=-1.0/64.0*sqrt(19635.0/AIR_PI)*stheta2*stheta*ctheta*(39*ctheta4-26*ctheta2+3); \ 97 res[39]*=sin(3*phi); res[33]*=cos(3*phi); \ 98 res[38]=res[34]=3.0/128.0*sqrt(2*595.0/AIR_PI)*stheta2*(143*ctheta4*ctheta2-143*ctheta4+33*ctheta2-1); \ 99 res[38]*=sin(2*phi); res[34]*=cos(2*phi); \ 100 res[37]=res[35]=-3.0/64.0*sqrt(17.0/AIR_PI)*stheta*ctheta*(715*ctheta4*ctheta2-1001*ctheta4+385*ctheta2-35); \ 101 res[37]*=sin(phi); res[35]*=cos(phi); \ 102 res[36]=1.0/256.0*sqrt(17.0/AIR_PI)*(6435*ctheta4*ctheta4-12012*ctheta2*ctheta4+6930*ctheta4-1260*ctheta2+35); \ 103 return 8; /* higher orders currently not implemented */ \ 104 } 105 106 TIJK_EVAL_ESH_BASIS(double, d) 107 TIJK_EVAL_ESH_BASIS(float, f) 108 109 #define TIJK_EVAL_ESH(TYPE, SUF) \ 110 TYPE \ 111 tijk_eval_esh_##SUF(TYPE *coeffs, unsigned int order, \ 112 TYPE theta, TYPE phi) { \ 113 TYPE basis[_TIJK_MAX_ESH_LEN]; \ 114 TYPE res=0.0; \ 115 unsigned int i; \ 116 if (order!=tijk_eval_esh_basis_##SUF(basis, order, theta, phi)) \ 117 return 0; /* there has been an error. */ \ 118 for (i=0; i<tijk_esh_len[order/2]; i++) \ 119 res+=basis[i]*coeffs[i]; \ 120 return res; \ 121 } 122 123 TIJK_EVAL_ESH(double, d) 124 TIJK_EVAL_ESH(float, f) 125 126 #define TIJK_ESH_SP(TYPE, SUF) \ 127 TYPE \ 128 tijk_esh_sp_##SUF(TYPE *A, TYPE *B, unsigned int order) { \ 129 TYPE res=0.0; \ 130 if (order<=tijk_max_esh_order) { \ 131 unsigned int i; \ 132 for (i=0; i<tijk_esh_len[order/2]; i++) { \ 133 res+=A[i]*B[i]; \ 134 } \ 135 } \ 136 return res; \ 137 } 138 139 TIJK_ESH_SP(double, d) 140 TIJK_ESH_SP(float, f) 141 142 /* DOES NOT work in-place (with res==ten) */ 143 #define TIJK_3D_SYM_TO_ESH(TYPE, SUF) \ 144 int \ 145 tijk_3d_sym_to_esh_##SUF(TYPE *res, const TYPE *ten, \ 146 const tijk_type *type) { \ 147 const TYPE *m; /* conversion matrix */ \ 148 unsigned int i, j, n; \ 149 if (res==ten) return -1; \ 150 n=type->num; \ 151 if (type==tijk_2o3d_sym) { \ 152 m=_tijk_sym2esh_o2_##SUF; \ 153 } else if (type==tijk_4o3d_sym) { \ 154 m=_tijk_sym2esh_o4_##SUF; \ 155 } else if (type==tijk_6o3d_sym) { \ 156 m=_tijk_sym2esh_o6_##SUF; \ 157 } else if (type==tijk_8o3d_sym) { \ 158 m=_tijk_sym2esh_o8_##SUF; \ 159 } else { \ 160 return -1; /* cannot do conversion */ \ 161 } \ 162 for (i=0; i<n; i++) { \ 163 res[i]=0; \ 164 for (j=0; j<n; j++) { \ 165 res[i]+=m[n*i+j]*ten[j]; \ 166 } \ 167 } \ 168 return type->order; \ 169 } 170 171 TIJK_3D_SYM_TO_ESH(double, d) 172 TIJK_3D_SYM_TO_ESH(float, f) 173 174 /* DOES NOT work in-place (with res==sh) */ 175 #define TIJK_ESH_TO_3D_SYM(TYPE, SUF) \ 176 const tijk_type* \ 177 tijk_esh_to_3d_sym_##SUF(TYPE *res, const TYPE *sh, unsigned int order) { \ 178 const TYPE *m; /* conversion matrix */ \ 179 const tijk_type *type; \ 180 unsigned int i, j, n; \ 181 if (res==sh) return NULL; \ 182 switch (order) { \ 183 case 2: \ 184 m=_tijk_esh2sym_o2_##SUF; \ 185 type=tijk_2o3d_sym; \ 186 break; \ 187 case 4: \ 188 m=_tijk_esh2sym_o4_##SUF; \ 189 type=tijk_4o3d_sym; \ 190 break; \ 191 case 6: \ 192 m=_tijk_esh2sym_o6_##SUF; \ 193 type=tijk_6o3d_sym; \ 194 break; \ 195 case 8: \ 196 m=_tijk_esh2sym_o8_##SUF; \ 197 type=tijk_8o3d_sym; \ 198 break; \ 199 default: \ 200 return NULL; /* cannot do the conversion */ \ 201 } \ 202 n=tijk_esh_len[order/2]; \ 203 for (i=0; i<n; i++) { \ 204 res[i]=0; \ 205 for (j=0; j<n; j++) { \ 206 res[i]+=m[n*i+j]*sh[j]; \ 207 } \ 208 } \ 209 return type; \ 210 } 211 212 TIJK_ESH_TO_3D_SYM(double, d) 213 TIJK_ESH_TO_3D_SYM(float, f) 214 215 /* Convolve in with kernel and write to out (in==out is permitted) */ 216 #define TIJK_ESH_CONVOLVE(TYPE, SUF) \ 217 void tijk_esh_convolve_##SUF(TYPE *out, const TYPE *in, \ 218 const TYPE *kernel, unsigned int order) { \ 219 unsigned int o, idx=0; \ 220 for (o=0; o<=order/2; o++) { \ 221 while (idx<tijk_esh_len[o]) { \ 222 *(out++)=*(in++)*kernel[o]; \ 223 idx++; \ 224 } \ 225 } \ 226 } 227 228 TIJK_ESH_CONVOLVE(double, d) 229 TIJK_ESH_CONVOLVE(float, f) 230 231 /* Deconvolve in with kernel and write to out (in==out is permitted) */ 232 #define TIJK_ESH_DECONVOLVE(TYPE, SUF) \ 233 void tijk_esh_deconvolve_##SUF(TYPE *out, const TYPE *in, \ 234 const TYPE *kernel, unsigned int order) { \ 235 unsigned int o, idx=0; \ 236 for (o=0; o<=order/2; o++) { \ 237 while (idx<tijk_esh_len[o]) { \ 238 *(out++)=*(in++)/kernel[o]; \ 239 idx++; \ 240 } \ 241 } \ 242 } 243 244 TIJK_ESH_DECONVOLVE(double, d) 245 TIJK_ESH_DECONVOLVE(float, f) 246 247 /* Make a deconvolution kernel that turns a given signal (rotationally 248 * symmetric in "compressed" form, i.e., one coefficient per SH order) 249 * into a rank-1 term of given order. 250 * Return 0 upon success 251 * 1 if order is unsupported 252 * 2 if any of the given signal values is zero 253 */ 254 #define TIJK_ESH_MAKE_KERNEL_RANK1(TYPE, SUF) \ 255 int tijk_esh_make_kernel_rank1_##SUF(TYPE *kernel, const TYPE *signl, \ 256 unsigned int order) { \ 257 unsigned int i; \ 258 TYPE rank1[TIJK_TYPE_MAX_NUM], rank1sh[TIJK_TYPE_MAX_NUM]; \ 259 TYPE v[3]={0.0,0.0,1.0}; \ 260 /* Need to determine SH coefficients of a z-aligned rank-1 tensor */ \ 261 const tijk_type *type = NULL; \ 262 switch (order) { \ 263 case 2: type=tijk_2o3d_sym; break; \ 264 case 4: type=tijk_4o3d_sym; break; \ 265 case 6: type=tijk_6o3d_sym; break; \ 266 case 8: type=tijk_8o3d_sym; break; \ 267 default: return 1; \ 268 } \ 269 (*type->sym->make_rank1_##SUF)(rank1, 1.0, v); \ 270 tijk_3d_sym_to_esh_##SUF(rank1sh, rank1, type); \ 271 for (i=0; i<1+order/2; i++) \ 272 if (signl[i]==0) { \ 273 return 2; \ 274 } \ 275 kernel[0]=signl[0]/rank1sh[0]; \ 276 if (order>=2) { \ 277 kernel[1]=signl[1]/rank1sh[3]; \ 278 if (order>=4) { \ 279 kernel[2]=signl[2]/rank1sh[10]; \ 280 if (order>=6) { \ 281 kernel[3]=signl[3]/rank1sh[21]; \ 282 if (order>=8) { \ 283 kernel[4]=signl[4]/rank1sh[36]; \ 284 } \ 285 } \ 286 } \ 287 } \ 288 return 0; \ 289 } 290 291 TIJK_ESH_MAKE_KERNEL_RANK1(double, d) 292 TIJK_ESH_MAKE_KERNEL_RANK1(float, f) 293 294 /* Make a deconvolution kernel that turns a given signal (rotationally 295 * symmetric in "compressed" form, i.e., one coefficient per SH order) 296 * into a truncated delta peak of given order. 297 * Return 0 upon success 298 * 1 if order is unsupported 299 * 2 if any of the given signal values is zero 300 */ 301 #define TIJK_ESH_MAKE_KERNEL_DELTA(TYPE, SUF) \ 302 int tijk_esh_make_kernel_delta_##SUF(TYPE *kernel, const TYPE *signl, \ 303 unsigned int order) { \ 304 /* we need a truncated delta peak of given order */ \ 305 TYPE deltash[TIJK_TYPE_MAX_NUM]; \ 306 unsigned int i; \ 307 if (order>tijk_max_esh_order || order%2!=0) \ 308 return 1; \ 309 for (i=0; i<1+order/2; i++) \ 310 if (signl[i]==0) { \ 311 return 2; \ 312 } \ 313 tijk_eval_esh_basis_##SUF(deltash, order, 0, 0); \ 314 kernel[0]=signl[0]/deltash[0]; \ 315 if (order>=2) { \ 316 kernel[1]=signl[1]/deltash[3]; \ 317 if (order>=4) { \ 318 kernel[2]=signl[2]/deltash[10]; \ 319 if (order>=6) { \ 320 kernel[3]=signl[3]/deltash[21]; \ 321 if (order>=8) { \ 322 kernel[4]=signl[4]/deltash[36]; \ 323 } \ 324 } \ 325 } \ 326 } \ 327 return 0; \ 328 } 329 330 TIJK_ESH_MAKE_KERNEL_DELTA(double, d) 331 TIJK_ESH_MAKE_KERNEL_DELTA(float, f) 332 333 334 #include "convertQuietPop.h" 335