15 template <
class InIteratorT,
class OutIteratorT>
17 SpinImage(
const Vec3f& axisPos,
23 for (; begin != end; ++begin)
25 Vec3f s = *begin - axisPos;
27 spin[1] =
s.dot(axisDir);
28 spin[0] = (
s - spin[1] * axisDir).length();
33 template <
class InIteratorT>
43 a =
std::sqrt(std::pow(i[1][0] - i[0][0], 2) + std::pow(i[1][1] - i[0][1], 2));
44 b =
std::sqrt(std::pow(i[2][0] - i[1][0], 2) + std::pow(i[2][1] - i[1][1], 2));
45 c =
std::sqrt(std::pow(i[0][0] - i[2][0], 2) + std::pow(i[0][1] - i[2][1], 2));
47 bot = (
a + b +
c) * (-
a + b +
c) * (
a - b +
c) * (
a + b -
c);
56 top1 = (i[1][1] - i[0][1]) *
c *
c - (i[2][1] - i[0][1]) *
a *
a;
57 top2 = (i[1][0] - i[0][0]) *
c *
c - (i[2][0] - i[0][0]) *
a *
a;
58 bot = (i[1][1] - i[0][1]) * (i[2][0] - i[0][0]) - (i[2][1] - i[0][1]) * (i[1][0] - i[0][0]);
60 (*center)[0] = i[0][0] + 0.5f * top1 / bot;
61 (*center)[1] = i[0][1] - 0.5f * top2 / bot;
69 if (samples.
size() < 8)
83 size_t k = samples.
size() >> 1;
86 for (
size_t i = 0; i < samples.
size(); ++i)
89 double a01 = n0xn1 * dsamples[k + 2];
90 double b01 = n0xn1 * dsamples[k + 3];
91 double a0 = ((dsamples[2] - dsamples[1]) % dsamples[k]) * dsamples[k + 2];
92 double a1 = ((dsamples[0] - dsamples[2]) % dsamples[k + 1]) * dsamples[k + 2];
93 double a = ((dsamples[0] - dsamples[2]) % (dsamples[1] - dsamples[0])) * dsamples[k + 2];
94 double b0 = ((dsamples[3] - dsamples[1]) % dsamples[k]) * dsamples[k + 3];
95 double b1 = ((dsamples[0] - dsamples[3]) % dsamples[k + 1]) * dsamples[k + 3];
96 double b = ((dsamples[0] - dsamples[3]) % (dsamples[1] - dsamples[0])) * dsamples[k + 3];
97 double cc = b01 / a01;
98 double ccc = b0 - a0 * cc;
99 double c = -(b1 - a1 * cc) / ccc;
100 double d = (-b +
a * cc) / ccc;
101 double p = (a0 *
c + a1 + a01 * d) / (2 * a01 *
c);
102 double q = (
a + a0 * d) / (a01 *
c);
103 double rt = p * p -
q;
114 double s1 =
c * t1 + d;
115 double s2 =
c * t2 + d;
117 Vec3f pos1 = samples[0] + s1 * samples[k];
118 Vec3f normal1 = pos1 - (samples[1] + t1 * samples[k + 1]);
120 Vec3f pos2 = samples[0] + s2 * samples[k];
121 Vec3f normal2 = pos2 - (samples[1] + t2 * samples[k + 1]);
126 SpinImage(pos1, normal1, samples.begin(), samples.begin() + k, std::back_inserter(spin1));
127 SpinImage(pos2, normal2, samples.begin(), samples.begin() + k, std::back_inserter(spin2));
128 float minorRadius1, majorRadius1, minorRadius2, majorRadius2,
129 distSum1 = std::numeric_limits<float>::infinity(),
130 distSum2 = std::numeric_limits<float>::infinity(), tmp;
132 if (CircleFrom3Points(spin1.
begin(), &minorRadius1, &minorCenter1))
134 majorRadius1 = minorCenter1[0];
137 for (
size_t i = 3; i < spin1.
size(); ++i)
138 distSum1 += (tmp = ((spin1[i] - minorCenter1).Length() - minorRadius1)) * tmp;
140 if (CircleFrom3Points(spin2.
begin(), &minorRadius2, &minorCenter2))
142 majorRadius2 = minorCenter2[0];
145 for (
size_t i = 3; i < spin2.
size(); ++i)
146 distSum2 += (tmp = ((spin2[i] - minorCenter2).Length() - minorRadius2)) * tmp;
148 if (distSum1 != std::numeric_limits<float>::infinity() && distSum1 < distSum2)
151 m_rminor = minorRadius1;
152 m_rmajor = majorRadius1;
153 m_center = pos1 + minorCenter1[1] * m_normal;
155 else if (distSum2 != std::numeric_limits<float>::infinity())
158 m_rminor = minorRadius2;
159 m_rmajor = majorRadius2;
160 m_center = pos2 + minorCenter2[1] * m_normal;
166 m_appleShaped = m_rmajor < m_rminor;
167 ComputeAppleParams();
174 if (samples.
size() < 8)
188 size_t k = samples.
size() >> 1;
191 for (
size_t i = 0; i < samples.
size(); ++i)
195 Vec3f pos1, normal1, pos2, normal2;
196 for (
size_t w = 0; w < k ; ++w)
198 for (
size_t x = 0 ;
x < k ; ++
x)
201 for (
size_t y = 0 ; y < k ; ++y)
203 for (
size_t z = 0 ; z < k; ++z)
205 double a01 = n0xn1 * dsamples[k + y];
206 double b01 = n0xn1 * dsamples[k + z];
212 double a0 = ((dsamples[y] - dsamples[
x]) % dsamples[k + w]) * dsamples[k + y];
213 double a1 = ((dsamples[w] - dsamples[y]) % dsamples[k +
x]) * dsamples[k + y];
214 double a = ((dsamples[w] - dsamples[y]) % (dsamples[
x] - dsamples[w])) *
216 double b0 = ((dsamples[z] - dsamples[
x]) % dsamples[k + w]) * dsamples[k + z];
217 double b1 = ((dsamples[w] - dsamples[z]) % dsamples[k +
x]) * dsamples[k + z];
218 double b = ((dsamples[w] - dsamples[z]) % (dsamples[
x] - dsamples[w])) *
220 double cc = b01 / a01;
221 double ccc = b0 - a0 * cc;
222 double c = -(b1 - a1 * cc) / ccc;
223 double d = (-b +
a * cc) / ccc;
224 double p = (a0 *
c + a1 + a01 * d) / (2 * a01 *
c);
225 double q = (
a + a0 * d) / (a01 *
c);
226 double rt = p * p -
q;
237 double s1 =
c * t1 + d;
238 double s2 =
c * t2 + d;
239 pos1 = samples[w] + s1 * samples[k + w];
240 normal1 = pos1 - (samples[
x] + t1 * samples[k +
x]);
242 pos2 = samples[w] + s2 * samples[k + w];
243 normal2 = pos2 - (samples[
x] + t2 * samples[k +
x]);
255 SpinImage(pos1, normal1, samples.
begin(), samples.
begin() + k, std::back_inserter(spin1));
256 SpinImage(pos2, normal2, samples.
begin(), samples.
begin() + k, std::back_inserter(spin2));
257 float minorRadius1, majorRadius1, minorRadius2, majorRadius2,
258 distSum1 = std::numeric_limits<float>::infinity(),
259 distSum2 = std::numeric_limits<float>::infinity(), tmp;
261 if (CircleFrom3Points(spin1.
begin(), &minorRadius1, &minorCenter1))
263 majorRadius1 = minorCenter1[0];
266 for (
size_t i = 3; i < spin1.
size(); ++i)
267 distSum1 += (tmp = ((spin1[i] - minorCenter1).Length() - minorRadius1)) * tmp;
269 if (CircleFrom3Points(spin2.
begin(), &minorRadius2, &minorCenter2))
271 majorRadius2 = minorCenter2[0];
274 for (
size_t i = 3; i < spin2.
size(); ++i)
275 distSum2 += (tmp = ((spin2[i] - minorCenter2).Length() - minorRadius2)) * tmp;
277 if (distSum1 != std::numeric_limits<float>::infinity() && distSum1 < distSum2)
280 m_rminor = minorRadius1;
281 m_rmajor = majorRadius1;
282 m_center = pos1 + minorCenter1[1] * m_normal;
284 else if (distSum2 != std::numeric_limits<float>::infinity())
287 m_rminor = minorRadius2;
288 m_rmajor = majorRadius2;
289 m_center = pos2 + minorCenter2[1] * m_normal;
295 m_appleShaped = m_rmajor < m_rminor;
296 ComputeAppleParams();
305 i->read((
char*)&m_normal,
sizeof(m_center));
306 i->read((
char*)&m_center,
sizeof(m_center));
307 i->read((
char*)&m_rminor,
sizeof(m_rminor));
308 i->read((
char*)&m_rmajor,
sizeof(m_rmajor));
312 for (
size_t j = 0; j < 3; ++j)
316 for (
size_t j = 0; j < 3; ++j)
323 m_appleShaped = m_rmajor < m_rminor;
324 ComputeAppleParams();
332 fread(&m_normal,
sizeof(m_normal), 1, i);
333 fread(&m_center,
sizeof(m_center), 1, i);
334 fread(&m_rminor,
sizeof(m_rminor), 1, i);
335 fread(&m_rmajor,
sizeof(m_rmajor), 1, i);
336 fread(&rot,
sizeof(rot), 1, i);
337 m_appleShaped = m_rmajor < m_rminor;
338 ComputeAppleParams();
344 for (
int i = 0; i < 3; i++)
346 m_normal[i] = array[i];
347 m_center[i] = array[i + 3];
351 m_appleShaped = m_rmajor < m_rminor;
352 ComputeAppleParams();
360 m_center += translate;
367 s[0] =
x[0] - param[0];
368 s[1] =
x[1] - param[1];
369 s[2] =
x[2] - param[2];
370 float g =
s[0] * param[3] +
s[1] * param[4] +
s[2] * param[5];
371 float f = param[5] *
s[1] - param[4] *
s[2];
373 float v = param[3] *
s[2] - param[5] *
s[0];
375 v = param[4] *
s[0] - param[3] *
s[1];
379 return std::sqrt(g * g + ((tmp = (f - param[6])) * tmp)) - param[7];
386 s[0] =
x[0] - param[0];
387 s[1] =
x[1] - param[1];
388 s[2] =
x[2] - param[2];
389 float g =
s[0] * param[3] +
s[1] * param[4] +
s[2] * param[5];
390 float f = param[5] *
s[1] - param[4] *
s[2];
392 float v = param[3] *
s[2] - param[5] *
s[0];
394 v = param[4] *
s[0] - param[3] *
s[1];
401 dg[3] =
s[0] - param[3] * g;
402 dg[4] =
s[1] - param[4] * g;
403 dg[5] =
s[2] - param[5] * g;
405 df[0] = (param[3] * g -
s[0]) / f;
406 df[1] = (param[4] * g -
s[1]) / f;
407 df[2] = (param[5] * g -
s[2]) / f;
412 float d =
std::sqrt(g * g + ((tmp = (f - param[6])) * tmp)) - param[7];
413 float dr = d + param[7];
414 float fr = f - param[6];
415 for (
unsigned int i = 0; i < 6; ++i)
417 gradient[i] = (g * dg[i] + fr * df[i]) / dr;
419 gradient[6] = -fr * dr;
426 float l =
std::sqrt(param[3] * param[3] + param[4] * param[4] + param[5] * param[5]);
427 for (
unsigned int i = 3; i < 6; ++i)
433 template <
class WeightT>
444 template <
class IteratorT>
453 int size = end - begin;
454 #pragma omp parallel for schedule(static) reduction(+ : chi)
455 for (
int idx = 0; idx < size; ++idx)
458 s[0] = begin[idx][0] - params[0];
459 s[1] = begin[idx][1] - params[1];
460 s[2] = begin[idx][2] - params[2];
461 ScalarType g =
s[0] * params[3] +
s[1] * params[4] +
s[2] * params[5];
466 v = params[4] *
s[0] - params[3] *
s[1];
471 chi += (
values[idx] = WeightT::Weigh(
472 std::sqrt(g * g + ((tmp = (f - params[6])) * tmp)) - params[7])) *
479 template <
class IteratorT>
488 int size = end - begin;
489 #pragma omp parallel for schedule(static)
490 for (
int idx = 0; idx < size; ++idx)
493 s[0] = begin[idx][0] - params[0];
494 s[1] = begin[idx][1] - params[1];
495 s[2] = begin[idx][2] - params[2];
496 ScalarType g =
s[0] * params[3] +
s[1] * params[4] +
s[2] * params[5];
501 dg[3] =
s[0] - params[3] * g;
502 dg[4] =
s[1] - params[4] * g;
503 dg[5] =
s[2] - params[5] * g;
505 df[0] = (params[3] * g -
s[0]) / temp[idx];
506 df[1] = (params[4] * g -
s[1]) / temp[idx];
507 df[2] = (params[5] * g -
s[2]) / temp[idx];
515 for (
unsigned int j = 0; j < 6; ++j)
517 matrix[idx *
NumParams + j] = (g * dg[j] + fr * df[j]) / dr;
521 WeightT::template DerivWeigh<NumParams>(d, matrix + idx *
NumParams);
529 std::sqrt(params[3] * params[3] + params[4] * params[4] + params[5] * params[5]);
530 for (
unsigned int i = 3; i < 6; ++i)
543 for (
size_t i = 0; i < 3; ++i)
545 param[i] = m_center[i];
547 for (
size_t i = 0; i < 3; ++i)
549 param[i + 3] = m_normal[i];
561 for (
size_t i = 0; i < 3; ++i)
563 m_center[i] = param[i];
565 for (
size_t i = 0; i < 3; ++i)
567 m_normal[i] = param[i + 3];
571 m_appleShaped = m_rmajor < m_rminor;
572 ComputeAppleParams();
581 o->write((
const char*)&m_normal,
sizeof(m_normal));
582 o->write((
const char*)&m_center,
sizeof(m_center));
583 o->write((
const char*)&m_rminor,
sizeof(m_rminor));
584 o->write((
const char*)&m_rmajor,
sizeof(m_rmajor));
588 (*o) << m_normal[0] <<
" " << m_normal[1] <<
" " << m_normal[2] <<
" " << m_center[0] <<
" "
589 << m_center[1] <<
" " << m_center[2] <<
" " << m_rminor <<
" " << m_rmajor <<
" ";
602 fwrite(&m_normal,
sizeof(m_normal), 1, o);
603 fwrite(&m_center,
sizeof(m_center), 1, o);
604 fwrite(&m_rminor,
sizeof(m_rminor), 1, o);
605 fwrite(&m_rmajor,
sizeof(m_rmajor), 1, o);
617 for (
int i = 0; i < 3; i++)
619 array[i] = m_normal[i];
620 array[i + 3] = m_center[i];
627 Torus::ComputeAppleParams()
631 m_cutOffAngle =
M_PI;
635 m_cutOffAngle = std::acos((2 * m_rmajor - m_rminor) / m_rminor) +
M_PI / 2.f;
636 m_appleHeight = std::sin(m_cutOffAngle) * m_rminor;