1 /*
2 * Stellarium
3 * Copyright (C) Fabien Chereau
4 * Author 2006 Johannes Gajdosik
5 *
6 * This program is free software; you can redistribute it and/or
7 * modify it under the terms of the GNU General Public License
8 * as published by the Free Software Foundation; either version 2
9 * of the License, or (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
19 */
20
21 #include "SphericMirrorCalculator.hpp"
22 #include <QSettings>
23 #include <QDebug>
24 #include "StelUtils.hpp"
25
SphericMirrorCalculator(const QSettings & conf)26 SphericMirrorCalculator::SphericMirrorCalculator(const QSettings& conf)
27 {
28 const Vec3f mirror_position(
29 conf.value("spheric_mirror/mirror_position_x",0.0).toFloat(),
30 conf.value("spheric_mirror/mirror_position_y",2.0).toFloat(),
31 conf.value("spheric_mirror/mirror_position_z",0.0).toFloat());
32 const float mirror_radius(conf.value("spheric_mirror/mirror_radius",0.25).toFloat());
33 DomeCenter = mirror_position * (-1.0f/mirror_radius);
34 const float dome_radius(conf.value("spheric_mirror/dome_radius",2.5).toFloat());
35 DomeRadius = dome_radius / mirror_radius;
36 const Vec3f projector_position(
37 conf.value("spheric_mirror/projector_position_x",0.0).toFloat(),
38 conf.value("spheric_mirror/projector_position_y",1.0).toFloat(),
39 conf.value("spheric_mirror/projector_position_z",-0.2).toFloat());
40 P = (projector_position - mirror_position) * (1.0f/mirror_radius);
41 PP = P.lengthSquared();
42 lP = std::sqrt(PP);
43 p = P * (1.0f/lP);
44 float image_distance_div_height = conf.value("spheric_mirror/image_distance_div_height",-1e38f).toFloat();
45 if (image_distance_div_height <= -1e38f)
46 {
47 const float scaling_factor = conf.value("spheric_mirror/scaling_factor", 0.8).toFloat();
48 image_distance_div_height = std::sqrt(PP-1.0f) * scaling_factor;
49 qDebug() << "INFO: spheric_mirror:scaling_factor is deprecated and may be removed in future versions.";
50 qDebug() << " In order to keep your setup unchanged, please use spheric_mirror:image_distance_div_height = "
51 << image_distance_div_height << " instead";
52 }
53 horzZoomFactor = conf.value("spheric_mirror/flip_horz",true).toBool() ? (-image_distance_div_height) : image_distance_div_height;
54 vertZoomFactor = conf.value("spheric_mirror/flip_vert",false).toBool() ? (-image_distance_div_height) : image_distance_div_height;
55
56 const float alpha = conf.value("spheric_mirror/projector_alpha",0.0f).toFloat() * (M_PI_180f);
57 const float phi = conf.value("spheric_mirror/projector_phi",0.0f).toFloat() * (M_PI_180f);
58 float delta = conf.value("spheric_mirror/projector_delta",-1e38f).toFloat();
59 if (delta <= -1e38f)
60 {
61 float x,y;
62 // before calling transform() horzZoomFactor,vertZoomFactor,
63 // alphaDeltaPhi must already be initialized
64 initRotMatrix(0.0,0.0,0.0);
65 transform(Vec3f(0,0,1),x,y);
66 const float zenith_y(conf.value("spheric_mirror/zenith_y",0.125).toFloat());
67 delta = -atan(y/image_distance_div_height) + atan(zenith_y/image_distance_div_height);
68 qDebug() << "INFO: spheric_mirror:zenith_y is deprecated and may be removed in future versions.";
69 qDebug() << " In order to keep your setup unchanged, please use spheric_mirror:projector_delta = "
70 << (delta*M_180_PIf) << " instead";
71 }
72 else
73 delta *= M_PI_180f;
74
75 initRotMatrix(alpha,delta,phi);
76 }
77
initRotMatrix(float alpha,float delta,float phi)78 void SphericMirrorCalculator::initRotMatrix(float alpha, float delta, float phi)
79 {
80 const float ca = cos(alpha);
81 const float sa = sin(alpha);
82 const float cd = cos(delta);
83 const float sd = sin(delta);
84 const float cp = cos(phi);
85 const float sp = sin(phi);
86
87 alphaDeltaPhi[0] = ca*cp - sa*sd*sp;
88 alphaDeltaPhi[1] = - sa*cp - ca*sd*sp;
89 alphaDeltaPhi[2] = - cd*sp;
90
91 alphaDeltaPhi[3] = sa*cd;
92 alphaDeltaPhi[4] = ca*cd;
93 alphaDeltaPhi[5] = - sd;
94
95 alphaDeltaPhi[6] = - ca*sp - sa*sd*cp;
96 alphaDeltaPhi[7] = sa*sp - ca*sd*cp;
97 alphaDeltaPhi[8] = - cd*cp;
98 /*
99 // check if alphaDeltaPhi is an orthogonal matrix:
100 for (int i=0;i<3;i++) {
101 for (int j=0;j<3;j++) {
102 float prod0 = 0;
103 float prod1 = 0;
104 for (int k=0;k<3;k++) {
105 prod0 += alphaDeltaPhi[3*i+k]*alphaDeltaPhi[3*j+k];
106 prod1 += alphaDeltaPhi[i+3*k]*alphaDeltaPhi[j+3*k];
107 }
108 if (i==j) {
109 prod0 -= 1.0;
110 prod1 -= 1.0;
111 }
112 if (fabs(prod0)>1e-10) {
113 qDebug() << "i: " << i << ", j: " << j
114 << ", prod0: " << prod0 << ", prod1: " << prod1;
115 Q_ASSERT(0);
116 }
117 if (fabs(prod1)>1e-10) {
118 qDebug << "i: " << i << ", j: " << j
119 << ", prod0: " << prod0 << ", prod1: " << prod1;
120 Q_ASSERT(0);
121 }
122 }
123 }
124 */
125 }
126
127
transform(const Vec3f & v,float & xb,float & yb) const128 bool SphericMirrorCalculator::transform(const Vec3f &v, float &xb,float &yb) const
129 {
130 const Vec3f S = DomeCenter + (v * (DomeRadius/v.length()));
131 const Vec3f SmP = S - P;
132 const float P_SmP = P.dot(SmP);
133 const bool rval = ( (PP-1.0f)*SmP.dot(SmP) > P_SmP*P_SmP );
134
135 const float lS = S.length();
136 const Vec3f s(S/lS);
137 float t_min = 0;
138 float t_max = 1;
139 Vec3f Q;
140 // more iterations would be more accurate,
141 // but I keep this number of iterations for exact zenith_y compatibility:
142 for (int i=0;i<10;i++)
143 {
144 const float t = 0.5f * (t_min+t_max);
145 Q = p*t + s*(1.0f-t);
146 Q.normalize();
147 Vec3f Qp = P-Q;
148 Qp.normalize();
149 Vec3f Qs = S-Q;
150 Qs.normalize();
151 if ( (Qp-Qs).dot(Q) > 0.0f )
152 t_max = t;
153 else
154 t_min = t;
155 }
156 Vec3f x = Q-P;
157
158 // rotate
159 const float zb = alphaDeltaPhi[1]*x[0] + alphaDeltaPhi[4]*x[1] + alphaDeltaPhi[7]*x[2];
160 xb = (horzZoomFactor/zb) * (alphaDeltaPhi[0]*x[0] + alphaDeltaPhi[3]*x[1] + alphaDeltaPhi[6]*x[2]);
161 yb = (vertZoomFactor/zb) * (alphaDeltaPhi[2]*x[0] + alphaDeltaPhi[5]*x[1] + alphaDeltaPhi[8]*x[2]);
162
163 return rval;
164 }
165
166
retransform(float x,float y,Vec3f & v) const167 bool SphericMirrorCalculator::retransform(float x,float y, Vec3f &v) const
168 {
169 x /= horzZoomFactor;
170 y /= vertZoomFactor;
171 v[0] = alphaDeltaPhi[0]*x + alphaDeltaPhi[1] + alphaDeltaPhi[2]*y;
172 v[1] = alphaDeltaPhi[3]*x + alphaDeltaPhi[4] + alphaDeltaPhi[5]*y;
173 v[2] = alphaDeltaPhi[6]*x + alphaDeltaPhi[7] + alphaDeltaPhi[8]*y;
174 const float vv = v.dot(v);
175 const float Pv = P.dot(v);
176 const float discr = Pv*Pv-(P.dot(P)-1.0f)*vv;
177 if (discr < 0)
178 return false;
179
180 const Vec3f Q = P + v*((-Pv-std::sqrt(discr))/vv);
181 const Vec3f w = v - Q*(2*v.dot(Q));
182 const Vec3f MQ = Q - DomeCenter;
183 float f = -MQ.dot(w);
184 f += std::sqrt(f*f - (MQ.dot(MQ)-DomeRadius*DomeRadius)*vv);
185 const Vec3f S = Q + w*(f/vv);
186 v = S - DomeCenter;
187 v *= (1.0f/DomeRadius);
188 return true;
189 }
190
retransform(float x,float y,Vec3f & v,Vec3f & vX,Vec3f & vY) const191 bool SphericMirrorCalculator::retransform(float x,float y, Vec3f &v, Vec3f &vX,Vec3f &vY) const
192 {
193 x /= horzZoomFactor;
194 const float dx = 1.0f/horzZoomFactor;
195 y /= vertZoomFactor;
196 const float dy = 1.0f/vertZoomFactor;
197
198 v[0] = alphaDeltaPhi[0]*x + alphaDeltaPhi[1] + alphaDeltaPhi[2]*y;
199 v[1] = alphaDeltaPhi[3]*x + alphaDeltaPhi[4] + alphaDeltaPhi[5]*y;
200 v[2] = alphaDeltaPhi[6]*x + alphaDeltaPhi[7] + alphaDeltaPhi[8]*y;
201
202 vX[0] = alphaDeltaPhi[0]*dx;
203 vX[1] = alphaDeltaPhi[3]*dx;
204 vX[2] = alphaDeltaPhi[6]*dx;
205
206 vY[0] = alphaDeltaPhi[2]*dy;
207 vY[1] = alphaDeltaPhi[5]*dy;
208 vY[2] = alphaDeltaPhi[8]*dy;
209
210 const float vv = v.dot(v);
211 const float vvX = 2.0f*v.dot(vX);
212 const float vvY = 2.0f*v.dot(vY);
213
214 const float Pv = P.dot(v);
215 const float PvX = P.dot(vX);
216 const float PvY = P.dot(vY);
217
218 const float discr = Pv*Pv-(P.dot(P)-1.0f)*vv;
219 const float discr_x = 2.0f*Pv*PvX-(P.dot(P)-1.0f)*vvX;
220 const float discr_y = 2.0f*Pv*PvY-(P.dot(P)-1.0f)*vvY;
221
222 if (discr < 0)
223 return false;
224
225 const Vec3f Q = P + v*((-Pv-std::sqrt(discr))/vv);
226 const Vec3f Q_x = vX*((-Pv-std::sqrt(discr))/vv) + v*( (vv*(-PvX-0.5f*discr_x/std::sqrt(discr)) - vvX*(-Pv-std::sqrt(discr))) /(vv*vv));
227 const Vec3f Q_y = vY*((-Pv-std::sqrt(discr))/vv) + v*( (vv*(-PvY-0.5f*discr_y/std::sqrt(discr)) - vvY*(-Pv-std::sqrt(discr))) /(vv*vv));
228
229 const Vec3f w = v - Q*(2*v.dot(Q));
230 const Vec3f w_x = vX - Q_x*(2*v.dot(Q)) - Q*(2*(vX.dot(Q)+v.dot(Q_x)));
231 const Vec3f w_y = vY - Q_y*(2*v.dot(Q)) - Q*(2*(vY.dot(Q)+v.dot(Q_y)));
232
233 const Vec3f MQ = Q - DomeCenter;
234 // MQ_x = Q_x
235 // MQ_y = Q_y
236
237 float f = -MQ.dot(w);
238 float f_x = -Q_x.dot(w)-MQ.dot(w_x);
239 float f_y = -Q_y.dot(w)-MQ.dot(w_y);
240
241 float f1 = f + std::sqrt(f*f - (MQ.dot(MQ)-DomeRadius*DomeRadius)*vv);
242 float f1_x = f_x + 0.5f*(2*f*f_x - (MQ.dot(MQ)-DomeRadius*DomeRadius)*vvX - 2.f*MQ.dot(Q_x)*vv ) / std::sqrt(f*f - (MQ.dot(MQ)-DomeRadius*DomeRadius)*vv);
243 float f1_y = f_y + 0.5f*(2*f*f_y - (MQ.dot(MQ)-DomeRadius*DomeRadius)*vvY - 2.f*MQ.dot(Q_y)*vv ) / std::sqrt(f*f - (MQ.dot(MQ)-DomeRadius*DomeRadius)*vv);
244
245 const Vec3f S = Q + w*(f1/vv);
246 const Vec3f S_x = Q_x + w*((vv*f1_x-vvX*f1)/(vv*vv)) + w_x*(f1/vv);
247 const Vec3f S_y = Q_y + w*((vv*f1_y-vvY*f1)/(vv*vv)) + w_y*(f1/vv);
248
249 v = S - DomeCenter;
250 vX = S_x;
251 vY = S_y;
252
253 v *= (1.0f/DomeRadius);
254 vX *= (1.0f/DomeRadius);
255 vY *= (1.0f/DomeRadius);
256 return true;
257 }
258