eos  0.7.1
texture_extraction.hpp
1 /*
2  * Eos - A 3D Morphable Model fitting library written in modern C++11/14.
3  *
4  * File: include/eos/render/texture_extraction.hpp
5  *
6  * Copyright 2014, 2015 Patrik Huber
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 #pragma once
21 
22 #ifndef TEXTURE_EXTRACTION_HPP_
23 #define TEXTURE_EXTRACTION_HPP_
24 
25 #include "eos/render/detail/texture_extraction_detail.hpp"
26 #include "eos/render/Mesh.hpp"
27 #include "eos/render/render_affine.hpp"
28 #include "eos/render/detail/render_detail.hpp"
29 
30 #include "opencv2/core/core.hpp"
31 #include "opencv2/imgproc/imgproc.hpp"
32 
33 #include <tuple>
34 #include <cassert>
35 #include <future>
36 
37 namespace eos {
38  namespace render {
39 
45  NearestNeighbour,
46  Bilinear,
47  Area
48 };
49 
50 // Forward declarations:
51 inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, cv::Mat depthbuffer, bool compute_view_angle, TextureInterpolation mapping_type, int isomap_resolution);
52 namespace detail { cv::Mat interpolate_black_line(cv::Mat isomap); }
53 
76 inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, bool compute_view_angle = false, TextureInterpolation mapping_type = TextureInterpolation::NearestNeighbour, int isomap_resolution = 512)
77 {
78  // Render the model to get a depth buffer:
79  cv::Mat depthbuffer;
80  std::tie(std::ignore, depthbuffer) = render::render_affine(mesh, affine_camera_matrix, image.cols, image.rows);
81  // Note: There's potential for optimisation here - we don't need to do everything that is done in render_affine to just get the depthbuffer.
82 
83  // Now forward the call to the actual texture extraction function:
84  return extract_texture(mesh, affine_camera_matrix, image, depthbuffer, compute_view_angle, mapping_type, isomap_resolution);
85 };
86 
107 inline cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, cv::Mat depthbuffer, bool compute_view_angle = false, TextureInterpolation mapping_type = TextureInterpolation::NearestNeighbour, int isomap_resolution = 512)
108 {
109  assert(mesh.vertices.size() == mesh.texcoords.size());
110  assert(image.type() == CV_8UC3); // the other cases are not yet supported
111 
112  using cv::Mat;
113  using cv::Vec2f;
114  using cv::Vec3f;
115  using cv::Vec4f;
116  using cv::Vec3b;
117  using std::min;
118  using std::max;
119  using std::floor;
120  using std::ceil;
121 
122  affine_camera_matrix = detail::calculate_affine_z_direction(affine_camera_matrix);
123 
124  Mat isomap = Mat::zeros(isomap_resolution, isomap_resolution, CV_8UC4);
125  // #Todo: We should handle gray images, but output a 4-channel isomap nevertheless I think.
126 
127  std::vector<std::future<void>> results;
128  for (const auto& triangle_indices : mesh.tvi) {
129 
130  // Note: If there's a performance problem, there's no need to capture the whole mesh - we could capture only the three required vertices with their texcoords.
131  auto extract_triangle = [&mesh, &affine_camera_matrix, &triangle_indices, &depthbuffer, &isomap, &mapping_type, &image, &compute_view_angle]() {
132 
133  // Find out if the current triangle is visible:
134  // We do a second rendering-pass here. We use the depth-buffer of the final image, and then, here,
135  // check if each pixel in a triangle is visible. If the whole triangle is visible, we use it to extract
136  // the texture.
137  // Possible improvement: - If only part of the triangle is visible, split it
138 
139  // This could be optimized in 2 ways though:
140  // - Use render(), or as in render(...), transfer the vertices once, not in a loop over all triangles (vertices are getting transformed multiple times)
141  // - We transform them later (below) a second time. Only do it once.
142 
143  // Project the triangle vertices to screen coordinates, and use the depthbuffer to check whether the triangle is visible:
144  const Vec4f v0 = Mat(affine_camera_matrix * Mat(mesh.vertices[triangle_indices[0]]));
145  const Vec4f v1 = Mat(affine_camera_matrix * Mat(mesh.vertices[triangle_indices[1]]));
146  const Vec4f v2 = Mat(affine_camera_matrix * Mat(mesh.vertices[triangle_indices[2]]));
147 
148  if (!detail::is_triangle_visible(v0, v1, v2, depthbuffer))
149  {
150  //continue;
151  return;
152  }
153 
154  float alpha_value;
155  if (compute_view_angle)
156  {
157  // Calculate how well visible the current triangle is:
158  // (in essence, the dot product of the viewing direction (0, 0, 1) and the face normal)
159  const Vec3f face_normal = calculate_face_normal(Vec3f(Mat(mesh.vertices[triangle_indices[0]]).rowRange(0, 3)), Vec3f(Mat(mesh.vertices[triangle_indices[1]]).rowRange(0, 3)), Vec3f(Mat(mesh.vertices[triangle_indices[2]]).rowRange(0, 3)));
160  // Transform the normal to "screen" (kind of "eye") space using the upper 3x3 part of the affine camera matrix (=the translation can be ignored):
161  Vec3f face_normal_transformed = Mat(affine_camera_matrix.rowRange(0, 3).colRange(0, 3) * Mat(face_normal));
162  face_normal_transformed /= cv::norm(face_normal_transformed, cv::NORM_L2); // normalise to unit length
163  // Implementation notes regarding the affine camera matrix and the sign:
164  // If the matrix given were the model_view matrix, the sign would be correct.
165  // However, affine_camera_matrix includes glm::ortho, which includes a z-flip.
166  // So we need to flip one of the two signs.
167  // * viewing_direction(0.0f, 0.0f, 1.0f) is correct if affine_camera_matrix were only a model_view matrix
168  // * affine_camera_matrix includes glm::ortho, which flips z, so we flip the sign of viewing_direction.
169  // We don't need the dot product since viewing_direction.xy are 0 and .z is 1:
170  const float angle = -face_normal_transformed[2]; // flip sign, see above
171  assert(angle >= -1.f && angle <= 1.f);
172  // angle is [-1, 1].
173  // * +1 means 0° (same direction)
174  // * 0 means 90°
175  // * -1 means 180° (facing opposite directions)
176  // It's a linear relation, so +0.5 is 45° etc.
177  // An angle larger than 90° means the vertex won't be rendered anyway (because it's back-facing) so we encode 0° to 90°.
178  if (angle < 0.0f) {
179  alpha_value = 0.0f;
180  } else {
181  alpha_value = angle * 255.0f;
182  }
183  }
184  else {
185  // no visibility angle computation - if the triangle/pixel is visible, set the alpha chan to 255 (fully visible pixel).
186  alpha_value = 255.0f;
187  }
188 
189  // Todo: Documentation
190  cv::Point2f src_tri[3];
191  cv::Point2f dst_tri[3];
192 
193  Vec4f vec(mesh.vertices[triangle_indices[0]][0], mesh.vertices[triangle_indices[0]][1], mesh.vertices[triangle_indices[0]][2], 1.0f);
194  Vec4f res = Mat(affine_camera_matrix * Mat(vec));
195  src_tri[0] = Vec2f(res[0], res[1]);
196 
197  vec = Vec4f(mesh.vertices[triangle_indices[1]][0], mesh.vertices[triangle_indices[1]][1], mesh.vertices[triangle_indices[1]][2], 1.0f);
198  res = Mat(affine_camera_matrix * Mat(vec));
199  src_tri[1] = Vec2f(res[0], res[1]);
200 
201  vec = Vec4f(mesh.vertices[triangle_indices[2]][0], mesh.vertices[triangle_indices[2]][1], mesh.vertices[triangle_indices[2]][2], 1.0f);
202  res = Mat(affine_camera_matrix * Mat(vec));
203  src_tri[2] = Vec2f(res[0], res[1]);
204 
205  dst_tri[0] = cv::Point2f(isomap.cols*mesh.texcoords[triangle_indices[0]][0], isomap.rows*mesh.texcoords[triangle_indices[0]][1] - 1.0f);
206  dst_tri[1] = cv::Point2f(isomap.cols*mesh.texcoords[triangle_indices[1]][0], isomap.rows*mesh.texcoords[triangle_indices[1]][1] - 1.0f);
207  dst_tri[2] = cv::Point2f(isomap.cols*mesh.texcoords[triangle_indices[2]][0], isomap.rows*mesh.texcoords[triangle_indices[2]][1] - 1.0f);
208 
209  // Get the inverse Affine Transform from original image: from dst to src
210  Mat warp_mat_org_inv = cv::getAffineTransform(dst_tri, src_tri);
211  warp_mat_org_inv.convertTo(warp_mat_org_inv, CV_32FC1);
212 
213  // We now loop over all pixels in the triangle and select, depending on the mapping type, the corresponding texel(s) in the source image
214  for (int x = min(dst_tri[0].x, min(dst_tri[1].x, dst_tri[2].x)); x < max(dst_tri[0].x, max(dst_tri[1].x, dst_tri[2].x)); ++x) {
215  for (int y = min(dst_tri[0].y, min(dst_tri[1].y, dst_tri[2].y)); y < max(dst_tri[0].y, max(dst_tri[1].y, dst_tri[2].y)); ++y) {
216  if (detail::is_point_in_triangle(cv::Point2f(x, y), dst_tri[0], dst_tri[1], dst_tri[2])) {
217  if (mapping_type == TextureInterpolation::Area) {
218 
219  // calculate positions of 4 corners of pixel in image (src)
220  Vec3f homogenous_dst_upper_left(x - 0.5f, y - 0.5f, 1.0f);
221  Vec3f homogenous_dst_upper_right(x + 0.5f, y - 0.5f, 1.0f);
222  Vec3f homogenous_dst_lower_left(x - 0.5f, y + 0.5f, 1.0f);
223  Vec3f homogenous_dst_lower_right(x + 0.5f, y + 0.5f, 1.0f);
224 
225  Vec2f src_texel_upper_left = Mat(warp_mat_org_inv * Mat(homogenous_dst_upper_left));
226  Vec2f src_texel_upper_right = Mat(warp_mat_org_inv * Mat(homogenous_dst_upper_right));
227  Vec2f src_texel_lower_left = Mat(warp_mat_org_inv * Mat(homogenous_dst_lower_left));
228  Vec2f src_texel_lower_right = Mat(warp_mat_org_inv * Mat(homogenous_dst_lower_right));
229 
230  float min_a = min(min(src_texel_upper_left[0], src_texel_upper_right[0]), min(src_texel_lower_left[0], src_texel_lower_right[0]));
231  float max_a = max(max(src_texel_upper_left[0], src_texel_upper_right[0]), max(src_texel_lower_left[0], src_texel_lower_right[0]));
232  float min_b = min(min(src_texel_upper_left[1], src_texel_upper_right[1]), min(src_texel_lower_left[1], src_texel_lower_right[1]));
233  float max_b = max(max(src_texel_upper_left[1], src_texel_upper_right[1]), max(src_texel_lower_left[1], src_texel_lower_right[1]));
234 
235  cv::Vec3i color;
236  int num_texels = 0;
237 
238  for (int a = ceil(min_a); a <= floor(max_a); ++a)
239  {
240  for (int b = ceil(min_b); b <= floor(max_b); ++b)
241  {
242  if (detail::is_point_in_triangle(cv::Point2f(a, b), src_texel_upper_left, src_texel_lower_left, src_texel_upper_right) || detail::is_point_in_triangle(cv::Point2f(a, b), src_texel_lower_left, src_texel_upper_right, src_texel_lower_right)) {
243  if (a < image.cols && b < image.rows) { // if src_texel in triangle and in image
244  num_texels++;
245  color += image.at<Vec3b>(b, a);
246  }
247  }
248  }
249  }
250  if (num_texels > 0)
251  color = color / num_texels;
252  else { // if no corresponding texel found, nearest neighbour interpolation
253  // calculate corresponding position of dst_coord pixel center in image (src)
254  Vec3f homogenous_dst_coord = Vec3f(x, y, 1.0f);
255  Vec2f src_texel = Mat(warp_mat_org_inv * Mat(homogenous_dst_coord));
256 
257  if ((cvRound(src_texel[1]) < image.rows) && cvRound(src_texel[0]) < image.cols) {
258  color = image.at<Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]));
259  }
260  }
261  isomap.at<Vec3b>(y, x) = color;
262  }
263  else if (mapping_type == TextureInterpolation::Bilinear) {
264 
265  // calculate corresponding position of dst_coord pixel center in image (src)
266  Vec3f homogenous_dst_coord(x, y, 1.0f);
267  Vec2f src_texel = Mat(warp_mat_org_inv * Mat(homogenous_dst_coord));
268 
269  // calculate distances to next 4 pixels
270  using std::sqrt;
271  using std::pow;
272  float distance_upper_left = sqrt(pow(src_texel[0] - floor(src_texel[0]), 2) + pow(src_texel[1] - floor(src_texel[1]), 2));
273  float distance_upper_right = sqrt(pow(src_texel[0] - floor(src_texel[0]), 2) + pow(src_texel[1] - ceil(src_texel[1]), 2));
274  float distance_lower_left = sqrt(pow(src_texel[0] - ceil(src_texel[0]), 2) + pow(src_texel[1] - floor(src_texel[1]), 2));
275  float distance_lower_right = sqrt(pow(src_texel[0] - ceil(src_texel[0]), 2) + pow(src_texel[1] - ceil(src_texel[1]), 2));
276 
277  // normalise distances
278  float sum_distances = distance_lower_left + distance_lower_right + distance_upper_left + distance_upper_right;
279  distance_lower_left /= sum_distances;
280  distance_lower_right /= sum_distances;
281  distance_upper_left /= sum_distances;
282  distance_upper_right /= sum_distances;
283 
284  // set color depending on distance from next 4 pixels
285  for (int color = 0; color < 3; ++color) {
286  float color_upper_left = image.at<Vec3b>(floor(src_texel[1]), floor(src_texel[0]))[color] * distance_upper_left;
287  float color_upper_right = image.at<Vec3b>(floor(src_texel[1]), ceil(src_texel[0]))[color] * distance_upper_right;
288  float color_lower_left = image.at<Vec3b>(ceil(src_texel[1]), floor(src_texel[0]))[color] * distance_lower_left;
289  float color_lower_right = image.at<Vec3b>(ceil(src_texel[1]), ceil(src_texel[0]))[color] * distance_lower_right;
290 
291  isomap.at<Vec3b>(y, x)[color] = color_upper_left + color_upper_right + color_lower_left + color_lower_right;
292  }
293  }
294  else if (mapping_type == TextureInterpolation::NearestNeighbour) {
295 
296  // calculate corresponding position of dst_coord pixel center in image (src)
297  const Mat homogenous_dst_coord(Vec3f(x, y, 1.0f));
298  const Vec2f src_texel = Mat(warp_mat_org_inv * homogenous_dst_coord);
299 
300  if ((cvRound(src_texel[1]) < image.rows) && (cvRound(src_texel[0]) < image.cols) && cvRound(src_texel[0]) > 0 && cvRound(src_texel[1]) > 0)
301  {
302  cv::Vec4b isomap_pixel;
303  isomap.at<cv::Vec4b>(y, x)[0] = image.at<Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]))[0];
304  isomap.at<cv::Vec4b>(y, x)[1] = image.at<Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]))[1];
305  isomap.at<cv::Vec4b>(y, x)[2] = image.at<Vec3b>(cvRound(src_texel[1]), cvRound(src_texel[0]))[2];
306  isomap.at<cv::Vec4b>(y, x)[3] = static_cast<uchar>(alpha_value); // pixel is visible
307  }
308  }
309  }
310  }
311  }
312  }; // end lambda auto extract_triangle();
313  results.emplace_back(std::async(extract_triangle));
314  } // end for all mesh.tvi
315  // Collect all the launched tasks:
316  for (auto&& r : results) {
317  r.get();
318  }
319 
320  // Workaround for the black line in the isomap (see GitHub issue #4):
321  if (mesh.texcoords.size() <= 3448)
322  {
323  isomap = detail::interpolate_black_line(isomap);
324  }
325 
326  return isomap;
327 };
328 
329 namespace detail {
330 
331 // Workaround for the pixels that don't get filled in extract_texture().
332 // There's a vertical line of missing values in the middle of the isomap,
333 // as well as a few pixels on a horizontal line around the mouth. They
334 // manifest themselves as black lines in the final isomap. This function
335 // just fills these missing values by interpolating between two neighbouring
336 // pixels. See GitHub issue #4.
337 cv::Mat interpolate_black_line(cv::Mat isomap)
338 {
339  assert(isomap.type() == CV_8UC4);
340  // Replace the vertical black line ("missing data"):
341  int col = isomap.cols / 2;
342  for (int row = 0; row < isomap.rows; ++row)
343  {
344  if (isomap.at<cv::Vec4b>(row, col) == cv::Vec4b(0, 0, 0, 0))
345  {
346  isomap.at<cv::Vec4b>(row, col) = isomap.at<cv::Vec4b>(row, col - 1) * 0.5f + isomap.at<cv::Vec4b>(row, col + 1) * 0.5f;
347  }
348  }
349 
350  // Replace the horizontal line around the mouth that occurs in the
351  // isomaps of resolution 512x512 and higher:
352  if (isomap.rows == 512) // num cols is 512 as well
353  {
354  int r = 362;
355  for (int c = 206; c <= 306; ++c)
356  {
357  if (isomap.at<cv::Vec4b>(r, c) == cv::Vec4b(0, 0, 0, 0))
358  {
359  isomap.at<cv::Vec4b>(r, c) = isomap.at<cv::Vec4b>(r - 1, c) * 0.5f + isomap.at<cv::Vec4b>(r + 1, c) * 0.5f;
360  }
361  }
362  }
363  if (isomap.rows == 1024) // num cols is 1024 as well
364  {
365  int r = 724;
366  for (int c = 437; c <= 587; ++c)
367  {
368  if (isomap.at<cv::Vec4b>(r, c) == cv::Vec4b(0, 0, 0, 0))
369  {
370  isomap.at<cv::Vec4b>(r, c) = isomap.at<cv::Vec4b>(r - 1, c) * 0.5f + isomap.at<cv::Vec4b>(r + 1, c) * 0.5f;
371  }
372  }
373  r = 725;
374  for (int c = 411; c <= 613; ++c)
375  {
376  if (isomap.at<cv::Vec4b>(r, c) == cv::Vec4b(0, 0, 0, 0))
377  {
378  isomap.at<cv::Vec4b>(r, c) = isomap.at<cv::Vec4b>(r - 1, c) * 0.5f + isomap.at<cv::Vec4b>(r + 1, c) * 0.5f;
379  }
380  }
381  }
382  // Higher resolutions are probably affected as well but not used so far in practice.
383 
384  return isomap;
385 };
386 } /* namespace detail */
387 
388  } /* namespace render */
389 } /* namespace eos */
390 
391 #endif /* TEXTURE_EXTRACTION_HPP_ */
cv::Vec3f calculate_face_normal(const cv::Vec3f &v0, const cv::Vec3f &v1, const cv::Vec3f &v2)
Definition: utils.hpp:95
std::pair< cv::Mat, cv::Mat > render(Mesh mesh, cv::Mat model_view_matrix, cv::Mat projection_matrix, int viewport_width, int viewport_height, const boost::optional< Texture > &texture=boost::none, bool enable_backface_culling=false, bool enable_near_clipping=true, bool enable_far_clipping=true)
Definition: render.hpp:125
TextureInterpolation
Definition: texture_extraction.hpp:44
cv::Mat extract_texture(Mesh mesh, cv::Mat affine_camera_matrix, cv::Mat image, cv::Mat depthbuffer, bool compute_view_angle, TextureInterpolation mapping_type, int isomap_resolution)
Definition: texture_extraction.hpp:107
std::pair< cv::Mat, cv::Mat > render_affine(Mesh mesh, cv::Mat affine_camera_matrix, int viewport_width, int viewport_height, bool do_backface_culling=true)
Definition: render_affine.hpp:52
Namespace containing all of eos&#39;s 3D model fitting functionality.
std::vector< cv::Vec4f > vertices
3D vertex positions.
Definition: Mesh.hpp:47
This class represents a 3D mesh consisting of vertices, vertex colour information and texture coordin...
Definition: Mesh.hpp:45
std::vector< std::array< int, 3 > > tvi
Triangle vertex indices.
Definition: Mesh.hpp:51
std::vector< cv::Vec2f > texcoords
Texture coordinates for each vertex.
Definition: Mesh.hpp:49