asr_approx_mvbb: util.cpp Source File

Go to the documentation of this file.
00001 // ========================================================================================
00002 //  ApproxMVBB
00003 //  Copyright (C) 2014 by Gabriel Nützi <nuetzig (at) imes (d0t) mavt (d0t) ethz (døt) ch>
00004 //
00005 //  This Source Code Form is subject to the terms of the Mozilla Public
00006 //  License, v. 2.0. If a copy of the MPL was not distributed with this
00007 //  file, You can obtain one at http://mozilla.org/MPL/2.0/.
00008 // ========================================================================================
00009 
00010 #include <iostream>
00011 #include <limits>
00012 
00013 #include <ApproxMVBB/Diameter/Utils/util.hpp>
00014 
00015 namespace ApproxMVBB{
00016 namespace Diameter{
00017 
00018 /* partially sort a list of points
00019 
00020    given a "diameter", points which are 'outside'
00021    the sphere with a given threshold are put at the beginning
00022    of the list,
00023    the returned value is the index of the last point 'outside'
00024    ie
00025    points from #first    to #index are outside
00026      "     "    #index+1 to #last  are inside
00027 
00028    If there are no points outside, #first-1 is returned.
00029 
00030    Given a segment [AB], its centre C is (A+B)/2.
00031    The dot product MA.MB is equal to MC^2 - |AB|^2 / 4
00032 
00033    1. Being naive, ie to test if points are outside
00034       the ball of diameter [AB], yield to condition
00035       MA.MB > 0
00036 
00037    2. Being a little smarter, ie to  test if points are outside
00038       the ball of center (A+B)/2 and diameter D (assume D>|AB|)
00039       yield to condition
00040       MA.MB > ( D^2 - |AB|^2 ) / 4
00041 
00042 */
00043 
00044 int _LastPointOutsideSphereWithDiameter( TypeSegment *theSeg,
00045                                          const double squareDiameter,
00046                                          const double **theList,
00047                                          const int first,
00048                                          int *last,
00049                                          const int dim,
00050                                          const int _reduction_mode_ )
00051 {
00052   int i;
00053   int index = first - 1;
00054   int l = *last;
00055 
00056   double mamb, am2;
00057 
00058   double minThreshold;
00059   double medThreshold;
00060   double maxThreshold;
00061 
00062   double ab  = sqrt( theSeg->squareDiameter );
00063   double ab2 = theSeg->squareDiameter;
00064 
00065   double R  = sqrt( squareDiameter );
00066   double R2 = squareDiameter;
00067 
00068   if ( first > *last ) return( first - 1 );
00069 
00070   if ( squareDiameter <= theSeg->squareDiameter ) {
00071 
00072     maxThreshold = medThreshold = 0.0;
00073     minThreshold = -0.23205080756887729352 * theSeg->squareDiameter;
00074 
00075   } else {
00076 
00077     minThreshold = (R - .86602540378443864676*ab)
00078       * (R - .86602540378443864676*ab)
00079       - 0.25 * ab2;
00080 
00081     medThreshold = 0.5 * ab2 + 0.25 * R2
00082       - .43301270189221932338 * R * sqrt( 4 * ab2 - R2 );
00083 
00084     maxThreshold = 0.25 * (R2 - ab2);
00085   }
00086 
00087 
00088 
00089   /* 0 : no reduction
00090      1 : just inside the smallest sphere
00091      2 : complete reduction
00092   */
00093   switch ( _reduction_mode_ ) {
00094 
00095   case 0 :
00096 
00097     /* NO REDUCTION CASE
00098      */
00099     if ( dim == 2 ) {
00100       for ( i=first; i<=l; i++ ) {
00101         mamb = _ScalarProduct2D( theList[i], theSeg->extremity1,
00102                                  theList[i], theSeg->extremity2 );
00103         if ( mamb > maxThreshold ) {
00104           index ++;
00105           _SwapPoints( theList, index, i );
00106         }
00107       }
00108       return( index );
00109     }
00110 
00111     if ( dim == 3 ) {
00112       for ( i=first; i<=l; i++ ) {
00113         mamb = _ScalarProduct3D( theList[i], theSeg->extremity1,
00114                                  theList[i], theSeg->extremity2 );
00115         if ( mamb > maxThreshold ) {
00116           index ++;
00117           _SwapPoints( theList, index, i );
00118         }
00119       }
00120       return( index );
00121     }
00122 
00123     for ( i=first; i<=l; i++ ) {
00124       mamb = _ScalarProduct( theList[i], theSeg->extremity1,
00125                              theList[i], theSeg->extremity2, dim );
00126       if ( mamb > maxThreshold ) {
00127         index ++;
00128         _SwapPoints( theList, index, i );
00129       }
00130     }
00131     return( index );
00132 
00133     /* END
00134        NO REDUCTION CASE
00135      */
00136     break;
00137 
00138   default :
00139   case 1 :
00140 
00141     /* REDUCTION INSIDE THE SMALLEST SPHERE
00142      */
00143 
00144     if ( dim == 2 ) {
00145       for ( i=first; i<=l; i++ ) {
00146         mamb = _ScalarProduct2D( theList[i], theSeg->extremity1,
00147                                  theList[i], theSeg->extremity2 );
00148         if ( mamb > maxThreshold ) {
00149           index ++;
00150           _SwapPoints( theList, index, i );
00151         }
00152         else if ( mamb <= minThreshold ) {
00153           _SwapPoints( theList, i, l );
00154           i--;   l--;
00155         }
00156       }
00157       *last = l;
00158       return( index );
00159     }
00160 
00161     if ( dim == 3 ) {
00162       for ( i=first; i<=l; i++ ) {
00163         mamb = _ScalarProduct3D( theList[i], theSeg->extremity1,
00164                                  theList[i], theSeg->extremity2 );
00165         if ( mamb > maxThreshold ) {
00166           index ++;
00167           _SwapPoints( theList, index, i );
00168         }
00169         else if ( mamb <= minThreshold ) {
00170           _SwapPoints( theList, i, l );
00171           i--;   l--;
00172         }
00173       }
00174       *last = l;
00175       return( index );
00176     }
00177 
00178     for ( i=first; i<=l; i++ ) {
00179       mamb = _ScalarProduct( theList[i], theSeg->extremity1,
00180                              theList[i], theSeg->extremity2, dim );
00181       if ( mamb > maxThreshold ) {
00182         index ++;
00183         _SwapPoints( theList, index, i );
00184       }
00185       else if ( mamb <= minThreshold ) {
00186         _SwapPoints( theList, i, l );
00187         i--;   l--;
00188       }
00189     }
00190     *last = l;
00191     return( index );
00192 
00193     /* END
00194        REDUCTION INSIDE THE SMALLEST SPHERE
00195      */
00196     break;
00197 
00198   case 2 :
00199 
00200     /* COMPLETE REDUCTION
00201 
00202        On suppose implicitement
00203        que l'ensemble des points est
00204        compris dans l'intersection des spheres centrees
00205        sur A et B et de rayon |AB|.
00206        En effet, les conditions (pour un rayon dans
00207        l'intervalle [min,med]) pour l'elimination ne le
00208        verifient pas.
00209      */
00210     if ( dim == 2 ) {
00211 
00212       if ( squareDiameter <= theSeg->squareDiameter ) {
00213 
00214         for ( i=first; i<=l; i++ ) {
00215           mamb = _ScalarProduct2D( theList[i], theSeg->extremity1,
00216                                    theList[i], theSeg->extremity2 );
00217           if ( mamb > maxThreshold ) {
00218             index ++;
00219             _SwapPoints( theList, index, i );
00220           }
00221           else if (  mamb > minThreshold ) {
00222             am2 = _SquareDistance2D( theList[i], theSeg->extremity1 );
00223             if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) ) - mamb*mamb < 0 ) {
00224               _SwapPoints( theList, i, l );
00225               i--;   l--;
00226             }
00227           }
00228           else {
00229             /* if ( mamb <= minThreshold ) */
00230             _SwapPoints( theList, i, l );
00231             i--;   l--;
00232           }
00233         }
00234 
00235       } else {
00236 
00237         for ( i=first; i<=l; i++ ) {
00238           mamb = _ScalarProduct2D( theList[i], theSeg->extremity1,
00239                                    theList[i], theSeg->extremity2  );
00240           if ( mamb > maxThreshold ) {
00241             index ++;
00242             _SwapPoints( theList, index, i );
00243           }
00244           else if ( mamb > medThreshold ) {
00245             continue;
00246           }
00247           else if ( mamb > minThreshold ) {
00248             am2 = _SquareDistance2D( theList[i], theSeg->extremity1 );
00249             if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) )
00250                  - (R2 - ab2 - mamb)*(R2 - ab2 - mamb) < 0 ) {
00251               _SwapPoints( theList, i, l );
00252               i--;   l--;
00253             }
00254           }
00255           else {
00256             /* if ( mamb <= minThreshold ) */
00257             _SwapPoints( theList, i, l );
00258             i--;   l--;
00259           }
00260         }
00261 
00262       }
00263       *last = l;
00264       return( index );
00265     }
00266 
00267     if ( dim == 3 ) {
00268 
00269       if ( squareDiameter <= theSeg->squareDiameter ) {
00270 
00271         for ( i=first; i<=l; i++ ) {
00272           mamb = _ScalarProduct3D( theList[i], theSeg->extremity1,
00273                                    theList[i], theSeg->extremity2 );
00274           if ( mamb > maxThreshold ) {
00275             index ++;
00276             _SwapPoints( theList, index, i );
00277           }
00278           else if (  mamb > minThreshold ) {
00279             am2 = _SquareDistance3D( theList[i], theSeg->extremity1 );
00280             if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) ) - mamb*mamb < 0 ) {
00281               _SwapPoints( theList, i, l );
00282               i--;   l--;
00283             }
00284           }
00285           else {
00286             /* if ( mamb <= minThreshold ) */
00287             _SwapPoints( theList, i, l );
00288             i--;   l--;
00289           }
00290         }
00291 
00292       } else {
00293 
00294         for ( i=first; i<=l; i++ ) {
00295           mamb = _ScalarProduct3D( theList[i], theSeg->extremity1,
00296                                    theList[i], theSeg->extremity2  );
00297           if ( mamb > maxThreshold ) {
00298             index ++;
00299             _SwapPoints( theList, index, i );
00300           }
00301           else if ( mamb > medThreshold ) {
00302             continue;
00303           }
00304           else if ( mamb > minThreshold ) {
00305             am2 = _SquareDistance3D( theList[i], theSeg->extremity1 );
00306             if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) )
00307                  - (R2 - ab2 - mamb)*(R2 - ab2 - mamb) < 0 ) {
00308               _SwapPoints( theList, i, l );
00309               i--;   l--;
00310             }
00311           }
00312           else {
00313             /* if ( mamb <= minThreshold ) */
00314             _SwapPoints( theList, i, l );
00315             i--;   l--;
00316           }
00317         }
00318 
00319       }
00320       *last = l;
00321       return( index );
00322     }
00323 
00324     if ( squareDiameter <= theSeg->squareDiameter ) {
00325 
00326       for ( i=first; i<=l; i++ ) {
00327         mamb = _ScalarProduct( theList[i], theSeg->extremity1,
00328                                theList[i], theSeg->extremity2, dim );
00329         if ( mamb > maxThreshold ) {
00330           index ++;
00331           _SwapPoints( theList, index, i );
00332         }
00333         else if (  mamb > minThreshold ) {
00334           am2 = _SquareDistance( theList[i], theSeg->extremity1, dim );
00335           if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) ) - mamb*mamb < 0 ) {
00336             _SwapPoints( theList, i, l );
00337             i--;   l--;
00338           }
00339         }
00340         else {
00341           /* if ( mamb <= minThreshold ) */
00342           _SwapPoints( theList, i, l );
00343           i--;   l--;
00344         }
00345       }
00346 
00347     } else {
00348 
00349       for ( i=first; i<=l; i++ ) {
00350         mamb = _ScalarProduct( theList[i], theSeg->extremity1,
00351                                theList[i], theSeg->extremity2, dim );
00352         if ( mamb > maxThreshold ) {
00353           index ++;
00354           _SwapPoints( theList, index, i );
00355         }
00356         else if ( mamb > medThreshold ) {
00357           continue;
00358         }
00359         else if ( mamb > minThreshold ) {
00360           am2 = _SquareDistance( theList[i], theSeg->extremity1, dim );
00361           if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) )
00362                - (R2 - ab2 - mamb)*(R2 - ab2 - mamb) < 0 ) {
00363             _SwapPoints( theList, i, l );
00364             i--;   l--;
00365           }
00366         }
00367         else {
00368           /* if ( mamb <= minThreshold ) */
00369           _SwapPoints( theList, i, l );
00370           i--;   l--;
00371         }
00372       }
00373 
00374     }
00375     *last = l;
00376     return( index );
00377 
00378     /* END
00379        COMPLETE REDUCTION
00380     */
00381   }
00382   return( -1 );
00383 }
00384 
00385 
00386 int _LastPointOutsideSphereAndBoundWithDiameter( TypeSegment *theSeg,
00387                                                  const double squareDiameter,
00388                                                  const double **theList,
00389                                                  const int first,
00390                                                  int *last,
00391                                                  const int dim,
00392                                                  const int _reduction_mode_,
00393                                                  double *bound )
00394 {
00395   int i;
00396   int index = first - 1;
00397   int l = *last;
00398 
00399   double mamb, am2;
00400 
00401   double minThreshold;
00402   double medThreshold;
00403   double maxThreshold;
00404 
00405   double ab  = sqrt( theSeg->squareDiameter );
00406   double ab2 = theSeg->squareDiameter;
00407 
00408   double R  = sqrt( squareDiameter );
00409   double R2 = squareDiameter;
00410 
00411   /* bound
00412      c'est la plus grande valeur pour les points
00413      l'interieur de la sphere
00414      si squareDiameter > theSeg->squareDiameter
00415      cette valeur peut etre positive
00416   */
00417   double b = *bound = (- theSeg->squareDiameter * 0.25);
00418 
00419 
00420   if ( squareDiameter <= theSeg->squareDiameter ) {
00421 
00422     maxThreshold = medThreshold = 0.0;
00423     minThreshold = -0.23205080756887729352 * theSeg->squareDiameter;
00424 
00425     return( _LastPointOutsideSphereWithDiameter( theSeg, squareDiameter,
00426                                                  theList, first, last, dim,
00427                                                  _reduction_mode_ ) );
00428 
00429   } else {
00430 
00431     minThreshold = (R - .86602540378443864676*ab)
00432       * (R - .86602540378443864676*ab)
00433       - 0.25 * ab2;
00434 
00435     medThreshold = 0.5 * ab2 + 0.25 * R2
00436       - .43301270189221932338 * R * sqrt( 4 * ab2 - R2 );
00437 
00438     maxThreshold = 0.25 * (R2 - ab2);
00439   }
00440 
00441 
00442 
00443   /* 0 : no reduction
00444      1 : just inside the smallest sphere
00445      2 : complete reduction
00446   */
00447   switch ( _reduction_mode_ ) {
00448 
00449   case 0 :
00450 
00451     /* NO REDUCTION CASE
00452      */
00453     if ( dim == 2 ) {
00454       for ( i=first; i<=l; i++ ) {
00455         mamb = _ScalarProduct2D( theList[i], theSeg->extremity1,
00456                                  theList[i], theSeg->extremity2 );
00457         if ( mamb > maxThreshold ) {
00458           index ++;
00459           _SwapPoints( theList, index, i );
00460           continue;
00461         }
00462         if ( b < mamb ) b = mamb;
00463       }
00464       *bound = b;
00465       return( index );
00466     }
00467 
00468     if ( dim == 3 ) {
00469       for ( i=first; i<=l; i++ ) {
00470         mamb = _ScalarProduct3D( theList[i], theSeg->extremity1,
00471                                  theList[i], theSeg->extremity2 );
00472         if ( mamb > maxThreshold ) {
00473           index ++;
00474           _SwapPoints( theList, index, i );
00475           continue;
00476         }
00477         if ( b < mamb ) b = mamb;
00478       }
00479       *bound = b;
00480       return( index );
00481     }
00482 
00483     for ( i=first; i<=l; i++ ) {
00484       mamb = _ScalarProduct( theList[i], theSeg->extremity1,
00485                              theList[i], theSeg->extremity2, dim );
00486       if ( mamb > maxThreshold ) {
00487         index ++;
00488         _SwapPoints( theList, index, i );
00489         continue;
00490       }
00491       if ( b < mamb ) b = mamb;
00492     }
00493     *bound = b;
00494     return( index );
00495 
00496     /* END
00497        NO REDUCTION CASE
00498      */
00499     break;
00500 
00501   default :
00502   case 1 :
00503 
00504     /* REDUCTION INSIDE THE SMALLEST SPHERE
00505      */
00506 
00507     if ( dim == 2 ) {
00508       for ( i=first; i<=l; i++ ) {
00509         mamb = _ScalarProduct2D( theList[i], theSeg->extremity1,
00510                                  theList[i], theSeg->extremity2 );
00511         if ( mamb > maxThreshold ) {
00512           index ++;
00513           _SwapPoints( theList, index, i );
00514           continue;
00515         }
00516         if ( mamb <= minThreshold ) {
00517           _SwapPoints( theList, i, l );
00518           i--;   l--;
00519           continue;
00520         }
00521         if ( b < mamb ) b = mamb;
00522       }
00523       *last = l;
00524       *bound = b;
00525       return( index );
00526     }
00527 
00528     if ( dim == 3 ) {
00529       for ( i=first; i<=l; i++ ) {
00530         mamb = _ScalarProduct3D( theList[i], theSeg->extremity1,
00531                                  theList[i], theSeg->extremity2 );
00532         if ( mamb > maxThreshold ) {
00533           index ++;
00534           _SwapPoints( theList, index, i );
00535           continue;
00536         }
00537         if ( mamb <= minThreshold ) {
00538           _SwapPoints( theList, i, l );
00539           i--;   l--;
00540           continue;
00541         }
00542         if ( b < mamb ) b = mamb;
00543       }
00544       *last = l;
00545       *bound = b;
00546       return( index );
00547     }
00548 
00549     for ( i=first; i<=l; i++ ) {
00550       mamb = _ScalarProduct( theList[i], theSeg->extremity1,
00551                              theList[i], theSeg->extremity2, dim );
00552       if ( mamb > maxThreshold ) {
00553         index ++;
00554         _SwapPoints( theList, index, i );
00555         continue;
00556       }
00557       if ( mamb <= minThreshold ) {
00558         _SwapPoints( theList, i, l );
00559         i--;   l--;
00560         continue;
00561       }
00562       if ( b < mamb ) b = mamb;
00563     }
00564     *last = l;
00565     *bound = b;
00566     return( index );
00567 
00568     /* END
00569        REDUCTION INSIDE THE SMALLEST SPHERE
00570      */
00571     break;
00572 
00573   case 2 :
00574 
00575     /* COMPLETE REDUCTION
00576 
00577        On suppose implicitement
00578        que l'ensemble des points est
00579        compris dans l'intersection des spheres centrees
00580        sur A et B et de rayon |AB|.
00581        En effet, les conditions (pour un rayon dans
00582        l'intervalle [min,med]) pour l'elimination ne le
00583        verifient pas.
00584      */
00585     if ( dim == 2 ) {
00586 
00587       if ( squareDiameter <= theSeg->squareDiameter ) {
00588 
00589         for ( i=first; i<=l; i++ ) {
00590           mamb = _ScalarProduct2D( theList[i], theSeg->extremity1,
00591                                    theList[i], theSeg->extremity2 );
00592           if ( mamb > maxThreshold ) {
00593             index ++;
00594             _SwapPoints( theList, index, i );
00595             continue;
00596           }
00597           if (  mamb > minThreshold ) {
00598             am2 = _SquareDistance2D( theList[i], theSeg->extremity1 );
00599             if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) ) - mamb*mamb < 0 ) {
00600               _SwapPoints( theList, i, l );
00601               i--;   l--;
00602               continue;
00603             }
00604             if ( b < mamb ) b = mamb;
00605             continue;
00606           }
00607           /* if ( mamb <= minThreshold ) */
00608           _SwapPoints( theList, i, l );
00609           i--;   l--;
00610         }
00611 
00612       } else {
00613 
00614         for ( i=first; i<=l; i++ ) {
00615           mamb = _ScalarProduct2D( theList[i], theSeg->extremity1,
00616                                    theList[i], theSeg->extremity2  );
00617           if ( mamb > maxThreshold ) {
00618             index ++;
00619             _SwapPoints( theList, index, i );
00620             continue;
00621           }
00622           if ( mamb > medThreshold ) {
00623             if ( b < mamb ) b = mamb;
00624             continue;
00625           }
00626           if ( mamb > minThreshold ) {
00627             am2 = _SquareDistance2D( theList[i], theSeg->extremity1 );
00628             if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) )
00629                  - (R2 - ab2 - mamb)*(R2 - ab2 - mamb) < 0 ) {
00630               _SwapPoints( theList, i, l );
00631               i--;   l--;
00632               continue;
00633             }
00634             if ( b < mamb ) b = mamb;
00635             continue;
00636           }
00637           /* if ( mamb <= minThreshold ) */
00638           _SwapPoints( theList, i, l );
00639           i--;   l--;
00640         }
00641 
00642       }
00643       *last = l;
00644       *bound = b;
00645       return( index );
00646     }
00647 
00648     if ( dim == 3 ) {
00649 
00650       if ( squareDiameter <= theSeg->squareDiameter ) {
00651 
00652         for ( i=first; i<=l; i++ ) {
00653           mamb = _ScalarProduct3D( theList[i], theSeg->extremity1,
00654                                    theList[i], theSeg->extremity2 );
00655           if ( mamb > maxThreshold ) {
00656             index ++;
00657             _SwapPoints( theList, index, i );
00658             continue;
00659           }
00660           if (  mamb > minThreshold ) {
00661             am2 = _SquareDistance3D( theList[i], theSeg->extremity1 );
00662             if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) ) - mamb*mamb < 0 ) {
00663               _SwapPoints( theList, i, l );
00664               i--;   l--;
00665               continue;
00666             }
00667             if ( b < mamb ) b = mamb;
00668             continue;
00669           }
00670           /* if ( mamb <= minThreshold ) */
00671           _SwapPoints( theList, i, l );
00672           i--;   l--;
00673         }
00674 
00675       } else {
00676 
00677         for ( i=first; i<=l; i++ ) {
00678           mamb = _ScalarProduct3D( theList[i], theSeg->extremity1,
00679                                    theList[i], theSeg->extremity2  );
00680           if ( mamb > maxThreshold ) {
00681             index ++;
00682             _SwapPoints( theList, index, i );
00683             continue;
00684           }
00685           if ( mamb > medThreshold ) {
00686             if ( b < mamb ) b = mamb;
00687             continue;
00688           }
00689           if ( mamb > minThreshold ) {
00690             am2 = _SquareDistance3D( theList[i], theSeg->extremity1 );
00691             if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) )
00692                  - (R2 - ab2 - mamb)*(R2 - ab2 - mamb) < 0 ) {
00693               _SwapPoints( theList, i, l );
00694               i--;   l--;
00695               continue;
00696             }
00697             if ( b < mamb ) b = mamb;
00698             continue;
00699           }
00700           /* if ( mamb <= minThreshold ) */
00701           _SwapPoints( theList, i, l );
00702           i--;   l--;
00703         }
00704 
00705       }
00706       *last = l;
00707       *bound = b;
00708       return( index );
00709     }
00710 
00711     if ( squareDiameter <= theSeg->squareDiameter ) {
00712 
00713       for ( i=first; i<=l; i++ ) {
00714         mamb = _ScalarProduct( theList[i], theSeg->extremity1,
00715                                theList[i], theSeg->extremity2, dim );
00716         if ( mamb > maxThreshold ) {
00717           index ++;
00718           _SwapPoints( theList, index, i );
00719           continue;
00720         }
00721         if (  mamb > minThreshold ) {
00722           am2 = _SquareDistance( theList[i], theSeg->extremity1, dim );
00723           if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) ) - mamb*mamb < 0 ) {
00724             _SwapPoints( theList, i, l );
00725             i--;   l--;
00726             continue;
00727           }
00728           if ( b < mamb ) b = mamb;
00729           continue;
00730         }
00731         /* if ( mamb <= minThreshold ) */
00732         _SwapPoints( theList, i, l );
00733         i--;   l--;
00734       }
00735 
00736     } else {
00737 
00738       for ( i=first; i<=l; i++ ) {
00739         mamb = _ScalarProduct( theList[i], theSeg->extremity1,
00740                                theList[i], theSeg->extremity2, dim );
00741         if ( mamb > maxThreshold ) {
00742           index ++;
00743           _SwapPoints( theList, index, i );
00744           continue;
00745         }
00746         if ( mamb > medThreshold ) {
00747           if ( b < mamb ) b = mamb;
00748           continue;
00749         }
00750         if ( mamb > minThreshold ) {
00751           am2 = _SquareDistance( theList[i], theSeg->extremity1, dim );
00752           if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) )
00753                - (R2 - ab2 - mamb)*(R2 - ab2 - mamb) < 0 ) {
00754             _SwapPoints( theList, i, l );
00755             i--;   l--;
00756             continue;
00757           }
00758           if ( b < mamb ) b = mamb;
00759           continue;
00760         }
00761         /* if ( mamb <= minThreshold ) */
00762         _SwapPoints( theList, i, l );
00763         i--;   l--;
00764       }
00765 
00766     }
00767     *last = l;
00768     *bound = b;
00769     return( index );
00770 
00771     /* END
00772        COMPLETE REDUCTION
00773     */
00774   }
00775   return( -1 );
00776 }
00777 
00778 //void _CountPointsInSpheres( TypeSegment *theSeg,
00779 //                          const double squareDiameter,
00780 //                          const double **theList,
00781 //                          const int first,
00782 //                          const int last,
00783 //                          const int dim )
00784 //{
00785 //  double minThreshold;
00786 //  double medThreshold;
00787 //  double maxThreshold;
00788 //
00789 //  double ab  = sqrt( theSeg->squareDiameter );
00790 //  double ab2 = theSeg->squareDiameter;
00791 //
00792 //  double R  = sqrt( squareDiameter );
00793 //  double R2 = squareDiameter;
00794 //
00795 //  int n[5] = {0,0,0,0,0};
00796 //  int i;
00797 //
00798 //  double mamb, am2;
00799 //
00800 //  minThreshold = (R - .86602540378443864676*ab)
00801 //    * (R - .86602540378443864676*ab)
00802 //    - 0.25 * ab2;
00803 //
00804 //  medThreshold = 0.5 * ab2 + 0.25 * R2
00805 //    - .43301270189221932338 * R * sqrt( 4 * ab2 - R2 );
00806 //
00807 //  maxThreshold = 0.25 * (R2 - ab2);
00808 //
00809 //
00810 //  for ( i=first; i<=last; i++ ) {
00811 //
00812 //    mamb = _ScalarProduct( theList[i], theSeg->extremity1,
00813 //                         theList[i], theSeg->extremity2, dim );
00814 //
00815 //    if ( mamb > maxThreshold ) {
00816 //      n[0] ++ ;
00817 //      continue;
00818 //    }
00819 //
00820 //    if ( mamb > medThreshold ) {
00821 //      n[1] ++ ;
00822 //      continue;
00823 //    }
00824 //
00825 //    if ( mamb > minThreshold ) {
00826 //      am2 = _SquareDistance( theList[i], theSeg->extremity1, dim );
00827 //      if ( 3.0 * ( am2 * ab2 - (am2 - mamb)*(am2 - mamb) )
00828 //         - (R2 - ab2 - mamb)*(R2 - ab2 - mamb) < 0 ) {
00829 //      n[2] ++;
00830 //      continue;
00831 //      }
00832 //      n[3] ++;
00833 //      continue;
00834 //    }
00835 //
00836 //    n[4] ++;
00837 //  }
00838 //
00839 //  printf(" diametre courant = %g  -  double normale = %g\n",
00840 //       R, ab );
00841 //  printf(" %8d points dont\n", last-first+1 );
00842 //  printf(" - %6d : candidats extremites       R=%g\n", n[0], sqrt(maxThreshold+0.25*ab2) );
00843 //  printf(" - %6d : rien a faire               R=%g\n", n[1], sqrt(medThreshold+0.25*ab2)  );
00844 //  printf(" - %6d : a tester                   \n", n[2]+n[3] );
00845 //  printf("   + %6d : a eliminer\n", n[2] );
00846 //  printf("   + %6d : a conserver\n", n[3] );
00847 //  printf(" - %6d : elimines directement       R=%g\n", n[4], sqrt(minThreshold+0.25*ab2)  );
00848 //  printf("----------\n" );
00849 //  printf(" %8d\n", n[0]+n[1]+n[2]+n[3]+n[4] );
00850 //
00851 //}
00852 
00853 
00854 /* Find the farthest point from a sphere
00855 
00856    returned value :
00857    - its index if there are some points outside the sphere
00858    - else #(first index) - 1
00859 
00860    As the sphere diameter is an approximation of the set diameter,
00861    points verifying MA.MB <= (1.5 -sqrt(3)) AB^2
00862    can be removed
00863 
00864 */
00865 
00866 
00867 int _FarthestPointFromSphere( TypeSegment *theSeg,
00868                               const double **theList,
00869                               const int first,
00870                               int *last,
00871                               const int dim,
00872                               const int _reduction_mode_ )
00873 {
00874   int i, l = (*last);
00875   int index = first - 1;
00876   double diff, maxdiff = 0.0;
00877   /* threshold = 1.5 - sqrt(3)
00878    */
00879   double threshold = -0.23205080756887729352 * theSeg->squareDiameter;
00880 
00881   if ( l < first ) return( index );
00882 
00883 
00884   switch ( _reduction_mode_ ) {
00885 
00886   case 0 :
00887 
00888     /* NO REDUCTION CASE
00889      */
00890 
00891     if ( dim == 2 ) {
00892       for ( i=first; i<=l; i++ ) {
00893         diff = _ScalarProduct2D( theList[i], theSeg->extremity1,
00894                                  theList[i], theSeg->extremity2 );
00895         if ( maxdiff < diff ) {
00896           index   = i;
00897           maxdiff = diff;
00898         }
00899       }
00900       return( index );
00901     }
00902 
00903     if ( dim == 3 ) {
00904       for ( i=first; i<=l; i++ ) {
00905         diff = _ScalarProduct3D( theList[i], theSeg->extremity1,
00906                                  theList[i], theSeg->extremity2 );
00907         if ( maxdiff < diff ) {
00908           index   = i;
00909           maxdiff = diff;
00910         }
00911       }
00912     return( index );
00913     }
00914 
00915     for ( i=first; i<=l; i++ ) {
00916       diff = _ScalarProduct( theList[i], theSeg->extremity1,
00917                              theList[i], theSeg->extremity2, dim );
00918       if ( maxdiff < diff ) {
00919         index   = i;
00920         maxdiff = diff;
00921       }
00922     }
00923     return( index );
00924 
00925     /* END
00926        NO REDUCTION CASE
00927     */
00928     break;
00929 
00930   default :
00931   case 1 :
00932 
00933     /* REDUCTION INSIDE THE SMALLEST SPHERE
00934      */
00935 
00936     /* AB = diameter extremities
00937        MA.MB = (MC+CA).(MC+CB) = MC^2 + CA.CB + MC ( CB+CA )
00938        = MC^2 - R^2   + 0
00939     */
00940     if ( dim == 2 ) {
00941       for ( i=first; i<=l; i++ ) {
00942         diff = _ScalarProduct2D( theList[i], theSeg->extremity1,
00943                                  theList[i], theSeg->extremity2 );
00944         if ( diff > maxdiff ) {
00945           index   = i;
00946           maxdiff = diff;
00947         }
00948         else if ( diff <= threshold ) {
00949           _SwapPoints( theList, i, l );
00950           i --;   l --;
00951           continue;
00952         }
00953       }
00954       *last = l;
00955       return( index );
00956     }
00957 
00958 
00959     if ( dim == 3 ) {
00960       for ( i=first; i<=l; i++ ) {
00961         diff = _ScalarProduct3D( theList[i], theSeg->extremity1,
00962                                  theList[i], theSeg->extremity2 );
00963         if ( maxdiff < diff ) {
00964           index   = i;
00965           maxdiff = diff;
00966         }
00967         else if ( diff <= threshold ) {
00968           _SwapPoints( theList, i, l );
00969           i --;   l --;
00970           continue;
00971         }
00972       }
00973       *last = l;
00974       return( index );
00975     }
00976 
00977 
00978     for ( i=first; i<=l; i++ ) {
00979       diff = _ScalarProduct( theList[i], theSeg->extremity1,
00980                              theList[i], theSeg->extremity2, dim );
00981       if ( maxdiff < diff ) {
00982         index   = i;
00983         maxdiff = diff;
00984       }
00985       else if ( diff <= threshold ) {
00986         _SwapPoints( theList, i, l );
00987         i --;   l --;
00988         continue;
00989       }
00990     }
00991     *last = l;
00992     return( index );
00993     /* END
00994        REDUCTION INSIDE THE SMALLEST SPHERE
00995      */
00996   }
00997 
00998   return( -1 );
00999 }
01000 
01001 double _MaximalSegmentInTwoLists( TypeSegment *theSeg,
01002                                   const int index1,
01003                                   const double **theList1,
01004                                   int *first1,
01005                                   int *last1,
01006                                   const double **theList2,
01007                                   int *first2,
01008                                   int *last2,
01009                                   int dim )
01010 {
01011 
01012   int f1=*first1;
01013   int l1=*last1;
01014 
01015   int i1 = index1;
01016 
01017   int f2=*first2;
01018   int l2=*last2;
01019 
01020   int i2;
01021 
01022   const double *ref;
01023 
01024   double d, dprevious;
01025 
01026 
01027   theSeg->extremity1 = (double*)NULL;
01028   theSeg->extremity2 = (double*)NULL;
01029   theSeg->squareDiameter = std::numeric_limits<double>::lowest();
01030 
01031 
01032   if ( *first1 < 0 || *last1 < 0 ) return( -1.0 );
01033   if ( *first1 > *last1 ) return( 0.0 );
01034   if ( *first2 < 0 || *last2 < 0 ) return( -1.0 );
01035   if ( *first2 > *last2 ) return( 0.0 );
01036   if ( *first2 > *last2 ) {
01037     l2 = *first2;
01038     f2 = *last2;
01039   }
01040   if ( index1 < f1 || index1 > l1 ) return( -1.0 );
01041 
01042 
01043   do {
01044 
01045     dprevious = theSeg->squareDiameter;
01046 
01047     /* reference point in list #1
01048      */
01049     ref = theList1[i1];
01050     /* point #i1 will be compared against all other points
01051        in list #2
01052        reject it at the beginning of the list
01053        do not consider it in the future
01054     */
01055     _SwapPoints( theList1, f1, i1 );
01056     f1 ++;
01057 
01058     /* find the furthest point from 'ref'
01059      */
01060     d = _MaximalDistanceFromPoint( &i2, ref, theList2, f2, l2, dim );
01061 
01062     /* better estimate
01063      */
01064     if ( d > theSeg->squareDiameter ) {
01065 
01066       theSeg->extremity1 = ref;
01067       theSeg->extremity2 = theList2[i2];
01068       theSeg->squareDiameter = d;
01069 
01070       if ( f1 <= l1 ) {
01071 
01072         dprevious = theSeg->squareDiameter;
01073 
01074         /* reference point in list #1
01075          */
01076         ref = theList2[i2];
01077         /* point #i2 will be compared against all other points
01078            in list #1
01079            reject it at the beginning of the list
01080            do not consider it in the future
01081         */
01082         _SwapPoints( theList2, f2, i2 );
01083         f2 ++;
01084 
01085         /* find the furthest point from 'ref'
01086          */
01087         d = _MaximalDistanceFromPoint( &i1, ref, theList1, f1, l1, dim );
01088 
01089         /* better estimate
01090          */
01091         if ( d > theSeg->squareDiameter ) {
01092 
01093           theSeg->extremity1 = theList1[i1];
01094           theSeg->extremity2 = ref;
01095           theSeg->squareDiameter = d;
01096 
01097         }
01098 
01099       }
01100 
01101     }
01102 
01103 
01104   } while( theSeg->squareDiameter > dprevious && f1 <= l1 && f2 <= l2 );
01105 
01106 
01107   *first1 = f1;
01108   *last1  = l1;
01109   *first2 = f2;
01110   *last2  = l2;
01111   return( theSeg->squareDiameter );
01112 
01113 }
01114 
01115 double _MaximalSegmentInOneList( TypeSegment *theSeg,
01116                                  const int index,
01117                                  const double **theList,
01118                                  int *first,
01119                                  int *last,
01120                                  const int dim )
01121 {
01122   int f=*first;
01123   int l=*last;
01124 
01125   int i = index;
01126   const double *ref;
01127 
01128   double d, dprevious;
01129 
01130   theSeg->extremity1 = (double*)NULL;
01131   theSeg->extremity2 = (double*)NULL;
01132   theSeg->squareDiameter = std::numeric_limits<double>::lowest();
01133 
01134 
01135 
01136   if ( *first < 0 || *last < 0 ) return( -1.0 );
01137   if ( *first > *last ) return( 0.0 );
01138   if ( index < f || index > l ) return( -1.0 );
01139   if ( f == l ) {
01140     theSeg->extremity1 = theList[i];
01141     theSeg->extremity2 = theList[i];
01142     return( 0.0 );
01143   }
01144 
01145 
01146 
01147   do {
01148 
01149     dprevious = theSeg->squareDiameter;
01150 
01151 
01152     ref = theList[i];
01153     /* point #i will be compared against all other points
01154        reject it at the beginning of the list
01155        do not consider it in the future
01156     */
01157     _SwapPoints( theList, f, i );
01158     f ++;
01159 
01160 
01161     /* find the furthest point from 'ref'
01162      */
01163     d = _MaximalDistanceFromPoint( &i, ref, theList, f, l, dim );
01164     if ( d > theSeg->squareDiameter ) {
01165       theSeg->extremity1 = ref;
01166       theSeg->extremity2 = theList[i];
01167       theSeg->squareDiameter = d;
01168     }
01169 
01170 
01171   } while ( theSeg->squareDiameter > dprevious && f <= l );
01172 
01173   *first = f;
01174   *last  = l;
01175   return( theSeg->squareDiameter );
01176 }
01177 
01178 double _MaximalDistanceFromPoint( int *index,
01179                                   const double *ref,
01180                                   const double **theList,
01181                                   const int first,
01182                                   const int last,
01183                                   const int dim )
01184 {
01185   int f=first;
01186   int l=last;
01187   int i;
01188   double dmax, d;
01189 
01190   *index = -1;
01191   if ( first < 0 || last < 0 ) return( -1.0 );
01192   if ( first > last ) return( 0.0 );
01193 
01194 
01195 
01196   if ( dim == 2 ) {
01197 
01198     dmax   = _SquareDistance2D( theList[f], ref );
01199     *index = f;
01200 
01201     for ( i=f+1; i<=l; i++ ) {
01202       d = _SquareDistance2D( theList[i], ref );
01203       if ( d > dmax ) {
01204         dmax   = d;
01205         *index = i;
01206       }
01207     }
01208     return( dmax );
01209 
01210   }
01211 
01212   if ( dim == 3 ) {
01213 
01214     dmax   = _SquareDistance3D( theList[f], ref );
01215     *index = f;
01216 
01217     for ( i=f+1; i<=l; i++ ) {
01218       d = _SquareDistance3D( theList[i], ref );
01219       if ( d > dmax ) {
01220         dmax   = d;
01221         *index = i;
01222       }
01223     }
01224     return( dmax );
01225 
01226   }
01227 
01228 
01229 
01230 
01231   dmax   = _SquareDistance( theList[f], ref, dim );
01232   *index = f;
01233 
01234   for ( i=f+1; i<=l; i++ ) {
01235     d = _SquareDistance( theList[i], ref, dim );
01236     if ( d > dmax ) {
01237         dmax   = d;
01238         *index = i;
01239     }
01240   }
01241   return( dmax );
01242 
01243 }
01244 
01245 
01246 /* Swap two points
01247  */
01248 
01249 void _SwapPoints( const double **theList, const int i, const int j )
01250 {
01251   const double *tmp;
01252   tmp = theList[i];
01253   theList[i] = theList[j];
01254   theList[j] = tmp;
01255 }
01256 
01257 static int _base_ =  100000000;
01258 
01259 void _InitCounter( TypeCounter *c )
01260 {
01261   c->c1 = c->c2 = 0;
01262 }
01263 void _AddToCounter( TypeCounter *c, const int i )
01264 {
01265   c->c2 += i;
01266   while ( c->c2 >= _base_ ) {
01267     c->c2 -= _base_;
01268     c->c1 ++;
01269   }
01270   while ( c->c2 < 0 ) {
01271     c->c2 += _base_;
01272     c->c1 --;
01273   }
01274 }
01275 double _GetCounterAverage( TypeCounter *c, const int i )
01276 {
01277   return( (_base_ / (double)i ) * c->c1 + c->c2 / (double)i );
01278 }
01279 
01280 
01281 #ifdef _STATS_
01282 TypeCounter scalarProducts;
01283 
01284 void _InitScalarProductCounter()
01285 {
01286   _InitCounter( &scalarProducts );
01287 }
01288 static void _IncScalarProductCounter()
01289 {
01290   _AddToCounter( &scalarProducts, 1 );
01291 }
01292 double _GetScalarProductAverage( int n )
01293 {
01294   return( _GetCounterAverage( &scalarProducts, n ) );
01295 }
01296 #endif
01297 
01298 
01299 
01300 
01301 /* square distance
01302  */
01303 double _SquareDistance( const double *a, const double *b, const int dim )
01304 {
01305   register int i;
01306   register double d = 0.0;
01307   register double ba;
01308   for ( i=0; i<dim; i++ ) {
01309     ba = b[i] - a[i];
01310     d += ba * ba;
01311   }
01312 
01313 #ifdef _STATS_
01314   _IncScalarProductCounter();
01315 #endif
01316 
01317   return( d );
01318 }
01319 
01320 
01321 
01322 
01323 
01324 double _SquareDistance3D( const double *a, const double *b )
01325 {
01326 
01327 #ifdef _STATS_
01328   _IncScalarProductCounter();
01329 #endif
01330 
01331   return( (b[0]-a[0])*(b[0]-a[0]) +
01332           (b[1]-a[1])*(b[1]-a[1]) +
01333           (b[2]-a[2])*(b[2]-a[2]) );
01334 }
01335 
01336 
01337 double _SquareDistance2D( const double *a, const double *b )
01338 {
01339 
01340 #ifdef _STATS_
01341   _IncScalarProductCounter();
01342 #endif
01343 
01344   return( (b[0]-a[0])*(b[0]-a[0]) +
01345           (b[1]-a[1])*(b[1]-a[1]) );
01346 }
01347 
01348 
01349 
01350 
01351 
01352 /* dot product
01353  */
01354 double _ScalarProduct( const double *a, const double *b,
01355                        const double *c, const double *d, const int dim )
01356 {
01357   register int i;
01358   register double scalar = 0.0;
01359   register double ab, cd;
01360   for ( i=0; i<dim; i++ ) {
01361     ab = b[i] - a[i];
01362     cd = d[i] - c[i];
01363     scalar += ab * cd;
01364   }
01365 
01366 #ifdef _STATS_
01367   _IncScalarProductCounter();
01368 #endif
01369 
01370   return( scalar );
01371 }
01372 
01373 
01374 
01375 double _ScalarProduct3D( const double *a, const double *b,
01376                          const double *c, const double *d )
01377 {
01378 
01379 #ifdef _STATS_
01380   _IncScalarProductCounter();
01381 #endif
01382 
01383   return( (b[0]-a[0])*(d[0]-c[0]) +
01384           (b[1]-a[1])*(d[1]-c[1]) +
01385           (b[2]-a[2])*(d[2]-c[2]) );
01386 }
01387 
01388 double _ScalarProduct2D( const double *a, const double *b,
01389                          const double *c, const double *d )
01390 {
01391 
01392 #ifdef _STATS_
01393   _IncScalarProductCounter();
01394 #endif
01395 
01396   return( (b[0]-a[0])*(d[0]-c[0]) +
01397           (b[1]-a[1])*(d[1]-c[1]) );
01398 }
01399 
01400 int _FindPointInList( const double **theList,
01401                       const int first,
01402                       const int last,
01403                       double x0,
01404                       double x1 )
01405 {
01406   int i, j=-1;
01407   double e=1e-5;
01408   for ( i=first; i<=last; i++ ) {
01409     if ( theList[i][0]-e < x0 && theList[i][0]+e > x0 &&
01410          theList[i][1]-e < x1 && theList[i][1]+e > x1 ) {
01411       if ( j == -1 ) {
01412         j = i;
01413       } else {
01414         fprintf( stderr, "found again at #%d\n", i );
01415       }
01416     }
01417   }
01418   return( j );
01419 }
01420 
01421 
01422 
01423 
01424 
01425 
01426 
01427 
01428 
01429 
01430 
01431 
01432 
01433 
01434 
01435 
01436 
01437 
01438 
01439 
01440 
01441 
01442 double _QuadraticDiameterInOneList( TypeSegment *theDiam,
01443                                     const double **theList,
01444                                     const int first,
01445                                     const int last,
01446                                     const int dim )
01447 {
01448   int i, j;
01449   double d;
01450 
01451   theDiam->extremity1 = (double*)NULL;
01452   theDiam->extremity2 = (double*)NULL;
01453   theDiam->squareDiameter = std::numeric_limits<double>::lowest();
01454 
01455   d = 0.0;
01456 
01457   if ( dim == 2 ) {
01458 
01459     for ( i=first; i<=last-1; i++ )
01460     for ( j=i+1;   j<=last;   j++ ) {
01461       d = _SquareDistance2D( theList[i], theList[j] );
01462       if ( d > theDiam->squareDiameter ) {
01463         theDiam->squareDiameter = d;
01464         theDiam->extremity1 = theList[i];
01465         theDiam->extremity2 = theList[j];
01466       }
01467     }
01468     return( theDiam->squareDiameter );
01469 
01470   }
01471 
01472 
01473 
01474   if ( dim == 3 ) {
01475 
01476     for ( i=first; i<=last-1; i++ )
01477     for ( j=i+1;   j<=last;   j++ ) {
01478       d = _SquareDistance3D( theList[i], theList[j] );
01479       if ( d > theDiam->squareDiameter ) {
01480         theDiam->squareDiameter = d;
01481         theDiam->extremity1 = theList[i];
01482         theDiam->extremity2 = theList[j];
01483       }
01484     }
01485     return( theDiam->squareDiameter );
01486 
01487   }
01488 
01489 
01490 
01491   for ( i=first; i<=last-1; i++ )
01492   for ( j=i+1;   j<=last;   j++ ) {
01493     d = _SquareDistance( theList[i], theList[j], dim );
01494     if ( d > theDiam->squareDiameter ) {
01495       theDiam->squareDiameter = d;
01496       theDiam->extremity1 = theList[i];
01497       theDiam->extremity2 = theList[j];
01498     }
01499   }
01500 
01501   return( theDiam->squareDiameter );
01502 }
01503 
01504 
01505 
01506 
01507 
01508 
01509 
01510 
01511 double _QuadraticDiameterInTwoLists( TypeSegment *theDiam,
01512                                      int   *index1,
01513                                      int   *index2,
01514                                      const double **theList1,
01515                                      const int first1,
01516                                      const int last1,
01517                                      const double **theList2,
01518                                      const int first2,
01519                                      const int last2,
01520                                      const int dim )
01521 {
01522   int i, j;
01523   double d;
01524 
01525   /*
01526   theDiam->extremity1 = (double*)NULL;
01527   theDiam->extremity2 = (double*)NULL;
01528   theDiam->squareDiameter = std::numeric_limits<double>::lowest();
01529   */
01530 
01531   if ( index1 != NULL && index2 != NULL ) {
01532     *index1 = first1-1;
01533     *index2 = first2-1;
01534   }
01535 
01536   d = 0.0;
01537 
01538   if ( dim == 2 ) {
01539 
01540     for ( i=first1; i<=last1; i++ )
01541     for ( j=first2; j<=last2; j++ ) {
01542       d = _SquareDistance2D( theList1[i], theList2[j] );
01543       if ( d > theDiam->squareDiameter ) {
01544         theDiam->squareDiameter = d;
01545         theDiam->extremity1 = theList1[i];
01546         theDiam->extremity2 = theList2[j];
01547         if ( index1 != NULL && index2 != NULL ) {
01548           *index1 = i;
01549           *index2 = j;
01550         }
01551       }
01552     }
01553     return( theDiam->squareDiameter );
01554 
01555   }
01556 
01557 
01558 
01559   if ( dim == 3 ) {
01560 
01561     for ( i=first1; i<=last1; i++ )
01562     for ( j=first2; j<=last2; j++ ) {
01563       d = _SquareDistance3D( theList1[i], theList2[j] );
01564       if ( d > theDiam->squareDiameter ) {
01565         theDiam->squareDiameter = d;
01566         theDiam->extremity1 = theList1[i];
01567         theDiam->extremity2 = theList2[j];
01568         if ( index1 != NULL && index2 != NULL ) {
01569           *index1 = i;
01570           *index2 = j;
01571         }
01572       }
01573     }
01574     return( theDiam->squareDiameter );
01575 
01576   }
01577 
01578 
01579 
01580   for ( i=first1; i<=last1; i++ )
01581   for ( j=first2; j<=last2; j++ ) {
01582     d = _SquareDistance( theList1[i], theList2[j], dim );
01583     if ( d > theDiam->squareDiameter ) {
01584       theDiam->squareDiameter = d;
01585       theDiam->extremity1 = theList1[i];
01586       theDiam->extremity2 = theList2[j];
01587       if ( index1 != NULL && index2 != NULL ) {
01588         *index1 = i;
01589         *index2 = j;
01590       }
01591     }
01592   }
01593 
01594   return( theDiam->squareDiameter );
01595 }
01596 
01597 }
01598 }