Перейти к содержимому
Compvision.ru
mrgloom

ICP(Iterative Closest Point)

Recommended Posts

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

библиотека посерьёзней и более приятная. (нужен буст, но собралась без проблем, только в дебаг версии 1 функция почему то не была определена check_query_in_bound)

http://www.cvlibs.net/software/libicp.html

прямая ссылка

http://www.cvlibs.net/downloads/libicp.zip

пока правда не понял какие ограничения в демо, только поворот и сдвиг.

  • Like 1

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

попробовал, почему то не всегда работает нормально при малом кол-ве точек.

и опять же только rigid регистрация(т.е. только поворот и сдвиг)

и не всегда нормально сходится, наверно надо как то давать приближение по другому.

пример работы.

LibICP_example.png

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

я вот подумал, а что если ransac приспособить под это дело.

т.е. берем рендомно из 1 сета точек 4 точки, так же из второго, считаем матрицу перехода, затем трансформируем точки первого сета, считаем метрику(какую не уверен ну например сумма мин расстояний для каждой точки) повторяем всё это k раз.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

попробовал прикрутить ранзак для этих дел. но почеу то работает неправильно, видимо выбрана неправильная метрика.

я для каждой точки нахожу ближайшую точку, потом суммирую и делю на кол-во точек.

хотя возможно надо просто какое то ограничение на матрицу трансформации давать.

еще есть идея, что если мы не имеем пар точек от дескрипторов, то можно попытаться найти пары точек, определяя дескриптор как shape context.

начальное состояние

init_state.png

расстрел из 1000 попыток

ransac_result.png

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

похоже я понял почему ничего не получилось, т.к. когда я беру рэндомно по 4 точки и считаю матрицу перехода, то если мы имеем лишь похожие фигуры(т.е. ни одна пара не образует нормального преобразования), то перебрав даже все пары мы не получим того преобразования которое требуется.

т.е. получается, что надо не точки выбирать рэндомно, а саму матрицу, а это никаких ограничений, хотя возможно можно вывести примерные ограничения, например взять ограничивающие прямоугольники и из этого вывести ограничение на матрицу.

http://www.cse.psu.edu/~rcollins/CSE486/lecture15_6pp.pdf

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

всё таки непонятно по обычному ICP (где только поворот и сдвиг)

там написано для каждой точки первого сета найти ближайшую точку из второго сета и на основе этого сделать трансформацию.

http://www.cs.technion.ac.il/~cs236329/tutorials/ICP.pdf

так вот как делать эту трансформацию?

вроде как можно как СЛАУ через МНК как я пробовал тут

http://www.compvision.ru/forum/index.php?showtopic=999

но может надо подругому?

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

vector<Point> vec_pair;

    for(int i=0;i<vec_M.size();++i)

    {

        double min_dist=INT_MAX;

        int id=-1;

        for(int j=0;j<vec_T.size();++j)

        {

            double metric= sqrt(double(vec_T[j].x-vec_M[i].x)*(vec_T[j].x-vec_M[i].x)+

                (vec_T[j].y-vec_M[i].y)*(vec_T[j].y-vec_M[i].y));

            if(min_dist>metric)

            {

                min_dist=metric;

                id=j;

            }

        }

        line(img,vec_M[i],vec_T[id],cvScalar(0,0,255));

        vec_pair.push_back(vec_T[id]);

    }

    Mat m= estimateRigidTransform(vec_M,vec_pair,1);

    cout<<m;

попробовал так, но не дает нормальный результат, возможно надо отсекать какие либо пары?

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

А здесь не то что Вам нужно?:

http://code.google.com/p/gmmreg/

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

да я это уже видел,и даже оно завелось, но там алгоритмы сложные и не очень понятно как управлять параметрами.

но похоже придётся всё равно к этому вернутся.

там экзешник, где все алгоритмы требует, какую то довольно большую стороннюю либу.(какую не помню давно собирал)

впрочем там экзешник есть, можно не пересобирать.

а чтобы завести под питон мне пришлось ставить еще 4 сторонние библиотеки для питона.

причем scipy почему то через pip ставится отказался.

на матлабе тоже пробовал, но тоже давно.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

вообщем попробовал на тестовом примере там где рыбу из двух частей хотят составить.

EM_TPS

EM_GRBF

TPS_L2

GRBF_L2

TPS_KC

GRBF_KC

rigid

в итоге если не указывать ctrl_pts. (хотя я не понял, что это такое то ли начальное приближение, то ли их можно как то сформировать из имеющихся точек)

# 'ctrl_pts' serve as the control points when thin plate splines (TPS) or

# Gaussian radial basis functions (GRBF) are used in the nonrigid registration

# if 'ctrl_pts' is not provided, model will be used as ctrl_pts

# the program 'gmmreg_aux' can be used to generate ctrl pts from regular grid pts

сработали нормально только TPS_L2,TPS_KC.

но на моем сете 6к и 8к точек, алгоритм выделил слишком много памяти и благополучно вылетел.

вроде память выделяется когда как раз не указаны ctrl_pts.

вот так работает rigid - видимо обычный ICP

rigid_no_guess.png

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

забавно, но эти оказались просто сеткой с каким то шагом, причем если их использовать память не кушается, странно почему алгоритм не может сам генерировать эту сетку, а надо ее делать отдельно.

но все равно ничего хорошего не получилось,т.к. надо видимо примерное приближение всё таки.

TPS_L2_no_guess.png

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

http://en.wikipedia.org/wiki/Weierstrass_transform

начал смотреть rigid версию в gmmreg, похоже код там поддерживает еще и скейл.

втретил там некий GaussTransform и похоже это не размытие, как то связано с методом сопряженного градиента(conjugate gradient method) или каким то еще итеративным методом, но причём тут гауссиан?

double energy = GaussTransform(gmmreg->transformed_model,

gmmreg->scene, scale, gradient);

вообщем похоже все хвосты ведут на какой то минимайзер, у которого виртуальные функции определенные в самом классе rigid gmmreg.

только всё таки непонятно, что это зам метод.

//:   An object that represents a function from R^n -> R.

//    It is commonly used to express the

//    interface of a minimizer.

class vnl_cost_function : public vnl_unary_function<double, vnl_vector<double> >

{

 public:


  //! Default constructor   

  vnl_cost_function():dim(0) {}


  //! Construct with a specified number of unknowns

  vnl_cost_function(int number_of_unknowns):dim(number_of_unknowns) {}


  virtual ~vnl_cost_function() {}


  //:  The main function.  Given the parameter vector x, compute the value of f(x).

  virtual double f(vnl_vector<double> const& x);


  //:  Calculate the gradient of f at parameter vector x.

  virtual void gradf(vnl_vector<double> const& x, vnl_vector<double>& gradient);


  //:  Compute one or both of f and g.

  //   Normally implemented in terms of the above two, but may be faster if specialized. f != 0 => compute f

  virtual void compute(vnl_vector<double> const& x, double *f, vnl_vector<double>* g);


  //:  Return the number of unknowns

  int get_number_of_unknowns() const { return dim; }


  //:  Compute finite-difference gradient

  void fdgradf(vnl_vector<double> const& x, vnl_vector<double>& gradient, double stepsize = 1e-5);


  //:  Called when error is printed for user.

  virtual double reported_error(double f_value) { return f_value; }


  //:  Conveniences for printing grad, fdgrad

  vnl_vector<double> gradf(vnl_vector<double> const& x);

  vnl_vector<double> fdgradf(vnl_vector<double> const& x);


protected:


    //! Set number of unknowns. 

    void set_number_of_unknowns(int number_of_unknowns) { dim=number_of_unknowns; }


public:

    int dim;

};

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Контрольные точки, насколько я понял, это точки по которым считается сплайн-деформация.

Например для лица это будут глаза, уголки губ и т.д.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Есть еще ресурс с OpenCV-шным кодом:

http://www.morethantechnical.com/2010/06/06/iterative-closest-point-icp-with-opencv-w-code/

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

а можете пояснить для какого рода задач вы используете алгоритм?

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Контрольные точки, насколько я понял, это точки по которым считается сплайн-деформация.

Например для лица это будут глаза, уголки губ и т.д.

какие контрольные точки?

последний код отсюда http://www.morethantechnical.com/2010/06/06/iterative-closest-point-icp-with-opencv-w-code/

какой то мутный я его пробовал, вроде даже он не заработал на моих примерах. вообщем что то там мне не понравилось, во всяком случае чтобы его исправить надо понимать что там делается.

по поводу conjugate gradient method там это помоему задача energy minimization

типа есть transformed_model, а есть scene и видимо итеративно смотрят как хорошо они наложились при каком либо преобразовании.

написал еще китайцу который делал gmmreg, но видимо это бесполезно.а код довольно жёсткий, минимум комментариев, иногда вставлены куски-аналоги из матлаба.

а можете пояснить для какого рода задач вы используете алгоритм?

регистрация 2-х наборов точек.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

http://www.morethantechnical.com/2012/08/12/registering-point-clouds-rigidly-with-scale-using-pcl-wcode/#more-1129

а тут вот еще rigid через PCL

и как то скейл определяют непонятно через PCA

double s = sqrt(ev_B[0])/sqrt(ev_A[0]);
через первые собственные значения.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

вообщем похоже там по коду используется этот алгоритм minimization with constraints

http://en.wikipedia.org/wiki/Limited-memory_BFGS

//: Limited memory Broyden Fletcher Goldfarb Shannon minimization with constraints.

// Lower and upper bounds may be specified for the variables to be optimized.

// The algorithm miminizes a nonlinear function f(x) of n variables

// subject to simple bound constraints of l <= x <= u.

осталось понять какая задаётся функция для минимизации.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

вообщем вся соль похоже в этих функциях, но почему они так считаются не понятно.

возможно это связан ос темой "robust point set registration using mixture of Gaussians"

//непонятно назначение, по сути просто меняет знак для энергии и градиента

double gmmreg_rigid_func::eval(double& f, vnl_matrix<double>& g)

{

  /* L2 version and KC version are equivalent in the rigid case */

  g =  -g;

  return -f;

}


//основная функция которая вычисляет функцию от 3-х переменных(х,у смещение+ поворот)

//непонятно каким боком тут GaussTransform(вроде бы используется градиент с предыдущей итерации)

double gmmreg_rigid_func::f(const vnl_vector<double>& x)

{

  // std::cout << "f begin" << std::endl;

  gmmreg->perform_transform(x);  //x это вроде параметры которые мы пытаемся определить

  // std::cout << "transform done" << std::endl;

  double energy = GaussTransform(gmmreg->transformed_model,

      gmmreg->scene, scale, gradient);

  // std::cout << "gauss transform done" << std::endl;

  return eval(energy, gradient);

  // std::cout << "f end" << std::endl;

}


//вычисление градиента

void gmmreg_rigid_func::gradf(const vnl_vector<double>& x,

    vnl_vector<double>& g)

{

//код из матлаба как называется этот метод?

//похоже как то связано с "robust point set registration using mixture of Gaussians" т.к. rigid_costfunc указывает на это

/*     case 'rigid2d'

 *         [f, grad] = rigid_costfunc(transformed_model, scene, scale);

 *         grad = grad';

 *         g(1) = sum(grad(1,);

 *         g(2) = sum(grad(2,);

 *         grad = grad*model;

 *         theta = param(3);

 *         r = [-sin(theta) -cos(theta);

 *              cos(theta)  -sin(theta)];

 *         g(3) = sum(sum(grad.*r));

 *

 */

	int i = 0;

	vnl_matrix<double> r;

	vnl_matrix<double> gm;

	if (d == 2) 

	{ //rigid2d

		g[0] = gradient.get_column(0).sum();

		g[1] = gradient.get_column(1).sum();

		r.set_size(2,2);

		double theta = x[2];

		r[0][0] = -sin(theta);

		r[0][1] = -cos(theta);

		r[1][0] = cos(theta);

		r[1][1] = -sin(theta);

		gm = gradient.transpose()*gmmreg->model; //что дает перемножение градиента на модель?

		g[2] = 0;

		for (i=0;i<2;++i)

		{

		  for (int j=0;j<2;++j)

		  {

			g[2] += r[i][j]*gm[i][j];

		  }

		}

	}

	else if (d == 3)

	{ //rigid3d

/*

*     case 'rigid3d'

*        [f,grad] = rigid_costfunc(transformed_model, scene, scale);

*         [r,gq] = quaternion2rotation(param(1:4));

*         grad = grad';

*         gm = grad*model;

*         g(1) = sum(sum(gm.*gq{1}));

*         g(2) = sum(sum(gm.*gq{2}));

*         g(3) = sum(sum(gm.*gq{3}));

*         g(4) = sum(sum(gm.*gq{4}));

*         g(5) = sum(grad(1,);

*         g(6) = sum(grad(2,);

*         g(7) = sum(grad(3,);

*/

    g[4] = gradient.get_column(0).sum();

    g[5] = gradient.get_column(1).sum();

    g[6] = gradient.get_column(2).sum();

    vnl_vector<double> q;

    q.set_size(4);

    for (i=0;i<4;++i) q[i] = x[i];

    r.set_size(3,3);

    vnl_matrix<double> g1,g2,g3,g4;

    g1.set_size(3,3);

    g2.set_size(3,3);

    g3.set_size(3,3);

    g4.set_size(3,3);

    quaternion2rotation(q,r,g1,g2,g3,g4);

    gm = gradient.transpose()*gmmreg->model;

    g[0] = 0;

    for (i=0;i<3;++i) {

      for (int j=0;j<3;++j) {

        g[0] += g1[i][j]*gm[i][j];

      }

    }

    g[1] = 0;

    for (i=0;i<3;++i) {

      for (int j=0;j<3;++j) {

        g[1] += g2[i][j]*gm[i][j];

      }

    }

    g[2] = 0;

    for (i=0;i<3;++i) {

      for (int j=0;j<3;++j) {

        g[2] += g3[i][j]*gm[i][j];

      }

    }

    g[3] = 0;

    for (i=0;i<3;++i) {

      for (int j=0;j<3;++j) {

        g[3] += g4[i][j]*gm[i][j];

      }

    }

  }

}


[/code]

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

GaussTransform - похоже нужен для получения более менее гладкого градиента, чтобы минимизатор работал лучше.

На кусочно-гладкой функции производная будет "вылетать" в бесконечность на стыках гладких кусков.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

там используется метод

Robust Point Set Registration Using Gaussian Mixture Models

т.е. сначала наборы точек представляются как Gaussian Mixture Models ,а потом уже ка кто совмещаются.

но как это относится к коду пока непонятно.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

При совмещении распределений обычно минимизируется дивергенция Кульбака — Лейблера.

http://ru.wikipedia.org/wiki/%D0%A0%D0%B0%D1%81%D1%81%D1%82%D0%BE%D1%8F%D0%BD%D0%B8%D0%B5_%D0%9A%D1%83%D0%BB%D1%8C%D0%B1%D0%B0%D0%BA%D0%B0_%E2%80%94_%D0%9B%D0%B5%D0%B9%D0%B1%D0%BB%D0%B5%D1%80%D0%B0

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

они там предлагают в пейпэре минимизировать L2 distance между микстурами гауссиан.

и там это формулируется как inner product of gaussian mixture, только непонятно как выглядят сами микстуры гауссиан, как векторы?

и опять же я не очень понимаю, ну допустим мы сформировали положительно определенную cost function и пытаемся каким то там методом найти её минимум, допустим она зависит от 3-х параметров, как мы определяем куда двигать эти параметры чтобы сошлось?

(хотя можно минимизатор рассматривать как черный ящик, а задать только целевую функцию и ограничения)

тут вот нашел разные методы\метрики для того чтобы посчитать разницу между микстурами гауссиан. (правда тулбокс для аудио)

http://kom.aau.dk/project/isound/isp_gmmdistance.html

что еще означает affinity matrix которая часто фигурирует.

там аж 3 метода с этим связанных во введении описывается

Instead of assuming the one-to-one correspondence

based on the nearest neighbor criterion, one-to-many

relaxations have been proposed to allow for fuzzy correspondences.

One of the most notable contributions

along these lines is the work in Chui and Rangarajan

[2], wherein a soft assignment technique and the deterministic

annealing are combined to alternatively estimate

the transformation and update the correspondence. The

key idea of their work is to assume that each model

point corresponds to a weighted sum of the scene points,

instead of the closest scene point alone. These weights

are taken from an affinity matrix whose entries are

proportional to a Gaussian function of the pairwise distances

between the moving model and the fixed scene.

In other words, this affinity matrix assigns each model

point to all points in the scene. The closer the scene point

is to the model point, the larger is the assigned weight.

Clearly, the closest point correspondence in ICP can be

viewed as a binary version of this graduated assignment

scheme. In Chui and Rangarajan [2], the registration

problem is expressed as an joint optimization over the

transformation parameter and correspondence matrix.

However, in practice, it is implemented as an alternating

update strategy: at each iteration, a transformation is estimated

from the virtual correspondence between model

points and the linear combinations of scene points, then

the correspondence matrix gets modified as the model

points move towards the scene and the scale parameter

in the Gaussian function for computing the weights

decreases.

Gaussian affinity matrix

Recently Myronenko et al. [6] proposed yet

another robust non-rigid point set registration algorithm.

Like in Chui and Rangarajan [2], Myronenko et al. [6]

maintain the same Gaussian affinity matrix and also

adopt a similar alternating update strategy interpreted in

an expectation maximization (EM) [20] framework. The

major difference between these two algorithms is that

in Chui and Rangarajan [2] the non-rigid deformation

is modeled by thin-plate splines while in Myronenko

et al. [6] they use Gaussian radial basis functions.

Another trick is to

model each of the two point sets by a probability distribution

and then a distance measure between the two

distributions is minimized over the transformation space

to yield the desired transformation. For instance, Tsin

and Kanade [4] proposed a kernel correlation (KC) based

point set registration approach where the cost function

is proportional to the correlation of two kernel density

estimates. Glaunes et al. [24] formulated the problem

of aligning two unlabeled point sets by first modeling

the point sets as weighted sums of Dirac measures and

then finding the optimal diffeomorphic transformation

between the two discrete distributions.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Создайте учётную запись или войдите для комментирования

Вы должны быть пользователем, чтобы оставить комментарий

Создать учётную запись

Зарегистрируйтесь для создания учётной записи. Это просто!

Зарегистрировать учётную запись

Войти

Уже зарегистрированы? Войдите здесь.

Войти сейчас


  • Сейчас на странице   0 пользователей

    Нет пользователей, просматривающих эту страницу

×