123

#define sqrt5 2.236067977499789696

static int Stopping_Rule(double x0, double x1, double tolerance)
{
   double xm = 0.5 * fabs( x1 + x0 );

   if ( xm <= 1.0 ) return ( fabs( x1 - x0 ) < tolerance ) ? 1 : 0;
   return ( fabs( x1 - x0 ) < tolerance * xm ) ? 1 : 0;
}
void ScanViewer::Min_Search_Golden_Section
(
        qreal* a, qreal *fa , /// угол и значение площади для минимальноо угла + результат работы    ( 0 радиан    )
        qreal* b, qreal* fb , /// угол и значение площади для максимального угла + результат работы  ( Pi/2 радиан )
        qreal tolerance, QVector < int > x , QVector < int > y   /// Точность + массивы данных от линеек и энкодера (  vector x ( encoder in mm ) , vector y ( beams in mm ) )
        /// in case of light grid vector < int > y => QVector <QPair < int , int > limit;
)
{


    /// auto metrics = [this,x,y](qreal A)->qreal  ->  функция , принимающая на вход вектора x , y скрытно и угол поворота явно
    /// результат - double функционал - площадь

   auto metrics = [this,x,y](qreal A)->qreal   /// После поворота системы координат на А радиан найдем функционал
                                               /// Функционал - площадь или периметр построенный на [ Xmax-Xmin , Ymax - Ymin ]
    {
        //// Выделим переменные под предельные значения координат
        qreal xmax  , xmin;
        qreal ymax  , ymin;
        QPoint AS = rotate(QPoint(x[0],y[0]),A);  /// повернем первую точку из множества чтобы инициализировать максимумы и минимумы

        xmax = xmin = AS.x();
        ymax = ymin = AS.y();

        /// пройдемся по всем точкам для определения реальных максимумов и минимумов
            for ( int m = 0 ; m < x.size() ; m++ )
            {
                QPoint AS = rotate(QPoint(x[m],y[m]),A);
                if (AS.x()>xmax) xmax = AS.x();
                if (AS.x()<xmin) xmin = AS.x();
                if (AS.y()>ymax) ymax = AS.y();
                if (AS.y()<ymin) ymin = AS.y();
                // #form first optimization factor
            }

    return fabs(xmax-xmin)*(ymax-ymin); /// функционал - площадь
    };

   static const double lambda = 0.5 * (sqrt5 - 1.0);     /// Константа формулы золотого сечения
   static const double mu = 0.5 * (3.0 - sqrt5);         /// Константа формулы золотого сечения

   double x1;  /// внутренняя левая точка поискового интервала ( угол )
   double x2;  /// внутренняя правая точка поискового интервала ( угол )
   double fx1; /// внутренняя правая точка поискового интервала ( площадь или периметр )
   double fx2; /// внутренняя правая точка поискового интервала ( плозадь или периметр )

   /// Поиск угла и плозади для внутренних точек отрезка

   /// x1 делит отрезок [A,B] в золотом сечении
   /// x2 делит отрезок [B,A] в золотом сечении

   x1 = *b - lambda * (*b - *a);
   x2 = *a + lambda * (*b - *a);

   /// Определим функционал внутренних точек ( площадь и периметр )

   fx1 = metrics(x1);
   fx2 = metrics(x2);

   while (  ! Stopping_Rule( *a, *b, tolerance) ) {   /// Пока не достигнута заданная точность


      if (fx1 > fx2) { /// если площадь прямоугольника со сторонами \\ осям в точке 1 > 2
         *a = x1;      /// сужаем поисковый интервал слева до точки X1 + захват результата на выход
         *fa = fx1;    /// пересчет левого значения площади не нужен - он известен и так
         if ( Stopping_Rule( *a, *b, tolerance) ) break;   /// Если поисковый интервал стал подходить по точности выходим
         x1 = x2;       /// левая внутренняя точка стала правой
         fx1 = fx2;     /// соответствующая запись значения площади
         x2 = *b - mu * (*b - *a);   /// пересчет правой точки ( ее значение и плозадь заранее неизвестны )
         fx2 = metrics(x2);
      }


      else {          /// если площадь прямоугольника со сторонами \\ осям в точке 1 < 2
         *b = x2;       /// сужаем поисковый интервал справа до точки X2 + захват результата на выход
         *fb = fx2;     /// пересчет правого значения площади не нужен - он известен и так
         if ( Stopping_Rule( *a, *b, tolerance) ) break;   /// Если поисковый интервал стал подходить по точности выходим
         x2 = x1;       /// правая внутренняя точка стала левой
         fx2 = fx1;     /// соответствующая запись значения площади
         x1 = *a + mu * (*b - *a); /// пересчет левой точки ( ее значение и плозадь заранее неизвестны )
         fx1 = metrics(x1);
      }
   }



   return;
}











/// Вектора с данными + переменные под размеры и угол

void ScanViewer::getSizes(QVector <int> x , QVector < int > y , QSize * return_size  ,qreal *Angular)
{

    /// Аналог функции из Min_Search_Golden_Section
    auto metrics = [this,x,y](qreal A)->qreal
    {
        qreal xmax,xmin;
        qreal ymax , ymin;
        QPoint AS = rotate(QPoint(x[0],y[0]),A);
        xmax = xmin = AS.x();
        ymax = ymin = AS.y();

            for ( int m = 0 ; m < x.size() ; m++ )
            {
                QPoint AS = rotate(QPoint(x[m],y[m]),A);
                if (AS.x()>xmax) xmax = AS.x();
                if (AS.x()<xmin) xmin = AS.x();
                if (AS.y()>ymax) ymax = AS.y();
                if (AS.y()<ymin) ymin = AS.y();
                // #form first optimization factor
            }

    return fabs(xmax-xmin)*(ymax-ymin);
    };

    /// начальный и конечный угол , начальное и конечное значение функционала ( см . Min_Search_Golden_Section )
    /// 0.005 радиана - точность поиска угла
    /// x , y -> данные с линеек

    qreal eps = 0.005;
    qreal a = 0-3*eps; /// Расширенный интервал из-за особенностей метода поиска
    qreal b = M_PI/2.0+3*eps;
    qreal fa = metrics(a);
    qreal fb = metrics(b);

    Min_Search_Golden_Section(&a,&fa,&b,&fb,0.005,x,y);

    /// Получив значение угла , минимизирующего площадь , осталось определить размеры

    /// выделим переменные под максимумы и минимумы и инициализируем их первой точкой после поворота на угол

    qreal xmax,xmin;
    qreal ymax , ymin;
    QPoint AT = rotate(QPoint(x[0],y[0]),a);
    xmax = xmin = AT.x();
    ymax = ymin = AT.y();


    /// пройдемся по всем данным чтобы определить габариты

    for ( int m = 0 ; m < x.size() ; m++ )
    {
        QPoint AT = rotate(QPoint(x[m],y[m]),a);
        if (AT.x()>xmax) xmax = AT.x();
        if (AT.x()<xmin) xmin = AT.x();
        if (AT.y()>ymax) ymax = AT.y();
        if (AT.y()<ymin) ymin = AT.y();
        // #form first optimization factor
    }

/// вернем щначение размеров , углов и точки для построения .
*return_size = QSize(abs(xmax-xmin),abs(ymax-ymin));
*Angular = a;


    return ;



}

/// Поворот точки вокруг начала координат на угол Ang радиан

QPoint ScanViewer::rotate(QPoint Z,qreal Ang)
{
    /// сохраним на будущее синус и косинус чтобы лишний раз не считать

    qreal C = cos(Ang);
    qreal S = sin(Ang);

    /// домножение на матрицу поворота в раскрытом виде
    /// приведение типов необходимо чтобы перемножать координаты в целых мм на вещественные синусы и косинусы
    /// не теряя информации

    qreal X = qreal(Z.x())*C-qreal(Z.y())*S;
    qreal Y = qreal(Z.x())*S+qreal(Z.y())*C;

    /// Возвращаем точку .
    return QPoint(X,Y);
}