Jump to content
Compvision.ru
Sign in to follow this  
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

Share this post


Link to post
Share on other sites

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

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

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

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

LibICP_example.png

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

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

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

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

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

init_state.png

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

ransac_result.png

Share this post


Link to post
Share on other sites

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

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

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

Share this post


Link to post
Share on other sites

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

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

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

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

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

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

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

Share this post


Link to post
Share on other sites

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;

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

Share this post


Link to post
Share on other sites

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

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

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

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

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

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

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

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

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

TPS_L2_no_guess.png

Share this post


Link to post
Share on other sites

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;

};

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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

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

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

последний код отсюда 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-х наборов точек.

Share this post


Link to post
Share on other sites

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]);
через первые собственные значения.

Share this post


Link to post
Share on other sites

вообщем похоже там по коду используется этот алгоритм 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.

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

Share this post


Link to post
Share on other sites

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

возможно это связан ос темой "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]

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

Robust Point Set Registration Using Gaussian Mixture Models

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

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

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

они там предлагают в пейпэре минимизировать 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.

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this  

  • Recently Browsing   0 members

    No registered users viewing this page.

×