Soutenez-nous
Publicité
+ Répondre à la discussion
Page 4 sur 5 PremièrePremière 12345 DernièreDernière
Affichage des résultats 61 à 80 sur 81
  1. #61
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    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 :
    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;
    	}
    
    }

  2. #62
    Expert Confirmé
    Homme Profil pro
    Inscrit en
    septembre 2006
    Messages
    2 384
    Détails du profil
    Informations personnelles :
    Sexe : Homme

    Informations forums :
    Inscription : septembre 2006
    Messages : 2 384
    Points : 2 888
    Points
    2 888

    Par défaut

    Citation Envoyé par Nehmé Voir le message
    … 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 :
    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…)

  3. #63
    Membre chevronné

    Inscrit en
    septembre 2006
    Messages
    717
    Détails du profil
    Informations forums :
    Inscription : septembre 2006
    Messages : 717
    Points : 752
    Points
    752

    Par défaut

    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 :
    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);
    }

  4. #64
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    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.

  5. #65
    Membre Expert Avatar de Nemerle
    Inscrit en
    octobre 2003
    Messages
    1 106
    Détails du profil
    Informations personnelles :
    Âge : 43

    Informations forums :
    Inscription : octobre 2003
    Messages : 1 106
    Points : 1 071
    Points
    1 071

    Par défaut

    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).
    Nemerle, mathématicopilier de bars, membre du triumvirat du CSTM, 3/4 centre

  6. #66
    Membre chevronné

    Inscrit en
    septembre 2006
    Messages
    717
    Détails du profil
    Informations forums :
    Inscription : septembre 2006
    Messages : 717
    Points : 752
    Points
    752

    Par défaut

    Citation Envoyé par Nehmé Voir le message
    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.

  7. #67
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    Citation Envoyé par Sylvain Togni Voir le message
    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 !
    Images attachées Images attachées

  8. #68
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    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 !

  9. #69
    Membre chevronné

    Inscrit en
    septembre 2006
    Messages
    717
    Détails du profil
    Informations forums :
    Inscription : septembre 2006
    Messages : 717
    Points : 752
    Points
    752

    Par défaut

    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 :
    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);
    }

  10. #70
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    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.

  11. #71
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    Salut Sylvain,

    j'arrive pas à identifier l'algorithme mathématique que t'utilises pour calculer l'intersection de 2 droites avec ce code:
    Code :
    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.

  12. #72
    Expert Confirmé Sénior

    Inscrit en
    janvier 2007
    Messages
    10 134
    Détails du profil
    Informations personnelles :
    Âge : 56

    Informations forums :
    Inscription : janvier 2007
    Messages : 10 134
    Points : 12 719
    Points
    12 719

    Par défaut

    Citation Envoyé par Nehmé Voir le message
    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.


    Citation Envoyé par Sylvain Togni Voir le message
    Car comme souviron34 je pense qu'ils nous manque une information capitale pour pouvoir efficacement t'aider : la finalité de la chose.
    "Un homme sage ne croit que la moitié de ce qu’il lit. Plus sage encore, il sait laquelle".

    Consultant indépendant.
    Architecture systèmes complexes. Programmation grosses applications critiques. Ergonomie.
    C, Fortran, XWindow/Motif, Java

    Je ne réponds pas aux MP techniques

  13. #73
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    Citation Envoyé par souviron34 Voir le message




    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.

  14. #74
    Expert Confirmé
    Homme Profil pro
    Inscrit en
    septembre 2006
    Messages
    2 384
    Détails du profil
    Informations personnelles :
    Sexe : Homme

    Informations forums :
    Inscription : septembre 2006
    Messages : 2 384
    Points : 2 888
    Points
    2 888

    Par défaut

    Citation Envoyé par Nehmé Voir le message
    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…
    Images attachées Images attachées

  15. #75
    Expert Confirmé Sénior

    Inscrit en
    janvier 2007
    Messages
    10 134
    Détails du profil
    Informations personnelles :
    Âge : 56

    Informations forums :
    Inscription : janvier 2007
    Messages : 10 134
    Points : 12 719
    Points
    12 719

    Par défaut

    Citation Envoyé par Nehmé Voir le message
    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 ??
    "Un homme sage ne croit que la moitié de ce qu’il lit. Plus sage encore, il sait laquelle".

    Consultant indépendant.
    Architecture systèmes complexes. Programmation grosses applications critiques. Ergonomie.
    C, Fortran, XWindow/Motif, Java

    Je ne réponds pas aux MP techniques

  16. #76
    Membre chevronné

    Inscrit en
    septembre 2006
    Messages
    717
    Détails du profil
    Informations forums :
    Inscription : septembre 2006
    Messages : 717
    Points : 752
    Points
    752

    Par défaut

    Citation Envoyé par Nehmé Voir le message
    j'arrive pas à identifier l'algorithme mathématique que t'utilises pour calculer l'intersection de 2 droites avec ce code:
    Code :
    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 :
    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 :
    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.

  17. #77
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    Citation Envoyé par souviron34 Voir le message
    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.

  18. #78
    Expert Confirmé Sénior

    Inscrit en
    janvier 2007
    Messages
    10 134
    Détails du profil
    Informations personnelles :
    Âge : 56

    Informations forums :
    Inscription : janvier 2007
    Messages : 10 134
    Points : 12 719
    Points
    12 719

    Par défaut

    on y revient toujours

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

    Subject 2.09
    "Un homme sage ne croit que la moitié de ce qu’il lit. Plus sage encore, il sait laquelle".

    Consultant indépendant.
    Architecture systèmes complexes. Programmation grosses applications critiques. Ergonomie.
    C, Fortran, XWindow/Motif, Java

    Je ne réponds pas aux MP techniques

  19. #79
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    Citation Envoyé par souviron34 Voir le message
    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

  20. #80
    Membre régulier
    Profil pro Nehmé Bilal
    Inscrit en
    septembre 2008
    Messages
    248
    Détails du profil
    Informations personnelles :
    Nom : Nehmé Bilal
    Âge : 30
    Localisation : Canada

    Informations forums :
    Inscription : septembre 2008
    Messages : 248
    Points : 79
    Points
    79

    Par défaut

    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

Liens sociaux

Règles de messages

  • Vous ne pouvez pas créer de nouvelles discussions
  • Vous ne pouvez pas envoyer des réponses
  • Vous ne pouvez pas envoyer des pièces jointes
  • Vous ne pouvez pas modifier vos messages
  •