Discussion :

[Nehmé]

au cas où quelqu'un veut s'en servir, voici le code que j'ai utilisé pour l'algorithme du plus petit quadrilataire englobant.
on commence par trouver l'enveloppe convexe, ensuite on continue par la reduction des segments de l'enveloppe jusqu'à 4 ou un autre nombre de segments, en utilisant le prolongement des segments adjacents de JeitEmgie mais sans la droite oscillante. Le résultat n'est pas encore optimale à cause du fait que j'ai pas encore integrer la droite oscillante et/ou d'autres lacunes possibles dans l'algorithme.
Les commentaires sont en anglais parce que la partie enveloppe convexe je l'ai pris quelque part (Graham scan) et parce que je travail au Canada/Ontario, et donc je dois mettre les commentaires en anglais. vous verrez que mon anglais n'est pas génial !
Si quelqu'un ajoutes ou optimise l'algo , svp reposter le code ici ou envoyez moi une copie sur <nehmebilal@gmail.com>

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
 
void build_hull(deque<pair<double,double>> raw_points,
								deque<pair<double,double>> &output)
{
	if(raw_points.size()<=4)
	{
		output = raw_points;
		return;
	}
	output.clear();
//
// The initial array of points is stored in deque raw_points. I first
// sort it, which gives me the far left and far right points of the hull.
// These are special values, and they are stored off separately in the left
// and right members.
//
// I then go through the list of raw_points, and one by one determine whether
// each point is above or below the line formed by the right and left points.
// If it is above, the point is moved into the upper_partition_points sequence. If it
// is below, the point is moved into the lower_partition_points sequence. So the output
// of this routine is the left and right points, and the sorted points that are in the
// upper and lower partitions.
//
 
  std::pair<double,double> left;
  std::pair<double,double> right;
  std::deque< std::pair<double,double> > upper_partition_points;
  std::deque< std::pair<double,double> > lower_partition_points;
 
 
	//
	// Step one in partitioning the points is to sort the raw data
	//
	std::sort( raw_points.begin(), raw_points.end() );
	//
	// The the far left and far right points, remove them from the
	// sorted sequence and store them in special members
	//
	left = raw_points.front();
	raw_points.pop_front();
	//raw_points.erase(raw_points.begin());
	right = raw_points.back();
	raw_points.pop_back();
	//
	// Now put the remaining points in one of the two output sequences
	//
	for ( size_t i = 0 ; i < raw_points.size() ; i++ )
	{
			int dir = direction( left, right, raw_points[ i ] );
			if ( dir < 0 )
					upper_partition_points.push_back( raw_points[ i ] );
			else
					lower_partition_points.push_back( raw_points[ i ] );
	}
 
  //std::deque< std::pair<double,double> > lower_hull;
  std::deque< std::pair<double,double> > upper_hull;
 
  build_half_hull(lower_partition_points, output, left, right, 1 );
  build_half_hull(upper_partition_points, upper_hull, left, right, -1 );
 
	for(deque< pair<double,double> >::reverse_iterator it = upper_hull.rbegin()+1;
		it != upper_hull.rend()-1; it++)
	{
		output.push_back(*it);
	}
 
	if(output.size()<=4)
		return;
 
	simplifyHull(output,4);
}
//
// Building the hull consists of two procedures: building the lower and
// then the upper hull. The two procedures are nearly identical - the main
// difference between the two is the test for convexity. When building the upper
// hull, our rull is that the middle point must always be *above* the line formed
// by its two closest neighbors. When building the lower hull, the rule is that point
// must be *below* its two closest neighbors. We pass this information to the 
// building routine as the last parameter, which is either -1 or 1.
//
 
 
 
 
 
 
// This is the method that builds either the upper or the lower half convex
// hull. It takes as its input a copy of the input array, which will be the
// sorted list of points in one of the two halfs. It produces as output a list
// of the points in the corresponding convex hull.
//
// The factor should be 1 for the lower hull, and -1 for the upper hull.
 
void build_half_hull(std::deque< std::pair<double,double> > input,
                      std::deque< std::pair<double,double> > &output,
											std::pair<double,double> left, 
											std::pair<double,double> right,
                      int factor )
{
  //
  // The hull will always start with the left point, and end with the right
  // point. According, we start by adding the left point as the first point
  // in the output sequence, and make sure the right point is the last point 
  // in the input sequence.
  //
  input.push_back( right );
  output.push_back( left );
  //
  // The construction loop runs until the input is exhausted
  //
  while ( input.size() != 0 ) {
      //
      // Repeatedly add the leftmost point to the null, then test to see 
      // if a convexity violation has occured. If it has, fix things up
      // by removing the next-to-last point in the output suqeence until 
      // convexity is restored.
      //
      output.push_back( input.front() );
      //plot_hull( f, "adding a new point" );
			input.pop_front();
      //input.erase( input.begin() );
      while ( output.size() >= 3 ) {
          size_t end = output.size() - 1;
          if ( factor * direction( output[ end - 2 ], 
                                   output[ end ], 
                                   output[ end - 1 ] ) <= 0 ) {
              output.erase( output.begin() + end - 1 );
              //plot_hull( f, "backtracking" );
          }
          else
              break;
      }
  }
}
 
 
// We can do this by by translating the points so that p1 is at the origin,
// then taking the cross product of p0 and p2. The result will be positive,
// negative, or 0, meaning respectively that p2 has turned right, left, or
// is on a straight line.
//
static int direction( std::pair<double,double> p0,
                      std::pair<double,double> p1,
                      std::pair<double,double> p2 )
{
    return ( (p0.first - p1.first ) * (p2.second - p1.second ) )
           - ( (p2.first - p1.first ) * (p0.second - p1.second ) );
}
 
 
 
 
 
void simplifyHull(std::deque< std::pair<double,double> > &input,int numberOfEdges)
{
	int size = input.size();
	int i=0;
	int j;
	double a1,a2,b1,b2;
	double x, y, x1, y1, x2, y2;
	map<double,int> air;
	map<int,pair<double,double>> newPoint;
	double surface;
 
	while(size>numberOfEdges)
	{
		if(findIntersectionPoint(input[i],input[i+1],input[i+2],input[i+3],x,y,surface))
		{
			air[surface] = i;
			newPoint[i] = pair<double,double>(x,y);
		}
		i++;
 
		if(i == size-3)
		{
			if(findIntersectionPoint(input[i],input[i+1],input[i+2],input[0],x,y,surface))
			{
				air[surface] = i;
				newPoint[i] = pair<double,double>(x,y);
			}
			i++;
 
			if(findIntersectionPoint(input[i],input[i+1],input[0],input[1],x,y,surface))
			{
				air[surface] = i;
				newPoint[i] = pair<double,double>(x,y);
			}
			i++;
 
			if(findIntersectionPoint(input[i],input[0],input[1],input[2],x,y,surface))
			{
				air[surface] = i;
				newPoint[i] = pair<double,double>(x,y);
			}
 
			map<double,int>::iterator it = air.begin();
			j = it->second;
			if(j > size-4)
			{
				if(j == size-1)
				{
					input.pop_front();
					input.pop_front();
					input.push_front(newPoint[j]);
				}
				else if(j == size-2)
				{
					input.pop_back();
					input.pop_front();
					input.push_back(newPoint[j]);
				}
				else if(j == size-3)
				{
					input.pop_back();
					input.pop_back();
					input.push_back(newPoint[j]);
				}
			}
			else
			{
				input.erase(input.begin()+j+1, input.begin() + j+3);
				input.insert(input.begin()+j+1, newPoint[j]);
			}
			size += -1;
			i=0;
			air.clear();
		}
 
 
 
	}
 
}
 
 
 
 
bool findIntersectionPoint( pair<double,double> p0, 
													  pair<double,double> p1,
														pair<double,double> p2,
														pair<double,double> p3,
														double &x, double &y, double &surface)
 
{
	double deltax1;
	double deltay1;
	double deltax2;
	double deltay2;
	double x1,x2,y1,y2;
	double a1,a2,b1,b2;
 
	x1 = p1.first;
	y1 = p1.second;
	x2 = p2.first;
	y2 = p2.second;
 
	deltax1 = x1 - p0.first;
	deltay1 = y1-p0.second;
	deltax2 = x2 - p3.first;
	deltay2 = y2-p3.second;
 
	if(deltax1*deltay1*deltax2*deltay2 == 0)
	{
		if(deltax1 == 0 && deltax2 == 0) // the 2 segments are vertical
			return false;
 
		if(deltay1 == 0 && deltay2 == 0)	// the 2 segments are horizontal
			return false;
 
		if(deltax1 == 0 && deltay2 == 0)
		{
			x = x1;
			y = y2;
			surface = 0.5*abs( (x - x2)*(y - y1) ); // right triangle
			return true;
		}
 
		if(deltay1 == 0 && deltax2 == 0)
		{
			x = x2;
			y = y1;
			surface = 0.5*abs( (x - x1)*(y - y2) ); // right triangle
			return true;
		}
 
		if(deltax1 == 0)
		{
			x = x1;
			a2 = (y2 - p3.second)/(x2 - p3.first);
			b2 = y2-(a2*x2);
			y = a2*x + b2;
			surface = 0.5*abs( (x2 - x)*(y1 - y) );
			return true;
		}
		if(deltay1 == 0)
		{
			y = y1;
			a2 = (y2 - p3.second)/(x2 - p3.first);
			b2 = y2-(a2*x2);
			x = (y - b2)/a2;
			surface = 0.5*abs( (x1 - x)*(y2 - y) );
			return true;
		}
		if(deltax2 == 0)
		{
			x = x2;
			a1 = (y1 - p0.second)/(x1 - p0.first);
			b1 = y1-(a1*x1);
			y = a1*x + b1;
			surface = 0.5*abs( (x1 - x)*(y2 - y) );
			return true;
		}
		if(deltay2 == 0)
		{
			y = y2;
			a1 = (y1 - p0.second)/(x1 - p0.first);
			b1 = y1-(a1*x1);
			x = (y - b1)/a1;
			surface = 0.5*abs( (x2 - x)*(y1 - y) );
			return true;
		}
	}
	else
	{
		a1 = (y1-p0.second)/(x1 - p0.first);
		a2 = (y2-p3.second)/(x2 - p3.first);
 
		b1 = y1-(a1*x1);
		b2 = y2-(a2*x2);
		x = (b2-b1)/(a1-a2);
		y = (a1*x) + b1;
		surface = 0.5*abs( ( (x1 - x)*(y2 - y) ) - ( (x2 - x)*(y1 - y) ) );
		return true;
	}
 
}

[JeitEmgie]

… Le résultat n'est pas encore optimale à cause du fait que j'ai pas encore integrer la droite oscillante et/ou d'autres lacunes possibles dans l'algorithme.
…

pour cette optimisation, je vais regarder quand j'aurai un peu de temps si l'on peut affirmer qu'il n'y a qu'un maximum ou un minimum à la somme des 2 surfaces quand l'angle b de la droite oscillante au point S1 varie de 0 à PI - a1…

si c'est le cas, alors on pourra mettre en place un deuxième test (après le test a0+a1+a2 < 360°) pour savoir s'il vaut la peine de calculer ce minimum…

en effet si l'hypothèse est vérifiée cela signifie que l'angle b en allant de 0 à Pi - a1 ou de Pi - a1 vers 0, les surfaces doivent décroître pour que l'on aie un minimum… (si elles croissent on a un maximum donc >= aux 2 triangles quand b=0 et b=Pi-a1)

donc en calculant 4 surfaces (donc 6 surfaces de triangle) pour les valeurs de b = { 0, 0 + delta, PI-a1, PI-a1-delta} et en comparant le signe de leurs deltas on a notre réponse… (…SI… l'on peut garantir qu'il n'y a qu'un seul min ou max…)

(6 triangles car b=0 et b=PI-a1 ont un des 2 triangles = 0 … donc 4x2 - les 2 que l'on sait nulles par avance = 6)

Pour rappel la surface d'un triangle dont on connaît un côté A et les 2 angles adjacents b, c est :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
1/2 * |A|^2 + (sin(a) * sin(b) / sin (PI - b - c))

côté optimisation, maintenir un tableau contenant les sin des angles intérieurs et des sin(PI - a[i]) peut faire économiser pas mal de calcul "caché" (éviter de calculer X fois les mêmes sinus…)

[Sylvain Togni]

Moi je partirais plus vers un algorithme de type EM (Expectation / Minimisation) utilisant les moindres carrés. On part d'une approximation grossière, puis on l'affine petit à petit en minimisant les moindres carrés. Ca donnerait quelque chose comme ça (vite fait, à améliorer) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
// Computes the squared distance between segment [a,b] and point p
double squaredDistance(std::pair<double,double> const& a,
                       std::pair<double,double> const& b,
                       std::pair<double,double> const& p)
{
  double abx = b.first - a.first;
  double aby = b.second - a.second;
  double apx = p.first - a.first;
  double apy = p.second - a.second;
  double abap = abx*apx + aby*apy;
  if (abap <= 0)
    return apx*apx + apy*apy;
  double abab = abx*abx + aby*aby;
  if (abab <= abap)
  {
    double dx = b.first - p.first;
    double dy = b.second - p.second;
    return dx*dx + dy*dy;
  }
  else
  {
    double abh = apx*aby - apy*abx;
    return abh*abh/abab;
  }
}
 
// Finds the index of the nearest polygon's edge to the point p
std::size_t findNearesEdge(std::deque< std::pair<double,double> > const& polygon,
                           std::pair<double,double> const& p)
{
  std::size_t result = 0;
  double smallerD2 = DBL_MAX;
  for(std::size_t i = 0; i < polygon.size(); ++i)
  {
    double d2 = squaredDistance(i == 0 ? polygon.back() : polygon[i - 1], polygon[i], p);
    if (d2 < smallerD2)
    {
      smallerD2 = d2;
      result = i;
    }
  }
  return result;
}
 
// Mean least squares fitter of a line segment
struct LineFitter
{
  double s1;
  double sx;
  double sy;
  double sxx;
  double sxy;
  double syy;
 
  LineFitter()
  {
    s1 = sx = sy = sxx = sxy = syy = 0.0;
  }
 
  void addPoint(std::pair<double,double> const& p)
  {
    s1 += 1.0;
    sx += p.first;
    sy += p.second;
    sxx += p.first*p.first;
    sxy += p.first*p.second;
    syy += p.second*p.second;
  }
 
  void getCoefs(double& a, double& b, double& c) const
  {
    double y = (sxy*s1 - sx*sy)*2.0;
    double x = sy*sy - sx*sx + (sxx - syy)*s1;
    double theta = std::atan2(y, x)/2.0;
    a = -std::sin(theta);
    b = std::cos(theta);
    c = (sx*a + sy*b)/-s1;
  }
};
 
// Computes the intersection point between line 'a' and line 'b'
bool intersection(LineFitter const& a, LineFitter const& b, std::pair<double,double>& p)
{
  double a1, b1, c1;
  a.getCoefs(a1, b1, c1);
  double a2, b2, c2;
  b.getCoefs(a2, b2, c2);
  double det = a2*b1 - b2*a1;
  if (det == 0)
    return false;
  p.first = (b2*c1 - b1*c2)/det;
  p.second = (a1*c2 - a2*c1)/det;
  return true;
}
 
// Improves polygon approximation
bool adjustPolygonApproximation(std::deque< std::pair<double,double> > const& points,
                                std::deque< std::pair<double,double> >& polygon)
{
  // Update edge fitters
  std::vector<LineFitter> edgeFitters(polygon.size());
  for(std::size_t i = 0; i < points.size(); ++i)
    edgeFitters[findNearesEdge(polygon, points[i])].addPoint(points[i]);
 
  // Update polygon points
  for(std::size_t j = 0; j < polygon.size(); ++j)
  {
    if (!intersection(edgeFitters[j == 0 ? polygon.size() - 1 : j - 1], edgeFitters[j], polygon[j]))
      return false;
  }
  return true;
}
 
void findBestQuadrangle(std::deque< std::pair<double,double> > const& points,
                        std::deque< std::pair<double,double> >& result)
{
  // Finds an initial polygon approximation (the bounding box, for example)
  double xMin = DBL_MAX;
  double xMax = -DBL_MAX;
  double yMin = DBL_MAX;
  double yMax = -DBL_MAX;
  for(std::size_t i = 0; i < points.size(); ++i)
  {
    double x = points[i].first;
    if (x < xMin) xMin = x;
    if (x > xMax) xMax = x;
    double y = points[i].second;
    if (y < yMin) yMin = y;
    if (y > yMax) yMax = y;
  }
  result.clear();
  result.push_back(std::make_pair(xMin, yMin));
  result.push_back(std::make_pair(xMin, yMax));
  result.push_back(std::make_pair(xMax, yMax));
  result.push_back(std::make_pair(xMax, yMin));
 
  // Iterative improvements
  for(int i = 0; i < 20; ++i) // 20 iterations, for example
    adjustPolygonApproximation(points, result);
}

[Nehmé]

Le problème avec ta méthode Sylvain c'est que t'as pas vraiment dit grand chose sur qu'est ce que fait ton algorithme. tous qu'on sait c'est:
EM (Expectation / Minimisation) et moindres carrés

Je doute que quelqu'un voudra lire tout ton code pour savoir qu'est ce que ton algorithme fait.

[Nemerle]

Franchement, sur ce type de problème, je conseille un algo génétique qui solutionne le truc en 2 coups de cuillère à pot.

C'est exactement le même problème que de trouver le plus grand cercle inscrit dans une région ou des parcelles sont interdites (voir par exemple le site d'interstice, qui propose une applet Java).

[Sylvain Togni]

Le problème avec ta méthode Sylvain c'est que t'as pas vraiment dit grand chose sur qu'est ce que fait ton algorithme. tous qu'on sait c'est:
EM (Expectation / Minimisation) et moindres carrés

Je doute que quelqu'un voudra lire tout ton code pour savoir qu'est ce que ton algorithme fait.

Well, en fait j'ai écrit ce bout de code en pensant que tu l'essayerai, il n'y a pratiquement qu'un copier/coller à faire vu que j'ai repris les structures de données du code que tu as posté juste avant. Et s'il donne de meilleurs résultats j'étais prêt à l'expliquer plus en détails.

Car comme souviron34 je pense qu'ils nous manque une information capitale pour pouvoir efficacement t'aider : la finalité de la chose.

[Nehmé]

Well, en fait j'ai écrit ce bout de code en pensant que tu l'essayerai, il n'y a pratiquement qu'un copier/coller à faire vu que j'ai repris les structures de données du code que tu as posté juste avant. Et s'il donne de meilleurs résultats j'étais prêt à l'expliquer plus en détails.

Car comme souviron34 je pense qu'ils nous manque une information capitale pour pouvoir efficacement t'aider : la finalité de la chose.

Salut Sylvain,

Desolé je ne savais pas que tu avais fais ton code pour etre compatible avec le mien... Je viens de tester ton algorithme et il semble donner des resultats pas mal du tout. Il donne des bons resultats pour les surfaces complexes mais panique dans certain cas. L'autre chose, pour des cas simple comme une surface en forme rectangulaire par exemple, si la surface n'est pas orienté dans la direction des axes, l'approximation et un peu grossière.
La photo jointe montre les résultats.

Je n'ai pas encore lu ton code, mais les resultats que j ai obtenu avec sont assez bon. Je vais le lire dès que j'aurai le temps pour comprendre ton algo et essayer de eliminer les quelques erreurs et voir qu'est ce que ca donnerai.

Merci !

[Nehmé]

si t'as le temps, ca aidera beaucoup d'avoir quelque explications sur ton algo.
Je crois qu'il donnera des meilleurs resultats que ce que j'ai actuellement.
t'es pas obligé de detailler, juste les grande lignes.

merci encore !

[Sylvain Togni]

Good !

Alors pour l'explication commençons par les moindres carrés, tu connais peut-être déjà, ça consiste à ajuster (fitting) un modèle mathématique (dans notre cas une ligne droite) à un ensemble d'observations (dans notre cas des points) de sorte à minimiser une erreur quadratique (dans notre cas la distance euclidienne droite/points).

Si cela fonctionne très bien pour des modèles simples comme une ligne droite, c'est bien plus compliqué voir impossible pour des modèles plus complexes comme un quadrilatère. Mais un quadrilatère étant composé de quatre segments de droites, il suffit de segmenter l'ensemble des points en quatre groupes et appliquer les moindres carrés sur chaque groupe pour obtenir les quatre côtés du quadrilatère. Le problème étant comment segmenter l'ensemble des points en quatre groupes : l'idéal est d'associer chaque point au côté du quadrilatère le plus proche. Mais comme on ne connait pas encore le quadrilatère, on tourne en rond.

C'est là qu'intervient l'algorithme EM (Expectation Maximisation). L'idée est de partir d'une approximation, même grossière, du résultat (dans mon code j'ai pris le bounding box mais il y a peut-être mieux à faire), en déduire une segmentation des points en quatre groupes, appliquer les moindres carrés sur chaque groupe, calculer le nouveau quadrilatère et ainsi de suite un certain nombre de fois.

Petit à petit le quadrilatère obtenu converge vers la solution optimale. Sauf dans certains cas particuliers ou il peut y avoir divergence. Un moyen de minimiser le risque de divergence c'est d'utiliser une segmentation floue : les points ne sont plus associés à un unique côté mais aux quatre côtés à la fois avec des probabilités différentes, inversement proportionnelles à la distance (std::exp(-d2)). Voici ce que ça donne (le code est plus simple, en plus) :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
// Computes the squared distance between segment [a,b] and point p
double squaredDistance(std::pair<double,double> const& a,
                       std::pair<double,double> const& b,
                       std::pair<double,double> const& p)
{
  double abx = b.first - a.first;
  double aby = b.second - a.second;
  double apx = p.first - a.first;
  double apy = p.second - a.second;
  double abap = abx*apx + aby*apy;
  if (abap <= 0)
    return apx*apx + apy*apy;
  double abab = abx*abx + aby*aby;
  if (abab <= abap)
  {
    double dx = b.first - p.first;
    double dy = b.second - p.second;
    return dx*dx + dy*dy;
  }
  else
  {
    double abh = apx*aby - apy*abx;
    return abh*abh/abab;
  }
}
 
// Mean least squares fitter of a line segment
struct LineFitter
{
  double s1;
  double sx;
  double sy;
  double sxx;
  double sxy;
  double syy;
 
  LineFitter()
  {
    s1 = sx = sy = sxx = sxy = syy = 0.0;
  }
 
  void addPoint(std::pair<double,double> const& p, double w)
  {
    s1 += w;
    sx += w*p.first;
    sy += w*p.second;
    sxx += w*p.first*p.first;
    sxy += w*p.first*p.second;
    syy += w*p.second*p.second;
  }
 
  void getCoefs(double& a, double& b, double& c) const
  {
    double y = (sxy*s1 - sx*sy)*2.0;
    double x = sy*sy - sx*sx + (sxx - syy)*s1;
    double theta = std::atan2(y, x)/2.0;
    a = -std::sin(theta);
    b = std::cos(theta);
    c = (sx*a + sy*b)/-s1;
  }
};
 
// Computes the intersection point between line 'a' and line 'b'
bool intersection(LineFitter const& a, LineFitter const& b, std::pair<double,double>& p)
{
  double a1, b1, c1;
  a.getCoefs(a1, b1, c1);
  double a2, b2, c2;
  b.getCoefs(a2, b2, c2);
  double det = a2*b1 - b2*a1;
  if (det == 0)
    return false;
  p.first = (b2*c1 - b1*c2)/det;
  p.second = (a1*c2 - a2*c1)/det;
  return true;
}
 
// Improves polygon approximation
bool adjustPolygonApproximation(std::deque< std::pair<double,double> > const& points,
                                std::deque< std::pair<double,double> >& polygon)
{
  // Update edge fitters
  std::vector<LineFitter> edgeFitters(polygon.size());
  for(std::size_t i = 0; i < points.size(); ++i)
  {
    for(std::size_t j = 0; j < polygon.size(); ++j)
    {
      double d2 = squaredDistance(polygon[j], polygon[(j + 1)%polygon.size()], points[i]);
      edgeFitters[j].addPoint(points[i], std::exp(-d2));
    }
  }
 
  // Update polygon points
  for(std::size_t j = 0; j < polygon.size(); ++j)
  {
    if (!intersection(edgeFitters[j == 0 ? polygon.size() - 1 : j - 1], edgeFitters[j], polygon[j]))
      return false;
  }
  return true;
}
 
void findBestQuadrangle(std::deque< std::pair<double,double> > const& points,
                        std::deque< std::pair<double,double> >& result)
{
  // Finds an initial polygon approximation (the bounding box, for example)
  double xMin = DBL_MAX;
  double xMax = -DBL_MAX;
  double yMin = DBL_MAX;
  double yMax = -DBL_MAX;
  for(std::size_t i = 0; i < points.size(); ++i)
  {
    double x = points[i].first;
    if (x < xMin) xMin = x;
    if (x > xMax) xMax = x;
    double y = points[i].second;
    if (y < yMin) yMin = y;
    if (y > yMax) yMax = y;
  }
  result.clear();
  result.push_back(std::make_pair(xMin, yMin));
  result.push_back(std::make_pair(xMin, yMax));
  result.push_back(std::make_pair(xMax, yMax));
  result.push_back(std::make_pair(xMax, yMin));
 
  // Iterative improvements
  for(int i = 0; i < 20; ++i) // 20 iterations, for example
    adjustPolygonApproximation(points, result);
}

[Nehmé]

Salut,

j'ai essayé l'algorithme avec la logique floue, mais les cas de divergence existe toujours. J'ai toujours pas lu ton code comme il faut, à cause de certaine contrainte de temps au travail. Je crois que je vais l'essayé en partant d'un bounding box orienté dont j'ai dejà programmé l'algorithme et j'imagine que les résultats seront pas mal mieux.

[Nehmé]

Salut Sylvain,

j'arrive pas à identifier l'algorithme mathématique que t'utilises pour calculer l'intersection de 2 droites avec ce code:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
  void addPoint(std::pair<double,double> const& p)
  {
    s1 += 1.0;		// number of points
    sx += p.first;
    sy += p.second;
    sxx += p.first*p.first;
    sxy += p.first*p.second;
    syy += p.second*p.second;
  }
 
  void getCoefs(double& a, double& b, double& c) const
  {
    double y = (sxy*s1 - sx*sy)*2.0;
    double x = sy*sy - sx*sx + (sxx - syy)*s1;
    double theta = atan2(y, x)/2.0;
		a = -sin(theta);
		b = cos(theta);
    c = (sx*a + sy*b)/-s1;
  }
};
 
 
 
// Computes the intersection point between line 'a' and line 'b'
bool intersection(LineFitter const& a, LineFitter const& b, std::pair<double,double>& p)
{
  double a1, b1, c1;
  a.getCoefs(a1, b1, c1);
  double a2, b2, c2;
  b.getCoefs(a2, b2, c2);
  double det = a2*b1 - b2*a1;
  if (det == 0)
    return false;
  p.first = (b2*c1 - b1*c2)/det;
  p.second = (a1*c2 - a2*c1)/det;
  return true;
}

peux tu m'expliquer d'où ca viens ?

J'ai essayé ton algorithme en partant d'un bounding rectangle orienté et une seul iteration suffit pour donner un bon resultat, mais la divergence existe toujours dans certain cas. c'est pourquoi j'essais de voir de ou ca viens dans le code.

[souviron34]

J'ai essayé ton algorithme en partant d'un bounding rectangle orienté et une seul iteration suffit pour donner un bon resultat, mais la divergence existe toujours dans certain cas. c'est pourquoi j'essais de voir de ou ca viens dans le code.

Car comme souviron34 je pense qu'ils nous manque une information capitale pour pouvoir efficacement t'aider : la finalité de la chose.

[Nehmé]

La finalité de la chose est d'approximer une surface par un quadrilatère. tous ce que je veux c'est que le quadrilatère couvre le mieux la surface. S'il dépasse un peu vers l'exterieur ou vers l'intérieur ce n'est pas grave. Cette surface est une face d'une roche et j'aimerai que le quadrilatère représente le mieux cette surface en terme de périmètre ou surface et en terme de position.

[JeitEmgie]

La finalité de la chose est d'approximer une surface par un quadrilatère. tous ce que je veux c'est que le quadrilatère couvre le mieux la surface. S'il dépasse un peu vers l'exterieur ou vers l'intérieur ce n'est pas grave. Cette surface est une face d'une roche et j'aimerai que le quadrilatère représente le mieux cette surface en terme de périmètre ou surface et en terme de position.

alors, quand vous détectez que l'on est dans la situation où l'optimum est la droite "oscillante" :
soit vous prenez comme position celle où P02-P13 est perpendiculaire à S1-S'03 (voir le schéma ci-avant…)
soit vous prenez la parallèle à S'02; S'13 passsant par S1

vous obtiendrez un compromis raisonnable… aussi bien vis-à-vis de la forme que de la surface…

et pour les autres cas vous gardez la comparaison des surfaces des triangles ajoutés…

OU

dans tous les cas, vous prenez l'une des 2 droites citées ci-dessus,
et le critère de choix doit alors être un tri par la surface des paires de triangles consécutifs…

[souviron34]

S'il dépasse un peu vers l'exterieur ou vers l'intérieur ce n'est pas grave.

je croyais qu'on parlait du plus petit quadrilatère englobant ???

La solution de Jeitmgie est correcte, mais il me semble qu'il y en a d'autres, assez simples..

Je réfléchirais là-dessus...

Mais j'aurais quand même aimé savoir : ce qu'on a au départ (sans doute le maillage 3D ??). Ce qu'on veut avoir à la fin : pourquoi découper par tranche ? une raison particulière ? pourquoi des quadrilatères ? une raison particulière ?

Est-ce qu'en fait le fond du problème n'est pas d'avoir l'enveloppe du maillage 3D ??

[Sylvain Togni]

j'arrive pas à identifier l'algorithme mathématique que t'utilises pour calculer l'intersection de 2 droites avec ce code:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
  void addPoint(std::pair<double,double> const& p)
  {
    s1 += 1.0;		// number of points
    sx += p.first;
    sy += p.second;
    sxx += p.first*p.first;
    sxy += p.first*p.second;
    syy += p.second*p.second;
  }

Ca c'est les moindres carrés classiques : le calcul des moments (sommes) d'ordre 0, 1 et 2 en X et Y

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
  void getCoefs(double& a, double& b, double& c) const
  {
    double y = (sxy*s1 - sx*sy)*2.0;
    double x = sy*sy - sx*sx + (sxx - syy)*s1;
    double theta = atan2(y, x)/2.0;
		a = -sin(theta);
		b = cos(theta);
    c = (sx*a + sy*b)/-s1;
  }

Toujours les moindres carrés : on calcule les coefficients (a, b, c) de la droite définie par l'équation ax + by + c = 0 qui minimisent l'erreur quadratique.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
// Computes the intersection point between line 'a' and line 'b'
bool intersection(LineFitter const& a, LineFitter const& b, std::pair<double,double>& p)
{
  double a1, b1, c1;
  a.getCoefs(a1, b1, c1);
  double a2, b2, c2;
  b.getCoefs(a2, b2, c2);
  double det = a2*b1 - b2*a1;
  if (det == 0)
    return false;
  p.first = (b2*c1 - b1*c2)/det;
  p.second = (a1*c2 - a2*c1)/det;
  return true;
}

Et là c'est un calcul classique d'intersection entre 2 droites a1*x + b1*y + c1 = 0 et a2*x + b2*y + c2 = 0 (système de 2 équations à 2 inconnues). Si det == 0 c'est que les 2 droites sont parallèles, c'est une cause de divergence, il faudrait étudier les cas où cela arrive pour trouver une solution de remplacement.

J'ai essayé ton algorithme en partant d'un bounding rectangle orienté et une seul iteration suffit pour donner un bon resultat, mais la divergence existe toujours dans certain cas. c'est pourquoi j'essais de voir de ou ca viens dans le code.

Tu utilise bien le second code (segmentation floue) ? A part le problème de côtés parallèles, je ne vois pas d'autres cas de divergence.

[Nehmé]

je croyais qu'on parlait du plus petit quadrilatère englobant ???

La solution de Jeitmgie est correcte, mais il me semble qu'il y en a d'autres, assez simples..

Je réfléchirais là-dessus...

Mais j'aurais quand même aimé savoir : ce qu'on a au départ (sans doute le maillage 3D ??). Ce qu'on veut avoir à la fin : pourquoi découper par tranche ? une raison particulière ? pourquoi des quadrilatères ? une raison particulière ?

Est-ce qu'en fait le fond du problème n'est pas d'avoir l'enveloppe du maillage 3D ??

Je veux te dire pourquoi vu que tu insistes
Les données d'entrée de mon algorithme sont effectivement des maillages 3D, pourquoi des quadrilatère : parce que les données de sorties de mon algorithme doivent être écrit dans un format d'un logiciel dont les élements possible sont soit un hexaèdre soit un quadrilatère, ce qui fait en sorte que mon élements le plus simple est un quadrilatère. Pourquoi je découpe en tranche: parce que selon la nature de mes données, je sais que les tranche sont toujours perpendiculaire à Z et donc je fais juste profiter de la situation pour simplifier mon algorithme. Le but c'est de passer du maillage 3D triangulaire au maillage en quadrilatère ou hexahères.

[souviron34]

on y revient toujours

http://www.faqs.org/faqs/graphics/algorithms-faq/

Subject 2.09

[Nehmé]

on y revient toujours

http://www.faqs.org/faqs/graphics/algorithms-faq/

Subject 2.09

lol .. Il me semble que t'as pas lu la discussion surviron.

2.09: How can I find the minimum area rectangle enclosing a set of points?

Ceci n'est pas ce qu'on cherche et d'ailleurs on a dejà la solution

[Nehmé]

Salut,

Vue que j'ai été affecté à un autre projet et que je pourrai pas revenir à celui ci avant un bout de temps, est ce que je devrais déclarer cette solution comme résolu selon la politique du Forum ?
(JeitEmgie Je n'ai pas eu le temps de tester ton dernier critère, je le ferai dès que je reviendrai sur ce projet )

Merci à tous