From 5e6fa616623cf5e704082e6c45d659e566f3afd0 Mon Sep 17 00:00:00 2001 From: Alex Yu Date: Sat, 24 Aug 2019 23:47:07 -0700 Subject: [PATCH] Initial commit, working synthetic dataset generator (dependencies are ugly) --- CMakeLists.txt | 130 +++++++++++++ LICENSE.md | 13 ++ README.md | 50 +++++ intrin.txt | 12 ++ partlabels.cpp | 497 +++++++++++++++++++++++++++++++++++++++++++++++++ sample.cpp | 482 +++++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 1184 insertions(+) create mode 100644 CMakeLists.txt create mode 100644 LICENSE.md create mode 100644 README.md create mode 100644 intrin.txt create mode 100644 partlabels.cpp create mode 100644 sample.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..f9a2dec --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,130 @@ +cmake_minimum_required( VERSION 2.8 ) + +if(NOT CMAKE_BUILD_TYPE) + set(CMAKE_BUILD_TYPE Release) +endif() + +project( smplsynth ) + +set( CMAKE_CXX_STACK_SIZE "10000000" ) +set( CMAKE_CXX_STANDARD 11 ) +set( CMAKE_MODULE_PATH "${PROJECT_SOURCE_DIR}/cmake_modules" ) +if (WIN32) + set( COLOR_TERM OFF ) +else() + set( COLOR_TERM ON ) +endif() + +if( NOT CMAKE_BUILD_TYPE ) + set( CMAKE_BUILD_TYPE Release ) +endif() + +set( PROJ_NAME "smplsynth" ) +set( PART_LABELS_NAME "partlabels" ) +set( OUTPUT_NAME "smplsynth" ) + +include( CheckCXXCompilerFlag ) +CHECK_CXX_COMPILER_FLAG( "-std=c++11" COMPILER_SUPPORTS_CXX11 ) +CHECK_CXX_COMPILER_FLAG( "-std=c++0x" COMPILER_SUPPORTS_CXX0X ) + +if( COMPILER_SUPPORTS_CXX11 ) + set( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11" ) +elseif( COMPILER_SUPPORTS_CXX0X ) + set( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++0x" ) +else() + message( STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support. Please use a different C++ compiler." ) +endif() + +if ( CMAKE_COMPILER_IS_GNUCXX ) + set( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-deprecated -Wno-deprecated-declarations -O3 -g" ) +endif ( CMAKE_COMPILER_IS_GNUCXX ) + +set( Boost_USE_STATIC_LIBS ON ) +set( Boost_USE_STATIC ON ) + +if( NOT DEFINED EIGEN_INCLUDE_DIRS ) + if ( DEFINED Eigen_DIR ) + set( EIGEN_INCLUDE_DIRS ${Eigen_DIR} ) + endif ( DEFINED Eigen_DIR ) +endif( NOT DEFINED EIGEN_INCLUDE_DIRS ) +if ( NOT DEFINED Eigen_DIR ) + find_package( Eigen QUIET ) +endif ( NOT DEFINED Eigen_DIR ) +set( EIGEN_INCLUDE_DIR EIGEN_INCLUDE_DIRS ) +include_directories(${EIGEN_INCLUDE_DIRS}) +message(STATUS "Found Eigen: ${EIGEN_INCLUDE_DIRS}") + +# require OpenCV +find_package( OpenCV REQUIRED ) +if( OpenCV_FOUND ) + message( STATUS "Found OpenCV: ${OpenCV_INCLUDE_DIRS}" ) +endif( OpenCV_FOUND ) + +# require PCL +find_package( PCL REQUIRED ) +if( PCL_FOUND ) + message( STATUS "Found PCL: ${PCL_INCLUDE_DIRS}" ) +endif( PCL_FOUND ) + +# require OpenARK +find_package( OpenARK REQUIRED ) +if( OpenARK_FOUND ) + message( STATUS "Found OpenARK: ${OpenARK_INCLUDE_DIR}" ) +endif( OpenARK_FOUND ) + +# require K4A +find_package( k4a REQUIRED ) +if( k4a_FOUND) + message( STATUS "Found k4a (Azure Kinect SDK)" ) +endif( k4a_FOUND ) + +find_package( Boost REQUIRED COMPONENTS system program_options ) + +# require Ceres +find_package( Ceres REQUIRED ) +IF(Ceres_FOUND) + MESSAGE(STATUS "Found Ceres: ${CERES_INCLUDE_DIRS} ${CERES_LIBRARIES}") +ENDIF(Ceres_FOUND) + +include_directories( + ${OpenCV_INCLUDE_DIRS} + ${PCL_INCLUDE_DIRS} + ${CERES_INCLUDE_DIRS} + ${OpenARK_INCLUDE_DIR} + ${Boost_INCLUDE_DIRS} +) + +set( + DEPENDENCIES + ${OpenCV_LIBRARIES} + ${PCL_LIBRARIES} + ${OpenARK_LIB} + ${CERES_LIBRARIES} + ${Boost_LIBRARIES} + k4a::k4a +) + +add_definitions( + ${PCL_DEFINITIONS} +) +foreach( DEFINITION ${PCL_DEFINITIONS} ) + set( TARGET_COMPILE_FLAGS "${TARGET_COMPILE_FLAGS} ${DEFINITION}" ) +endforeach() + +add_executable( ${PROJ_NAME} sample.cpp ) +target_include_directories( ${PROJ_NAME} PRIVATE ${INCLUDE_DIR} ) +set_target_properties( ${PROJ_NAME} PROPERTIES OUTPUT_NAME ${OUTPUT_NAME} ) +target_link_libraries( ${PROJ_NAME} ${DEPENDENCIES} ) +set_target_properties( ${PROJ_NAME} PROPERTIES COMPILE_FLAGS ${TARGET_COMPILE_FLAGS} ) + +add_executable( ${PART_LABELS_NAME} partlabels.cpp ) +target_include_directories( ${PART_LABELS_NAME} PRIVATE ${INCLUDE_DIR} ) +set_target_properties( ${PART_LABELS_NAME} PROPERTIES OUTPUT_NAME ${PART_LABELS_NAME} ) +target_link_libraries( ${PART_LABELS_NAME} ${DEPENDENCIES} ) +set_target_properties( ${PART_LABELS_NAME} PROPERTIES COMPILE_FLAGS ${TARGET_COMPILE_FLAGS} ) + +if ( MSVC ) + set_property(TARGET ${PROJ_NAME} APPEND PROPERTY LINK_FLAGS /DEBUG) +else () + target_link_libraries( ${PROJ_NAME} -pthread ) +endif ( MSVC ) diff --git a/LICENSE.md b/LICENSE.md new file mode 100644 index 0000000..d5795ec --- /dev/null +++ b/LICENSE.md @@ -0,0 +1,13 @@ +Copyright 2019 Alex Yu + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. diff --git a/README.md b/README.md new file mode 100644 index 0000000..a7d9333 --- /dev/null +++ b/README.md @@ -0,0 +1,50 @@ +## Fast marching method + +- To use in a project, include fmm.hpp and use `fmm::fmm(image, seeds, weight_map_type[, image_segmentation_threshold])` + +- This is header-only and you do not need to build this, unless you want to build the sample program. + +- All functions are templated to be OpenCV Mat compatible, also includes bare-bones Image struct which allows you to use any row-major contiguous matrix type: `Image(rows, cols, data_ptr)` + +Comment on top of fmm.hpp copied here: +```cpp +/** Fast marching method implementation + * usage: fmm::fmm(image, seeds, weight_map_type = IDENTITY, + * segmentation_threshold = disabled, + * normalize_output_geodesic_distances = true, + * output = nullptr) + * > returns an image, either geodesic distance map or, if + * segmentation_threshold is given, a segmentation mask + * image: input image (can be OpenCV Mat or fmm::Image) + * seeds: std::vector of points (each can be OpenCV Point or fmm::Point) + * weight_map_type: transformation to apply to input image to use as FMM +* weight function. Can be one of: + * fmm::weight::IDENTITY no transformation (avoids a copy) + * fmm::weight::GRADIENT gradient magnitude (using Sobel) + * fmm::weight::ABSDIFF absolute difference from average + * grayscale value of seeds + * fmm::weight::LAPLACIAN image Laplacian magnitude + * segmentation_threshold: if specified, sets pixels with geodesic value less + * than or equal to this threshold to 1 and others to 0 + * in the output image. If not given,the geodesic + * distances map will be returned. + * normalize_output_geodesic_distances: if true, normalizes geodesic distances + * values to be in the interval [0, 1]. + * If segmentation_threshold is specified, + * this occurs prior to segmentation. + * Default true. + * output: optionally, a pointer to an already-allocated output image. + * This allows you to avoid a copy if you already have one + * allocated. By default a new image is created, and this + * is not necessary. + * + * fmm::Image usage (optional) + * - To make owning image fmm::Image(rows, cols) + * - To make non-owning image that maps to row-major data (of same type, or char/uchar): + * fmm::Image(rows, cols, data_ptr) + * */ +``` + +## License + +Apache 2.0 diff --git a/intrin.txt b/intrin.txt new file mode 100644 index 0000000..cadafe4 --- /dev/null +++ b/intrin.txt @@ -0,0 +1,12 @@ +fx 606.438 +cx 637.294 +fy 606.351 +cy 366.992 +k0 0.777798 +k1 -2.93384 +k2 1.6463 +k3 0.655163 +k4 -2.76696 +k5 1.57894 +p0 0.000662754 +p1 6.69302e-05 diff --git a/partlabels.cpp b/partlabels.cpp new file mode 100644 index 0000000..de2a832 --- /dev/null +++ b/partlabels.cpp @@ -0,0 +1,497 @@ +#include +#include +#include +#include +#include +#include +#include +#include "opencv2/core.hpp" +#include "opencv2/imgproc.hpp" +#include "opencv2/highgui.hpp" +//#include +#include +#include +#include +#include + +#include "Avatar.h" +#include "Calibration.h" +#include "HumanDetector.h" +#include "Util.h" +//#include "concaveman.h" + +#define BEGIN_PROFILE auto start = std::chrono::high_resolution_clock::now() +#define PROFILE(x) do{printf("%s: %f ms\n", #x, std::chrono::duration(std::chrono::high_resolution_clock::now() - start).count()); start = std::chrono::high_resolution_clock::now(); }while(false) + +//#define DEBUG_PROPOSE +//#define DEBUG + +constexpr char WIND_NAME[] = "Result"; + +namespace { +using namespace ark; + +namespace random_util { + +template +/** mersenne twister PRNG */ +inline T randint_mt(T lo, T hi) { + thread_local static std::mt19937 rg(std::random_device{}()); + std::uniform_int_distribution pick(lo, hi); + return pick(rg); +} + +template +/** xorshift-based PRNG */ +inline T randint(T lo, T hi) { + if (hi <= lo) return lo; + static unsigned long x = std::random_device{}(), y = std::random_device{}(), z = std::random_device{}(); + unsigned long t; + x ^= x << 16; + x ^= x >> 5; + x ^= x << 1; + t = x; + x = y; + y = z; + z = t ^ x ^ y; + return z % (hi - lo + 1) + lo; +} + +inline double uniform(double min_inc = 0., double max_exc = 1.) { + thread_local static std::mt19937 rg(std::random_device{}()); + std::uniform_real_distribution normal(min_inc, max_exc); + return normal(rg); +} + +inline double randn(double mean = 0, double variance = 1) { + thread_local static std::mt19937 rg(std::random_device{}()); + std::normal_distribution normal(mean, variance); + return normal(rg); +} + +} // random_util + +typedef nanoflann::KDTreeEigenMatrixAdaptor kdTree; + + + +void getJoints(HumanAvatar& ava, + const CameraIntrin& intrin, + std::vector& joints_out) { + for (auto i = 0; i < HumanAvatar::NUM_JOINTS; ++i) { + auto pt = ava.getJointPosition(i); + joints_out.emplace_back( + pt.x() * intrin.fx / pt.z() + intrin.cx, + -pt.y() * intrin.fy / pt.z() + intrin.cy + ); + } +} + +inline void _fromSpherical(double rho, double theta, + double phi, Eigen::Vector3d& out) { + out[0] = rho * sin(phi) * cos(theta); + out[1] = rho * cos(phi); + out[2] = rho * sin(phi) * sin(theta); +} + +void randomizeParams(HumanAvatar& ava, double shape_sigma = 1.0) { + // Shape keys + Eigen::template Map > w(ava.w()); + for (int i = 0; i < HumanAvatar::NUM_SHAPEKEYS; ++i) { + w(i) = random_util::randn() * shape_sigma; + } + + // Pose + for (int c = 0; c < ava.posePrior.nComps; ++c) { + Eigen::Matrix r; + for (int i = 0; i < (HumanAvatar::NUM_JOINTS-1) * 3; ++i) { + r(i) = random_util::randn(); + } + r *= ava.posePrior.cov_cho[c]; + r += ava.posePrior.mean.row(c); + for (int i = 0; i < (HumanAvatar::NUM_JOINTS-1); ++i) { + auto joint = ava.getJoint(i + 1); + Eigen::AngleAxisd angle_axis; + angle_axis.angle() = r.segment<3>(i*3).norm(); + angle_axis.axis() = r.segment<3>(i*3)/angle_axis.angle(); + joint->rotation = angle_axis; + } + } + for (int c = 0; c < ava.posePrior.nComps; ++c) { + Eigen::Matrix r; + for (int i = 0; i < (HumanAvatar::NUM_JOINTS-1) * 3; ++i) { + r(i) = random_util::randn(); + } + r *= ava.posePrior.cov_cho[c]; + r += ava.posePrior.mean.row(c); + for (int i = 0; i < (HumanAvatar::NUM_JOINTS-1); ++i) { + auto joint = ava.getJoint(i + 1); + Eigen::AngleAxisd angle_axis; + angle_axis.angle() = r.segment<3>(i*3).norm(); + angle_axis.axis() = r.segment<3>(i*3)/angle_axis.angle(); + joint->rotation = angle_axis; + } + } + + // Root position + Eigen::Vector3d pos; + pos.x() = random_util::uniform(-1.0, 1.0); + pos.y() = random_util::uniform(-0.5, 0.5); + pos.z() = random_util::uniform(2.2, 4.5); + ava.setCenterPosition(pos); + + // Root rotation + const Eigen::Vector3d axis_up(0., 1., 0.); + double angle_up = random_util::uniform(0.0, std::sqrt(2 * PI)); + angle_up *= angle_up; + Eigen::AngleAxisd aa_up(angle_up, axis_up); + + double theta = random_util::uniform(0, 2 * PI); + double phi = random_util::uniform(-PI/2, PI/2); + Eigen::Vector3d axis_perturb; + _fromSpherical(1.0, theta, phi, axis_perturb); + double angle_perturb = random_util::randn(0.0, 0.2); + Eigen::AngleAxisd aa_perturb(angle_perturb, axis_perturb); + + ava.setCenterRotation(aa_perturb * aa_up); +} + +inline void paintTriangleNN( + cv::Mat& output_assigned_joint_mask, + const cv::Size& image_size, + const std::vector& projected, + const std::vector& assigned_joint, + const cv::Vec3i& face) { + std::pair xf[3] = + { + {projected[face[0]].x, 0}, + {projected[face[1]].x, 1}, + {projected[face[2]].x, 2} + }; + std::sort(xf, xf+3); + + // reorder points for convenience + auto a = projected[face[xf[0].second]], + b = projected[face[xf[1].second]], + c = projected[face[xf[2].second]]; + a.x = std::floor(a.x); + c.x = std::ceil(c.x); + if (a.x == c.x) return; + + auto assigned_a = assigned_joint[face[xf[0].second]], + assigned_b = assigned_joint[face[xf[1].second]], + assigned_c = assigned_joint[face[xf[2].second]]; + /* + const auto az = model_points[face[xf[0].second]].z, + bz = model_points[face[xf[1].second]].z, + cz = model_points[face[xf[2].second]].z; + */ + + int minxi = std::max(a.x, 0), + maxxi = std::min(c.x, image_size.width-1), + midxi = std::floor(b.x); + + if (a.x != b.x) { + double mhi = (c.y-a.y)/(c.x-a.x); + double bhi = a.y - a.x * mhi; + double mlo = (b.y-a.y)/(b.x-a.x); + double blo = a.y - a.x * mlo; + if (b.y > c.y) { + std::swap(mlo, mhi); + std::swap(blo, bhi); + } + for (int i = minxi; i <= std::min(midxi, image_size.width-1); ++i) { + int minyi = std::max(std::floor(mlo * i + blo), 0), + maxyi = std::min(std::ceil(mhi * i + bhi), image_size.height-1); + if (minyi > maxyi) continue; + + for (int j = minyi; j <= maxyi; ++j) { + auto& out = output_assigned_joint_mask.at(j, i); + int dista = (a.x - i) * (a.x - i) + (a.y - j) * (a.y - j); + int distb = (b.x - i) * (b.x - i) + (b.y - j) * (b.y - j); + int distc = (c.x - i) * (c.x - i) + (c.y - j) * (c.y - j); + if (dista < distb && dista < distc) { + out = assigned_a; + } else if (distb < distc) { + out = assigned_b; + } else { + out = assigned_c; + } + } + } + } + if (b.x != c.x) { + double mhi = (c.y-a.y)/(c.x-a.x); + double bhi = a.y - a.x * mhi; + double mlo = (c.y-b.y)/(c.x-b.x); + double blo = b.y - b.x * mlo; + if (b.y > a.y) { + std::swap(mlo, mhi); + std::swap(blo, bhi); + } + for (int i = std::max(midxi, 0)+1; i <= maxxi; ++i) { + int minyi = std::max(std::floor(mlo * i + blo), 0), + maxyi = std::min(std::ceil(mhi * i + bhi), image_size.height-1); + if (minyi > maxyi) continue; + + double w1v = (b.y - c.y) * (i - c.x); + double w2v = (c.y - a.y) * (i - c.x); + for (int j = minyi; j <= maxyi; ++j) { + auto& out = output_assigned_joint_mask.at(j, i); + int dista = (a.x - i) * (a.x - i) + (a.y - j) * (a.y - j); + int distb = (b.x - i) * (b.x - i) + (b.y - j) * (b.y - j); + int distc = (c.x - i) * (c.x - i) + (c.y - j) * (c.y - j); + if (dista < distb && dista < distc) { + out = assigned_a; + } else if (distb < distc) { + out = assigned_b; + } else { + out = assigned_c; + } + } + } + } +} + +void run(int num_threads, int max_num_to_gen, std::string out_path, const cv::Size& image_size, const CameraIntrin& intrin, int starting_number, bool overwrite) +{ + // Loadm mesh + std::ifstream meshIfs(util::resolveRootPath("data/avatar-model/mesh.txt")); + int nFaces; + meshIfs >> nFaces; + std::vector mesh; + mesh.reserve(nFaces); + for (int i = 0; i < nFaces;++i) { + mesh.emplace_back(); + auto& face = mesh.back(); + meshIfs >> face[0] >> face[1] >> face[2]; + } + + // Load joint assignments + std::ifstream skelIfs(util::resolveRootPath("data/avatar-model/skeleton.txt")); + int nJoints, nPoints; + skelIfs >> nJoints >> nPoints; + skelIfs.ignore(1000, '\n'); + std::string _; + for (int i = 0; i < nJoints; ++i) + std::getline(skelIfs, _); + std::vector assignedJoint(nPoints); + for (int i = 0; i < nPoints; ++i) { + int nAssign; skelIfs >> nAssign; + double largestWeight = -DBL_MAX; + for (int j = 0; j < nAssign; ++j) { + int id; double weight; + skelIfs >> id >> weight; + if (weight > largestWeight) { + largestWeight = weight; + assignedJoint[i] = id; + } + } + } + + using boost::filesystem::path; + path outPath (out_path); + path partMaskPath = outPath / "part_mask"; + path jointsPath = outPath / "joint"; + + if (!boost::filesystem::exists(outPath)) { + boost::filesystem::create_directories(outPath); + } + if (!boost::filesystem::exists(jointsPath)) { + boost::filesystem::create_directories(jointsPath); + } + if (!boost::filesystem::exists(partMaskPath)) { + boost::filesystem::create_directories(partMaskPath); + } + + boost::lockfree::queue que; + for (int i = starting_number; i < starting_number+max_num_to_gen; ++i) { + std::stringstream ss_img_id; + ss_img_id << std::setw(8) << std::setfill('0') << std::to_string(i); + auto jointFilePath = jointsPath / ("joint_" + ss_img_id.str() + ".yml"); + if (!boost::filesystem::exists(jointFilePath.string())) break; + if (!overwrite) { + auto partMaskFilePath = + partMaskPath / ("part_mask_" + ss_img_id.str() + ".png"); + if (boost::filesystem::exists(partMaskFilePath.string())) { + continue; + } + } + que.push(i); + } + + typedef std::pair FaceType; + auto faceComp = [](const FaceType& a, const FaceType& b) { + return a.first > b.first; + }; + + auto worker = [&]() { + HumanAvatar ava(HumanDetector::HUMAN_MODEL_PATH, HumanDetector::HUMAN_MODEL_SHAPE_KEYS); + std::vector projected(ava.getCloud()->size()); + + std::vector faces; + faces.reserve(nFaces); + for (int i = 0; i < nFaces;++i) { + faces.emplace_back(0.f, mesh[i]); + } + + while(true) { + int i; + if (!que.pop(i)) break; + std::stringstream ss_img_id; + ss_img_id << std::setw(8) << std::setfill('0') << std::to_string(i); + + std::string jointFilePath = (jointsPath / ("joint_" + ss_img_id.str() + ".yml")).string(); + + cv::FileStorage fs2(jointFilePath, cv::FileStorage::READ); + std::vector p; + fs2["pos"] >> p; + std::vector w; + fs2["shape"] >> w; + std::vector r; + fs2["rots"] >> r; + fs2.release(); + + std::copy(r.begin(), r.end(), ava.r()); + std::copy(w.begin(), w.end(), ava.w()); + std::copy(p.begin(), p.end(), ava.p()); + + ava.update(); + + auto modelCloud = ava.getCloud(); + auto& modelPoints = modelCloud->points; + + // Compute part labels + for (size_t i = 0; i < modelCloud->size(); ++i) { + const auto& pt = modelPoints[i]; + projected[i].x = static_cast(pt.x) + * intrin.fx / pt.z + intrin.cx; + projected[i].y = -static_cast(pt.y) * intrin.fy / pt.z + intrin.cy; + } + + // Sort faces by decreasing center depth + // so that when painted front faces will cover back faces + for (int i = 0; i < nFaces;++i) { + auto& face = faces[i].second; + faces[i].first = + (modelPoints[face[0]].z + modelPoints[face[1]].z + modelPoints[face[2]].z) / 3.f; + } + std::sort(faces.begin(), faces.end(), faceComp); + + // Paint the faces using nearest neighbors + cv::Mat partMaskMap(image_size, CV_8U); + partMaskMap.setTo(255); + for (int i = 0; i < nFaces;++i) { + paintTriangleNN(partMaskMap, image_size, projected, assignedJoint, faces[i].second); + } + + const std::string partMaskImgPath = (partMaskPath / ("part_mask_" + ss_img_id.str() + ".png")).string(); + cv::imwrite(partMaskImgPath, partMaskMap); + cout << "Wrote " << partMaskImgPath << endl; + } + }; + std::vector threads; + for (int i = 0; i < num_threads; ++i) { + threads.emplace_back(worker); + } + for (auto& t : threads) { + t.join(); + } +} +} + +int main(int argc, char** argv) { + namespace po = boost::program_options; + + std::string outPath, intrinPath; + int startingNumber, numToGen, numThreads; + bool overwrite; + cv::Size size; + + po::options_description desc("Option arguments"); + po::options_description descPositional("OpenARK Synthetic Avatar Dataset Part Label Mask Generator (c) Alex Yu 2019\nPosition arguments"); + po::options_description descCombined(""); + desc.add_options() + ("help", "produce help message") + ("overwrite,o", po::bool_switch(&overwrite), "If specified, overwrites existing files. Else, skips over them.") + (",j", po::value(&numThreads)->default_value(boost::thread::hardware_concurrency()), "Number of threads") + ; + + descPositional.add_options() + ("output_path", po::value(&outPath)->required(), "Input/Output DataSet Path") + ("num_to_gen", po::value(&numToGen)->default_value(INT_MAX), "Maximum number of images to generate") + ("starting_number", po::value(&startingNumber)->default_value(0), "Number of images to generate") + ("intrin_path", po::value(&intrinPath)->default_value(""), "Path to camera intrinsics file (default: uses hardcoded K4A intrinsics)") + ("width", po::value(&size.width)->default_value(1280), "Width of generated images") + ("height", po::value(&size.height)->default_value(720), "Height of generated imaes") + ; + descCombined.add(descPositional); + descCombined.add(desc); + po::variables_map vm; + + po::positional_options_description posopt; + posopt.add("output_path", 1); + posopt.add("num_to_gen", 1); + posopt.add("starting_number", 1); + posopt.add("intrin_path", 1); + posopt.add("width", 1); + posopt.add("height", 1); + try { + po::store(po::command_line_parser(argc, argv).options(descCombined) + .positional(posopt).run(), + vm); + } catch (std::exception& e) { + std::cerr << "Error: " << e.what() << "\n"; + std::cerr << descPositional << "\n" << desc << "\n"; + return 1; + } + + if ( vm.count("help") ) + { + std::cout << descPositional << "\n" << desc << "\n"; + return 0; + } + + try { + po::notify(vm); + } catch (std::exception& e) { + std::cerr << "Error: " << e.what() << "\n"; + std::cerr << descPositional << "\n" << desc << "\n"; + return 1; + } + + CameraIntrin intrin; + if (!intrinPath.empty()) intrin.readFile(intrinPath); + else { + intrin.clear(); + intrin.fx = 606.438; + intrin.fy = 606.351; + intrin.cx = 637.294; + intrin.cy = 366.992; + } + /* + auto viewer = Visualizer::getPCLVisualizer(); + HumanAvatar ava(HumanDetector::HUMAN_MODEL_PATH, HumanDetector::HUMAN_MODEL_SHAPE_KEYS); + pcl::PolygonMesh mesh; + auto cloud = ava.getCloud(); + pcl::toPCLPointCloud2(*cloud, mesh.cloud); + std::ifstream meshIfs(util::resolveRootPath("data/avatar-model/mesh.txt")); + int nFaces; + meshIfs >> nFaces; + for (int i = 0; i < nFaces;++i) { + mesh.polygons.emplace_back(); + auto& face = mesh.polygons.back().vertices; + face.resize(3); + meshIfs >> face[0] >> face[1] >> face[2]; + } + viewer->setBackgroundColor(0, 0, 0); + viewer->addPolygonMesh(mesh,"meshes",0); + viewer->addCoordinateSystem(1.0); + viewer->initCameraParameters(); + viewer->spin(); + */ + + run(numThreads, numToGen, outPath, size, intrin, + startingNumber, overwrite); + return 0; +} diff --git a/sample.cpp b/sample.cpp new file mode 100644 index 0000000..b0cfcc9 --- /dev/null +++ b/sample.cpp @@ -0,0 +1,482 @@ +#include +#include +#include +#include +#include +#include +#include +#include "opencv2/core.hpp" +#include "opencv2/imgproc.hpp" +#include "opencv2/highgui.hpp" +//#include +#include +#include +#include +#include + +#include "Avatar.h" +#include "Calibration.h" +#include "HumanDetector.h" +#include "Util.h" + +//#define DEBUG_PROPOSE +//#define DEBUG + +constexpr char WIND_NAME[] = "Result"; + +namespace { +using namespace ark; + +namespace random_util { + +template +/** mersenne twister PRNG */ +inline T randint_mt(T lo, T hi) { + thread_local static std::mt19937 rg(std::random_device{}()); + std::uniform_int_distribution pick(lo, hi); + return pick(rg); +} + +template +/** xorshift-based PRNG */ +inline T randint(T lo, T hi) { + if (hi <= lo) return lo; + static unsigned long x = std::random_device{}(), y = std::random_device{}(), z = std::random_device{}(); + unsigned long t; + x ^= x << 16; + x ^= x >> 5; + x ^= x << 1; + t = x; + x = y; + y = z; + z = t ^ x ^ y; + return z % (hi - lo + 1) + lo; +} + +inline double uniform(double min_inc = 0., double max_exc = 1.) { + thread_local static std::mt19937 rg(std::random_device{}()); + std::uniform_real_distribution normal(min_inc, max_exc); + return normal(rg); +} + +inline double randn(double mean = 0, double variance = 1) { + thread_local static std::mt19937 rg(std::random_device{}()); + std::normal_distribution normal(mean, variance); + return normal(rg); +} + +} // random_util + +typedef nanoflann::KDTreeEigenMatrixAdaptor kdTree; + +void getJoints(HumanAvatar& ava, + const CameraIntrin& intrin, + std::vector& joints_out) { + for (auto i = 0; i < HumanAvatar::NUM_JOINTS; ++i) { + auto pt = ava.getJointPosition(i); + joints_out.emplace_back( + pt.x() * intrin.fx / pt.z() + intrin.cx, + -pt.y() * intrin.fy / pt.z() + intrin.cy + ); + } +} + +inline void _fromSpherical(double rho, double theta, + double phi, Eigen::Vector3d& out) { + out[0] = rho * sin(phi) * cos(theta); + out[1] = rho * cos(phi); + out[2] = rho * sin(phi) * sin(theta); +} + +void randomizeParams(HumanAvatar& ava, double shape_sigma = 1.0) { + // Shape keys + Eigen::template Map > w(ava.w()); + for (int i = 0; i < HumanAvatar::NUM_SHAPEKEYS; ++i) { + w(i) = random_util::randn() * shape_sigma; + } + + // Pose + for (int c = 0; c < ava.posePrior.nComps; ++c) { + Eigen::Matrix r; + for (int i = 0; i < (HumanAvatar::NUM_JOINTS-1) * 3; ++i) { + r(i) = random_util::randn(); + } + r *= ava.posePrior.cov_cho[c]; + r += ava.posePrior.mean.row(c); + for (int i = 0; i < (HumanAvatar::NUM_JOINTS-1); ++i) { + auto joint = ava.getJoint(i + 1); + Eigen::AngleAxisd angle_axis; + angle_axis.angle() = r.segment<3>(i*3).norm(); + angle_axis.axis() = r.segment<3>(i*3)/angle_axis.angle(); + joint->rotation = angle_axis; + } + } + for (int c = 0; c < ava.posePrior.nComps; ++c) { + Eigen::Matrix r; + for (int i = 0; i < (HumanAvatar::NUM_JOINTS-1) * 3; ++i) { + r(i) = random_util::randn(); + } + r *= ava.posePrior.cov_cho[c]; + r += ava.posePrior.mean.row(c); + for (int i = 0; i < (HumanAvatar::NUM_JOINTS-1); ++i) { + auto joint = ava.getJoint(i + 1); + Eigen::AngleAxisd angle_axis; + angle_axis.angle() = r.segment<3>(i*3).norm(); + angle_axis.axis() = r.segment<3>(i*3)/angle_axis.angle(); + joint->rotation = angle_axis; + } + } + + // Root position + Eigen::Vector3d pos; + pos.x() = random_util::uniform(-1.0, 1.0); + pos.y() = random_util::uniform(-0.5, 0.5); + pos.z() = random_util::uniform(2.2, 4.5); + ava.setCenterPosition(pos); + + // Root rotation + const Eigen::Vector3d axis_up(0., 1., 0.); + double angle_up = random_util::uniform(0.0, std::sqrt(2 * PI)); + angle_up *= angle_up; + Eigen::AngleAxisd aa_up(angle_up, axis_up); + + double theta = random_util::uniform(0, 2 * PI); + double phi = random_util::uniform(-PI/2, PI/2); + Eigen::Vector3d axis_perturb; + _fromSpherical(1.0, theta, phi, axis_perturb); + double angle_perturb = random_util::randn(0.0, 0.2); + Eigen::AngleAxisd aa_perturb(angle_perturb, axis_perturb); + + ava.setCenterRotation(aa_perturb * aa_up); +} + +inline void paintTriangleBary( + cv::Mat& output_depth, + const cv::Size& image_size, + const std::vector& projected, + const std::vector>& + model_points, + const cv::Vec3i& face) { + std::pair xf[3] = + { + {projected[face[0]].x, 0}, + {projected[face[1]].x, 1}, + {projected[face[2]].x, 2} + }; + std::sort(xf, xf+3); + + // reorder points for convenience + auto a = projected[face[xf[0].second]], + b = projected[face[xf[1].second]], + c = projected[face[xf[2].second]]; + a.x = std::floor(a.x); + c.x = std::ceil(c.x); + if (a.x == c.x) return; + const auto az = model_points[face[xf[0].second]].z, + bz = model_points[face[xf[1].second]].z, + cz = model_points[face[xf[2].second]].z; + + int minxi = std::max(a.x, 0), + maxxi = std::min(c.x, image_size.width-1), + midxi = std::floor(b.x); + + double denom = 1.0 / ((b.y - c.y) * (a.x - c.x) + (c.x - b.x) * (a.y - c.y)); + if (a.x != b.x) { + double mhi = (c.y-a.y)/(c.x-a.x); + double bhi = a.y - a.x * mhi; + double mlo = (b.y-a.y)/(b.x-a.x); + double blo = a.y - a.x * mlo; + if (b.y > c.y) { + std::swap(mlo, mhi); + std::swap(blo, bhi); + } + for (int i = minxi; i <= std::min(midxi, image_size.width-1); ++i) { + int minyi = std::max(std::floor(mlo * i + blo), 0), + maxyi = std::min(std::ceil(mhi * i + bhi), image_size.height-1); + if (minyi > maxyi) continue; + + double w1v = (b.y - c.y) * (i - c.x); + double w2v = (c.y - a.y) * (i - c.x); + for (int j = minyi; j <= maxyi; ++j) { + float w1 = (w1v + (c.x - b.x) * (j - c.y)) * denom; + float w2 = (w2v + (a.x - c.x) * (j - c.y)) * denom; + output_depth.at(j, i) = + w1 * az + w2 * bz + (1. - w1 - w2) * cz; + } + } + } + if (b.x != c.x) { + double mhi = (c.y-a.y)/(c.x-a.x); + double bhi = a.y - a.x * mhi; + double mlo = (c.y-b.y)/(c.x-b.x); + double blo = b.y - b.x * mlo; + if (b.y > a.y) { + std::swap(mlo, mhi); + std::swap(blo, bhi); + } + for (int i = std::max(midxi, 0)+1; i <= maxxi; ++i) { + int minyi = std::max(std::floor(mlo * i + blo), 0), + maxyi = std::min(std::ceil(mhi * i + bhi), image_size.height-1); + if (minyi > maxyi) continue; + + double w1v = (b.y - c.y) * (i - c.x); + double w2v = (c.y - a.y) * (i - c.x); + for (int j = minyi; j <= maxyi; ++j) { + float w1 = (w1v + (c.x - b.x) * (j - c.y)) * denom; + float w2 = (w2v + (a.x - c.x) * (j - c.y)) * denom; + output_depth.at(j, i) = + w1 * az + w2 * bz + (1. - w1 - w2) * cz; + } + } + } +} + +void run(int num_threads, int num_to_gen, std::string out_path, const cv::Size& image_size, const CameraIntrin& intrin, int starting_number, bool overwrite) +{ + std::ifstream meshIfs(util::resolveRootPath("data/avatar-model/mesh.txt")); + int nFaces; + meshIfs >> nFaces; + std::vector mesh; + mesh.reserve(nFaces); + for (int i = 0; i < nFaces;++i) { + mesh.emplace_back(); + auto& face = mesh.back(); + meshIfs >> face[0] >> face[1] >> face[2]; + } + + using boost::filesystem::path; + path outPath (out_path); + path intrinPath = outPath / "intrin.txt"; + path depthPath = outPath / "depth_exr"; + path jointsPath = outPath / "joint"; + + if (!boost::filesystem::exists(outPath)) { + boost::filesystem::create_directories(outPath); + } + if (!boost::filesystem::exists(depthPath)) { + boost::filesystem::create_directories(depthPath); + } if (!boost::filesystem::exists(jointsPath)) { + boost::filesystem::create_directories(jointsPath); + } + + intrin.writeFile(intrinPath.string()); + + boost::lockfree::queue que; + for (int i = starting_number; i < starting_number+num_to_gen; ++i) { + if (!overwrite) { + std::stringstream ss_img_id; + ss_img_id << std::setw(8) << std::setfill('0') << std::to_string(i); + auto depthFilePath = + depthPath / ("depth_" + ss_img_id.str() + ".exr"); + std::ifstream testIfs(depthFilePath.string()); + if (testIfs) { + continue; + } + } + que.push(i); + } + + typedef std::pair FaceType; + auto faceComp = [](const FaceType& a, const FaceType& b) { + return a.first > b.first; + }; + + auto worker = [&]() { + HumanAvatar ava(HumanDetector::HUMAN_MODEL_PATH, HumanDetector::HUMAN_MODEL_SHAPE_KEYS); + std::vector projected(ava.getCloud()->size()); + + std::vector faces; + faces.reserve(nFaces); + for (int i = 0; i < nFaces;++i) { + faces.emplace_back(0.f, mesh[i]); + } + + while(true) { + int i; + if (!que.pop(i)) break; + std::stringstream ss_img_id; + ss_img_id << std::setw(8) << std::setfill('0') << std::to_string(i); + + randomizeParams(ava); + + ava.update(); + + auto modelCloud = ava.getCloud(); + auto& modelPoints = modelCloud->points; + + // Render depth + for (size_t i = 0; i < modelCloud->size(); ++i) { + const auto& pt = modelPoints[i]; + projected[i].x = static_cast(pt.x) + * intrin.fx / pt.z + intrin.cx; + projected[i].y = -static_cast(pt.y) * intrin.fy / pt.z + intrin.cy; + } + + // Sort faces by decreasing center depth + // so that when painted front faces will cover back faces + for (int i = 0; i < nFaces;++i) { + auto& face = faces[i].second; + faces[i].first = + (modelPoints[face[0]].z + modelPoints[face[1]].z + modelPoints[face[2]].z) / 3.f; + } + std::sort(faces.begin(), faces.end(), faceComp); + + // Paint the faces using Barycentric coordinate interpolation + cv::Mat renderedDepth = cv::Mat::zeros(image_size, CV_32F); + for (int i = 0; i < nFaces;++i) { + paintTriangleBary(renderedDepth, image_size, projected, modelPoints, faces[i].second); + } + //cv::Mat visual; + //cv::normalize(renderedDepth, visual, 0.0, 1.0, cv::NORM_MINMAX); + //visual.convertTo(visual, CV_8UC1, 255.0); + //cv::applyColorMap(visual, visual, cv::COLORMAP_HOT); + //cv::imshow(WIND_NAME, visual); + // cv::waitKey(0); + + const std::string depthImgPath = (depthPath / ("depth_" + ss_img_id.str() + ".exr")).string(); + cv::imwrite(depthImgPath, renderedDepth); + cout << "Wrote " << depthImgPath << endl; + + // Output labels + std::vector joints; + getJoints(ava, intrin, joints); + + const std::string jointFilePath = (jointsPath / ("joint_" + ss_img_id.str() + ".yml")).string(); + cv::FileStorage fs3(jointFilePath, cv::FileStorage::WRITE); + fs3 << "joints" << joints; + + // Also write xyz positions + std::vector jointsXYZ; + for (auto i = 0; i < HumanAvatar::NUM_JOINTS; ++i) { + auto pt = ava.getJointPosition(i); + jointsXYZ.emplace_back(pt.x(), pt.y(), pt.z()); + } + fs3 << "joints_xyz" << jointsXYZ; + + // Also write OpenARK avatar parameters + cv::Point3d p(ava.p()[0], ava.p()[1], ava.p()[2]); + fs3 << "pos" << p; + + std::vector w(HumanAvatar::NUM_SHAPEKEYS); + std::copy(ava.w(), ava.w() + w.size(), w.begin()); + fs3 << "shape" << w; + + std::vector r(HumanAvatar::NUM_JOINTS * HumanAvatar::NUM_ROT_PARAMS); + std::copy(ava.r(), ava.r() + r.size(), r.begin()); + fs3 << "rots" << r; + + // Convert to SMPL parameters + Eigen::VectorXd smplParams = ava.smplParams(); + std::vector smplParamsVec(smplParams.rows()); + std::copy(smplParams.data(), smplParams.data() + smplParams.rows(), smplParamsVec.begin()); + fs3 << "smpl_params" << smplParamsVec; + + fs3.release(); + } + }; + std::vector threads; + for (int i = 0; i < num_threads; ++i) { + threads.emplace_back(worker); + } + for (auto& t : threads) { + t.join(); + } +} +} + +int main(int argc, char** argv) { + namespace po = boost::program_options; + + std::string outPath, intrinPath; + int startingNumber, numToGen, numThreads; + bool overwrite; + cv::Size size; + + po::options_description desc("Option arguments"); + po::options_description descPositional("OpenARK Synthetic Avatar Depth Image Dataset Generator (c) Alex Yu 2019\nPosition arguments"); + po::options_description descCombined(""); + desc.add_options() + ("help", "produce help message") + ("overwrite,o", po::bool_switch(&overwrite), "If specified, overwrites existing files. Else, skips over them.") + (",j", po::value(&numThreads)->default_value(boost::thread::hardware_concurrency()), "Number of threads") + ; + + descPositional.add_options() + ("output_path", po::value(&outPath)->required(), "Output Path") + ("num_to_gen", po::value(&numToGen)->default_value(1), "Number of images to generate") + ("starting_number", po::value(&startingNumber)->default_value(0), "Number of images to generate") + ("intrin_path", po::value(&intrinPath)->default_value(""), "Path to camera intrinsics file (default: uses hardcoded K4A intrinsics)") + ("width", po::value(&size.width)->default_value(1280), "Width of generated images") + ("height", po::value(&size.height)->default_value(720), "Height of generated imaes") + ; + descCombined.add(descPositional); + descCombined.add(desc); + po::variables_map vm; + + po::positional_options_description posopt; + posopt.add("output_path", 1); + posopt.add("num_to_gen", 1); + posopt.add("starting_number", 1); + posopt.add("intrin_path", 1); + posopt.add("width", 1); + posopt.add("height", 1); + try { + po::store(po::command_line_parser(argc, argv).options(descCombined) + .positional(posopt).run(), + vm); + } catch (std::exception& e) { + std::cerr << "Error: " << e.what() << "\n"; + std::cerr << descPositional << "\n" << desc << "\n"; + return 1; + } + + if ( vm.count("help") ) + { + std::cout << descPositional << "\n" << desc << "\n"; + return 0; + } + + try { + po::notify(vm); + } catch (std::exception& e) { + std::cerr << "Error: " << e.what() << "\n"; + std::cerr << descPositional << "\n" << desc << "\n"; + return 1; + } + + CameraIntrin intrin; + if (!intrinPath.empty()) intrin.readFile(intrinPath); + else { + intrin.clear(); + intrin.fx = 606.438; + intrin.fy = 606.351; + intrin.cx = 637.294; + intrin.cy = 366.992; + } + /* + auto viewer = Visualizer::getPCLVisualizer(); + HumanAvatar ava(HumanDetector::HUMAN_MODEL_PATH, HumanDetector::HUMAN_MODEL_SHAPE_KEYS); + pcl::PolygonMesh mesh; + auto cloud = ava.getCloud(); + pcl::toPCLPointCloud2(*cloud, mesh.cloud); + std::ifstream meshIfs(util::resolveRootPath("data/avatar-model/mesh.txt")); + int nFaces; + meshIfs >> nFaces; + for (int i = 0; i < nFaces;++i) { + mesh.polygons.emplace_back(); + auto& face = mesh.polygons.back().vertices; + face.resize(3); + meshIfs >> face[0] >> face[1] >> face[2]; + } + viewer->setBackgroundColor(0, 0, 0); + viewer->addPolygonMesh(mesh,"meshes",0); + viewer->addCoordinateSystem(1.0); + viewer->initCameraParameters(); + viewer->spin(); + */ + + run(numThreads, numToGen, outPath, size, intrin, + startingNumber, overwrite); + return 0; +}