I have been using ceres-solver for a long time now and it is an amazing tool. My usage up until now was not based on reusable code and I am trying to improve this. Ceres uses a specific structure with a specific templated method as an interface to its automatic differentiation. In the problem that I am trying to solve, inheritance makes sense because the different cost functions that I need to have are very similar to each other. I have created an example that is similar (but it makes no sense, sorry). Imagine we want to be able to find the polygon that has a given area. In my example, polygons can be either triangles or rectangles. With this in mind, it makes sense to have a base class that implements everything and specific classes that implement, in this case, the area computation for each specific polygon:
ShapeCostFunction
class shapeAreaCostFunction
{
public:
shapeAreaCostFunction(double desired_area): desired_area_(desired_area){}
template<typename T>
bool operator()(const T* shape, T* residual) const{
residual[0] = T(desired_area_) - area(shape);
return true;
}
template<typename T>
virtual T area(const T* shape) const = 0;
protected:
double desired_area_;
};
RectangleCostFunction
#include "shapeAreaCostFunction.h"
#include "areaLibrary.h"
class rectangleAreaCostFunction : public shapeAreaCostFunction
{
public:
rectangleAreaCostFunction(double desired_area): shapeAreaCostFunction(desired_area){}
template<typename T>
T area(const T* triangle) const
{
return rectangleArea(triangle);
}
};
TriangleCostFunction
#include "shapeAreaCostFunction.h"
#include "areaLibrary.h"
class triangleAreaCostFunction : public shapeAreaCostFunction
{
public:
triangleAreaCostFunction(double desired_area): shapeAreaCostFunction(desired_area){}
template<typename T>
T area(const T* triangle) const
{
return triangleArea(triangle);
}
};
AreaLibrary
template<typename T>
T rectangleArea(const T* rectangle)
{
return rectangle[0]*rectangle[1];
}
template<typename T>
T triangleArea(const T* triangle)
{
return rectangleArea(triangle)/T(2);
}
Main
#include <ceres/ceres.h>
#include <iostream>
#include "rectangleAreaCostFunction.h"
#include "triangleAreaCostFunction.h"
#include "areaLibrary.h"
int main(int argc, char** argv){
// Initialize glogging
//google::InitGoogleLogging(argv[0]);
// Get values
/// Get total area
double total_area;
std::cout<<"Enter the desired area: ";
std::cin>>total_area;
/// Get initial rectangle
double rect[2];
std::cout<<"Enter initial rectangle base: ";
std::cin>>rect[0];
std::cout<<"Enter initial rectangle height: ";
std::cin>>rect[1];
/// Get initial triagnle
double tri[2];
std::cout<<"Enter initial triangle base: ";
std::cin>>tri[0];
std::cout<<"Enter initial triangle height: ";
std::cin>>tri[1];
// Copy initial values
double rect_ini[] = {rect[0],rect[1]};
double tri_ini[] = {tri[0],tri[1]};
// Create problem
ceres::Problem problem;
ceres::CostFunction* cost_function_rectangle = new ceres::AutoDiffCostFunction<rectangleAreaCostFunction, 1, 2>(
new rectangleAreaCostFunction(total_area));
ceres::CostFunction* cost_function_triangle = new ceres::AutoDiffCostFunction<triangleAreaCostFunction, 1, 2>(
new triangleAreaCostFunction(total_area));
problem.AddResidualBlock(cost_function_rectangle, NULL, rect);
problem.AddResidualBlock(cost_function_triangle, NULL, tri);
// Solve
ceres::Solver::Options options;
options.linear_solver_type = ceres::DENSE_QR;
options.minimizer_progress_to_stdout = true;
options.max_num_iterations = 10;
ceres::Solver::Summary summary;
ceres::Solve(options, &problem, &summary);
// Compute final areas
double rect_area = rectangleArea(rect);
double tri_area = triangleArea(tri);
// Display results
std::cout << summary.FullReport() << std::endl;
std::cout<<"Rectangle: ("<<rect_ini[0]<<","<<rect_ini[1]<<") -> ("<<rect[0]<<","<<rect[1]<<") total area: "<<rect_area<<"("<< rect_area - total_area<<")"<<std::endl;
std::cout<<"Triangle: ("<<tri_ini[0]<<","<<tri_ini[1]<<") -> ("<<tri[0]<<","<<tri[1]<<") total area: "<<tri_area<<"("<< tri_area - total_area<<")"<<std::endl;
// Exit
return 0;
}
The problem with this is that templated functions cannot be virtual as explained several times in stackoverflow (here and here). However, it seems there are some workarounds using boost::any
. I have tried to use this in my example with no success. I have also tried to move the template from the class method to the class, similar to here, but ceres do not accept it as a cost function.
My questions are (and please, keep in mind that I am restricted to have the method template<typename T> bool operator()(...)const
otherwise I cannot interact with ceres):
template<typename T> bool operator()(...)const
class method? Thank you in advance.
I can think of two approaches.
First, compose lambdas. Second, use CRTP.
This is best done with c++14.
template<class Area>
auto cost_function(Area area, double desired){
return [=](auto const* shape, auto* residual){
using T=std::decay_t<decltype(*shape)>;
residual[0] = T(desired_area_) - area(shape);
return true;
};
}
auto triangle = [](auto* shape){return triangleArea(shape);};
To create a triangle area cost function:
auto tri_cost = cost_function(triangle, 3.14159);
and to get the type, decltype(tri_cost)
.
So:
auto tri_cost = cost_function(triangle, 3.14159);
ceres::CostFunction* cost_function_triangle = new ceres::AutoDiffCostFunction<decltype(tri_cost), 1, 2>(
new decltype(tri_cost)(tri_cost));
You can do a similar composition technique without lambdas, but it is more tedious. You can also wrap some of these naked new's up in helper functions.
template<class D>
class shapeAreaCostFunction {
public:
shapeAreaCostFunction(double desired_area): desired_area_(desired_area){}
template<typename T>
bool operator()(const T* shape, T* residual) const{
residual[0] = T(desired_area_) - static_cast<D const*>(this)->area(shape);
return true;
}
protected:
double desired_area_;
};
modify derived types like this:
class triangleAreaCostFunction :
public shapeAreaCostFunction<triangleAreaCostFunction>
{
using base=shapeAreaCostFunction<triangleAreaCostFunction>;
public:
triangleAreaCostFunction(double desired_area): base(desired_area){}
template<typename T>
T area(const T* triangle) const
{
return triangleArea(triangle);
}
};
this is known as using the CRTP to implement static polymorphism.