1 // Copyright (c) 2018, ETH Zurich and UNC Chapel Hill.
2 // All rights reserved.
3 //
4 // Redistribution and use in source and binary forms, with or without
5 // modification, are permitted provided that the following conditions are met:
6 //
7 //     * Redistributions of source code must retain the above copyright
8 //       notice, this list of conditions and the following disclaimer.
9 //
10 //     * Redistributions in binary form must reproduce the above copyright
11 //       notice, this list of conditions and the following disclaimer in the
12 //       documentation and/or other materials provided with the distribution.
13 //
14 //     * Neither the name of ETH Zurich and UNC Chapel Hill nor the names of
15 //       its contributors may be used to endorse or promote products derived
16 //       from this software without specific prior written permission.
17 //
18 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
19 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
20 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
21 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
22 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
23 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
24 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
26 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
27 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
28 // POSSIBILITY OF SUCH DAMAGE.
29 //
30 // Author: Johannes L. Schoenberger (jsch-at-demuc-dot-de)
31 
32 #include "base/undistortion.h"
33 
34 #include <fstream>
35 
36 #include "base/camera_models.h"
37 #include "base/pose.h"
38 #include "base/warp.h"
39 #include "util/misc.h"
40 
41 namespace colmap {
42 namespace {
43 
44 template <typename Derived>
WriteMatrix(const Eigen::MatrixBase<Derived> & matrix,std::ofstream * file)45 void WriteMatrix(const Eigen::MatrixBase<Derived>& matrix,
46                  std::ofstream* file) {
47   typedef typename Eigen::MatrixBase<Derived>::Index index_t;
48   for (index_t r = 0; r < matrix.rows(); ++r) {
49     for (index_t c = 0; c < matrix.cols() - 1; ++c) {
50       *file << matrix(r, c) << " ";
51     }
52     *file << matrix(r, matrix.cols() - 1) << std::endl;
53   }
54 }
55 
56 // Write projection matrix P = K * [R t] to file and prepend given header.
WriteProjectionMatrix(const std::string & path,const Camera & camera,const Image & image,const std::string & header)57 void WriteProjectionMatrix(const std::string& path, const Camera& camera,
58                            const Image& image, const std::string& header) {
59   CHECK_EQ(camera.ModelId(), PinholeCameraModel::model_id);
60 
61   std::ofstream file(path, std::ios::trunc);
62   CHECK(file.is_open()) << path;
63 
64   Eigen::Matrix3d calib_matrix = Eigen::Matrix3d::Identity();
65   calib_matrix(0, 0) = camera.FocalLengthX();
66   calib_matrix(1, 1) = camera.FocalLengthY();
67   calib_matrix(0, 2) = camera.PrincipalPointX();
68   calib_matrix(1, 2) = camera.PrincipalPointY();
69 
70   const Eigen::Matrix3x4d proj_matrix = calib_matrix * image.ProjectionMatrix();
71 
72   if (!header.empty()) {
73     file << header << std::endl;
74   }
75 
76   WriteMatrix(proj_matrix, &file);
77 }
78 
WriteCOLMAPCommands(const bool geometric,const std::string & workspace_path,const std::string & workspace_format,const std::string & pmvs_option_name,const std::string & output_prefix,const std::string & indent,std::ofstream * file)79 void WriteCOLMAPCommands(const bool geometric,
80                          const std::string& workspace_path,
81                          const std::string& workspace_format,
82                          const std::string& pmvs_option_name,
83                          const std::string& output_prefix,
84                          const std::string& indent, std::ofstream* file) {
85   if (geometric) {
86     *file << indent << "$COLMAP_EXE_PATH/colmap patch_match_stereo \\" << std::endl;
87     *file << indent << "  --workspace_path " << workspace_path << " \\"
88           << std::endl;
89     *file << indent << "  --workspace_format " << workspace_format << " \\"
90           << std::endl;
91     if (workspace_format == "PMVS") {
92       *file << indent << "  --pmvs_option_name " << pmvs_option_name << " \\"
93             << std::endl;
94     }
95     *file << indent << "  --PatchMatchStereo.max_image_size 2000 \\"
96           << std::endl;
97     *file << indent << "  --PatchMatchStereo.geom_consistency true"
98           << std::endl;
99   } else {
100     *file << indent << "$COLMAP_EXE_PATH/colmap patch_match_stereo \\" << std::endl;
101     *file << indent << "  --workspace_path " << workspace_path << " \\"
102           << std::endl;
103     *file << indent << "  --workspace_format " << workspace_format << " \\"
104           << std::endl;
105     if (workspace_format == "PMVS") {
106       *file << indent << "  --pmvs_option_name " << pmvs_option_name << " \\"
107             << std::endl;
108     }
109     *file << indent << "  --PatchMatchStereo.max_image_size 2000 \\"
110           << std::endl;
111     *file << indent << "  --PatchMatchStereo.geom_consistency false"
112           << std::endl;
113   }
114 
115   *file << indent << "$COLMAP_EXE_PATH/colmap stereo_fusion \\" << std::endl;
116   *file << indent << "  --workspace_path " << workspace_path << " \\"
117         << std::endl;
118   *file << indent << "  --workspace_format " << workspace_format << " \\"
119         << std::endl;
120   if (workspace_format == "PMVS") {
121     *file << indent << "  --pmvs_option_name " << pmvs_option_name << " \\"
122           << std::endl;
123   }
124   if (geometric) {
125     *file << indent << "  --input_type geometric \\" << std::endl;
126   } else {
127     *file << indent << "  --input_type photometric \\" << std::endl;
128   }
129   *file << indent << "  --output_path "
130         << JoinPaths(workspace_path, output_prefix + "fused.ply") << std::endl;
131 
132   *file << indent << "$COLMAP_EXE_PATH/colmap poisson_mesher \\" << std::endl;
133   *file << indent << "  --input_path "
134         << JoinPaths(workspace_path, output_prefix + "fused.ply") << " \\"
135         << std::endl;
136   *file << indent << "  --output_path "
137         << JoinPaths(workspace_path, output_prefix + "meshed-poisson.ply")
138         << std::endl;
139 
140   *file << indent << "$COLMAP_EXE_PATH/colmap delaunay_mesher \\" << std::endl;
141   *file << indent << "  --input_path "
142         << JoinPaths(workspace_path, output_prefix) << " \\" << std::endl;
143   *file << indent << "  --input_type dense " << std::endl;
144   *file << indent << "  --output_path "
145         << JoinPaths(workspace_path, output_prefix + "meshed-delaunay.ply")
146         << std::endl;
147 }
148 
149 }  // namespace
150 
COLMAPUndistorter(const UndistortCameraOptions & options,const Reconstruction & reconstruction,const std::string & image_path,const std::string & output_path)151 COLMAPUndistorter::COLMAPUndistorter(const UndistortCameraOptions& options,
152                                      const Reconstruction& reconstruction,
153                                      const std::string& image_path,
154                                      const std::string& output_path)
155     : options_(options),
156       image_path_(image_path),
157       output_path_(output_path),
158       reconstruction_(reconstruction) {}
159 
Run()160 void COLMAPUndistorter::Run() {
161   PrintHeading1("Image undistortion");
162 
163   CreateDirIfNotExists(JoinPaths(output_path_, "images"));
164   CreateDirIfNotExists(JoinPaths(output_path_, "sparse"));
165   CreateDirIfNotExists(JoinPaths(output_path_, "stereo"));
166   CreateDirIfNotExists(JoinPaths(output_path_, "stereo/depth_maps"));
167   CreateDirIfNotExists(JoinPaths(output_path_, "stereo/normal_maps"));
168   CreateDirIfNotExists(JoinPaths(output_path_, "stereo/consistency_graphs"));
169   reconstruction_.CreateImageDirs(JoinPaths(output_path_, "images"));
170   reconstruction_.CreateImageDirs(JoinPaths(output_path_, "stereo/depth_maps"));
171   reconstruction_.CreateImageDirs(
172       JoinPaths(output_path_, "stereo/normal_maps"));
173   reconstruction_.CreateImageDirs(
174       JoinPaths(output_path_, "stereo/consistency_graphs"));
175 
176   ThreadPool thread_pool;
177   std::vector<std::future<void>> futures;
178   futures.reserve(reconstruction_.NumRegImages());
179   for (size_t i = 0; i < reconstruction_.NumRegImages(); ++i) {
180     futures.push_back(
181         thread_pool.AddTask(&COLMAPUndistorter::Undistort, this, i));
182   }
183 
184   for (size_t i = 0; i < futures.size(); ++i) {
185     if (IsStopped()) {
186       break;
187     }
188 
189     std::cout << StringPrintf("Undistorting image [%d/%d]", i + 1,
190                               futures.size())
191               << std::endl;
192 
193     futures[i].get();
194   }
195 
196   std::cout << "Writing reconstruction..." << std::endl;
197   Reconstruction undistorted_reconstruction = reconstruction_;
198   UndistortReconstruction(options_, &undistorted_reconstruction);
199   undistorted_reconstruction.Write(JoinPaths(output_path_, "sparse"));
200 
201   std::cout << "Writing configuration..." << std::endl;
202   WritePatchMatchConfig();
203   WriteFusionConfig();
204 
205   std::cout << "Writing scripts..." << std::endl;
206   WriteScript(false);
207   WriteScript(true);
208 
209   GetTimer().PrintMinutes();
210 }
211 
Undistort(const size_t reg_image_idx) const212 void COLMAPUndistorter::Undistort(const size_t reg_image_idx) const {
213   const image_t image_id = reconstruction_.RegImageIds().at(reg_image_idx);
214   const Image& image = reconstruction_.Image(image_id);
215   const Camera& camera = reconstruction_.Camera(image.CameraId());
216 
217   const std::string output_image_path =
218       JoinPaths(output_path_, "images", image.Name());
219 
220   Bitmap distorted_bitmap;
221   const std::string input_image_path = JoinPaths(image_path_, image.Name());
222   if (!distorted_bitmap.Read(input_image_path)) {
223     std::cerr << "ERROR: Cannot read image at path " << input_image_path
224               << std::endl;
225     return;
226   }
227 
228   Bitmap undistorted_bitmap;
229   Camera undistorted_camera;
230   UndistortImage(options_, distorted_bitmap, camera, &undistorted_bitmap,
231                  &undistorted_camera);
232 
233   undistorted_bitmap.Write(output_image_path);
234 }
235 
WritePatchMatchConfig() const236 void COLMAPUndistorter::WritePatchMatchConfig() const {
237   const auto path = JoinPaths(output_path_, "stereo/patch-match.cfg");
238   std::ofstream file(path, std::ios::trunc);
239   CHECK(file.is_open()) << path;
240   for (const auto image_id : reconstruction_.RegImageIds()) {
241     const auto& image = reconstruction_.Image(image_id);
242     file << image.Name() << std::endl;
243     file << "__auto__, 20" << std::endl;
244   }
245 }
246 
WriteFusionConfig() const247 void COLMAPUndistorter::WriteFusionConfig() const {
248   const auto path = JoinPaths(output_path_, "stereo/fusion.cfg");
249   std::ofstream file(path, std::ios::trunc);
250   CHECK(file.is_open()) << path;
251   for (const auto image_id : reconstruction_.RegImageIds()) {
252     const auto& image = reconstruction_.Image(image_id);
253     file << image.Name() << std::endl;
254   }
255 }
256 
WriteScript(const bool geometric) const257 void COLMAPUndistorter::WriteScript(const bool geometric) const {
258   const std::string path =
259       JoinPaths(output_path_, geometric ? "run-colmap-geometric.sh"
260                                         : "run-colmap-photometric.sh");
261   std::ofstream file(path, std::ios::trunc);
262   CHECK(file.is_open()) << path;
263 
264   file << "# You must set $COLMAP_EXE_PATH to " << std::endl
265        << "# the directory containing the COLMAP executables." << std::endl;
266   WriteCOLMAPCommands(geometric, ".", "COLMAP", "option-all", "", "", &file);
267 }
268 
PMVSUndistorter(const UndistortCameraOptions & options,const Reconstruction & reconstruction,const std::string & image_path,const std::string & output_path)269 PMVSUndistorter::PMVSUndistorter(const UndistortCameraOptions& options,
270                                  const Reconstruction& reconstruction,
271                                  const std::string& image_path,
272                                  const std::string& output_path)
273     : options_(options),
274       image_path_(image_path),
275       output_path_(output_path),
276       reconstruction_(reconstruction) {}
277 
Run()278 void PMVSUndistorter::Run() {
279   PrintHeading1("Image undistortion (CMVS/PMVS)");
280 
281   CreateDirIfNotExists(JoinPaths(output_path_, "pmvs"));
282   CreateDirIfNotExists(JoinPaths(output_path_, "pmvs/txt"));
283   CreateDirIfNotExists(JoinPaths(output_path_, "pmvs/visualize"));
284   CreateDirIfNotExists(JoinPaths(output_path_, "pmvs/models"));
285 
286   ThreadPool thread_pool;
287   std::vector<std::future<void>> futures;
288   futures.reserve(reconstruction_.NumRegImages());
289   for (size_t i = 0; i < reconstruction_.NumRegImages(); ++i) {
290     futures.push_back(
291         thread_pool.AddTask(&PMVSUndistorter::Undistort, this, i));
292   }
293 
294   for (size_t i = 0; i < futures.size(); ++i) {
295     if (IsStopped()) {
296       thread_pool.Stop();
297       std::cout << "WARNING: Stopped the undistortion process. Image point "
298                    "locations and camera parameters for not yet processed "
299                    "images in the Bundler output file is probably wrong."
300                 << std::endl;
301       break;
302     }
303 
304     std::cout << StringPrintf("Undistorting image [%d/%d]", i + 1,
305                               futures.size())
306               << std::endl;
307 
308     futures[i].get();
309   }
310 
311   std::cout << "Writing bundle file..." << std::endl;
312   Reconstruction undistorted_reconstruction = reconstruction_;
313   UndistortReconstruction(options_, &undistorted_reconstruction);
314   const std::string bundle_path = JoinPaths(output_path_, "pmvs/bundle.rd.out");
315   undistorted_reconstruction.ExportBundler(bundle_path,
316                                            bundle_path + ".list.txt");
317 
318   std::cout << "Writing visibility file..." << std::endl;
319   WriteVisibilityData();
320 
321   std::cout << "Writing option file..." << std::endl;
322   WriteOptionFile();
323 
324   std::cout << "Writing scripts..." << std::endl;
325   WritePMVSScript();
326   WriteCMVSPMVSScript();
327   WriteCOLMAPScript(false);
328   WriteCOLMAPScript(true);
329   WriteCMVSCOLMAPScript(false);
330   WriteCMVSCOLMAPScript(true);
331 
332   GetTimer().PrintMinutes();
333 }
334 
Undistort(const size_t reg_image_idx) const335 void PMVSUndistorter::Undistort(const size_t reg_image_idx) const {
336   const std::string output_image_path = JoinPaths(
337       output_path_, StringPrintf("pmvs/visualize/%08d.jpg", reg_image_idx));
338   const std::string proj_matrix_path =
339       JoinPaths(output_path_, StringPrintf("pmvs/txt/%08d.txt", reg_image_idx));
340 
341   const image_t image_id = reconstruction_.RegImageIds().at(reg_image_idx);
342   const Image& image = reconstruction_.Image(image_id);
343   const Camera& camera = reconstruction_.Camera(image.CameraId());
344 
345   Bitmap distorted_bitmap;
346   const std::string input_image_path = JoinPaths(image_path_, image.Name());
347   if (!distorted_bitmap.Read(input_image_path)) {
348     std::cerr << StringPrintf("ERROR: Cannot read image at path %s",
349                               input_image_path.c_str())
350               << std::endl;
351     return;
352   }
353 
354   Bitmap undistorted_bitmap;
355   Camera undistorted_camera;
356   UndistortImage(options_, distorted_bitmap, camera, &undistorted_bitmap,
357                  &undistorted_camera);
358 
359   undistorted_bitmap.Write(output_image_path);
360   WriteProjectionMatrix(proj_matrix_path, undistorted_camera, image, "CONTOUR");
361 }
362 
WriteVisibilityData() const363 void PMVSUndistorter::WriteVisibilityData() const {
364   const auto path = JoinPaths(output_path_, "pmvs/vis.dat");
365   std::ofstream file(path, std::ios::trunc);
366   CHECK(file.is_open()) << path;
367 
368   file << "VISDATA" << std::endl;
369   file << reconstruction_.NumRegImages() << std::endl;
370 
371   const std::vector<image_t>& reg_image_ids = reconstruction_.RegImageIds();
372 
373   for (size_t i = 0; i < reg_image_ids.size(); ++i) {
374     const image_t image_id = reg_image_ids[i];
375     const Image& image = reconstruction_.Image(image_id);
376     std::unordered_set<image_t> visible_image_ids;
377     for (point2D_t point2D_idx = 0; point2D_idx < image.NumPoints2D();
378          ++point2D_idx) {
379       const Point2D& point2D = image.Point2D(point2D_idx);
380       if (point2D.HasPoint3D()) {
381         const Point3D& point3D = reconstruction_.Point3D(point2D.Point3DId());
382         for (const TrackElement track_el : point3D.Track().Elements()) {
383           if (track_el.image_id != image_id) {
384             visible_image_ids.insert(track_el.image_id);
385           }
386         }
387       }
388     }
389 
390     std::vector<image_t> sorted_visible_image_ids(visible_image_ids.begin(),
391                                                   visible_image_ids.end());
392     std::sort(sorted_visible_image_ids.begin(), sorted_visible_image_ids.end());
393 
394     file << i << " " << visible_image_ids.size();
395     for (const image_t visible_image_id : sorted_visible_image_ids) {
396       file << " " << visible_image_id;
397     }
398     file << std::endl;
399   }
400 }
401 
WritePMVSScript() const402 void PMVSUndistorter::WritePMVSScript() const {
403   const auto path = JoinPaths(output_path_, "run-pmvs.sh");
404   std::ofstream file(path, std::ios::trunc);
405   CHECK(file.is_open()) << path;
406 
407   file << "# You must set $PMVS_EXE_PATH to " << std::endl
408        << "# the directory containing the CMVS-PMVS executables." << std::endl;
409   file << "$PMVS_EXE_PATH/pmvs2 pmvs/ option-all" << std::endl;
410 }
411 
WriteCMVSPMVSScript() const412 void PMVSUndistorter::WriteCMVSPMVSScript() const {
413   const auto path = JoinPaths(output_path_, "run-cmvs-pmvs.sh");
414   std::ofstream file(path, std::ios::trunc);
415   CHECK(file.is_open()) << path;
416 
417   file << "# You must set $PMVS_EXE_PATH to " << std::endl
418        << "# the directory containing the CMVS-PMVS executables." << std::endl;
419   file << "$PMVS_EXE_PATH/cmvs pmvs/" << std::endl;
420   file << "$PMVS_EXE_PATH/genOption pmvs/" << std::endl;
421   file << "find pmvs/ -iname \"option-*\" | sort | while read file_name"
422        << std::endl;
423   file << "do" << std::endl;
424   file << "    option_name=$(basename \"$file_name\")" << std::endl;
425   file << "    if [ \"$option_name\" = \"option-all\" ]; then" << std::endl;
426   file << "        continue" << std::endl;
427   file << "    fi" << std::endl;
428   file << "    $PMVS_EXE_PATH/pmvs2 pmvs/ $option_name" << std::endl;
429   file << "done" << std::endl;
430 }
431 
WriteCOLMAPScript(const bool geometric) const432 void PMVSUndistorter::WriteCOLMAPScript(const bool geometric) const {
433   const std::string path =
434       JoinPaths(output_path_, geometric ? "run-colmap-geometric.sh"
435                                         : "run-colmap-photometric.sh");
436   std::ofstream file(path, std::ios::trunc);
437   CHECK(file.is_open()) << path;
438 
439   file << "# You must set $COLMAP_EXE_PATH to " << std::endl
440        << "# the directory containing the COLMAP executables." << std::endl;
441   WriteCOLMAPCommands(geometric, "pmvs", "PMVS", "option-all", "option-all-",
442                       "", &file);
443 }
444 
WriteCMVSCOLMAPScript(const bool geometric) const445 void PMVSUndistorter::WriteCMVSCOLMAPScript(const bool geometric) const {
446   const std::string path =
447       JoinPaths(output_path_, geometric ? "run-cmvs-colmap-geometric.sh"
448                                         : "run-cmvs-colmap-photometric.sh");
449   std::ofstream file(path, std::ios::trunc);
450   CHECK(file.is_open()) << path;
451 
452   file << "# You must set $PMVS_EXE_PATH to " << std::endl
453        << "# the directory containing the CMVS-PMVS executables" << std::endl;
454   file << "# and you must set $COLMAP_EXE_PATH to " << std::endl
455        << "# the directory containing the COLMAP executables." << std::endl;
456   file << "$PMVS_EXE_PATH/cmvs pmvs/" << std::endl;
457   file << "$PMVS_EXE_PATH/genOption pmvs/" << std::endl;
458   file << "find pmvs/ -iname \"option-*\" | sort | while read file_name"
459        << std::endl;
460   file << "do" << std::endl;
461   file << "    workspace_path=$(dirname \"$file_name\")" << std::endl;
462   file << "    option_name=$(basename \"$file_name\")" << std::endl;
463   file << "    if [ \"$option_name\" = \"option-all\" ]; then" << std::endl;
464   file << "        continue" << std::endl;
465   file << "    fi" << std::endl;
466   file << "    rm -rf \"$workspace_path/stereo\"" << std::endl;
467   WriteCOLMAPCommands(geometric, "pmvs", "PMVS", "$option_name",
468                       "$option_name-", "    ", &file);
469   file << "done" << std::endl;
470 }
471 
WriteOptionFile() const472 void PMVSUndistorter::WriteOptionFile() const {
473   const auto path = JoinPaths(output_path_, "pmvs/option-all");
474   std::ofstream file(path, std::ios::trunc);
475   CHECK(file.is_open()) << path;
476 
477   file << "# Generated by COLMAP - all images, no clustering." << std::endl;
478 
479   file << "level 1" << std::endl;
480   file << "csize 2" << std::endl;
481   file << "threshold 0.7" << std::endl;
482   file << "wsize 7" << std::endl;
483   file << "minImageNum 3" << std::endl;
484   file << "CPU " << std::thread::hardware_concurrency() << std::endl;
485   file << "setEdge 0" << std::endl;
486   file << "useBound 0" << std::endl;
487   file << "useVisData 1" << std::endl;
488   file << "sequence -1" << std::endl;
489   file << "maxAngle 10" << std::endl;
490   file << "quad 2.0" << std::endl;
491 
492   file << "timages " << reconstruction_.NumRegImages();
493   for (size_t i = 0; i < reconstruction_.NumRegImages(); ++i) {
494     file << " " << i;
495   }
496   file << std::endl;
497 
498   file << "oimages 0" << std::endl;
499 }
500 
CMPMVSUndistorter(const UndistortCameraOptions & options,const Reconstruction & reconstruction,const std::string & image_path,const std::string & output_path)501 CMPMVSUndistorter::CMPMVSUndistorter(const UndistortCameraOptions& options,
502                                      const Reconstruction& reconstruction,
503                                      const std::string& image_path,
504                                      const std::string& output_path)
505     : options_(options),
506       image_path_(image_path),
507       output_path_(output_path),
508       reconstruction_(reconstruction) {}
509 
Run()510 void CMPMVSUndistorter::Run() {
511   PrintHeading1("Image undistortion (CMP-MVS)");
512 
513   ThreadPool thread_pool;
514   std::vector<std::future<void>> futures;
515   futures.reserve(reconstruction_.NumRegImages());
516   for (size_t i = 0; i < reconstruction_.NumRegImages(); ++i) {
517     futures.push_back(
518         thread_pool.AddTask(&CMPMVSUndistorter::Undistort, this, i));
519   }
520 
521   for (size_t i = 0; i < futures.size(); ++i) {
522     if (IsStopped()) {
523       break;
524     }
525 
526     std::cout << StringPrintf("Undistorting image [%d/%d]", i + 1,
527                               futures.size())
528               << std::endl;
529 
530     futures[i].get();
531   }
532 
533   GetTimer().PrintMinutes();
534 }
535 
Undistort(const size_t reg_image_idx) const536 void CMPMVSUndistorter::Undistort(const size_t reg_image_idx) const {
537   const std::string output_image_path =
538       JoinPaths(output_path_, StringPrintf("%05d.jpg", reg_image_idx + 1));
539   const std::string proj_matrix_path =
540       JoinPaths(output_path_, StringPrintf("%05d_P.txt", reg_image_idx + 1));
541 
542   const image_t image_id = reconstruction_.RegImageIds().at(reg_image_idx);
543   const Image& image = reconstruction_.Image(image_id);
544   const Camera& camera = reconstruction_.Camera(image.CameraId());
545 
546   Bitmap distorted_bitmap;
547   const std::string input_image_path = JoinPaths(image_path_, image.Name());
548   if (!distorted_bitmap.Read(input_image_path)) {
549     std::cerr << "ERROR: Cannot read image at path " << input_image_path
550               << std::endl;
551     return;
552   }
553 
554   Bitmap undistorted_bitmap;
555   Camera undistorted_camera;
556   UndistortImage(options_, distorted_bitmap, camera, &undistorted_bitmap,
557                  &undistorted_camera);
558 
559   undistorted_bitmap.Write(output_image_path);
560   WriteProjectionMatrix(proj_matrix_path, undistorted_camera, image, "CONTOUR");
561 }
562 
PureImageUndistorter(const UndistortCameraOptions & options,const std::string & image_path,const std::string & output_path,const std::vector<std::pair<std::string,Camera>> & image_names_and_cameras)563 PureImageUndistorter::PureImageUndistorter(
564   const UndistortCameraOptions& options, const std::string& image_path,
565   const std::string& output_path,
566   const std::vector<std::pair<std::string, Camera>>& image_names_and_cameras)
567     : options_(options),
568       image_path_(image_path),
569       output_path_(output_path),
570       image_names_and_cameras_(image_names_and_cameras) {}
571 
Run()572 void PureImageUndistorter::Run() {
573   PrintHeading1("Image undistortion");
574 
575   CreateDirIfNotExists(output_path_);
576 
577   ThreadPool thread_pool;
578   std::vector<std::future<void>> futures;
579   size_t num_images = image_names_and_cameras_.size();
580   futures.reserve(num_images);
581   for (size_t i = 0; i < num_images; ++i) {
582     futures.push_back(thread_pool.AddTask(&PureImageUndistorter::Undistort,
583                                           this, i));
584   }
585 
586   for (size_t i = 0; i < futures.size(); ++i) {
587     if (IsStopped()) {
588       break;
589     }
590 
591     std::cout << StringPrintf("Undistorting image [%d/%d]", i + 1,
592                               futures.size())
593               << std::endl;
594 
595     futures[i].get();
596   }
597 
598   GetTimer().PrintMinutes();
599 }
600 
Undistort(const size_t image_idx) const601 void PureImageUndistorter::Undistort(const size_t image_idx) const {
602   const std::string& image_name = image_names_and_cameras_[image_idx].first;
603   const Camera& camera = image_names_and_cameras_[image_idx].second;
604 
605   const std::string output_image_path =
606   JoinPaths(output_path_, image_name);
607 
608   Bitmap distorted_bitmap;
609   const std::string input_image_path = JoinPaths(image_path_, image_name);
610   if (!distorted_bitmap.Read(input_image_path)) {
611     std::cerr << "ERROR: Cannot read image at path " << input_image_path
612               << std::endl;
613     return;
614   }
615 
616   Bitmap undistorted_bitmap;
617   Camera undistorted_camera;
618   UndistortImage(options_, distorted_bitmap, camera, &undistorted_bitmap,
619                  &undistorted_camera);
620 
621   undistorted_bitmap.Write(output_image_path);
622 }
623 
StereoImageRectifier(const UndistortCameraOptions & options,const Reconstruction & reconstruction,const std::string & image_path,const std::string & output_path,const std::vector<std::pair<image_t,image_t>> & stereo_pairs)624 StereoImageRectifier::StereoImageRectifier(
625     const UndistortCameraOptions& options, const Reconstruction& reconstruction,
626     const std::string& image_path, const std::string& output_path,
627     const std::vector<std::pair<image_t, image_t>>& stereo_pairs)
628     : options_(options),
629       image_path_(image_path),
630       output_path_(output_path),
631       stereo_pairs_(stereo_pairs),
632       reconstruction_(reconstruction) {}
633 
Run()634 void StereoImageRectifier::Run() {
635   PrintHeading1("Stereo rectification");
636 
637   ThreadPool thread_pool;
638   std::vector<std::future<void>> futures;
639   futures.reserve(stereo_pairs_.size());
640   for (const auto& stereo_pair : stereo_pairs_) {
641     futures.push_back(thread_pool.AddTask(&StereoImageRectifier::Rectify, this,
642                                           stereo_pair.first,
643                                           stereo_pair.second));
644   }
645 
646   for (size_t i = 0; i < futures.size(); ++i) {
647     if (IsStopped()) {
648       break;
649     }
650 
651     std::cout << StringPrintf("Rectifying image pair [%d/%d]", i + 1,
652                               futures.size())
653               << std::endl;
654 
655     futures[i].get();
656   }
657 
658   GetTimer().PrintMinutes();
659 }
660 
Rectify(const image_t image_id1,const image_t image_id2) const661 void StereoImageRectifier::Rectify(const image_t image_id1,
662                                    const image_t image_id2) const {
663   const Image& image1 = reconstruction_.Image(image_id1);
664   const Image& image2 = reconstruction_.Image(image_id2);
665   const Camera& camera1 = reconstruction_.Camera(image1.CameraId());
666   const Camera& camera2 = reconstruction_.Camera(image2.CameraId());
667 
668   const std::string image_name1 = StringReplace(image1.Name(), "/", "-");
669   const std::string image_name2 = StringReplace(image2.Name(), "/", "-");
670 
671   const std::string stereo_pair_name =
672       StringPrintf("%s-%s", image_name1.c_str(), image_name2.c_str());
673 
674   CreateDirIfNotExists(JoinPaths(output_path_, stereo_pair_name));
675 
676   const std::string output_image1_path =
677       JoinPaths(output_path_, stereo_pair_name, image_name1);
678   const std::string output_image2_path =
679       JoinPaths(output_path_, stereo_pair_name, image_name2);
680 
681   Bitmap distorted_bitmap1;
682   const std::string input_image1_path = JoinPaths(image_path_, image1.Name());
683   if (!distorted_bitmap1.Read(input_image1_path)) {
684     std::cerr << "ERROR: Cannot read image at path " << input_image1_path
685               << std::endl;
686     return;
687   }
688 
689   Bitmap distorted_bitmap2;
690   const std::string input_image2_path = JoinPaths(image_path_, image2.Name());
691   if (!distorted_bitmap2.Read(input_image2_path)) {
692     std::cerr << "ERROR: Cannot read image at path " << input_image2_path
693               << std::endl;
694     return;
695   }
696 
697   Eigen::Vector4d qvec;
698   Eigen::Vector3d tvec;
699   ComputeRelativePose(image1.Qvec(), image1.Tvec(), image2.Qvec(),
700                       image2.Tvec(), &qvec, &tvec);
701 
702   Bitmap undistorted_bitmap1;
703   Bitmap undistorted_bitmap2;
704   Camera undistorted_camera;
705   Eigen::Matrix4d Q;
706   RectifyAndUndistortStereoImages(
707       options_, distorted_bitmap1, distorted_bitmap2, camera1, camera2, qvec,
708       tvec, &undistorted_bitmap1, &undistorted_bitmap2, &undistorted_camera,
709       &Q);
710 
711   undistorted_bitmap1.Write(output_image1_path);
712   undistorted_bitmap2.Write(output_image2_path);
713 
714   const auto Q_path = JoinPaths(output_path_, stereo_pair_name, "Q.txt");
715   std::ofstream Q_file(Q_path, std::ios::trunc);
716   CHECK(Q_file.is_open()) << Q_path;
717   WriteMatrix(Q, &Q_file);
718 }
719 
UndistortCamera(const UndistortCameraOptions & options,const Camera & camera)720 Camera UndistortCamera(const UndistortCameraOptions& options,
721                        const Camera& camera) {
722   CHECK_GE(options.blank_pixels, 0);
723   CHECK_LE(options.blank_pixels, 1);
724   CHECK_GT(options.min_scale, 0.0);
725   CHECK_LE(options.min_scale, options.max_scale);
726   CHECK_NE(options.max_image_size, 0);
727   CHECK_GE(options.roi_min_x, 0.0);
728   CHECK_GE(options.roi_min_y, 0.0);
729   CHECK_LE(options.roi_max_x, 1.0);
730   CHECK_LE(options.roi_max_y, 1.0);
731   CHECK_LT(options.roi_min_x, options.roi_max_x);
732   CHECK_LT(options.roi_min_y, options.roi_max_y);
733 
734   Camera undistorted_camera;
735   undistorted_camera.SetModelId(PinholeCameraModel::model_id);
736   undistorted_camera.SetWidth(camera.Width());
737   undistorted_camera.SetHeight(camera.Height());
738 
739   // Copy focal length parameters.
740   const std::vector<size_t>& focal_length_idxs = camera.FocalLengthIdxs();
741   CHECK_LE(focal_length_idxs.size(), 2)
742       << "Not more than two focal length parameters supported.";
743   if (focal_length_idxs.size() == 1) {
744     undistorted_camera.SetFocalLengthX(camera.FocalLength());
745     undistorted_camera.SetFocalLengthY(camera.FocalLength());
746   } else if (focal_length_idxs.size() == 2) {
747     undistorted_camera.SetFocalLengthX(camera.FocalLengthX());
748     undistorted_camera.SetFocalLengthY(camera.FocalLengthY());
749   }
750 
751   // Copy principal point parameters.
752   undistorted_camera.SetPrincipalPointX(camera.PrincipalPointX());
753   undistorted_camera.SetPrincipalPointY(camera.PrincipalPointY());
754 
755   // Modify undistorted camera parameters based on ROI if enabled
756   size_t roi_min_x = 0;
757   size_t roi_min_y = 0;
758   size_t roi_max_x = camera.Width();
759   size_t roi_max_y = camera.Height();
760 
761   const bool roi_enabled = options.roi_min_x > 0.0 || options.roi_min_y > 0.0 ||
762                            options.roi_max_x < 1.0 || options.roi_max_y < 1.0;
763 
764   if (roi_enabled) {
765     roi_min_x = static_cast<size_t>(
766         std::round(options.roi_min_x * static_cast<double>(camera.Width())));
767     roi_min_y = static_cast<size_t>(
768         std::round(options.roi_min_y * static_cast<double>(camera.Height())));
769     roi_max_x = static_cast<size_t>(
770         std::round(options.roi_max_x * static_cast<double>(camera.Width())));
771     roi_max_y = static_cast<size_t>(
772         std::round(options.roi_max_y * static_cast<double>(camera.Height())));
773 
774     // Make sure that the roi is valid.
775     roi_min_x = std::min(roi_min_x, camera.Width() - 1);
776     roi_min_y = std::min(roi_min_y, camera.Height() - 1);
777     roi_max_x = std::max(roi_max_x, roi_min_x + 1);
778     roi_max_y = std::max(roi_max_y, roi_min_y + 1);
779 
780     undistorted_camera.SetWidth(roi_max_x - roi_min_x);
781     undistorted_camera.SetHeight(roi_max_y - roi_min_y);
782 
783     undistorted_camera.SetPrincipalPointX(camera.PrincipalPointX() -
784                                           static_cast<double>(roi_min_x));
785     undistorted_camera.SetPrincipalPointY(camera.PrincipalPointY() -
786                                           static_cast<double>(roi_min_y));
787   }
788 
789   // Scale the image such the the boundary of the undistorted image.
790   if (roi_enabled || (camera.ModelId() != SimplePinholeCameraModel::model_id &&
791                       camera.ModelId() != PinholeCameraModel::model_id)) {
792     // Determine min/max coordinates along top / bottom image border.
793 
794     double left_min_x = std::numeric_limits<double>::max();
795     double left_max_x = std::numeric_limits<double>::lowest();
796     double right_min_x = std::numeric_limits<double>::max();
797     double right_max_x = std::numeric_limits<double>::lowest();
798 
799     for (size_t y = roi_min_y; y < roi_max_y; ++y) {
800       // Left border.
801       const Eigen::Vector2d world_point1 =
802           camera.ImageToWorld(Eigen::Vector2d(0.5, y + 0.5));
803       const Eigen::Vector2d undistorted_point1 =
804           undistorted_camera.WorldToImage(world_point1);
805       left_min_x = std::min(left_min_x, undistorted_point1(0));
806       left_max_x = std::max(left_max_x, undistorted_point1(0));
807       // Right border.
808       const Eigen::Vector2d world_point2 =
809           camera.ImageToWorld(Eigen::Vector2d(camera.Width() - 0.5, y + 0.5));
810       const Eigen::Vector2d undistorted_point2 =
811           undistorted_camera.WorldToImage(world_point2);
812       right_min_x = std::min(right_min_x, undistorted_point2(0));
813       right_max_x = std::max(right_max_x, undistorted_point2(0));
814     }
815 
816     // Determine min, max coordinates along left / right image border.
817 
818     double top_min_y = std::numeric_limits<double>::max();
819     double top_max_y = std::numeric_limits<double>::lowest();
820     double bottom_min_y = std::numeric_limits<double>::max();
821     double bottom_max_y = std::numeric_limits<double>::lowest();
822 
823     for (size_t x = roi_min_x; x < roi_max_x; ++x) {
824       // Top border.
825       const Eigen::Vector2d world_point1 =
826           camera.ImageToWorld(Eigen::Vector2d(x + 0.5, 0.5));
827       const Eigen::Vector2d undistorted_point1 =
828           undistorted_camera.WorldToImage(world_point1);
829       top_min_y = std::min(top_min_y, undistorted_point1(1));
830       top_max_y = std::max(top_max_y, undistorted_point1(1));
831       // Bottom border.
832       const Eigen::Vector2d world_point2 =
833           camera.ImageToWorld(Eigen::Vector2d(x + 0.5, camera.Height() - 0.5));
834       const Eigen::Vector2d undistorted_point2 =
835           undistorted_camera.WorldToImage(world_point2);
836       bottom_min_y = std::min(bottom_min_y, undistorted_point2(1));
837       bottom_max_y = std::max(bottom_max_y, undistorted_point2(1));
838     }
839 
840     const double cx = undistorted_camera.PrincipalPointX();
841     const double cy = undistorted_camera.PrincipalPointY();
842 
843     // Scale such that undistorted image contains all pixels of distorted image.
844     const double min_scale_x =
845         std::min(cx / (cx - left_min_x),
846                  (undistorted_camera.Width() - 0.5 - cx) / (right_max_x - cx));
847     const double min_scale_y = std::min(
848         cy / (cy - top_min_y),
849         (undistorted_camera.Height() - 0.5 - cy) / (bottom_max_y - cy));
850 
851     // Scale such that there are no blank pixels in undistorted image.
852     const double max_scale_x =
853         std::max(cx / (cx - left_max_x),
854                  (undistorted_camera.Width() - 0.5 - cx) / (right_min_x - cx));
855     const double max_scale_y = std::max(
856         cy / (cy - top_max_y),
857         (undistorted_camera.Height() - 0.5 - cy) / (bottom_min_y - cy));
858 
859     // Interpolate scale according to blank_pixels.
860     double scale_x = 1.0 / (min_scale_x * options.blank_pixels +
861                             max_scale_x * (1.0 - options.blank_pixels));
862     double scale_y = 1.0 / (min_scale_y * options.blank_pixels +
863                             max_scale_y * (1.0 - options.blank_pixels));
864 
865     // Clip the scaling factors.
866     scale_x = Clip(scale_x, options.min_scale, options.max_scale);
867     scale_y = Clip(scale_y, options.min_scale, options.max_scale);
868 
869     // Scale undistorted camera dimensions.
870     const size_t orig_undistorted_camera_width = undistorted_camera.Width();
871     const size_t orig_undistorted_camera_height = undistorted_camera.Height();
872     undistorted_camera.SetWidth(static_cast<size_t>(
873         std::max(1.0, scale_x * undistorted_camera.Width())));
874     undistorted_camera.SetHeight(static_cast<size_t>(
875         std::max(1.0, scale_y * undistorted_camera.Height())));
876 
877     // Scale the principal point according to the new dimensions of the camera.
878     undistorted_camera.SetPrincipalPointX(
879         undistorted_camera.PrincipalPointX() *
880         static_cast<double>(undistorted_camera.Width()) /
881         static_cast<double>(orig_undistorted_camera_width));
882     undistorted_camera.SetPrincipalPointY(
883         undistorted_camera.PrincipalPointY() *
884         static_cast<double>(undistorted_camera.Height()) /
885         static_cast<double>(orig_undistorted_camera_height));
886   }
887 
888   if (options.max_image_size > 0) {
889     const double max_image_scale_x =
890         options.max_image_size /
891         static_cast<double>(undistorted_camera.Width());
892     const double max_image_scale_y =
893         options.max_image_size /
894         static_cast<double>(undistorted_camera.Height());
895     const double max_image_scale =
896         std::min(max_image_scale_x, max_image_scale_y);
897     if (max_image_scale < 1.0) {
898       undistorted_camera.Rescale(max_image_scale);
899     }
900   }
901 
902   return undistorted_camera;
903 }
904 
UndistortImage(const UndistortCameraOptions & options,const Bitmap & distorted_bitmap,const Camera & distorted_camera,Bitmap * undistorted_bitmap,Camera * undistorted_camera)905 void UndistortImage(const UndistortCameraOptions& options,
906                     const Bitmap& distorted_bitmap,
907                     const Camera& distorted_camera, Bitmap* undistorted_bitmap,
908                     Camera* undistorted_camera) {
909   CHECK_EQ(distorted_camera.Width(), distorted_bitmap.Width());
910   CHECK_EQ(distorted_camera.Height(), distorted_bitmap.Height());
911 
912   *undistorted_camera = UndistortCamera(options, distorted_camera);
913   undistorted_bitmap->Allocate(static_cast<int>(undistorted_camera->Width()),
914                                static_cast<int>(undistorted_camera->Height()),
915                                distorted_bitmap.IsRGB());
916   distorted_bitmap.CloneMetadata(undistorted_bitmap);
917 
918   WarpImageBetweenCameras(distorted_camera, *undistorted_camera,
919                           distorted_bitmap, undistorted_bitmap);
920 }
921 
UndistortReconstruction(const UndistortCameraOptions & options,Reconstruction * reconstruction)922 void UndistortReconstruction(const UndistortCameraOptions& options,
923                              Reconstruction* reconstruction) {
924   const auto distorted_cameras = reconstruction->Cameras();
925   for (auto& camera : distorted_cameras) {
926     reconstruction->Camera(camera.first) =
927         UndistortCamera(options, camera.second);
928   }
929 
930   for (const auto& distorted_image : reconstruction->Images()) {
931     auto& image = reconstruction->Image(distorted_image.first);
932     const auto& distorted_camera = distorted_cameras.at(image.CameraId());
933     const auto& undistorted_camera = reconstruction->Camera(image.CameraId());
934     for (point2D_t point2D_idx = 0; point2D_idx < image.NumPoints2D();
935          ++point2D_idx) {
936       auto& point2D = image.Point2D(point2D_idx);
937       point2D.SetXY(undistorted_camera.WorldToImage(
938           distorted_camera.ImageToWorld(point2D.XY())));
939     }
940   }
941 }
942 
RectifyStereoCameras(const Camera & camera1,const Camera & camera2,const Eigen::Vector4d & qvec,const Eigen::Vector3d & tvec,Eigen::Matrix3d * H1,Eigen::Matrix3d * H2,Eigen::Matrix4d * Q)943 void RectifyStereoCameras(const Camera& camera1, const Camera& camera2,
944                           const Eigen::Vector4d& qvec,
945                           const Eigen::Vector3d& tvec, Eigen::Matrix3d* H1,
946                           Eigen::Matrix3d* H2, Eigen::Matrix4d* Q) {
947   CHECK(camera1.ModelId() == SimplePinholeCameraModel::model_id ||
948         camera1.ModelId() == PinholeCameraModel::model_id);
949   CHECK(camera2.ModelId() == SimplePinholeCameraModel::model_id ||
950         camera2.ModelId() == PinholeCameraModel::model_id);
951 
952   // Compute the average rotation between the first and the second camera.
953   Eigen::AngleAxisd rvec(
954       Eigen::Quaterniond(qvec(0), qvec(1), qvec(2), qvec(3)));
955   rvec.angle() *= -0.5;
956 
957   Eigen::Matrix3d R2 = rvec.toRotationMatrix();
958   Eigen::Matrix3d R1 = R2.transpose();
959 
960   // Determine the translation, such that it coincides with the X-axis.
961   Eigen::Vector3d t = R2 * tvec;
962 
963   Eigen::Vector3d x_unit_vector(1, 0, 0);
964   if (t.transpose() * x_unit_vector < 0) {
965     x_unit_vector *= -1;
966   }
967 
968   const Eigen::Vector3d rotation_axis = t.cross(x_unit_vector);
969 
970   Eigen::Matrix3d R_x;
971   if (rotation_axis.norm() < std::numeric_limits<double>::epsilon()) {
972     R_x = Eigen::Matrix3d::Identity();
973   } else {
974     const double angle = std::acos(std::abs(t.transpose() * x_unit_vector) /
975                                    (t.norm() * x_unit_vector.norm()));
976     R_x = Eigen::AngleAxisd(angle, rotation_axis.normalized());
977   }
978 
979   // Apply the X-axis correction.
980   R1 = R_x * R1;
981   R2 = R_x * R2;
982   t = R_x * t;
983 
984   // Determine the intrinsic calibration matrix.
985   Eigen::Matrix3d K = Eigen::Matrix3d::Identity();
986   K(0, 0) = std::min(camera1.MeanFocalLength(), camera2.MeanFocalLength());
987   K(1, 1) = K(0, 0);
988   K(0, 2) = camera1.PrincipalPointX();
989   K(1, 2) = (camera1.PrincipalPointY() + camera2.PrincipalPointY()) / 2;
990 
991   // Compose the homographies.
992   *H1 = K * R1 * camera1.CalibrationMatrix().inverse();
993   *H2 = K * R2 * camera2.CalibrationMatrix().inverse();
994 
995   // Determine the inverse projection matrix that transforms disparity values
996   // to 3D world coordinates: [x, y, disparity, 1] * Q = [X, Y, Z, 1] * w.
997   *Q = Eigen::Matrix4d::Identity();
998   (*Q)(3, 0) = -K(1, 2);
999   (*Q)(3, 1) = -K(0, 2);
1000   (*Q)(3, 2) = K(0, 0);
1001   (*Q)(2, 3) = -1 / t(0);
1002   (*Q)(3, 3) = 0;
1003 }
1004 
RectifyAndUndistortStereoImages(const UndistortCameraOptions & options,const Bitmap & distorted_image1,const Bitmap & distorted_image2,const Camera & distorted_camera1,const Camera & distorted_camera2,const Eigen::Vector4d & qvec,const Eigen::Vector3d & tvec,Bitmap * undistorted_image1,Bitmap * undistorted_image2,Camera * undistorted_camera,Eigen::Matrix4d * Q)1005 void RectifyAndUndistortStereoImages(
1006     const UndistortCameraOptions& options, const Bitmap& distorted_image1,
1007     const Bitmap& distorted_image2, const Camera& distorted_camera1,
1008     const Camera& distorted_camera2, const Eigen::Vector4d& qvec,
1009     const Eigen::Vector3d& tvec, Bitmap* undistorted_image1,
1010     Bitmap* undistorted_image2, Camera* undistorted_camera,
1011     Eigen::Matrix4d* Q) {
1012   CHECK_EQ(distorted_camera1.Width(), distorted_image1.Width());
1013   CHECK_EQ(distorted_camera1.Height(), distorted_image1.Height());
1014   CHECK_EQ(distorted_camera2.Width(), distorted_image2.Width());
1015   CHECK_EQ(distorted_camera2.Height(), distorted_image2.Height());
1016 
1017   *undistorted_camera = UndistortCamera(options, distorted_camera1);
1018   undistorted_image1->Allocate(static_cast<int>(undistorted_camera->Width()),
1019                                static_cast<int>(undistorted_camera->Height()),
1020                                distorted_image1.IsRGB());
1021   distorted_image1.CloneMetadata(undistorted_image1);
1022 
1023   undistorted_image2->Allocate(static_cast<int>(undistorted_camera->Width()),
1024                                static_cast<int>(undistorted_camera->Height()),
1025                                distorted_image2.IsRGB());
1026   distorted_image2.CloneMetadata(undistorted_image2);
1027 
1028   Eigen::Matrix3d H1;
1029   Eigen::Matrix3d H2;
1030   RectifyStereoCameras(*undistorted_camera, *undistorted_camera, qvec, tvec,
1031                        &H1, &H2, Q);
1032 
1033   WarpImageWithHomographyBetweenCameras(H1.inverse(), distorted_camera1,
1034                                         *undistorted_camera, distorted_image1,
1035                                         undistorted_image1);
1036   WarpImageWithHomographyBetweenCameras(H2.inverse(), distorted_camera2,
1037                                         *undistorted_camera, distorted_image2,
1038                                         undistorted_image2);
1039 }
1040 
1041 }  // namespace colmap
1042