1 #include <string>
2 #include <sstream>
3 #include <iostream>
4 #include <fstream>
5 #include <cstdlib>
6 #include "breg3d_gdbicp_homography_generator.h"
7 
8 
9 #ifdef _MSC_VER
10 #  include "vcl_msvc_warnings.h"
11 #endif
12 
13 #include "vul/vul_file.h"
14 #include "vil/vil_image_view.h"
15 #include "vil/vil_save.h"
16 #include "vil/vil_convert.h"
17 #include <ihog/ihog_transform_2d.h>
18 
19 #include "vnl/vnl_matrix.h"
20 #include "vnl/vnl_matrix_fixed.h"
21 
compute_homography()22 ihog_transform_2d breg3d_gdbicp_homography_generator::compute_homography()
23 {
24   // this command must be in path.
25   std::string gdbicp_command("gdbicp");
26 
27   // create local dir for i/o
28   std::string io_dirname("./temp_gdbicp_io");
29   if (!vul_file::exists(io_dirname)) {
30     std::cout << "creating gdbicp io directory " << io_dirname << std::endl;
31     if (!vul_file::make_directory(io_dirname)) {
32       std::cerr << "error creating directory!\n";
33       return ihog_transform_2d();
34     }
35   }
36   // move to io directory
37   std::string og_dir = vul_file::get_cwd();
38   vul_file::change_directory(io_dirname);
39 
40   // create input files
41   std::string img0_fname("img0.tiff");
42   std::string img1_fname("img1.tiff");
43   std::string mask0_fname("mask0.tiff");
44   std::string mask1_fname("mask1.tiff");
45 
46   // save image0
47   vil_image_view<vxl_byte> img0_byte(img0_->ni(),img0_->nj());
48   vil_convert_stretch_range_limited<float>(*img0_,img0_byte,0.0f,255.0f);
49   vil_save(img0_byte,img0_fname.c_str());
50 
51   // save image1
52   vil_image_view<vxl_byte> img1_byte(img1_->ni(),img1_->nj());
53   vil_convert_stretch_range_limited<float>(*img1_,img1_byte,0.0f,255.0f);
54   vil_save(img1_byte,img1_fname.c_str());
55 
56   if (use_mask0_) {
57     // save mask0
58     vil_image_view<vxl_byte> mask0_byte(mask0_->ni(),mask0_->nj());
59     vil_convert_stretch_range_limited<float>(*mask0_,mask0_byte,0.0f,1.0f);
60     vil_save(mask0_byte,mask0_fname.c_str());
61   }
62   if (use_mask1_) {
63     // save mask1
64     vil_image_view<vxl_byte> mask1_byte(mask1_->ni(),mask1_->nj());
65     vil_convert_stretch_range_limited<float>(*mask1_,mask1_byte,0.0f,1.0f);
66     vil_save(mask1_byte,mask1_fname.c_str());
67   }
68   std::stringstream command;
69   command << gdbicp_command << ' ' << img0_fname << ' ' << img1_fname << ' ';
70   if (use_mask0_) {
71     command << "-mask_from " << mask0_fname << ' ';
72   }
73   if (use_mask1_) {
74     command << "-mask_to " << mask1_fname << ' ';
75   }
76   if (compute_projective_) {
77     command << "-model 1 ";
78   }
79   else {
80     command << "-model 0 ";
81   }
82   //command << "-no_render";
83 
84   std::cout << "running " << command.str() << std::endl;
85   int retval = std::system(command.str().c_str());
86   if (retval) {
87     std::cerr << "std::system(" << command.str()
88              << ") returned " << retval << ".\n"
89              << "  make sure " << gdbicp_command << " is in your path.\n";
90     vul_file::change_directory(og_dir);
91     return ihog_transform_2d();
92   }
93 
94   // read in output file
95   std::string output_fname("mosaic_" + vul_file::strip_extension(img0_fname) +"_to_" + vul_file::strip_extension(img1_fname) + ".xform");
96   if (!vul_file::exists(output_fname)) {
97     std::cerr << "error: output file " << output_fname << " does not exist!\n";
98     vul_file::change_directory(og_dir);
99     return ihog_transform_2d();
100   }
101   ihog_transform_2d xform = parse_gdbicp_output(output_fname);
102 
103   vul_file::change_directory(og_dir);
104   return xform;
105 }
106 
parse_gdbicp_output(const std::string & filename)107 ihog_transform_2d breg3d_gdbicp_homography_generator::parse_gdbicp_output(const std::string& filename)
108 {
109   std::ifstream ifs(filename.c_str());
110 
111   char curr_char = 0;
112   while (curr_char != '<' && ifs.good()) {
113     ifs.read(&curr_char,1);
114   }
115   std::string htype;
116   ifs >> htype;
117   std::cout << "transformation type = " << htype << std::endl;
118   unsigned hdim = 0;
119   ifs >> hdim;
120   std::cout << "transformation dimension = " << hdim << std::endl;
121 
122   vnl_matrix_fixed<double,3,3> H_centered;
123   double cx,cy,dx,dy;
124   if (compute_projective_) {
125     // full projective transformation
126     ifs >> H_centered;
127     std::cout << "H_centered =\n" << H_centered << std::endl;
128     ifs >> cx >> cy >> dx >> dy;
129   }
130   else {
131     // affine transformation
132     vnl_matrix_fixed<double,2,2> A;
133     ifs >> A;
134     std::cout << "A =\n" << A << std::endl;
135     H_centered.fill(0.0);
136     H_centered.update(A,0,0);
137     H_centered(2,2) = 1.0;
138     std::cout << "H_centered =\n" << H_centered << std::endl;
139     ifs >> dx >> dy >> cx >> cy;
140   }
141   std::cout << "cx=" << cx << " cy=" << cy << " dx=" << dx << " dy=" << dy << std::endl;
142   vnl_matrix_fixed<double,3,3> C,D;
143   C.set_identity();
144   D.set_identity();
145   C(0,2) = -cx;
146   C(1,2) = -cy;
147   D(0,2) = dx;
148   D(1,2) = dy;
149   vnl_matrix_fixed<double,3,3> H = D*H_centered*C;
150   std::cout << "H = " << H << std::endl;
151 
152   ihog_transform_2d xform;
153   if (compute_projective_) {
154     xform.set_projective(H);
155   }
156   else {
157     vnl_matrix_fixed<double,2,3> H_affine;
158     for (unsigned i=0; i<2; ++i) {
159       for (unsigned j=0; j<3; ++j) {
160         H_affine(i,j) = H(i,j);
161       }
162     }
163     xform.set_affine(H_affine);
164   }
165   return xform;
166 }
167