38 #ifndef PCL_REGISTRATION_IMPL_IA_FPCS_H_
39 #define PCL_REGISTRATION_IMPL_IA_FPCS_H_
43 #include <pcl/common/utils.h>
44 #include <pcl/registration/ia_fpcs.h>
45 #include <pcl/registration/transformation_estimation_3point.h>
46 #include <pcl/sample_consensus/sac_model_plane.h>
51 template <
typename Po
intT>
57 const float max_dist_sqr = max_dist * max_dist;
58 const std::size_t s = cloud->
size();
63 float mean_dist = 0.f;
66 std::vector<float> dists_sqr(2);
69 #pragma omp parallel for \
72 firstprivate(ids, dists_sqr) \
73 reduction(+:mean_dist, num) \
74 firstprivate(s, max_dist_sqr) \
75 num_threads(nr_threads)
76 for (
int i = 0; i < 1000; i++) {
78 if (dists_sqr[1] < max_dist_sqr) {
79 mean_dist += std::sqrt(dists_sqr[1]);
84 return (mean_dist / num);
88 template <
typename Po
intT>
95 const float max_dist_sqr = max_dist * max_dist;
96 const std::size_t s = indices.size();
101 float mean_dist = 0.f;
104 std::vector<float> dists_sqr(2);
107 #if OPENMP_LEGACY_CONST_DATA_SHARING_RULE
108 #pragma omp parallel for \
110 shared(tree, cloud, indices) \
111 firstprivate(ids, dists_sqr) \
112 reduction(+:mean_dist, num) \
113 num_threads(nr_threads)
115 #pragma omp parallel for \
117 shared(tree, cloud, indices, s, max_dist_sqr) \
118 firstprivate(ids, dists_sqr) \
119 reduction(+:mean_dist, num) \
120 num_threads(nr_threads)
122 for (
int i = 0; i < 1000; i++) {
123 tree.
nearestKSearch((*cloud)[indices[rand() % s]], 2, ids, dists_sqr);
124 if (dists_sqr[1] < max_dist_sqr) {
125 mean_dist += std::sqrt(dists_sqr[1]);
130 return (mean_dist / num);
134 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
140 , approx_overlap_(0.5f)
142 , score_threshold_(std::numeric_limits<float>::max())
144 , max_norm_diff_(90.f)
146 , fitness_score_(std::numeric_limits<float>::max())
148 , max_base_diameter_sqr_()
149 , use_normals_(false)
150 , normalize_delta_(true)
153 , coincidation_limit_()
155 , max_inlier_dist_sqr_()
156 , small_error_(0.00001f)
158 reg_name_ =
"pcl::registration::FPCSInitialAlignment";
166 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
174 final_transformation_ = guess;
176 std::vector<MatchingCandidates> all_candidates(max_iterations_);
179 #pragma omp parallel default(none) shared(abort, all_candidates, timer) \
180 num_threads(nr_threads_)
183 const unsigned int seed =
184 static_cast<unsigned int>(std::time(NULL)) ^ omp_get_thread_num();
186 PCL_DEBUG(
"[%s::computeTransformation] Using seed=%u\n", reg_name_.c_str(), seed);
187 #pragma omp for schedule(dynamic)
189 for (
int i = 0; i < max_iterations_; i++) {
190 #pragma omp flush(abort)
194 all_candidates[i] = candidates;
199 if (selectBase(base_indices, ratio) == 0) {
202 if (bruteForceCorrespondences(base_indices[0], base_indices[1], pairs_a) ==
204 bruteForceCorrespondences(base_indices[2], base_indices[3], pairs_b) ==
208 std::vector<pcl::Indices> matches;
209 if (determineBaseMatches(base_indices, matches, pairs_a, pairs_b, ratio) ==
212 handleMatches(base_indices, matches, candidates);
213 if (!candidates.empty())
214 all_candidates[i] = candidates;
220 abort = (!candidates.empty() ? candidates[0].fitness_score < score_threshold_
224 #pragma omp flush(abort)
230 finalCompute(all_candidates);
239 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
244 const unsigned int seed = std::time(
nullptr);
246 PCL_DEBUG(
"[%s::initCompute] Using seed=%u\n", reg_name_.c_str(), seed);
253 if (!input_ || !target_) {
254 PCL_ERROR(
"[%s::initCompute] Source or target dataset not given!\n",
259 if (!target_indices_ || target_indices_->empty()) {
260 target_indices_.reset(
new pcl::Indices(target_->size()));
262 for (
auto& target_index : *target_indices_)
263 target_index = index++;
264 target_cloud_updated_ =
true;
269 if (nr_samples_ != 0) {
270 const int ss =
static_cast<int>(indices_->size());
271 const int sample_fraction_src = std::max(1,
static_cast<int>(ss / nr_samples_));
274 for (
int i = 0; i < ss; i++)
275 if (rand() % sample_fraction_src == 0)
276 source_indices_->push_back((*indices_)[i]);
279 source_indices_ = indices_;
282 if (source_normals_ && target_normals_ && source_normals_->size() == input_->size() &&
283 target_normals_->size() == target_->size())
287 if (target_cloud_updated_) {
288 tree_->setInputCloud(target_, target_indices_);
289 target_cloud_updated_ =
false;
293 constexpr
int min_iterations = 4;
294 constexpr
float diameter_fraction = 0.3f;
297 Eigen::Vector4f pt_min, pt_max;
299 diameter_ = (pt_max - pt_min).norm();
302 float max_base_diameter = diameter_ * approx_overlap_ * 2.f;
303 max_base_diameter_sqr_ = max_base_diameter * max_base_diameter;
306 if (normalize_delta_) {
307 float mean_dist = getMeanPointDensity<PointTarget>(
308 target_, *target_indices_, 0.05f * diameter_, nr_threads_);
314 if (max_iterations_ == 0) {
315 float first_est = std::log(small_error_) /
316 std::log(1.0 - std::pow(
static_cast<double>(approx_overlap_),
317 static_cast<double>(min_iterations)));
319 static_cast<int>(first_est / (diameter_fraction * approx_overlap_ * 2.f));
323 if (score_threshold_ == std::numeric_limits<float>::max())
324 score_threshold_ = 1.f - approx_overlap_;
326 if (max_iterations_ < 4)
329 if (max_runtime_ < 1)
330 max_runtime_ = std::numeric_limits<int>::max();
333 max_pair_diff_ = delta_ * 2.f;
334 max_edge_diff_ = delta_ * 4.f;
335 coincidation_limit_ = delta_ * 2.f;
336 max_mse_ = powf(delta_ * 2.f, 2.f);
337 max_inlier_dist_sqr_ = powf(delta_ * 2.f, 2.f);
340 fitness_score_ = std::numeric_limits<float>::max();
346 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
351 const float too_close_sqr = max_base_diameter_sqr_ * 0.01;
353 Eigen::VectorXf coefficients(4);
356 Eigen::Vector4f centre_pt;
357 float nearest_to_plane = std::numeric_limits<float>::max();
361 for (
int i = 0; i < ransac_iterations_; i++) {
363 if (selectBaseTriangle(base_indices) < 0)
366 pcl::Indices base_triple(base_indices.begin(), base_indices.end() - 1);
371 const PointTarget* pt1 = &((*target_)[base_indices[0]]);
372 const PointTarget* pt2 = &((*target_)[base_indices[1]]);
373 const PointTarget* pt3 = &((*target_)[base_indices[2]]);
375 for (
const auto& target_index : *target_indices_) {
376 const PointTarget* pt4 = &((*target_)[target_index]);
381 float d4 = (pt4->getVector3fMap() - centre_pt.head(3)).squaredNorm();
385 if (d1 < too_close_sqr || d2 < too_close_sqr || d3 < too_close_sqr ||
386 d4 < too_close_sqr || d1 > max_base_diameter_sqr_ ||
387 d2 > max_base_diameter_sqr_ || d3 > max_base_diameter_sqr_)
392 if (dist_to_plane < nearest_to_plane) {
393 base_indices[3] = target_index;
394 nearest_to_plane = dist_to_plane;
399 if (nearest_to_plane != std::numeric_limits<float>::max()) {
402 setupBase(base_indices, ratio);
412 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
417 const auto nr_points = target_indices_->size();
421 base_indices[0] = (*target_indices_)[rand() % nr_points];
422 auto* index1 = base_indices.data();
425 for (
int i = 0; i < ransac_iterations_; i++) {
426 auto* index2 = &(*target_indices_)[rand() % nr_points];
427 auto* index3 = &(*target_indices_)[rand() % nr_points];
430 (*target_)[*index2].getVector3fMap() - (*target_)[*index1].getVector3fMap();
432 (*target_)[*index3].getVector3fMap() - (*target_)[*index1].getVector3fMap();
434 u.cross(v).squaredNorm();
437 if (t > best_t && u.squaredNorm() < max_base_diameter_sqr_ &&
438 v.squaredNorm() < max_base_diameter_sqr_) {
440 base_indices[1] = *index2;
441 base_indices[2] = *index3;
446 return (best_t == 0.f ? -1 : 0);
450 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
455 float best_t = std::numeric_limits<float>::max();
456 const pcl::Indices copy(base_indices.begin(), base_indices.end());
457 pcl::Indices temp(base_indices.begin(), base_indices.end());
460 for (
auto i = copy.begin(), i_e = copy.end(); i != i_e; ++i)
461 for (
auto j = copy.begin(), j_e = copy.end(); j != j_e; ++j) {
465 for (
auto k = copy.begin(), k_e = copy.end(); k != k_e; ++k) {
466 if (k == j || k == i)
469 auto l = copy.begin();
470 while (l == i || l == j || l == k)
481 float t = segmentToSegmentDist(temp, ratio_temp);
484 ratio[0] = ratio_temp[0];
485 ratio[1] = ratio_temp[1];
493 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
499 Eigen::Vector3f u = (*target_)[base_indices[1]].getVector3fMap() -
500 (*target_)[base_indices[0]].getVector3fMap();
501 Eigen::Vector3f v = (*target_)[base_indices[3]].getVector3fMap() -
502 (*target_)[base_indices[2]].getVector3fMap();
503 Eigen::Vector3f w = (*target_)[base_indices[0]].getVector3fMap() -
504 (*target_)[base_indices[2]].getVector3fMap();
512 float D = a * c - b * b;
513 float sN = 0.f, sD = D;
514 float tN = 0.f, tD = D;
517 if (D < small_error_) {
524 sN = (b * e - c * d);
525 tN = (a * e - b * d);
560 else if ((-d + b) > a)
570 ratio[0] = (std::abs(sN) < small_error_) ? 0.f : sN / sD;
571 ratio[1] = (std::abs(tN) < small_error_) ? 0.f : tN / tD;
573 Eigen::Vector3f x = w + (ratio[0] * u) - (ratio[1] * v);
578 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
583 const float max_norm_diff = 0.5f * max_norm_diff_ *
M_PI / 180.f;
587 float ref_norm_angle =
588 (use_normals_ ? ((*target_normals_)[idx1].getNormalVector3fMap() -
589 (*target_normals_)[idx2].getNormalVector3fMap())
594 auto it_out = source_indices_->begin(), it_out_e = source_indices_->end() - 1;
595 auto it_in_e = source_indices_->end();
596 for (; it_out != it_out_e; it_out++) {
597 auto it_in = it_out + 1;
598 const PointSource* pt1 = &(*input_)[*it_out];
599 for (; it_in != it_in_e; it_in++) {
600 const PointSource* pt2 = &(*input_)[*it_in];
604 if (std::abs(dist - ref_dist) < max_pair_diff_) {
607 const NormalT* pt1_n = &((*source_normals_)[*it_out]);
608 const NormalT* pt2_n = &((*source_normals_)[*it_in]);
611 (pt1_n->getNormalVector3fMap() - pt2_n->getNormalVector3fMap()).norm();
613 (pt1_n->getNormalVector3fMap() + pt2_n->getNormalVector3fMap()).norm();
615 float norm_diff = std::min<float>(std::abs(norm_angle_1 - ref_norm_angle),
616 std::abs(norm_angle_2 - ref_norm_angle));
617 if (norm_diff > max_norm_diff)
628 return (pairs.empty() ? -1 : 0);
632 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
636 std::vector<pcl::Indices>& matches,
639 const float (&ratio)[2])
655 cloud_e->resize(pairs_a.size() * 2);
656 auto it_pt = cloud_e->begin();
657 for (
const auto& pair : pairs_a) {
658 const PointSource* pt1 = &((*input_)[pair.index_match]);
659 const PointSource* pt2 = &((*input_)[pair.index_query]);
662 for (
int i = 0; i < 2; i++, it_pt++) {
663 it_pt->x = pt1->x + ratio[i] * (pt2->x - pt1->x);
664 it_pt->y = pt1->y + ratio[i] * (pt2->y - pt1->y);
665 it_pt->z = pt1->z + ratio[i] * (pt2->z - pt1->z);
671 tree_e->setInputCloud(cloud_e);
674 std::vector<float> dists_sqr;
677 for (
const auto& pair : pairs_b) {
678 const PointTarget* pt1 = &((*input_)[pair.index_match]);
679 const PointTarget* pt2 = &((*input_)[pair.index_query]);
682 for (
const float& r : ratio) {
684 pt_e.x = pt1->x + r * (pt2->x - pt1->x);
685 pt_e.y = pt1->y + r * (pt2->y - pt1->y);
686 pt_e.z = pt1->z + r * (pt2->z - pt1->z);
689 tree_e->radiusSearch(pt_e, coincidation_limit_, ids, dists_sqr);
690 for (
const auto&
id : ids) {
694 pairs_a[
static_cast<int>(std::floor((
id / 2.f)))].index_match;
696 pairs_a[
static_cast<int>(std::floor((
id / 2.f)))].index_query;
697 match_indices[2] = pair.index_match;
698 match_indices[3] = pair.index_query;
701 if (checkBaseMatch(match_indices, dist_base) < 0)
704 matches.push_back(match_indices);
710 return (!matches.empty() ? 0 : -1);
714 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
729 return (std::abs(d0 - dist_ref[0]) < max_edge_diff_ &&
730 std::abs(d1 - dist_ref[1]) < max_edge_diff_ &&
731 std::abs(d2 - dist_ref[2]) < max_edge_diff_ &&
732 std::abs(d3 - dist_ref[3]) < max_edge_diff_)
738 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
742 std::vector<pcl::Indices>& matches,
745 candidates.resize(1);
746 float fitness_score = std::numeric_limits<float>::max();
749 for (
auto& match : matches) {
750 Eigen::Matrix4f transformation_temp;
755 linkMatchWithBase(base_indices, match, correspondences_temp);
758 if (validateMatch(base_indices, match, correspondences_temp, transformation_temp) <
764 if (validateTransformation(transformation_temp, fitness_score) < 0)
768 candidates[0].fitness_score = fitness_score;
769 candidates[0].transformation = transformation_temp;
770 correspondences_temp.erase(correspondences_temp.end() - 1);
771 candidates[0].correspondences = correspondences_temp;
776 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
784 Eigen::Vector4f centre_base{0, 0, 0, 0}, centre_match{0, 0, 0, 0};
788 PointTarget centre_pt_base;
789 centre_pt_base.x = centre_base[0];
790 centre_pt_base.y = centre_base[1];
791 centre_pt_base.z = centre_base[2];
793 PointSource centre_pt_match;
794 centre_pt_match.x = centre_match[0];
795 centre_pt_match.y = centre_match[1];
796 centre_pt_match.z = centre_match[2];
801 auto it_match_orig = match_indices.begin();
802 for (
auto it_base = base_indices.cbegin(), it_base_e = base_indices.cend();
803 it_base != it_base_e;
804 it_base++, it_match_orig++) {
807 float best_diff_sqr = std::numeric_limits<float>::max();
810 for (
const auto& match_index : copy) {
814 float diff_sqr = std::abs(dist_sqr_1 - dist_sqr_2);
816 if (diff_sqr < best_diff_sqr) {
817 best_diff_sqr = diff_sqr;
818 best_index = match_index;
824 *it_match_orig = best_index;
829 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
835 Eigen::Matrix4f& transformation)
839 correspondences_temp.erase(correspondences_temp.end() - 1);
842 transformation_estimation_->estimateRigidTransformation(
843 *input_, *target_, correspondences_temp, transformation);
846 PointCloudSource match_transformed;
850 std::size_t nr_points = correspondences_temp.size();
852 for (std::size_t i = 0; i < nr_points; i++)
854 target_->points[base_indices[i]]);
857 return (mse < max_mse_ ? 0 : -1);
861 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
867 PointCloudSource source_transformed;
869 *input_, *source_indices_, source_transformed, transformation);
871 std::size_t nr_points = source_transformed.size();
872 std::size_t terminate_value =
873 fitness_score > 1 ? 0
874 :
static_cast<std::size_t
>((1.f - fitness_score) * nr_points);
876 float inlier_score_temp = 0;
878 std::vector<float> dists_sqr;
879 auto it = source_transformed.begin();
881 for (std::size_t i = 0; i < nr_points; it++, i++) {
883 tree_->nearestKSearch(*it, 1, ids, dists_sqr);
884 inlier_score_temp += (dists_sqr[0] < max_inlier_dist_sqr_ ? 1 : 0);
887 if (nr_points - i + inlier_score_temp < terminate_value)
892 inlier_score_temp /=
static_cast<float>(nr_points);
893 float fitness_score_temp = 1.f - inlier_score_temp;
895 if (fitness_score_temp > fitness_score)
898 fitness_score = fitness_score_temp;
903 template <
typename Po
intSource,
typename Po
intTarget,
typename NormalT,
typename Scalar>
906 finalCompute(
const std::vector<MatchingCandidates>& candidates)
909 int nr_candidates =
static_cast<int>(candidates.size());
911 float best_score = std::numeric_limits<float>::max();
912 for (
int i = 0; i < nr_candidates; i++) {
913 const float& fitness_score = candidates[i][0].fitness_score;
914 if (fitness_score < best_score) {
915 best_score = fitness_score;
921 if (!(best_index < 0)) {
922 fitness_score_ = candidates[best_index][0].fitness_score;
923 final_transformation_ = candidates[best_index][0].transformation;
924 *correspondences_ = candidates[best_index][0].correspondences;
927 converged_ = fitness_score_ < score_threshold_;
shared_ptr< const PointCloud< PointT > > ConstPtr
typename KdTreeReciprocal::Ptr KdTreeReciprocalPtr
std::string reg_name_
The registration method name.
typename PointCloudSource::Ptr PointCloudSourcePtr
int ransac_iterations_
The number of iterations RANSAC should run for.
TransformationEstimationPtr transformation_estimation_
A TransformationEstimation object, used to calculate the 4x4 rigid transformation.
int max_iterations_
The maximum number of iterations the internal optimization should run for.
void setIndices(const IndicesPtr &indices)
Provide a pointer to the vector of indices that represents the input data.
SampleConsensusModelPlane defines a model for 3D plane segmentation.
bool computeModelCoefficients(const Indices &samples, Eigen::VectorXf &model_coefficients) const override
Check whether the given index samples can form a valid plane model, compute the model coefficients fr...
double getTimeSeconds() const
Retrieve the time in seconds spent since the last call to reset().
virtual void finalCompute(const std::vector< MatchingCandidates > &candidates)
Final computation of best match out of vector of best matches.
void setupBase(pcl::Indices &base_indices, float(&ratio)[2])
Setup the base (four coplanar points) by ordering the points and computing intersection ratios and se...
int selectBaseTriangle(pcl::Indices &base_indices)
Select randomly a triplet of points with large point-to-point distances.
virtual int bruteForceCorrespondences(int idx1, int idx2, pcl::Correspondences &pairs)
Search for corresponding point pairs given the distance between two base points.
int selectBase(pcl::Indices &base_indices, float(&ratio)[2])
Select an approximately coplanar set of four points from the source cloud.
virtual int validateTransformation(Eigen::Matrix4f &transformation, float &fitness_score)
Validate the transformation by calculating the number of inliers after transforming the source cloud.
virtual int determineBaseMatches(const pcl::Indices &base_indices, std::vector< pcl::Indices > &matches, const pcl::Correspondences &pairs_a, const pcl::Correspondences &pairs_b, const float(&ratio)[2])
Determine base matches by combining the point pair candidate and search for coinciding intersection p...
virtual void linkMatchWithBase(const pcl::Indices &base_indices, pcl::Indices &match_indices, pcl::Correspondences &correspondences)
Sets the correspondences between the base B and the match M by using the distance of each point to th...
virtual void handleMatches(const pcl::Indices &base_indices, std::vector< pcl::Indices > &matches, MatchingCandidates &candidates)
Method to handle current candidate matches.
int checkBaseMatch(const pcl::Indices &match_indices, const float(&ds)[4])
Check if outer rectangle distance of matched points fit with the base rectangle.
float segmentToSegmentDist(const pcl::Indices &base_indices, float(&ratio)[2])
Calculate intersection ratios and segment to segment distances of base diagonals.
virtual int validateMatch(const pcl::Indices &base_indices, const pcl::Indices &match_indices, const pcl::Correspondences &correspondences, Eigen::Matrix4f &transformation)
Validate the matching by computing the transformation between the source and target based on the four...
FPCSInitialAlignment()
Constructor.
void computeTransformation(PointCloudSource &output, const Eigen::Matrix4f &guess) override
Rigid transformation computation method.
virtual bool initCompute()
Internal computation initialization.
search::KdTree is a wrapper class which inherits the pcl::KdTree class for performing search function...
int nearestKSearch(const PointT &point, int k, Indices &k_indices, std::vector< float > &k_sqr_distances) const override
Search for the k-nearest neighbors for the given query point.
void setInputCloud(const PointCloudConstPtr &cloud, const IndicesConstPtr &indices=IndicesConstPtr()) override
Provide a pointer to the input dataset.
Define standard C methods to do distance calculations.
Define methods for measuring time spent in code blocks.
void getMinMax3D(const pcl::PointCloud< PointT > &cloud, PointT &min_pt, PointT &max_pt)
Get the minimum and maximum values on each of the 3 (x-y-z) dimensions in a given pointcloud.
void transformPointCloud(const pcl::PointCloud< PointT > &cloud_in, pcl::PointCloud< PointT > &cloud_out, const Eigen::Matrix< Scalar, 4, 4 > &transform, bool copy_all_fields)
Apply a rigid transform defined by a 4x4 matrix.
unsigned int compute3DCentroid(ConstCloudIterator< PointT > &cloud_iterator, Eigen::Matrix< Scalar, 4, 1 > ¢roid)
Compute the 3D (X-Y-Z) centroid of a set of points and return it as a 3D vector.
double pointToPlaneDistance(const Point &p, double a, double b, double c, double d)
Get the distance from a point to a plane (unsigned) defined by ax+by+cz+d=0.
std::vector< MatchingCandidate, Eigen::aligned_allocator< MatchingCandidate > > MatchingCandidates
void ignore(const T &...)
Utility function to eliminate unused variable warnings.
float squaredEuclideanDistance(const PointType1 &p1, const PointType2 &p2)
Calculate the squared euclidean distance between the two given points.
float getMeanPointDensity(const typename pcl::PointCloud< PointT >::ConstPtr &cloud, float max_dist, int nr_threads=1)
Compute the mean point density of a given point cloud.
std::vector< pcl::Correspondence, Eigen::aligned_allocator< pcl::Correspondence > > Correspondences
float euclideanDistance(const PointType1 &p1, const PointType2 &p2)
Calculate the euclidean distance between the two given points.
IndicesAllocator<> Indices
Type used for indices in PCL.
shared_ptr< Indices > IndicesPtr
Correspondence represents a match between two entities (e.g., points, descriptors,...
A point structure representing normal coordinates and the surface curvature estimate.