00001
00002 #include <iostream>
00003 #include <boost/thread/thread.hpp>
00004 #include <boost/program_options.hpp>
00005
00006 #include <ndt_feature_reg/ndt_frame.h>
00007 #include <ndt_feature_reg/ndt_frame_proc.h>
00008 #include <ndt_feature_reg/ndt_frame_viewer.h>
00009
00010 #include <ndt_matcher_d2d_feature.h>
00011 #include <pointcloud_utils.h>
00012
00013 #include <cv.h>
00014 #include <highgui.h>
00015
00016 #include <pcl/visualization/pcl_visualizer.h>
00017 #include <pcl/io/io.h>
00018 #include <pcl/io/pcd_io.h>
00019 #include <pcl/common/transforms.h>
00020
00021 using namespace ndt_feature_reg;
00022 using namespace pcl;
00023 using namespace cv;
00024 using namespace std;
00025 using namespace lslgeneric;
00026
00027 namespace po = boost::program_options;
00028
00029 inline
00030 std::vector<std::string> split_string(const std::string &str, const std::string &separator)
00031 {
00032 std::string::size_type it = 0;
00033 std::vector<std::string> values;
00034 bool stop = false;
00035 while (!stop)
00036 {
00037 std::string::size_type start = it;
00038 it = str.find(separator, it+1);
00039 if (it == std::string::npos)
00040 {
00041 values.push_back(str.substr(start, std::string::npos));
00042 stop = true;
00043 }
00044 else
00045 {
00046 it++;
00047 values.push_back(str.substr(start, it-start-1));
00048 }
00049 }
00050 return values;
00051 }
00052
00053
00054 vector<pair<string,string> > load_associations(const string &association_file, vector<string> ×tamps, int skip_nb)
00055 {
00056
00057 vector<pair<string, string> > ret;
00058 std::ifstream file(association_file.c_str());
00059 std::string line;
00060 bool skip = false;
00061 int skip_counter = 0;
00062 while(std::getline (file,line))
00063 {
00064 vector<string> splits = split_string(line, std::string(" "));
00065 if (splits.size() < 4)
00066 {
00067 cerr << "couldn't parse : " << line << endl;
00068 continue;
00069 }
00070 if (!skip)
00071 {
00072 ret.push_back(pair<string,string>(splits[1], splits[3]));
00073 timestamps.push_back(splits[0]);
00074 }
00075
00076 if (skip_counter > skip_nb)
00077 {
00078 skip = false;
00079 skip_counter = 0;
00080 }
00081 else
00082 {
00083 skip = true;
00084 skip_counter++;
00085 }
00086 }
00087 return ret;
00088 }
00089
00090
00091
00092
00093 int main(int argc, char** argv)
00094 {
00095 cout << "--------------------------------------------------" << endl;
00096 cout << "Evaluations of f2f registration" << endl;
00097 cout << "--------------------------------------------------" << endl;
00098
00099 string association_name;
00100 string output_name;
00101 string times_name;
00102 double scale, max_inldist_xy, max_inldist_z;
00103 Eigen::Matrix<double,6,1> pose_increment_v;
00104 int nb_ransac;
00105 int nb_points;
00106 size_t support_size;
00107 size_t max_nb_frames;
00108 double max_var;
00109 double current_res;
00110 int delay_vis;
00111 double img_scale;
00112 double detector_thresh;
00113 int skip_nb;
00114 int freiburg_camera_nb;
00115 double max_kp_dist, min_kp_dist;
00116 double trim_factor;
00117 po::options_description desc("Allowed options");
00118 desc.add_options()
00119 ("help", "produce help message")
00120 ("debug", "print debug output")
00121 ("visualize", "visualize the output")
00122 ("association", po::value<string>(&association_name), "association file name (generated by association.py script)")
00123 ("output", po::value<string>(&output_name), "file name of the estimation output (to be used with the evaluation python script)")
00124 ("time_output", po::value<string>(×_name), "file name for saving the per-frame processing time (to be used with the evaluation python script)")
00125 ("scale", po::value<double>(&scale)->default_value(0.0002), "depth scale (depth = scale * pixel value)")
00126 ("max_inldist_xy", po::value<double>(&max_inldist_xy)->default_value(0.02), "max inlier distance in xy related to the camera plane (in meters)")
00127 ("max_inldist_z", po::value<double>(&max_inldist_z)->default_value(0.05), "max inlier distance in z - depth (in meters)")
00128 ("nb_ransac", po::value<int>(&nb_ransac)->default_value(1000), "max number of RANSAC iterations")
00129 ("nb_points", po::value<int>(&nb_points)->default_value(21), "number of surrounding points")
00130 ("support_size", po::value<size_t>(&support_size)->default_value(10), "width of patch around each feature in pixels")
00131 ("max_var", po::value<double>(&max_var)->default_value(0.3), "max variance of ndt cells")
00132 ("current_res", po::value<double>(¤t_res)->default_value(0.2), "resolution of ndt cells if doing full matching")
00133 ("skip_matching", "skip the ndt matching step")
00134 ("match_full", "skip the visual features step and perform only ndt matching")
00135 ("match_no_association", "skip the feature association step and match directly")
00136 ("estimate_di", "don't compute the ndt explicitly but estimate from neighbourhood depth")
00137 ("max_nb_frames", po::value<size_t>(&max_nb_frames)->default_value(10), "max number of frames to be keept in the processing step")
00138 ("delay_vis", po::value<int>(&delay_vis)->default_value(10), "param to opencv function waitKey (if visualize is used)")
00139 ("img_scale", po::value<double>(&img_scale)->default_value(1.), "if the rgb image should be scaled when feature are detected and extracted - mainly to boost the speed")
00140 ("detector_thresh", po::value<double>(&detector_thresh)->default_value(400.), "playaround parameter, used to test different detectorvalues")
00141 ("skip_nb", po::value<int>(&skip_nb)->default_value(0), "how many intermediate frames should be skipped between the registration is runned")
00142 ("old_reg_code", "use the old registration code - for debugging only")
00143 ("freiburg_camera_nb", po::value<int>(&freiburg_camera_nb)->default_value(1), "the camera used for the evaluation 1-> freiburg camera #1, 2-> freiburg camera #2")
00144 ("max_kp_dist", po::value<double>(&max_kp_dist)->default_value(10.), "max distance of keypoints used")
00145 ("min_kp_dist", po::value<double>(&min_kp_dist)->default_value(0.), "min distance of keypoints used")
00146 ("windowed_matching", "limit the search of the keypoint matches to a local region (will not work that good if the skip_nb is large)")
00147 ("trim_factor", po::value<double>(&trim_factor)->default_value(1.), "only utilize the trimfactor number of keypoints with the closest NDT distance - similar idea as in trimmed ICP")
00148 ("non_mean", "if the mean used in NDT matching (including the RANSAC step) are computed direclty from the closest depth image pixel")
00149 ;
00150
00151
00152 po::variables_map vm;
00153 po::store(po::parse_command_line(argc, argv, desc), vm);
00154 po::notify(vm);
00155
00156 if (!vm.count("output") || !vm.count("association"))
00157 {
00158 cout << "Missing association / output file names.\n";
00159 cout << desc << "\n";
00160 return 1;
00161 }
00162 if (vm.count("help"))
00163 {
00164 cout << desc << "\n";
00165 return 1;
00166 }
00167 bool debug = vm.count("debug");
00168 bool visualize = vm.count("visualize");
00169 bool skip_matching = vm.count("skip_matching");
00170 bool match_full = vm.count("match_full");
00171 bool match_no_association = vm.count("match_no_association");
00172 bool estimate_di = vm.count("estimate_di");
00173 bool old_reg_code = vm.count("old_reg_code");
00174 bool windowed_matching = vm.count("windowed_matching");
00175 bool non_mean = vm.count("non_mean");
00176
00177 if (times_name == "")
00178 {
00179 times_name = output_name + ".times";
00180 }
00181
00182 vector<string> timestamps;
00183 vector<pair<string,string> > associations = load_associations(association_name, timestamps, skip_nb);
00184 if (debug)
00185 {
00186 for (size_t i = 0; i < associations.size(); i++)
00187 {
00188 cout << "associations : " << associations[i].first << " <-> " << associations[i].second << endl;
00189 }
00190 cout << "loaded : " << associations.size() << " associations" << endl;
00191 }
00192
00193 double fx, fy, cx, cy, ds;
00194 std::vector<double> dist(5);
00195 switch (freiburg_camera_nb)
00196 {
00197 case 1:
00198 fx = 517.3;
00199 fy = 516.5;
00200 cx = 318.6;
00201 cy = 255.3;
00202 dist[0] = 0.2624;
00203 dist[1] = -0.9531;
00204 dist[2] = -0.0054;
00205 dist[3] = 0.0026;
00206 dist[4] = 1.1633;
00207 ds = 1.035;
00208 break;
00209 case 2:
00210 fx = 520.9;
00211 fy = 521.0;
00212 cx = 325.1;
00213 cy = 249.7;
00214 dist[0] = 0.2312;
00215 dist[1] = -0.7849;
00216 dist[2] = -0.0033;
00217 dist[3] = -0.0001;
00218 dist[4] = 0.9172;
00219 ds = 1.031;
00220 break;
00221 default:
00222 cerr << "unknown camera number : " << freiburg_camera_nb << endl;
00223 return 1;
00224 }
00225
00226
00227 double t1, t2;
00228
00229 lslgeneric::DepthCamera<pcl::PointXYZ> cameraparams(fx,fy,cx,cy,dist,ds*scale,false);
00230
00231 typedef Eigen::Transform<double,3,Eigen::Affine,Eigen::ColMajor> EigenTransform;
00232 std::vector<EigenTransform, Eigen::aligned_allocator<EigenTransform> > transformVector;
00233
00234 std::ofstream file_times(times_name.c_str(), ios::out);
00235 if (!file_times)
00236 {
00237 cerr << "Failed to open times file : " << times_name << endl;
00238 return -1;
00239 }
00240
00241 file_times<<"times = [";
00242
00243
00244 if (old_reg_code)
00245 {
00246 for (size_t i = 1; i < associations.size(); i++)
00247 {
00248 std::cout << "iter " << i << "(" << associations.size() << ")" << std::endl;
00249
00250 NDTFrame<pcl::PointXYZ> *frame1 = new NDTFrame<pcl::PointXYZ>();
00251 frame1->img = cv::imread(associations[i-1].first, 0);
00252 frame1->depth_img = cv::imread(associations[i-1].second, CV_LOAD_IMAGE_ANYDEPTH);
00253 if (i == 1)
00254 cameraparams.setupDepthPointCloudLookUpTable(frame1->depth_img.size());
00255 frame1->cameraParams = cameraparams;
00256 frame1->supportSize = support_size;
00257
00258 NDTFrame<pcl::PointXYZ> *frame2 = new NDTFrame<pcl::PointXYZ>();
00259 frame2->img = cv::imread(associations[i].first, 0);
00260 frame2->depth_img = cv::imread(associations[i].second, CV_LOAD_IMAGE_ANYDEPTH);
00261 frame2->cameraParams = cameraparams;
00262 frame2->supportSize = support_size;
00263
00264 frame1->current_res = current_res;
00265 frame2->current_res = current_res;
00266 NDTFrameProc<pcl::PointXYZ> proc(nb_ransac, max_inldist_xy, max_inldist_z);
00267 proc.detector = cv::FeatureDetector::create("SURF");
00268
00269 proc.img_scale = img_scale;
00270
00271 double t1 = getDoubleTime();
00272 proc.addFrame(frame1);
00273 proc.addFrame(frame2);
00274 proc.processFrames(skip_matching,estimate_di,match_full,match_no_association);
00275 double t2 = getDoubleTime();
00276
00277 if (visualize)
00278 viewKeypointMatches(&proc, delay_vis);
00279
00280 if (i == 1)
00281 transformVector.push_back(proc.transformVector[0]);
00282 transformVector.push_back(proc.transformVector[1]);
00283
00284 file_times << t2-t1 <<", ";
00285
00286
00287
00288 }
00289 }
00290 else
00291 {
00292
00293
00294 NDTFrameProc<pcl::PointXYZ> proc(nb_ransac, max_inldist_xy, max_inldist_z);
00295 proc.detector = cv::FeatureDetector::create("SURF");
00296
00297 proc.img_scale = img_scale;
00298 proc.trim_factor = trim_factor;
00299 proc.non_mean = non_mean;
00300 proc.pe.windowed = windowed_matching;
00301 proc.pe.maxDist = max_kp_dist;
00302 proc.pe.minDist = min_kp_dist;
00303
00304 for (size_t i = 0; i < associations.size(); i++)
00305 {
00306 std::cout << "iter " << i << "(" << associations.size() << ")" << std::endl;
00307 NDTFrame<pcl::PointXYZ> *frame = new NDTFrame<pcl::PointXYZ>();
00308 frame->img = cv::imread(associations[i].first, 0);
00309 frame->depth_img = cv::imread(associations[i].second, CV_LOAD_IMAGE_ANYDEPTH);
00310 frame->supportSize = support_size;
00311 frame->maxVar = max_var;
00312 frame->current_res = current_res;
00313 if (i == 0)
00314 cameraparams.setupDepthPointCloudLookUpTable(frame->depth_img.size());
00315 frame->cameraParams = cameraparams;
00316
00317
00318 double t1 = getDoubleTime();
00319 proc.addFrameIncremental(frame, skip_matching, estimate_di,match_full,match_no_association);
00320 double t2 = getDoubleTime();
00321
00322 if (visualize)
00323 viewKeypointMatches(&proc, delay_vis);
00324
00325
00326 file_times << t2-t1 <<", ";
00327 }
00328 transformVector = proc.transformVector;
00329 }
00330 file_times<<"];\n";
00331
00332
00333 cerr<<"sizes: "<<associations.size()<<" "<<transformVector.size()<<std::endl;
00334 assert(associations.size() == transformVector.size());
00335 Eigen::Transform<double,3,Eigen::Affine,Eigen::ColMajor> global_transform;
00336 global_transform.setIdentity();
00337
00338
00339 std::ofstream file(output_name.c_str(), ios::out);
00340 if (!file)
00341 {
00342 cerr << "Failed to create file : " << output_name << endl;
00343 return -1;
00344 }
00345
00346 for (size_t i = 0; i < transformVector.size(); i++)
00347 {
00348 global_transform = global_transform * transformVector[i];
00349 Eigen::Quaternion<double> tmp(global_transform.rotation());
00350 cout << timestamps[i] << " " << global_transform.translation().transpose() << " " << tmp.x() << " " << tmp.y() << " " << tmp.z() << " " << tmp.w() << endl;
00351 file << timestamps[i] << " " << global_transform.translation().transpose() << " " << tmp.x() << " " << tmp.y() << " " << tmp.z() << " " << tmp.w() << endl;
00352 }
00353 file.close();
00354
00355 return 0;
00356 }