Help with Ceres & sampling OpenCV Mat in cost function

[Andrew Hogue]

I'm new to Ceres and trying to determine the appropriate way to use it to solve my minimization problem.  I was hoping someone could provide suggestions as I'm stuck at the moment.

Problem: 

- I have a 3D model of a subject, a number of cameras with estimates on extrinsics, and intrinsics.

- I also have an Signed Distance Function image that is ZERO at the boundary of the subject, negative inside and positive outside

- My current estimate on Extrinsics is correct up to a small offset in translation and potentially in rotation 

- My minimization problem is thus to estimate the translation offset from the current extrinsics given all of the above that minimizes the projected distances.

Currently, trying to use the AutoDifferentiation with a templated cost function but to compute the residuals, I need to first:

- project internal mesh using known intrinsics and estimated extrinsics 

- then sample the SDF image (OpenCV mat) at the projected locations to sum up the error residuals over the mesh.

Questions:

1. Is this an appropriate method or am I modeling it wrong?  Any suggestions would be helpful.  Should I have a residual per vertex or for the entire model?

2. How do I properly use an OpenCV mat in this instance?  The templated parameter is obviously a Jet, and all calculations internal to the cost function will be thus using this templated Jet.  How does one convert the projected point (u,v) into an index I can then use to grab the value from my OpenCV mat?  

There isn't much info or suggestions on using OpenCV Mat's internally in this context, one suggestion was to use the cubic-interpolator and Grid2D interface.  Does anyone have any examples on this for OpenCV ?

Thanks in Advance!

-- Andrew

[Sameer Agarwal]

Hi Andrew,

Let me understand your problem better first.

The signed distance image is just a representation for the contour of the object right?

And when you say you have a 3d model, you are saying you have 3d mesh of the object placed in space ? or do you just have the mesh and you have to also figure out the translation and rotation of the mesh in space for it to line up with the images?

and you want to refine the camera parameters so that the contours of the mesh when projected into the image line up?

Before we talk about automatic differentiation etc, what is the objective function that you are trying to minimize? can you give its mathematical form?

Sameer

--
You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/0cf8d67f-1391-4b2a-87e8-bcdbd8e20271n%40googlegroups.com.

[Andrew Hogue]

Hi Sameer, sorry I just accidentally sent the message privately.  Here is a recap for the group:

Here is what I have:

- A number of Kinect Azures (6-8) with pre-determined intrinsics/extrinsics but they are not necessarily perfectly accurate, I'd like to improve the extrinsics here

- a 3D mesh (triangular mesh) 

-RGB images and Depth images from the cameras as well as a B&W matte image for identifying the subject

- an SDF image denoting the contour (zero on the boundary of the matte image)

Essentially, when I was thinking of using ceres to refine the extrinsics of the cameras using the SDF in image space to help since the gradient will aid in the direction towards the contours

The cost function is thus:

- p = K*[R|t] * v

- [u,v] = p/p_z

- sdf = SDF(u,v)

SumSDF += sdf^2;

then minimize the SumSDF cost over the mesh vertices....

[Dmitriy Korchemkin]

Is 3d-mesh a "ground-truth" triangular mesh of the object? In this case it might make sense optimizing distance from the depth data to the mesh.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/838138a6-ffa4-4afe-8757-76cf46f9d905n%40googlegroups.com.

--

Dmitriy A Korchemkin

[Andrew Hogue]

Sort of.  Either way I need to be able to sample the (u,v) from an OpenCV mat in the cost function, if the (u,v) comes from a calculation using the free templated parameters then how would I convert it to something I can use as an index, or I suppose convert the OpenCV mat into something I can sample from the cost function?

I'm hesitant to use just the depth data here anyways, I want the problem to incorporate an appearance model later on in the cost estimate as well but baby steps first :) 

[Dmitriy Korchemkin]

If interpolation of the data you store in cv::Mat makes sense from the modelling perspective - there is a ceres::BicubicInterpolator class  http://www.ceres-solver.org/nnls_modeling.html?highlight=evaluate#_CPPv4N5ceres19BiCubicInterpolatorE for obtaining a differentiable (and auto-differentiable) access to an underlying constant array of data.

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/cb0f9020-3eb3-4e96-b5aa-de37583cc594n%40googlegroups.com.

--

Dmitriy A Korchemkin

[Sameer Agarwal]

As Dimitry pointed out, this is exactly what the bicubicinterpolator is for. 

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/CAN8XMixqsK-3UQEMsL3oCOEefdsPMhksdMuGY_0WHkeCZPhWvg%40mail.gmail.com.

[Andrew Hogue]

Hi again,

This doesn't solve my issue however.... I now have a class with a BiCubicInterpolator defined as a member,  I create this in the constructor and populate it with the image data.

In my templated operator() for the cost function, it now needs to use this cubic interpolator to sample the image values.

So, I do the appropriate steps and take compute the projection using the free parameter for extrinsics and my fixed intrinsics.

I project and get a (u,v) which is now using the templated parameter, T.

However, I can't use this to sample with the BiCubicInterpolator since it has originally been defined with double as the type in the main class:

error C2664: 'void ceres::BiCubicInterpolator<ceres::Grid2D<double,1,true,true>>::Evaluate(const double &,const double &,double *) const': cannot convert argument 1 from 'T' to 'double'

So, what is the process to use the BiCubicInterpolator in the cost function with auto-differentiation?  Do I have to create a new interpolator in the operator() cost function each time it is called (which would be inefficient) or am I missing something completely when I define the member in my class...?

If anyone has a substantial example that solves a similar issue I would love to see it.

Thanks for any insight,

[Sameer Agarwal]

Andrew,

The interpolator has two Evaluate methods.

One is as you point out defined over doubles, but there is another one which allows you to use Jets for input and output.

https://github.com/ceres-solver/ceres-solver/blob/master/include/ceres/cubic_interpolation.h#L335

This is the one you should calling.

Sameer

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/917e57ba-8971-4283-a453-d214bde87215n%40googlegroups.com.

[Andrew Hogue]

ah, cool.  Do you have (or know of) a good example of using that within a cost function? 

Just took a look at the evaluator examples and I can see how to wrap the R,C into 2 different Jets passed to the function.... still fighting with it though :) 

If you know of any example I'd love to see it

Thanks again

[Sameer Agarwal]

Andrew I do not have an example handy, but a quick search on github reveals some examples

https://github.com/HKUST-Aerial-Robotics/edge_alignment/blob/5e76bd601760e9a929463a4bdc3eb2bdd78d1eb3/standalone/utils.h#L52

Sameer

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/65870dd7-7f72-4fa2-872e-177569751254n%40googlegroups.com.

[Dmitriy Korchemkin]

Here is a small "end to end" sample: https://gist.github.com/DmitriyKorchemkin/93efd0de71d6fb234f9807e568e06b05

[Sameer Agarwal]

Dmitriy,

Can you convince you to contribute this as an example to Ceres? 

Sameer

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/368dfa32-897e-486f-a49c-9c4bdc9f395fn%40googlegroups.com.

[Andrew Hogue]

thank you so much!  That helps immensely

[Dmitriy Korchemkin]

I'll submit change request to Gerrit on weekend.

[Sameer Agarwal]

Thank you!

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/9c6f3cb3-6edc-49ff-9b1b-64e3e79b0aa2n%40googlegroups.com.

[Andrew Hogue]

Also, something to note, please verify that in the interpolator it is Row/Col vs Col/Row...... 

For instance, I wanted to verify that the interpolator was going to work for me as you had above and populated it with my OpenCV data and then afterwards I pulled it out of the interpolator into a new image I could display it to verify correctness.

In order for me to obtain correct values per pixel, I had to have the grid be addressed using (Cols, Rows) in a snippet below.  This works but it is not obvious why by default it seems to be using Column Major even though the default should be RowMajor...... something to verify with your own implementation.

const int cols = imSDF.cols;

        const int rows = imSDF.rows;

        Eigen::ArrayXXd img_data = Eigen::ArrayXXd(rows,cols);

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

            for (int i = 0; i < cols; i++)

            {

                float sdf = imSDF.at<float>(cv::Point(i,j));

                img_data(j, i) = (double)sdf;

            }

        }

Grid theGrid(img_data.data(), 0, img_data.cols(), 0, img_data.rows());

        theInterpolator = new Interpolator(theGrid);

//   pulling data per pixel out here to test:

   cv::Mat imTest = cv::Mat::zeros(cv::Size(cols,rows), CV_32F);

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

            for (int i = 0; i < cols; i++)

            {

                ceres::Jet<double, 4> r_jet, c_jet,ff;

                r_jet.a = (double)j;

                r_jet.v[0] = 1.0;

                r_jet.v[1] = 1.1;

                r_jet.v[2] = 1.2;

                r_jet.v[3] = 1.3;

                c_jet.a = (double)i;

                c_jet.v[0] = 2.0;

                c_jet.v[1] = 2.1;

                c_jet.v[2] = 2.2;

                c_jet.v[3] = 2.3;

               

                theInterpolator->Evaluate(c_jet,r_jet, &ff);

                imTest.at<float>(cv::Point(i, j)) = (float)ff.a;

            }

        }

        cv::imshow("test", imTest);

[Dmitriy Korchemkin]

Are you sure the problem is not in Eigen::ArrayXXd being column-major by default (and other differences in notation between 3 libraries you're using together)?

Note also that BiCubicInterpolator::Evaluate expects coordinates supplied in order {row, col} (and opencv - in {col, row}, when used via cv::Mat::at(cv::Point) method).

[Andrew Hogue]

possibly... there's a lot going on in there that I have to sort out I'm sure.... something to verify and check if you are ever trying something similar..... 

I bet its possibly something I'm doing.... but at the moment it will just be documented and re-verified later...

Thanks for all your help!

[Andrew Hogue]

aha, you were right of course.

Switched everything back to row-major as default and then changed my Eigen::ArrayXXd as:

using EigArray = Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;

just to force it to be rowmajor.... works now :) 

Thanks, I hope this discussion helps someone later!

[Andrew Hogue]

Just to finish up the discussion on this thread.  I was able to get it working with your help!  Thanks!

If you take a look at this image you can see the results of before and after the nonlinear minimization.

I parameterized the problem as follows:

- solving for a 3D translation vector offset for my camera extrinsics

- added a residual block/cost function for each vertex in the input mesh

- cost function performs a projection using the intrinsics and extrinsics + offset from the parameters

- residual is computed using a distance field lookup image computed from a B&W matte image 

Thanks for all your help!

[Sameer Agarwal]

Awesome!

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/ab4ad086-e396-4515-8c3a-39b8993f593bn%40googlegroups.com.