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 "mvs/patch_match.h"
33
34 #include <numeric>
35 #include <unordered_set>
36
37 #include "mvs/consistency_graph.h"
38 #include "mvs/patch_match_cuda.h"
39 #include "mvs/workspace.h"
40 #include "util/math.h"
41 #include "util/misc.h"
42
43 #define PrintOption(option) std::cout << #option ": " << option << std::endl
44
45 namespace colmap {
46 namespace mvs {
47
PatchMatch(const PatchMatchOptions & options,const Problem & problem)48 PatchMatch::PatchMatch(const PatchMatchOptions& options, const Problem& problem)
49 : options_(options), problem_(problem) {}
50
~PatchMatch()51 PatchMatch::~PatchMatch() {}
52
Print() const53 void PatchMatchOptions::Print() const {
54 PrintHeading2("PatchMatchOptions");
55 PrintOption(max_image_size);
56 PrintOption(gpu_index);
57 PrintOption(depth_min);
58 PrintOption(depth_max);
59 PrintOption(window_radius);
60 PrintOption(window_step);
61 PrintOption(sigma_spatial);
62 PrintOption(sigma_color);
63 PrintOption(num_samples);
64 PrintOption(ncc_sigma);
65 PrintOption(min_triangulation_angle);
66 PrintOption(incident_angle_sigma);
67 PrintOption(num_iterations);
68 PrintOption(geom_consistency);
69 PrintOption(geom_consistency_regularizer);
70 PrintOption(geom_consistency_max_cost);
71 PrintOption(filter);
72 PrintOption(filter_min_ncc);
73 PrintOption(filter_min_triangulation_angle);
74 PrintOption(filter_min_num_consistent);
75 PrintOption(filter_geom_consistency_max_cost);
76 PrintOption(write_consistency_graph);
77 }
78
Print() const79 void PatchMatch::Problem::Print() const {
80 PrintHeading2("PatchMatch::Problem");
81
82 PrintOption(ref_image_idx);
83
84 std::cout << "src_image_idxs: ";
85 if (!src_image_idxs.empty()) {
86 for (size_t i = 0; i < src_image_idxs.size() - 1; ++i) {
87 std::cout << src_image_idxs[i] << " ";
88 }
89 std::cout << src_image_idxs.back() << std::endl;
90 } else {
91 std::cout << std::endl;
92 }
93 }
94
Check() const95 void PatchMatch::Check() const {
96 CHECK(options_.Check());
97
98 CHECK(!options_.gpu_index.empty());
99 const std::vector<int> gpu_indices = CSVToVector<int>(options_.gpu_index);
100 CHECK_EQ(gpu_indices.size(), 1);
101 CHECK_GE(gpu_indices[0], -1);
102
103 CHECK_NOTNULL(problem_.images);
104 if (options_.geom_consistency) {
105 CHECK_NOTNULL(problem_.depth_maps);
106 CHECK_NOTNULL(problem_.normal_maps);
107 CHECK_EQ(problem_.depth_maps->size(), problem_.images->size());
108 CHECK_EQ(problem_.normal_maps->size(), problem_.images->size());
109 }
110
111 CHECK_GT(problem_.src_image_idxs.size(), 0);
112
113 // Check that there are no duplicate images and that the reference image
114 // is not defined as a source image.
115 std::set<int> unique_image_idxs(problem_.src_image_idxs.begin(),
116 problem_.src_image_idxs.end());
117 unique_image_idxs.insert(problem_.ref_image_idx);
118 CHECK_EQ(problem_.src_image_idxs.size() + 1, unique_image_idxs.size());
119
120 // Check that input data is well-formed.
121 for (const int image_idx : unique_image_idxs) {
122 CHECK_GE(image_idx, 0) << image_idx;
123 CHECK_LT(image_idx, problem_.images->size()) << image_idx;
124
125 const Image& image = problem_.images->at(image_idx);
126 CHECK_GT(image.GetBitmap().Width(), 0) << image_idx;
127 CHECK_GT(image.GetBitmap().Height(), 0) << image_idx;
128 CHECK(image.GetBitmap().IsGrey()) << image_idx;
129 CHECK_EQ(image.GetWidth(), image.GetBitmap().Width()) << image_idx;
130 CHECK_EQ(image.GetHeight(), image.GetBitmap().Height()) << image_idx;
131
132 // Make sure, the calibration matrix only contains fx, fy, cx, cy.
133 CHECK_LT(std::abs(image.GetK()[1] - 0.0f), 1e-6f) << image_idx;
134 CHECK_LT(std::abs(image.GetK()[3] - 0.0f), 1e-6f) << image_idx;
135 CHECK_LT(std::abs(image.GetK()[6] - 0.0f), 1e-6f) << image_idx;
136 CHECK_LT(std::abs(image.GetK()[7] - 0.0f), 1e-6f) << image_idx;
137 CHECK_LT(std::abs(image.GetK()[8] - 1.0f), 1e-6f) << image_idx;
138
139 if (options_.geom_consistency) {
140 CHECK_LT(image_idx, problem_.depth_maps->size()) << image_idx;
141 const DepthMap& depth_map = problem_.depth_maps->at(image_idx);
142 CHECK_EQ(image.GetWidth(), depth_map.GetWidth()) << image_idx;
143 CHECK_EQ(image.GetHeight(), depth_map.GetHeight()) << image_idx;
144 }
145 }
146
147 if (options_.geom_consistency) {
148 const Image& ref_image = problem_.images->at(problem_.ref_image_idx);
149 const NormalMap& ref_normal_map =
150 problem_.normal_maps->at(problem_.ref_image_idx);
151 CHECK_EQ(ref_image.GetWidth(), ref_normal_map.GetWidth());
152 CHECK_EQ(ref_image.GetHeight(), ref_normal_map.GetHeight());
153 }
154 }
155
Run()156 void PatchMatch::Run() {
157 PrintHeading2("PatchMatch::Run");
158
159 Check();
160
161 patch_match_cuda_.reset(new PatchMatchCuda(options_, problem_));
162 patch_match_cuda_->Run();
163 }
164
GetDepthMap() const165 DepthMap PatchMatch::GetDepthMap() const {
166 return patch_match_cuda_->GetDepthMap();
167 }
168
GetNormalMap() const169 NormalMap PatchMatch::GetNormalMap() const {
170 return patch_match_cuda_->GetNormalMap();
171 }
172
GetSelProbMap() const173 Mat<float> PatchMatch::GetSelProbMap() const {
174 return patch_match_cuda_->GetSelProbMap();
175 }
176
GetConsistencyGraph() const177 ConsistencyGraph PatchMatch::GetConsistencyGraph() const {
178 const auto& ref_image = problem_.images->at(problem_.ref_image_idx);
179 return ConsistencyGraph(ref_image.GetWidth(), ref_image.GetHeight(),
180 patch_match_cuda_->GetConsistentImageIdxs());
181 }
182
PatchMatchController(const PatchMatchOptions & options,const std::string & workspace_path,const std::string & workspace_format,const std::string & pmvs_option_name)183 PatchMatchController::PatchMatchController(const PatchMatchOptions& options,
184 const std::string& workspace_path,
185 const std::string& workspace_format,
186 const std::string& pmvs_option_name)
187 : options_(options),
188 workspace_path_(workspace_path),
189 workspace_format_(workspace_format),
190 pmvs_option_name_(pmvs_option_name) {
191 std::vector<int> gpu_indices = CSVToVector<int>(options_.gpu_index);
192 }
193
Run()194 void PatchMatchController::Run() {
195 ReadWorkspace();
196 ReadProblems();
197 ReadGpuIndices();
198
199 thread_pool_.reset(new ThreadPool(gpu_indices_.size()));
200
201 // If geometric consistency is enabled, then photometric output must be
202 // computed first for all images without filtering.
203 if (options_.geom_consistency) {
204 auto photometric_options = options_;
205 photometric_options.geom_consistency = false;
206 photometric_options.filter = false;
207
208 for (size_t problem_idx = 0; problem_idx < problems_.size();
209 ++problem_idx) {
210 thread_pool_->AddTask(&PatchMatchController::ProcessProblem, this,
211 photometric_options, problem_idx);
212 }
213
214 thread_pool_->Wait();
215 }
216
217 for (size_t problem_idx = 0; problem_idx < problems_.size(); ++problem_idx) {
218 thread_pool_->AddTask(&PatchMatchController::ProcessProblem, this, options_,
219 problem_idx);
220 }
221
222 thread_pool_->Wait();
223
224 GetTimer().PrintMinutes();
225 }
226
ReadWorkspace()227 void PatchMatchController::ReadWorkspace() {
228 std::cout << "Reading workspace..." << std::endl;
229
230 Workspace::Options workspace_options;
231
232 auto workspace_format_lower_case = workspace_format_;
233 StringToLower(&workspace_format_lower_case);
234 if (workspace_format_lower_case == "pmvs") {
235 workspace_options.stereo_folder =
236 StringPrintf("stereo-%s", pmvs_option_name_.c_str());
237 }
238
239 workspace_options.max_image_size = options_.max_image_size;
240 workspace_options.image_as_rgb = false;
241 workspace_options.cache_size = options_.cache_size;
242 workspace_options.workspace_path = workspace_path_;
243 workspace_options.workspace_format = workspace_format_;
244 workspace_options.input_type = options_.geom_consistency ? "photometric" : "";
245
246 workspace_.reset(new Workspace(workspace_options));
247
248 if (workspace_format_lower_case == "pmvs") {
249 std::cout << StringPrintf("Importing PMVS workspace (option %s)...",
250 pmvs_option_name_.c_str())
251 << std::endl;
252 ImportPMVSWorkspace(*workspace_, pmvs_option_name_);
253 }
254
255 depth_ranges_ = workspace_->GetModel().ComputeDepthRanges();
256 }
257
ReadProblems()258 void PatchMatchController::ReadProblems() {
259 std::cout << "Reading configuration..." << std::endl;
260
261 problems_.clear();
262
263 const auto& model = workspace_->GetModel();
264
265 std::vector<std::string> config = ReadTextFileLines(
266 JoinPaths(workspace_path_, workspace_->GetOptions().stereo_folder,
267 "patch-match.cfg"));
268
269 std::vector<std::map<int, int>> shared_num_points;
270 std::vector<std::map<int, float>> triangulation_angles;
271
272 const float min_triangulation_angle_rad =
273 DegToRad(options_.min_triangulation_angle);
274
275 std::string ref_image_name;
276 std::unordered_set<int> ref_image_idxs;
277
278 struct ProblemConfig {
279 std::string ref_image_name;
280 std::vector<std::string> src_image_names;
281 };
282 std::vector<ProblemConfig> problem_configs;
283
284 for (size_t i = 0; i < config.size(); ++i) {
285 std::string& config_line = config[i];
286 StringTrim(&config_line);
287
288 if (config_line.empty() || config_line[0] == '#') {
289 continue;
290 }
291
292 if (ref_image_name.empty()) {
293 ref_image_name = config_line;
294 continue;
295 }
296
297 ref_image_idxs.insert(model.GetImageIdx(ref_image_name));
298
299 ProblemConfig problem_config;
300 problem_config.ref_image_name = ref_image_name;
301 problem_config.src_image_names = CSVToVector<std::string>(config_line);
302 problem_configs.push_back(problem_config);
303
304 ref_image_name.clear();
305 }
306
307 for (const auto& problem_config : problem_configs) {
308 PatchMatch::Problem problem;
309
310 problem.ref_image_idx = model.GetImageIdx(problem_config.ref_image_name);
311
312 if (problem_config.src_image_names.size() == 1 &&
313 problem_config.src_image_names[0] == "__all__") {
314 // Use all images as source images.
315 problem.src_image_idxs.clear();
316 problem.src_image_idxs.reserve(model.images.size() - 1);
317 for (const int image_idx : ref_image_idxs) {
318 if (image_idx != problem.ref_image_idx) {
319 problem.src_image_idxs.push_back(image_idx);
320 }
321 }
322 } else if (problem_config.src_image_names.size() == 2 &&
323 problem_config.src_image_names[0] == "__auto__") {
324 // Use maximum number of overlapping images as source images. Overlapping
325 // will be sorted based on the number of shared points to the reference
326 // image and the top ranked images are selected. Note that images are only
327 // selected if some points have a sufficient triangulation angle.
328
329 if (shared_num_points.empty()) {
330 shared_num_points = model.ComputeSharedPoints();
331 }
332 if (triangulation_angles.empty()) {
333 const float kTriangulationAnglePercentile = 75;
334 triangulation_angles =
335 model.ComputeTriangulationAngles(kTriangulationAnglePercentile);
336 }
337
338 const size_t max_num_src_images =
339 std::stoll(problem_config.src_image_names[1]);
340
341 const auto& overlapping_images =
342 shared_num_points.at(problem.ref_image_idx);
343 const auto& overlapping_triangulation_angles =
344 triangulation_angles.at(problem.ref_image_idx);
345
346 std::vector<std::pair<int, int>> src_images;
347 src_images.reserve(overlapping_images.size());
348 for (const auto& image : overlapping_images) {
349 if (ref_image_idxs.count(image.first) &&
350 overlapping_triangulation_angles.at(image.first) >=
351 min_triangulation_angle_rad) {
352 src_images.emplace_back(image.first, image.second);
353 }
354 }
355
356 const size_t eff_max_num_src_images =
357 std::min(src_images.size(), max_num_src_images);
358
359 std::partial_sort(src_images.begin(),
360 src_images.begin() + eff_max_num_src_images,
361 src_images.end(),
362 [](const std::pair<int, int>& image1,
363 const std::pair<int, int>& image2) {
364 return image1.second > image2.second;
365 });
366
367 problem.src_image_idxs.reserve(eff_max_num_src_images);
368 for (size_t i = 0; i < eff_max_num_src_images; ++i) {
369 problem.src_image_idxs.push_back(src_images[i].first);
370 }
371 } else {
372 problem.src_image_idxs.reserve(problem_config.src_image_names.size());
373 for (const auto& src_image_name : problem_config.src_image_names) {
374 problem.src_image_idxs.push_back(model.GetImageIdx(src_image_name));
375 }
376 }
377
378 if (problem.src_image_idxs.empty()) {
379 std::cout
380 << StringPrintf(
381 "WARNING: Ignoring reference image %s, because it has no "
382 "source images.",
383 problem_config.ref_image_name.c_str())
384 << std::endl;
385 } else {
386 problems_.push_back(problem);
387 }
388 }
389
390 std::cout << StringPrintf("Configuration has %d problems...",
391 problems_.size())
392 << std::endl;
393 }
394
ReadGpuIndices()395 void PatchMatchController::ReadGpuIndices() {
396 gpu_indices_ = CSVToVector<int>(options_.gpu_index);
397 if (gpu_indices_.size() == 1 && gpu_indices_[0] == -1) {
398 const int num_cuda_devices = GetNumCudaDevices();
399 CHECK_GT(num_cuda_devices, 0);
400 gpu_indices_.resize(num_cuda_devices);
401 std::iota(gpu_indices_.begin(), gpu_indices_.end(), 0);
402 }
403 }
404
ProcessProblem(const PatchMatchOptions & options,const size_t problem_idx)405 void PatchMatchController::ProcessProblem(const PatchMatchOptions& options,
406 const size_t problem_idx) {
407 if (IsStopped()) {
408 return;
409 }
410
411 const auto& model = workspace_->GetModel();
412
413 auto& problem = problems_.at(problem_idx);
414 const int gpu_index = gpu_indices_.at(thread_pool_->GetThreadIndex());
415 CHECK_GE(gpu_index, -1);
416
417 const std::string& stereo_folder = workspace_->GetOptions().stereo_folder;
418 const std::string output_type =
419 options.geom_consistency ? "geometric" : "photometric";
420 const std::string image_name = model.GetImageName(problem.ref_image_idx);
421 const std::string file_name =
422 StringPrintf("%s.%s.bin", image_name.c_str(), output_type.c_str());
423 const std::string depth_map_path =
424 JoinPaths(workspace_path_, stereo_folder, "depth_maps", file_name);
425 const std::string normal_map_path =
426 JoinPaths(workspace_path_, stereo_folder, "normal_maps", file_name);
427 const std::string consistency_graph_path = JoinPaths(
428 workspace_path_, stereo_folder, "consistency_graphs", file_name);
429
430 if (ExistsFile(depth_map_path) && ExistsFile(normal_map_path) &&
431 (!options.write_consistency_graph ||
432 ExistsFile(consistency_graph_path))) {
433 return;
434 }
435
436 PrintHeading1(StringPrintf("Processing view %d / %d", problem_idx + 1,
437 problems_.size()));
438
439 auto patch_match_options = options;
440
441 if (patch_match_options.depth_min < 0 || patch_match_options.depth_max < 0) {
442 patch_match_options.depth_min =
443 depth_ranges_.at(problem.ref_image_idx).first;
444 patch_match_options.depth_max =
445 depth_ranges_.at(problem.ref_image_idx).second;
446 CHECK(patch_match_options.depth_min > 0 &&
447 patch_match_options.depth_max > 0)
448 << " - You must manually set the minimum and maximum depth, since no "
449 "sparse model is provided in the workspace.";
450 }
451
452 patch_match_options.gpu_index = std::to_string(gpu_index);
453
454 if (patch_match_options.sigma_spatial <= 0.0f) {
455 patch_match_options.sigma_spatial = patch_match_options.window_radius;
456 }
457
458 std::vector<Image> images = model.images;
459 std::vector<DepthMap> depth_maps;
460 std::vector<NormalMap> normal_maps;
461 if (options.geom_consistency) {
462 depth_maps.resize(model.images.size());
463 normal_maps.resize(model.images.size());
464 }
465
466 problem.images = &images;
467 problem.depth_maps = &depth_maps;
468 problem.normal_maps = &normal_maps;
469
470 {
471 // Collect all used images in current problem.
472 std::unordered_set<int> used_image_idxs(problem.src_image_idxs.begin(),
473 problem.src_image_idxs.end());
474 used_image_idxs.insert(problem.ref_image_idx);
475
476 patch_match_options.filter_min_num_consistent =
477 std::min(static_cast<int>(used_image_idxs.size()) - 1,
478 patch_match_options.filter_min_num_consistent);
479
480 // Only access workspace from one thread at a time and only spawn resample
481 // threads from one master thread at a time.
482 std::unique_lock<std::mutex> lock(workspace_mutex_);
483
484 std::cout << "Reading inputs..." << std::endl;
485 for (const auto image_idx : used_image_idxs) {
486 images.at(image_idx).SetBitmap(workspace_->GetBitmap(image_idx));
487 if (options.geom_consistency) {
488 depth_maps.at(image_idx) = workspace_->GetDepthMap(image_idx);
489 normal_maps.at(image_idx) = workspace_->GetNormalMap(image_idx);
490 }
491 }
492 }
493
494 problem.Print();
495 patch_match_options.Print();
496
497 PatchMatch patch_match(patch_match_options, problem);
498 patch_match.Run();
499
500 std::cout << std::endl
501 << StringPrintf("Writing %s output for %s", output_type.c_str(),
502 image_name.c_str())
503 << std::endl;
504
505 patch_match.GetDepthMap().Write(depth_map_path);
506 patch_match.GetNormalMap().Write(normal_map_path);
507 if (options.write_consistency_graph) {
508 patch_match.GetConsistencyGraph().Write(consistency_graph_path);
509 }
510 }
511
512 } // namespace mvs
513 } // namespace colmap
514