GCC Code Coverage Report

Directory: ./
File: lib/exploragram/hexdom/spherical_harmonics_l4.cpp
Date: 2026-09-07 02:28:19
Exec Total Coverage
Lines: 0 84 0.0%
Functions: 0 8 0.0%
Branches: 0 34 0.0%
Line Branch Exec Source
1 /*
2 * Copyright (c) 2000-2022 Inria
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * * Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 * * Redistributions in binary form must reproduce the above copyright notice,
11 * this list of conditions and the following disclaimer in the documentation
12 * and/or other materials provided with the distribution.
13 * * Neither the name of the ALICE Project-Team nor the names of its
14 * contributors may be used to endorse or promote products derived from this
15 * software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27 * POSSIBILITY OF SUCH DAMAGE.
28 *
29 * Contact: Bruno Levy
30 *
31 * https://www.inria.fr/fr/bruno-levy
32 *
33 * Inria,
34 * Domaine de Voluceau,
35 * 78150 Le Chesnay - Rocquencourt
36 * FRANCE
37 *
38 */
39
40 #include <exploragram/hexdom/spherical_harmonics_l4.h>
41 #include <exploragram/hexdom/frame.h>
42 #include <cmath>
43
44 namespace GEO {
45 // ____ _ _ _ _ _ _ _ _ _
46 // / ___| _ __ | |__ ___ _ __ (_) ___ __ _ | | | | | | __ _ _ __ _ __ ___ ___ _ __ (_) ___ | | | || |
47 // \___ \ | '_ \ | '_ \ / _ \ | '__| | | / __| / _` | | | | |_| | / _` | | '__| | '_ ` _ \ / _ \ | '_ \ | | / __| | | | || |_
48 // ___) | | |_) | | | | | | __/ | | | | | (__ | (_| | | | | _ | | (_| | | | | | | | | | | (_) | | | | | | | | (__ | |___ |__ _|
49 // |____/ | .__/ |_| |_| \___| |_| |_| \___| \__,_| |_| |_| |_| \__,_| |_| |_| |_| |_| \___/ |_| |_| |_| \___| |_____| |_|
50 // |_|
51
52 ✗ double SphericalHarmonicL4::basis(index_t id, const vec3& v) {
53 ✗ double x=v.x, x2=x*x;
54 ✗ double y=v.y, y2=y*y;
55 ✗ double z=v.z, z2=z*z;
56 ✗ if (id==0) return (3./4.)*std::sqrt(35./M_PI) * x*y*(x2-y2);
57 ✗ if (id==1) return (3./4.)*std::sqrt(35./(2.*M_PI)) * z*y*(3*x2-y2);
58 ✗ if (id==2) return (3./4.)*std::sqrt(5./M_PI) * x*y*(7.*z2-1);
59 ✗ if (id==3) return (3./4.)*std::sqrt(5./(2.*M_PI)) * z*y*(7.*z2-3);
60 ✗ if (id==4) return (3./16.)*std::sqrt(1./M_PI) * (35.*z2*z2-30.*z2+3);
61 ✗ if (id==5) return (3./4.)*std::sqrt(5./(2.*M_PI)) * z*x*(7.*z*z-3);
62 ✗ if (id==6) return (3./8.)*std::sqrt(5./M_PI) * (x2-y2)*(7.*z*z-1);
63 ✗ if (id==7) return (3./4.)*std::sqrt(35./(2.*M_PI)) * z*x*(x2-3.*y2);
64 ✗ if (id==8) return (3./16.)*std::sqrt(35./M_PI) * (x2*(x2-3.*y2)-y2*(3.*x2-y2));
65 ✗ geo_assert_not_reached;
66
67 }
68
69 ✗ void SphericalHarmonicL4::Rz(double a) {
70 ✗ SphericalHarmonicL4 copy(*this);
71 double s,c;
72 ✗ for (index_t i=0; i<4; i++) {
73 ✗ s = sin(a*(4-i));
74 ✗ c = cos(a*(4-i));
75 ✗ coeff[i] = copy[8-i]*s + copy[i]*c;
76 ✗ coeff[8-i] = copy[8-i]*c - copy[i]*s;
77 }
78 ✗ }
79
80 ✗ void SphericalHarmonicL4::Ry(double alpha) {
81 ✗ SphericalHarmonicL4 c(*this);
82 ✗ double sa = sin(alpha), ca = cos(alpha);
83 ✗ double s2a = sin(2.*alpha), c2a = cos(2.*alpha);
84 ✗ double s3a = sin(3.*alpha), c3a = cos(3.*alpha);
85 ✗ double s4a = sin(4.*alpha), c4a = cos(4.*alpha);
86 ✗ coeff[0] = (c3a+ca*7.)*.125*c[0] + (3.*s3a+7.*sa)*sqrt(.0078125)*c[1] + -(c3a-ca)*sqrt(.109375)*c[2] + -(s3a-3.*sa)*sqrt(.0546875)*c[3];
87 ✗ coeff[1] = -(3.*s3a+7.*sa)*sqrt(.0078125)*c[0] + (9.*c3a+7.*ca)*.0625*c[1] + (3.*s3a-sa)*sqrt(.0546875)*c[2] + -(c3a-ca)*sqrt(.24609375)*c[3];
88 ✗ coeff[2] = -(c3a-ca)*sqrt(.109375)*c[0] + -(3.*s3a-sa)*sqrt(.0546875)*c[1] + (7.*c3a+ca)*.125*c[2] + (7.*s3a+3.*sa)*sqrt(.0078125)*c[3];
89 ✗ coeff[3] = (s3a-3.*sa)*sqrt(.0546875)*c[0] + -(c3a-ca)*sqrt(.24609375)*c[1] + -(7.*s3a+3.*sa)*sqrt(.0078125)*c[2] + (7.*c3a+9.*ca)*.0625*c[3];
90 ✗ coeff[4] = (35.*c4a+20.*c2a+9.)*.015625*c[4] + -(7.*s4a+2.*s2a)*sqrt(.009765625)*c[5] + -(7.*c4a-4.*c2a-3.)*sqrt(.0048828125)*c[6] + (s4a-2.*s2a)*sqrt(.068359375)*c[7] + (c4a-4.*c2a+3.)*sqrt(.008544921875)*c[8];
91 ✗ coeff[5] = (7.*s4a+2.*s2a)*sqrt(.009765625)*c[4] + (7.*c4a+c2a)*.125*c[5] + -(7.*s4a-2.*s2a)*sqrt(.0078125)*c[6] + -(c4a-c2a)*sqrt(.109375)*c[7] + (s4a-2.*s2a)*sqrt(.013671875)*c[8];
92 ✗ coeff[6] = -(7.*c4a-4.*c2a-3.)*sqrt(.0048828125)*c[4] + (7.*s4a-2.*s2a)*sqrt(.0078125)*c[5] + (7.*c4a+4.*c2a+5.)*.0625*c[6] + -(s4a+2.*s2a)*sqrt(.0546875)*c[7] + -(c4a+4.*c2a-5.)*sqrt(.0068359375)*c[8];
93 ✗ coeff[7] = -(s4a-2.*s2a)*sqrt(.068359375)*c[4] + -(c4a-c2a)*sqrt(.109375)*c[5] + (s4a+2.*s2a)*sqrt(.0546875)*c[6] + (c4a+7.*c2a)*.125*c[7] + -(s4a+14.*s2a)*sqrt(.001953125)*c[8];
94 ✗ coeff[8] = (c4a-4.*c2a+3.)*sqrt(.008544921875)*c[4] + -(s4a-2.*s2a)*sqrt(.013671875)*c[5] + -(c4a+4.*c2a-5.)*sqrt(.0068359375)*c[6] + (s4a+14.*s2a)*sqrt(.001953125)*c[7] + (c4a + 28.*c2a+35.)*.015625*c[8];
95 ✗ }
96
97 ✗ void SphericalHarmonicL4::Rx(double alpha) {
98 ✗ SphericalHarmonicL4 c(*this);
99 ✗ double sa = sin(alpha), ca = cos(alpha);
100 ✗ double s2a = sin(2.*alpha), c2a = cos(2.*alpha);
101 ✗ double s3a = sin(3.*alpha), c3a = cos(3.*alpha);
102 ✗ double s4a = sin(4.*alpha), c4a = cos(4.*alpha);
103 ✗ coeff[0] = (c3a+ca*7.)*.125*c[0] + (c3a-ca)*sqrt(.109375)*c[2] + -(s3a-3.*sa)*sqrt(.0546875)*c[5] + -(3.*s3a+7.*sa)*sqrt(.0078125)*c[7];
104 ✗ coeff[1] = (c4a+7.*c2a)*.125*c[1] + (c4a-c2a)*sqrt(.109375)*c[3] + - (s4a-2.*s2a)*sqrt(.068359375)*c[4] + -(s4a+2.*s2a)*sqrt(.0546875)*c[6] + -(s4a+14.*s2a)*sqrt(.001953125)*c[8];
105 ✗ coeff[2] = (c3a-ca)*sqrt(.109375)*c[0] + (7.*c3a+ca)*.125*c[2] + -(7.*s3a+3.*sa)*sqrt(.0078125)*c[5] + -(3.*s3a-sa)*sqrt(.0546875)*c[7];
106 ✗ coeff[3] = (c4a-c2a)*sqrt(.109375)*c[1] + (7.*c4a+c2a)*.125*c[3] + -(7.*s4a+2.*s2a)*sqrt(.009765625)*c[4] + -(7.*s4a-2.*s2a)*sqrt(.0078125)*c[6] + -(s4a-2.*s2a)*sqrt(.013671875)*c[8];
107 ✗ coeff[4] = (s4a-2.*s2a)*sqrt(.068359375)*c[1] + (7.*s4a+2.*s2a)*sqrt(.009765625)*c[3] + (35.*c4a+20.*c2a+9.)*.015625*c[4] + (7.*c4a-4.*c2a-3.)*sqrt(.0048828125)*c[6] + (c4a-4.*c2a+3.)*sqrt(.008544921875)*c[8];
108 ✗ coeff[5] = (s3a-3.*sa)*sqrt(.0546875)*c[0] + (7.*s3a+3*sa)*sqrt(.0078125)*c[2] + (7.*c3a+9.*ca)*.0625*c[5] + (c3a-ca)*sqrt(.24609375)*c[7];
109 ✗ coeff[6] = (s4a+2.*s2a)*sqrt(.0546875)*c[1] + (7.*s4a-2.*s2a)*sqrt(.0078125)*c[3] + (7.*c4a-4.*c2a-3.)*sqrt(.0048828125)*c[4] + (7.*c4a+4.*c2a+5.)*.0625*c[6] + (c4a+4.*c2a-5.)*sqrt(.0068359375)*c[8];
110 ✗ coeff[7] = (3.*s3a+7.*sa)*sqrt(.0078125)*c[0] + (3*s3a-sa)*sqrt(.0546875)*c[2] + (c3a-ca)*sqrt(.24609375)*c[5] + (9.*c3a+7.*ca)*.0625*c[7];
111 ✗ coeff[8] = (s4a+14.*s2a)*sqrt(.001953125)*c[1] + (s4a-2.*s2a)*sqrt(.013671875)*c[3] + (c4a-4.*c2a+3.)*sqrt(.008544921875)*c[4] + (c4a+4.*c2a-5.)*sqrt(.0068359375)*c[6] + (c4a+28.*c2a+35.)*.015625*c[8];
112 ✗ }
113
114 ✗ SphericalHarmonicL4 SphericalHarmonicL4::Ex() const {
115 ✗ return SphericalHarmonicL4(-sqrt(2.)*coeff[7], -sqrt(2.)*coeff[8]-sqrt(3.5)*coeff[6], -sqrt(3.5)*coeff[7]-sqrt(4.5)*coeff[5], -sqrt(4.5)*coeff[6]-sqrt(10.)*coeff[4], sqrt(10.)*coeff[3], sqrt(4.5)*coeff[2], sqrt(3.5)*coeff[1]+sqrt(4.5)*coeff[3], sqrt(2.)*coeff[0]+sqrt(3.5)*coeff[2], sqrt(2.)*coeff[1] );
116 }
117
118 ✗ SphericalHarmonicL4 SphericalHarmonicL4::Ey() const {
119 ✗ return SphericalHarmonicL4(sqrt(2.)*coeff[1], -sqrt(2.)*coeff[0]+sqrt(3.5)*coeff[2], -sqrt(3.5)*coeff[1]+sqrt(4.5)*coeff[3], -sqrt(4.5)*coeff[2], -sqrt(10.)*coeff[5], -sqrt(4.5)*coeff[6] + sqrt(10.)*coeff[4], -sqrt(3.5)*coeff[7]+sqrt(4.5)*coeff[5], -sqrt(2.)*coeff[8]+sqrt(3.5)*coeff[6], sqrt(2.)*coeff[7]);
120 }
121
122 ✗ SphericalHarmonicL4 SphericalHarmonicL4::Ez() const {
123 ✗ return SphericalHarmonicL4(4*coeff[8], 3*coeff[7], 2*coeff[6], coeff[5], 0, -coeff[3], -2*coeff[2], -3*coeff[1], -4*coeff[0]);
124 }
125
126
127 ✗ mat3 SphericalHarmonicL4::project_mat3(double grad_threshold, double dot_threshold, vec3* euler_prev) {
128 SphericalHarmonicL4 init_harmonics[5] = { rest_frame(),rest_frame(),rest_frame(),rest_frame(),rest_frame() };
129 vec3 init_rot[5] = { vec3(0, 0, 0),vec3(M_PI / 4., 0, 0),vec3(0, M_PI / 4., 0),vec3(0, 0, M_PI / 4.),vec3(M_PI / 4., 0, M_PI / 4.)};
130 ✗ FOR(i,5) init_harmonics[i].euler_rot(init_rot[i]);
131
132 mat3 W;
133 SphericalHarmonicL4 v;
134 double dot = -1.;
135
136 ✗ SphericalHarmonicL4 query = *this;
137 ✗ query = query / query.norm();
138
139
140 ✗ if (euler_prev) {
141 ✗ SphericalHarmonicL4 prev_seed = init_harmonics[0];
142
143 ✗ prev_seed.euler_rot(*euler_prev);
144 double tdot = prev_seed*query;
145 ✗ if (tdot>dot) {
146 dot = tdot;
147
148 ✗ W = euler_to_mat3(*euler_prev);
149 ✗ v = prev_seed;
150 }
151 }
152
153
154 ✗ FOR(i,5) {
155 ✗ double tdot = init_harmonics[i] * query;
156 ✗ if (tdot>dot) {
157 dot = tdot;
158
159 ✗ W = euler_to_mat3(init_rot[i]);
160 ✗ v = init_harmonics[i];
161 }
162 }
163
164 int cnt = 0;
165 double olddot = dot;
166 ✗ while (cnt < 10000) {
167 ✗ vec3 grad = vec3(query*v.Ex(), query*v.Ey(), query*v.Ez());
168 ✗ if (length(grad) < grad_threshold) break;
169 grad /= 8.; // constante au pif trouvée de manière experimentale ; en dessous de ça pb de convergence, au-dessus plus lent
170 ✗ v.Rx(grad[0]); v.Ry(grad[1]); v.Rz(grad[2]);
171
172 ✗ W = rotz(grad[2])*roty(grad[1])*rotx(grad[0])*W;
173 ✗ cnt++;
174 dot = v*query;
175
176 ✗ if (dot - olddot < dot_threshold) break;
177 olddot = dot;
178 }
179 ✗ if (cnt == 10000) GEO::Logger::out("HexDom") << "[error] SH projection infinite loop protection" << std::endl;
180
181 ✗ return W;
182 }
183
184 }
185