Skip to content

Commit cc43f7d

Browse files
clairerocheclaireroche
authored
Merge pull request #237 from LIHPC-Computational-Geometry/issue-273
issue-273: fix bug beta law indirect direction
2 parents 1eb2502 + cbde00b commit cc43f7d

3 files changed

Lines changed: 195 additions & 212 deletions

File tree

‎src/Core/Topo/EdgeMeshingPropertyBeta.cpp‎

Lines changed: 160 additions & 155 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,4 @@
11
/*----------------------------------------------------------------------------*/
2-
/*
3-
* \file EdgeMeshingPropertyBeta.cpp
4-
*
5-
* \author Eric Brière de l'Isle
6-
*
7-
* \date 24 janv. 2012
8-
*/
9-
/*----------------------------------------------------------------------------*/
102
#include "Topo/EdgeMeshingPropertyBeta.h"
113
#include "Utils/Common.h"
124
#include "Utils/MgxNumeric.h"
@@ -17,222 +9,235 @@ namespace Mgx3D {
179
/*----------------------------------------------------------------------------*/
1810
namespace Topo {
1911
/*----------------------------------------------------------------------------*/
20-
//#define _DEBUG_beta
21-
/*----------------------------------------------------------------------------*/
2212
EdgeMeshingPropertyBeta::
23-
EdgeMeshingPropertyBeta(int nbBras, double beta, bool isDirect,
24-
bool initWithFirstEdge, double meshingEdgeLength)
25-
: CoEdgeMeshingProperty(nbBras, beta_resserrement, isDirect)
26-
, m_beta(beta)
27-
, m_arm1(meshingEdgeLength)
28-
, m_initWithArm1(initWithFirstEdge)
13+
EdgeMeshingPropertyBeta(const int nbBras, const double beta, const bool isDirect,
14+
const bool initWithFirstEdge, const double meshingEdgeLength)
15+
: CoEdgeMeshingProperty(nbBras, beta_resserrement, isDirect)
16+
, m_beta(beta)
17+
, m_arm1(meshingEdgeLength)
18+
, m_initWithArm1(initWithFirstEdge)
2919
{
30-
if (!m_initWithArm1 && (m_beta<1.00001 || m_beta>1.01)){
31-
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
32-
messErr << "EdgeMeshingPropertyBeta, la valeur de beta doit être comprise entre 1 et 1.01, actuellement : "<<m_beta;
20+
if (!m_initWithArm1 && (m_beta < 1.00001 || m_beta > 1.01))
21+
{
22+
TkUtil::UTF8String messErr(TkUtil::Charset::UTF_8);
23+
messErr << "EdgeMeshingPropertyBeta: ";
24+
messErr << "la valeur de beta doit être comprise entre 1 et 1.01, actuellement : ";
25+
messErr << m_beta;
3326
throw TkUtil::Exception(messErr);
3427
}
35-
if (m_initWithArm1 && m_arm1 < 0){
36-
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
37-
messErr << "EdgeMeshingPropertyGeometric, la longueur du premier bras ne doit pas être négative : "<<m_arm1;
38-
throw TkUtil::Exception(messErr);
28+
if (m_initWithArm1 && m_arm1 < 0)
29+
{
30+
TkUtil::UTF8String messErr(TkUtil::Charset::UTF_8);
31+
messErr << "EdgeMeshingPropertyBeta: ";
32+
messErr << "la longueur du premier bras ne doit pas être négative : ";
33+
messErr << m_arm1;
34+
throw TkUtil::Exception(messErr);
3935
}
4036
if (nbBras < 1)
41-
throw TkUtil::Exception (TkUtil::UTF8String ("EdgeMeshingPropertyBeta, le nombre de bras doit être au moins de 1", TkUtil::Charset::UTF_8));
37+
{
38+
throw TkUtil::Exception(TkUtil::UTF8String(
39+
"EdgeMeshingPropertyBeta, le nombre de bras doit être au moins de 1", TkUtil::Charset::UTF_8));
40+
}
4241
}
42+
4343
/*----------------------------------------------------------------------------*/
44-
EdgeMeshingPropertyBeta::EdgeMeshingPropertyBeta (const CoEdgeMeshingProperty& prop)
45-
: CoEdgeMeshingProperty (prop.getNbEdges(), beta_resserrement, prop.getDirect())
46-
, m_beta(0.)
47-
, m_arm1(0.)
48-
, m_initWithArm1(false)
44+
EdgeMeshingPropertyBeta::EdgeMeshingPropertyBeta(const CoEdgeMeshingProperty &prop)
45+
: CoEdgeMeshingProperty(prop.getNbEdges(), beta_resserrement, prop.getDirect())
46+
, m_beta(0.)
47+
, m_arm1(0.)
48+
, m_initWithArm1(false)
4949
{
50-
const EdgeMeshingPropertyBeta* p = reinterpret_cast<const EdgeMeshingPropertyBeta*>(&prop);
51-
if (0 == p)
52-
{
53-
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
54-
messErr << "EdgeMeshingPropertyBeta, la propriété transmise en argument n'est pas de type EdgeMeshingPropertyBeta";
55-
throw TkUtil::Exception(messErr);
56-
} // if (0 == p)
57-
58-
m_beta = p->m_beta;
59-
m_arm1 = p->m_arm1;
60-
m_initWithArm1 = p->m_initWithArm1;
50+
const EdgeMeshingPropertyBeta *p = reinterpret_cast<const EdgeMeshingPropertyBeta *>(&prop);
51+
if (nullptr == p)
52+
{
53+
TkUtil::UTF8String messErr(TkUtil::Charset::UTF_8);
54+
messErr << "EdgeMeshingPropertyBeta, la propriété transmise en argument n'est pas de type EdgeMeshingPropertyBeta";
55+
throw TkUtil::Exception(messErr);
56+
}
57+
58+
m_beta = p->m_beta;
59+
m_arm1 = p->m_arm1;
60+
m_initWithArm1 = p->m_initWithArm1;
6161
}
62+
6263
/*----------------------------------------------------------------------------*/
63-
bool EdgeMeshingPropertyBeta::operator == (const CoEdgeMeshingProperty& cedp) const
64+
bool EdgeMeshingPropertyBeta::operator ==(const CoEdgeMeshingProperty &cedp) const
6465
{
65-
const EdgeMeshingPropertyBeta* props =
66-
dynamic_cast<const EdgeMeshingPropertyBeta*> (&cedp);
67-
if (0 == props)
68-
return false;
66+
const EdgeMeshingPropertyBeta *props =
67+
dynamic_cast<const EdgeMeshingPropertyBeta *>(&cedp);
68+
if (nullptr == props)
69+
{
70+
return false;
71+
}
6972

70-
if (m_beta != props->m_beta)
71-
return false;
73+
if (m_beta != props->m_beta)
74+
{
75+
return false;
76+
}
7277

73-
if (m_arm1 != props->m_arm1)
74-
return false;
78+
if (m_arm1 != props->m_arm1)
79+
{
80+
return false;
81+
}
7582

76-
if (m_initWithArm1 != props->m_initWithArm1)
77-
return false;
83+
if (m_initWithArm1 != props->m_initWithArm1)
84+
{
85+
return false;
86+
}
7887

79-
return CoEdgeMeshingProperty::operator == (cedp);
88+
return CoEdgeMeshingProperty::operator ==(cedp);
8089
}
90+
8191
/*----------------------------------------------------------------------------*/
8292
void EdgeMeshingPropertyBeta::setNbEdges(const int nb)
8393
{
8494
CoEdgeMeshingProperty::setNbEdges(nb);
8595
}
96+
8697
/*----------------------------------------------------------------------------*/
8798
void EdgeMeshingPropertyBeta::initCoeff(double length)
8899
{
89-
if (m_initWithArm1){
90-
// on calcule le coefficient de resserement
91-
m_beta = computeBeta(m_arm1);
92-
}
93-
initCoeff();
100+
if (m_initWithArm1)
101+
{
102+
// on calcule le coefficient de resserement
103+
m_beta = computeBeta(m_arm1);
104+
}
105+
initCoeff();
94106
}
95107

108+
/*----------------------------------------------------------------------------*/
96109
void EdgeMeshingPropertyBeta::initCoeff()
97110
{
98-
#ifdef _DEBUG_beta
99-
std::cout<<"EdgeMeshingPropertyBeta::initCoeff() "<<std::endl;
100-
#endif
101-
if (m_beta==0.0){
102-
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
111+
if (m_beta == 0.0)
112+
{
113+
TkUtil::UTF8String messErr(TkUtil::Charset::UTF_8);
103114
messErr << "Erreur interne: EdgeMeshingPropertyBeta, le coefficient de resserrement n'a pas été initialisé";
104115
throw TkUtil::Exception(messErr);
105116
}
106117

107118
m_dernierIndice = 0;
108119
}
109120

121+
/*----------------------------------------------------------------------------*/
110122
double EdgeMeshingPropertyBeta::
111123
nextCoeff()
112124
{
113-
m_dernierIndice+=1;
125+
m_dernierIndice += 1;
114126

115127
double feta = 0.0;
116128
double eta = 0.0;
117-
if (m_sens){
118-
eta = ((double)m_dernierIndice)/((double)m_nb_edges);
119-
feta=resserre(eta, m_beta);
129+
if (m_sens)
130+
{
131+
eta = static_cast<double>(m_dernierIndice) / static_cast<double>(m_nb_edges);
132+
feta = resserre(eta, m_beta);
120133
}
121-
else {
122-
eta = ((double)m_nb_edges-m_dernierIndice)/((double)m_nb_edges);
123-
feta = 1.0-resserre(eta, m_beta);
134+
else
135+
{
136+
eta = (static_cast<double>(m_nb_edges) - m_dernierIndice) / static_cast<double>(m_nb_edges);
137+
feta = 1.0 - resserre(eta, m_beta);
124138
}
125139

126-
#ifdef _DEBUG_beta
127-
std::cout<<"EdgeMeshingPropertyBeta::nextCoeff (m_beta="<<m_beta<<") return "<<feta<<std::endl;
128-
#endif
129-
130140
return feta;
131141
}
142+
132143
/*----------------------------------------------------------------------------*/
133-
double EdgeMeshingPropertyBeta::resserre(double eta, double beta)
144+
double EdgeMeshingPropertyBeta::resserre(const double eta, const double beta)
134145
{
135-
double p = std::log((beta + 1.0)/(beta - 1.0)) * (1.0 - eta);
136-
return 1.0 + beta * (1.0 - std::exp(p)) / (1.0 + std::exp(p));
146+
double p = std::log((beta + 1.0) / (beta - 1.0)) * (1.0 - eta);
147+
return 1.0 + beta * (1.0 - std::exp(p)) / (1.0 + std::exp(p));
137148
}
149+
138150
/*----------------------------------------------------------------------------*/
139151
TkUtil::UTF8String EdgeMeshingPropertyBeta::
140152
getScriptCommand() const
141153
{
142-
TkUtil::UTF8String o (TkUtil::Charset::UTF_8);
143-
o << getMgx3DAlias() << ".EdgeMeshingPropertyBeta("
144-
<< static_cast<long>(m_nb_edges) << ","
145-
<< Utils::Math::MgxNumeric::userRepresentation(m_beta);
154+
TkUtil::UTF8String o(TkUtil::Charset::UTF_8);
155+
o << getMgx3DAlias() << ".EdgeMeshingPropertyBeta( ";
156+
o << static_cast<long>(m_nb_edges) << ", ";
157+
o << Utils::Math::MgxNumeric::userRepresentation(m_beta);
158+
o << (getDirect() ? std::string(", True") : std::string(", False"));
146159
if (m_initWithArm1)
147160
{
148-
o << (getDirect() ? std::string(", True") : std::string(", False"));
149161
o << ", True"; // this is for m_initWithArm1
150162
o << ", " << Utils::Math::MgxNumeric::userRepresentation(m_arm1);
151163
}
152-
else if (!getDirect())
153-
{
154-
o << ", False";
155-
}
156164
o << ")";
157165
return o;
158166
}
167+
159168
/*----------------------------------------------------------------------------*/
160169
void EdgeMeshingPropertyBeta::
161-
addDescription(Utils::SerializedRepresentation& topoProprietes) const
170+
addDescription(Utils::SerializedRepresentation &topoProprietes) const
162171
{
163-
topoProprietes.addProperty (
164-
Utils::SerializedRepresentation::Property ("Beta Resserrement", m_beta));
165-
166-
if (getDirect())
167-
topoProprietes.addProperty (
168-
Utils::SerializedRepresentation::Property ("Sens", std::string("direct")));
169-
else
170-
topoProprietes.addProperty (
171-
Utils::SerializedRepresentation::Property ("Sens", std::string("inverse")));
172+
topoProprietes.addProperty(
173+
Utils::SerializedRepresentation::Property("Beta Resserrement", m_beta));
174+
175+
if (getDirect())
176+
{
177+
topoProprietes.addProperty(
178+
Utils::SerializedRepresentation::Property("Sens", std::string("direct")));
179+
}
180+
else
181+
{
182+
topoProprietes.addProperty(
183+
Utils::SerializedRepresentation::Property("Sens", std::string("inverse")));
184+
}
172185
}
173186

187+
/*----------------------------------------------------------------------------*/
174188
double EdgeMeshingPropertyBeta::
175-
computeBeta(const double lg)
189+
computeBeta(const double lg) const
176190
{
177-
#ifdef _DEBUG_beta
178-
std::cout<<"EdgeMeshingPropertyBeta::computeBeta()"<<std::endl;
179-
#endif
180-
181-
// recherche de beta pour laquelle feta=m_arm1
182-
double beta = 0.0;
183-
double eta = 0.0;
184-
if (m_sens){
185-
eta = (1.0)/((double)m_nb_edges);
186-
std::cout << "eta = " << eta << std::endl;
187-
}
188-
else {
189-
eta = ((double)m_nb_edges-1.0)/((double)m_nb_edges);
190-
std::cout << "eta = " << eta << std::endl;
191-
}
192-
193-
// on limite la recherche
194-
double beta1 = 1.00001;
195-
double beta2 = 1.01;
196-
double lg1 = resserre(eta, beta1);
197-
double lg2 = resserre(eta, beta2);
198-
199-
std::cout << "lg1 = " << lg1 << ", lg2 = " << lg2 << ", lg = " << lg << std::endl;
200-
if ((lg1<lg && lg2<lg) || (lg1>lg && lg2>lg)){
201-
TkUtil::UTF8String messErr (TkUtil::Charset::UTF_8);
202-
messErr << "EdgeMeshingPropertyBeta, on ne peut pas trouver de beta pour lg "
203-
<<lg<<" et nombre de bras de "<<(short)m_nb_edges;
204-
throw TkUtil::Exception(messErr);
205-
}
206-
uint iter = 0;
207-
bool sens = (lg1<lg2);
208-
while (!Utils::Math::MgxNumeric::isNearlyZero(beta2-beta1) && iter < 100){
209-
iter ++;
210-
beta = (beta1+beta2)/2.0;
211-
212-
double lg_iter = resserre(eta, beta);
213-
#ifdef _DEBUG_beta
214-
//std::cout<<"iter "<<iter<<std::endl;
215-
//std::cout<<" beta "<<beta<<", lg_iter "<<lg_iter<<std::endl;
216-
#endif
217-
218-
if (sens){
219-
if (lg_iter>lg)
220-
beta2 = beta;
221-
else
222-
beta1 = beta;
223-
}
224-
else {
225-
if (lg_iter<lg)
226-
beta2 = beta;
227-
else
228-
beta1 = beta;
229-
}
230-
} // end while
231-
232-
#ifdef _DEBUG_beta
233-
std::cout<<"EdgeMeshingPropertyGeometric::computeBeta return "<<beta<<", iter "<<iter<<std::endl;
234-
#endif
235-
return beta;
191+
// recherche de beta pour laquelle feta=m_arm1
192+
double beta = 0.0;
193+
const double eta = (1.0) / static_cast<double>(m_nb_edges);
194+
195+
// on limite la recherche
196+
double beta1 = 1.00001;
197+
double beta2 = 1.01;
198+
const double lg1 = resserre(eta, beta1);
199+
const double lg2 = resserre(eta, beta2);
200+
201+
if ((lg1 < lg && lg2 < lg) || (lg1 > lg && lg2 > lg))
202+
{
203+
TkUtil::UTF8String messErr(TkUtil::Charset::UTF_8);
204+
messErr << "EdgeMeshingPropertyBeta, on ne peut pas trouver de beta pour lg ";
205+
messErr << lg << " et nombre de bras de " << static_cast<short>(m_nb_edges);
206+
throw TkUtil::Exception(messErr);
207+
}
208+
uint iter = 0;
209+
const bool sens = (lg1 < lg2);
210+
while (!Utils::Math::MgxNumeric::isNearlyZero(beta2 - beta1) && iter < 100)
211+
{
212+
iter++;
213+
beta = (beta1 + beta2) / 2.0;
214+
const double lg_iter = resserre(eta, beta);
215+
216+
if (sens)
217+
{
218+
if (lg_iter > lg)
219+
{
220+
beta2 = beta;
221+
}
222+
else
223+
{
224+
beta1 = beta;
225+
}
226+
}
227+
else
228+
{
229+
if (lg_iter < lg)
230+
{
231+
beta2 = beta;
232+
}
233+
else
234+
{
235+
beta1 = beta;
236+
}
237+
}
238+
} // end while
239+
240+
return beta;
236241
}
237242

238243
/*----------------------------------------------------------------------------*/

0 commit comments

Comments
 (0)