diff --git a/include/multigrid_pclomp/multigrid_ndt_omp.h b/include/multigrid_pclomp/multigrid_ndt_omp.h
index 44d8a94b..d30d3592 100644
--- a/include/multigrid_pclomp/multigrid_ndt_omp.h
+++ b/include/multigrid_pclomp/multigrid_ndt_omp.h
@@ -112,7 +112,7 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
 #endif
 
   /** \brief Constructor.
-   * Sets \ref outlier_ratio_ to 0.35, \ref step_size_ to 0.05 and \ref resolution_ to 1.0
+   * Sets \ref outlier_ratio_ to 0.35, \ref params_.step_size to 0.05 and \ref params_.resolution to 1.0
    */
   MultiGridNormalDistributionsTransform();
 
@@ -128,11 +128,11 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
   virtual ~MultiGridNormalDistributionsTransform() {}
 
   void setNumThreads(int n) {
-    num_threads_ = n;
+    params_.num_threads = n;
   }
 
   inline int getNumThreads() const {
-    return num_threads_;
+    return params_.num_threads;
   }
 
   inline void setInputSource(const PointCloudSourceConstPtr &input) {
@@ -157,7 +157,7 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
    */
   inline void addTarget(const PointCloudTargetConstPtr &cloud, const std::string &target_id) {
     pcl::Registration<PointSource, PointTarget>::setInputTarget(cloud);
-    target_cells_.setLeafSize(resolution_, resolution_, resolution_);
+    target_cells_.setLeafSize(params_.resolution, params_.resolution, params_.resolution);
     target_cells_.setInputCloudAndFilter(cloud, target_id);
   }
 
@@ -174,8 +174,8 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
    */
   inline void setResolution(float resolution) {
     // Prevents unnecessary voxel initiations
-    if(resolution_ != resolution) {
-      resolution_ = resolution;
+    if(params_.resolution != resolution) {
+      params_.resolution = resolution;
       if(input_) init();
     }
   }
@@ -184,21 +184,21 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
    * \return side length of voxels
    */
   inline float getResolution() const {
-    return (resolution_);
+    return (params_.resolution);
   }
 
   /** \brief Get the newton line search maximum step length.
    * \return maximum step length
    */
   inline double getStepSize() const {
-    return (step_size_);
+    return (params_.step_size);
   }
 
   /** \brief Set/change the newton line search maximum step length.
    * \param[in] step_size maximum step length
    */
   inline void setStepSize(double step_size) {
-    step_size_ = step_size;
+    params_.step_size = step_size;
   }
 
   /** \brief Get the point cloud outlier ratio.
@@ -268,7 +268,7 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
   double calculateNearestVoxelTransformationLikelihood(const PointCloudSource &cloud) const;
 
   inline void setRegularizationScaleFactor(float regularization_scale_factor) {
-    regularization_scale_factor_ = regularization_scale_factor;
+    params_.regularization_scale_factor = regularization_scale_factor;
   }
 
   inline void setRegularizationPose(Eigen::Matrix4f regularization_pose) {
@@ -290,24 +290,36 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
     return ndt_result;
   }
 
+  /** \brief Set the transformation epsilon (maximum allowable translation squared
+   * difference between two consecutive transformations) in order for an optimization to
+   * be considered as having converged to the final solution. \param[in] epsilon the
+   * transformation epsilon in order for an optimization to be considered as having
+   * converged to the final solution.
+   */
+  inline void setTransformationEpsilon(double epsilon) {
+    params_.trans_epsilon = epsilon;
+  }
+
+  /** \brief Get the transformation epsilon (maximum allowable translation squared
+   * difference between two consecutive transformations) as set by the user.
+   */
+  inline double getTransformationEpsilon() {
+    return (params_.trans_epsilon);
+  }
+
+  inline void setMaximumIterations(int max_iterations) {
+    params_.max_iterations = max_iterations;
+  }
+  inline int getMaxIterations() const {
+    return params_.max_iterations;
+  }
+
   void setParams(const NdtParams &ndt_params) {
-    this->setTransformationEpsilon(ndt_params.trans_epsilon);
-    this->setStepSize(ndt_params.step_size);
-    this->setResolution(ndt_params.resolution);
-    this->setMaximumIterations(ndt_params.max_iterations);
-    setRegularizationScaleFactor(ndt_params.regularization_scale_factor);
-    setNumThreads(ndt_params.num_threads);
+    params_ = ndt_params;
   }
 
   NdtParams getParams() const {
-    NdtParams ndt_params;
-    ndt_params.trans_epsilon = transformation_epsilon_;
-    ndt_params.step_size = getStepSize();
-    ndt_params.resolution = getResolution();
-    ndt_params.max_iterations = max_iterations_;
-    ndt_params.regularization_scale_factor = regularization_scale_factor_;
-    ndt_params.num_threads = num_threads_;
-    return ndt_params;
+    return params_;
   }
 
   pcl::PointCloud<PointTarget> getVoxelPCD() const {
@@ -323,11 +335,9 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
   using pcl::Registration<PointSource, PointTarget>::input_;
   using pcl::Registration<PointSource, PointTarget>::target_;
   using pcl::Registration<PointSource, PointTarget>::nr_iterations_;
-  using pcl::Registration<PointSource, PointTarget>::max_iterations_;
   using pcl::Registration<PointSource, PointTarget>::previous_transformation_;
   using pcl::Registration<PointSource, PointTarget>::final_transformation_;
   using pcl::Registration<PointSource, PointTarget>::transformation_;
-  using pcl::Registration<PointSource, PointTarget>::transformation_epsilon_;
   using pcl::Registration<PointSource, PointTarget>::converged_;
 
   using pcl::Registration<PointSource, PointTarget>::update_visualizer_;
@@ -478,12 +488,6 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
 
   // double fitness_epsilon_;
 
-  /** \brief The side length of voxels. */
-  float resolution_;
-
-  /** \brief The maximum step length. */
-  double step_size_;
-
   /** \brief The ratio of outliers of points w.r.t. a normal distribution, Equation 6.7 [Magnusson 2009]. */
   double outlier_ratio_;
 
@@ -508,16 +512,16 @@ class MultiGridNormalDistributionsTransform : public pcl::Registration<PointSour
   Eigen::Vector3d h_ang_a2_, h_ang_a3_, h_ang_b2_, h_ang_b3_, h_ang_c2_, h_ang_c3_, h_ang_d1_, h_ang_d2_, h_ang_d3_, h_ang_e1_, h_ang_e2_, h_ang_e3_, h_ang_f1_, h_ang_f2_, h_ang_f3_;
 
   Eigen::Matrix<float, 16, 4> h_ang;
-  int num_threads_;
 
   Eigen::Matrix<double, 6, 6> hessian_;
   std::vector<Eigen::Matrix4f, Eigen::aligned_allocator<Eigen::Matrix4f>> transformation_array_;
   double nearest_voxel_transformation_likelihood_;
 
-  float regularization_scale_factor_;
   boost::optional<Eigen::Matrix4f> regularization_pose_;
   Eigen::Vector3f regularization_pose_translation_;
 
+  NdtParams params_;
+
 public:
   EIGEN_MAKE_ALIGNED_OPERATOR_NEW
 };
diff --git a/include/multigrid_pclomp/multigrid_ndt_omp_impl.hpp b/include/multigrid_pclomp/multigrid_ndt_omp_impl.hpp
index 54716dd1..fa86fb10 100644
--- a/include/multigrid_pclomp/multigrid_ndt_omp_impl.hpp
+++ b/include/multigrid_pclomp/multigrid_ndt_omp_impl.hpp
@@ -58,8 +58,8 @@
 
 template<typename PointSource, typename PointTarget>
 pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::MultiGridNormalDistributionsTransform(const MultiGridNormalDistributionsTransform &other) : BaseRegType(other), target_cells_(other.target_cells_) {
-  resolution_ = other.resolution_;
-  step_size_ = other.step_size_;
+  params_.resolution = other.params_.resolution;
+  params_.step_size = other.params_.step_size;
   outlier_ratio_ = other.outlier_ratio_;
   gauss_d1_ = other.gauss_d1_;
   gauss_d2_ = other.gauss_d2_;
@@ -67,20 +67,20 @@ pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::MultiGr
   trans_probability_ = other.trans_probability_;
   // No need to copy j_ang and h_ang, as those matrices are re-computed on every computeDerivatives() call
 
-  num_threads_ = other.num_threads_;
+  params_.num_threads = other.params_.num_threads;
   hessian_ = other.hessian_;
   transformation_array_ = other.transformation_array_;
   nearest_voxel_transformation_likelihood_ = other.nearest_voxel_transformation_likelihood_;
 
-  regularization_scale_factor_ = other.regularization_scale_factor_;
+  params_.regularization_scale_factor = other.params_.regularization_scale_factor;
   regularization_pose_ = other.regularization_pose_;
   regularization_pose_translation_ = other.regularization_pose_translation_;
 }
 
 template<typename PointSource, typename PointTarget>
 pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::MultiGridNormalDistributionsTransform(MultiGridNormalDistributionsTransform &&other) : BaseRegType(std::move(other)), target_cells_(std::move(other.target_cells_)) {
-  resolution_ = other.resolution_;
-  step_size_ = other.step_size_;
+  params_.resolution = other.params_.resolution;
+  params_.step_size = other.params_.step_size;
   outlier_ratio_ = other.outlier_ratio_;
   gauss_d1_ = other.gauss_d1_;
   gauss_d2_ = other.gauss_d2_;
@@ -88,12 +88,12 @@ pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::MultiGr
   trans_probability_ = other.trans_probability_;
   // No need to copy j_ang and h_ang, as those matrices are re-computed on every computeDerivatives() call
 
-  num_threads_ = other.num_threads_;
+  params_.num_threads = other.params_.num_threads;
   hessian_ = other.hessian_;
   transformation_array_ = other.transformation_array_;
   nearest_voxel_transformation_likelihood_ = other.nearest_voxel_transformation_likelihood_;
 
-  regularization_scale_factor_ = other.regularization_scale_factor_;
+  params_.regularization_scale_factor = other.params_.regularization_scale_factor;
   regularization_pose_ = other.regularization_pose_;
   regularization_pose_translation_ = other.regularization_pose_translation_;
 }
@@ -104,8 +104,8 @@ pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget> &pclomp:
 
   target_cells_ = other.target_cells_;
 
-  resolution_ = other.resolution_;
-  step_size_ = other.step_size_;
+  params_.resolution = other.params_.resolution;
+  params_.step_size = other.params_.step_size;
   outlier_ratio_ = other.outlier_ratio_;
   gauss_d1_ = other.gauss_d1_;
   gauss_d2_ = other.gauss_d2_;
@@ -113,12 +113,12 @@ pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget> &pclomp:
   trans_probability_ = other.trans_probability_;
   // No need to copy j_ang and h_ang, as those matrices are re-computed on every computeDerivatives() call
 
-  num_threads_ = other.num_threads_;
+  params_.num_threads = other.params_.num_threads;
   hessian_ = other.hessian_;
   transformation_array_ = other.transformation_array_;
   nearest_voxel_transformation_likelihood_ = other.nearest_voxel_transformation_likelihood_;
 
-  regularization_scale_factor_ = other.regularization_scale_factor_;
+  params_.regularization_scale_factor = other.params_.regularization_scale_factor;
   regularization_pose_ = other.regularization_pose_;
   regularization_pose_translation_ = other.regularization_pose_translation_;
 
@@ -130,8 +130,8 @@ pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget> &pclomp:
   BaseRegType::operator=(std::move(other));
 
   target_cells_ = std::move(other.target_cells_);
-  resolution_ = other.resolution_;
-  step_size_ = other.step_size_;
+  params_.resolution = other.params_.resolution;
+  params_.step_size = other.params_.step_size;
   outlier_ratio_ = other.outlier_ratio_;
   gauss_d1_ = other.gauss_d1_;
   gauss_d2_ = other.gauss_d2_;
@@ -139,12 +139,12 @@ pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget> &pclomp:
   trans_probability_ = other.trans_probability_;
   // No need to copy j_ang and h_ang, as those matrices are re-computed on every computeDerivatives() call
 
-  num_threads_ = other.num_threads_;
+  params_.num_threads = other.params_.num_threads;
   hessian_ = other.hessian_;
   transformation_array_ = other.transformation_array_;
   nearest_voxel_transformation_likelihood_ = other.nearest_voxel_transformation_likelihood_;
 
-  regularization_scale_factor_ = other.regularization_scale_factor_;
+  params_.regularization_scale_factor = other.params_.regularization_scale_factor;
   regularization_pose_ = other.regularization_pose_;
   regularization_pose_translation_ = other.regularization_pose_translation_;
 
@@ -154,22 +154,25 @@ pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget> &pclomp:
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 template<typename PointSource, typename PointTarget>
 pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::MultiGridNormalDistributionsTransform()
-    : target_cells_(), resolution_(1.0f), step_size_(0.1), outlier_ratio_(0.55), gauss_d1_(), gauss_d2_(), gauss_d3_(), trans_probability_(), regularization_pose_(boost::none) {
+    : target_cells_(), outlier_ratio_(0.55), gauss_d1_(), gauss_d2_(), gauss_d3_(), trans_probability_(), regularization_pose_(boost::none) {
   reg_name_ = "MultiGridNormalDistributionsTransform";
 
+  params_.trans_epsilon = 0.1;
+  params_.step_size = 0.1;
+  params_.resolution = 1.0f;
+  params_.max_iterations = 35;
+  params_.search_method = DIRECT7;
+  params_.num_threads = omp_get_max_threads();
+  params_.regularization_scale_factor = 0.0f;
+
   double gauss_c1, gauss_c2;
 
   // Initializes the gaussian fitting parameters (eq. 6.8) [Magnusson 2009]
   gauss_c1 = 10.0 * (1 - outlier_ratio_);
-  gauss_c2 = outlier_ratio_ / pow(resolution_, 3);
+  gauss_c2 = outlier_ratio_ / pow(params_.resolution, 3);
   gauss_d3_ = -log(gauss_c2);
   gauss_d1_ = -log(gauss_c1 + gauss_c2) - gauss_d3_;
   gauss_d2_ = -2 * log((-log(gauss_c1 * exp(-0.5) + gauss_c2) - gauss_d3_) / gauss_d1_);
-
-  transformation_epsilon_ = 0.1;
-  max_iterations_ = 35;
-
-  num_threads_ = omp_get_max_threads();
 }
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -182,7 +185,7 @@ void pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::co
 
   // Initializes the gaussian fitting parameters (eq. 6.8) [Magnusson 2009]
   gauss_c1 = 10 * (1 - outlier_ratio_);
-  gauss_c2 = outlier_ratio_ / pow(resolution_, 3);
+  gauss_c2 = outlier_ratio_ / pow(params_.resolution, 3);
   gauss_d3_ = -log(gauss_c2);
   gauss_d1_ = -log(gauss_c1 + gauss_c2) - gauss_d3_;
   gauss_d2_ = -2 * log((-log(gauss_c1 * exp(-0.5) + gauss_c2) - gauss_d3_) / gauss_d1_);
@@ -243,7 +246,7 @@ void pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::co
     }
 
     delta_p.normalize();
-    delta_p_norm = computeStepLengthMT(p, delta_p, delta_p_norm, step_size_, transformation_epsilon_ / 2, score, score_gradient, hessian, output);
+    delta_p_norm = computeStepLengthMT(p, delta_p, delta_p_norm, params_.step_size, params_.trans_epsilon / 2, score, score_gradient, hessian, output);
     delta_p *= delta_p_norm;
 
     transformation_ = (Eigen::Translation<float, 3>(static_cast<float>(delta_p(0)), static_cast<float>(delta_p(1)), static_cast<float>(delta_p(2))) * Eigen::AngleAxis<float>(static_cast<float>(delta_p(3)), Eigen::Vector3f::UnitX()) * Eigen::AngleAxis<float>(static_cast<float>(delta_p(4)), Eigen::Vector3f::UnitY()) * Eigen::AngleAxis<float>(static_cast<float>(delta_p(5)), Eigen::Vector3f::UnitZ()))
@@ -258,7 +261,7 @@ void pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::co
 
     nr_iterations_++;
 
-    if(nr_iterations_ >= max_iterations_ || (nr_iterations_ && (std::fabs(delta_p_norm) < transformation_epsilon_))) {
+    if(nr_iterations_ >= params_.max_iterations || (nr_iterations_ && (std::fabs(delta_p_norm) < params_.trans_epsilon))) {
       converged_ = true;
     }
   }
@@ -312,11 +315,11 @@ double pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::
   // Precompute Angular Derivatives (eq. 6.19 and 6.21)[Magnusson 2009]
   computeAngleDerivatives(p);
 
-  std::vector<std::vector<TargetGridLeafConstPtr>> neighborhoods(num_threads_);
-  std::vector<std::vector<float>> distancess(num_threads_);
+  std::vector<std::vector<TargetGridLeafConstPtr>> neighborhoods(params_.num_threads);
+  std::vector<std::vector<float>> distancess(params_.num_threads);
 
   // Update gradient and hessian for each point, line 17 in Algorithm 2 [Magnusson 2009]
-#pragma omp parallel for num_threads(num_threads_) schedule(guided, 8)
+#pragma omp parallel for num_threads(params_.num_threads) schedule(guided, 8)
   for(size_t idx = 0; idx < input_size; ++idx) {
     int thread_n = omp_get_thread_num();
 
@@ -342,7 +345,7 @@ double pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::
     auto &distances = distancess[thread_n];
 
     // Neighborhood search method other than kdtree is disabled in multigrid_ndt_omp
-    target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+    target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
 
     double sum_score_pt = 0;
     double nearest_voxel_score_pt = 0;
@@ -405,14 +408,14 @@ double pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::
     const float longitudinal_distance = dy * sin_yaw + dx * cos_yaw;
     const auto neighborhood_count_weight = static_cast<float>(total_neighborhood_count);
 
-    regularization_score = -regularization_scale_factor_ * neighborhood_count_weight * longitudinal_distance * longitudinal_distance;
+    regularization_score = -params_.regularization_scale_factor * neighborhood_count_weight * longitudinal_distance * longitudinal_distance;
 
-    regularization_gradient(0, 0) = regularization_scale_factor_ * neighborhood_count_weight * 2.0f * cos_yaw * longitudinal_distance;
-    regularization_gradient(1, 0) = regularization_scale_factor_ * neighborhood_count_weight * 2.0f * sin_yaw * longitudinal_distance;
+    regularization_gradient(0, 0) = params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * cos_yaw * longitudinal_distance;
+    regularization_gradient(1, 0) = params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * sin_yaw * longitudinal_distance;
 
-    regularization_hessian(0, 0) = -regularization_scale_factor_ * neighborhood_count_weight * 2.0f * cos_yaw * cos_yaw;
-    regularization_hessian(0, 1) = -regularization_scale_factor_ * neighborhood_count_weight * 2.0f * cos_yaw * sin_yaw;
-    regularization_hessian(1, 1) = -regularization_scale_factor_ * neighborhood_count_weight * 2.0f * sin_yaw * sin_yaw;
+    regularization_hessian(0, 0) = -params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * cos_yaw * cos_yaw;
+    regularization_hessian(0, 1) = -params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * cos_yaw * sin_yaw;
+    regularization_hessian(1, 1) = -params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * sin_yaw * sin_yaw;
     regularization_hessian(1, 0) = regularization_hessian(0, 1);
 
     score += regularization_score;
@@ -658,7 +661,7 @@ void pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::co
     std::vector<float> distances;
 
     // Neighborhood search method other than kdtree is disabled in multigrid_ndt_omp
-    target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+    target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
 
     for(auto &cell : neighborhood) {
       x_pt = (*input_)[idx];
@@ -985,7 +988,7 @@ double pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::
     std::vector<float> distances;
 
     // Neighborhood search method other than kdtree is disabled in multigrid_ndt_omp
-    target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+    target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
 
     for(auto &cell : neighborhood) {
       Eigen::Vector3d x_trans = Eigen::Vector3d(x_trans_pt.x, x_trans_pt.y, x_trans_pt.z);
@@ -1023,7 +1026,7 @@ double pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::
     std::vector<float> distances;
 
     // Neighborhood search method other than kdtree is disabled in multigrid_ndt_omp
-    target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+    target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
 
     for(auto &cell : neighborhood) {
       Eigen::Vector3d x_trans = Eigen::Vector3d(x_trans_pt.x, x_trans_pt.y, x_trans_pt.z);
@@ -1063,7 +1066,7 @@ double pclomp::MultiGridNormalDistributionsTransform<PointSource, PointTarget>::
     std::vector<float> distances;
 
     // Neighborhood search method other than kdtree is disabled in multigrid_ndt_omp
-    target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+    target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
 
     for(auto &cell : neighborhood) {
       Eigen::Vector3d x_trans = Eigen::Vector3d(x_trans_pt.x, x_trans_pt.y, x_trans_pt.z);
diff --git a/include/pclomp/ndt_omp.h b/include/pclomp/ndt_omp.h
index c337b83f..d29b49eb 100644
--- a/include/pclomp/ndt_omp.h
+++ b/include/pclomp/ndt_omp.h
@@ -96,7 +96,7 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
 #endif
 
   /** \brief Constructor.
-   * Sets \ref outlier_ratio_ to 0.35, \ref step_size_ to 0.05 and \ref resolution_ to 1.0
+   * Sets \ref outlier_ratio_ to 0.35, \ref step_size_ to 0.05 and \ref params_.resolution to 1.0
    */
   NormalDistributionsTransform();
 
@@ -104,11 +104,11 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
   virtual ~NormalDistributionsTransform() {}
 
   void setNumThreads(int n) {
-    num_threads_ = n;
+    params_.num_threads = n;
   }
 
   inline int getNumThreads() const {
-    return num_threads_;
+    return params_.num_threads;
   }
 
   /** \brief Provide a pointer to the input target (e.g., the point cloud that we want to align the input source to).
@@ -124,8 +124,8 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
    */
   inline void setResolution(float resolution) {
     // Prevents unnecessary voxel initiations
-    if(resolution_ != resolution) {
-      resolution_ = resolution;
+    if(params_.resolution != resolution) {
+      params_.resolution = resolution;
       if(input_) init();
     }
   }
@@ -134,21 +134,21 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
    * \return side length of voxels
    */
   inline float getResolution() const {
-    return (resolution_);
+    return (params_.resolution);
   }
 
   /** \brief Get the newton line search maximum step length.
    * \return maximum step length
    */
   inline double getStepSize() const {
-    return (step_size_);
+    return (params_.step_size);
   }
 
   /** \brief Set/change the newton line search maximum step length.
    * \param[in] step_size maximum step length
    */
   inline void setStepSize(double step_size) {
-    step_size_ = step_size;
+    params_.step_size = step_size;
   }
 
   /** \brief Get the point cloud outlier ratio.
@@ -166,11 +166,11 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
   }
 
   inline void setNeighborhoodSearchMethod(NeighborSearchMethod method) {
-    search_method = method;
+    params_.search_method = method;
   }
 
   inline NeighborSearchMethod getNeighborhoodSearchMethod() const {
-    return search_method;
+    return params_.search_method;
   }
 
   /** \brief Get the registration alignment probability.
@@ -226,7 +226,7 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
   double calculateNearestVoxelTransformationLikelihood(const PointCloudSource &cloud) const;
 
   inline void setRegularizationScaleFactor(float regularization_scale_factor) {
-    regularization_scale_factor_ = regularization_scale_factor;
+    params_.regularization_scale_factor = regularization_scale_factor;
   }
 
   inline void setRegularizationPose(Eigen::Matrix4f regularization_pose) {
@@ -248,26 +248,36 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
     return ndt_result;
   }
 
+  /** \brief Set the transformation epsilon (maximum allowable translation squared
+   * difference between two consecutive transformations) in order for an optimization to
+   * be considered as having converged to the final solution. \param[in] epsilon the
+   * transformation epsilon in order for an optimization to be considered as having
+   * converged to the final solution.
+   */
+  inline void setTransformationEpsilon(double epsilon) {
+    params_.trans_epsilon = epsilon;
+  }
+
+  /** \brief Get the transformation epsilon (maximum allowable translation squared
+   * difference between two consecutive transformations) as set by the user.
+   */
+  inline double getTransformationEpsilon() {
+    return (params_.trans_epsilon);
+  }
+
+  inline void setMaximumIterations(int max_iterations) {
+    params_.max_iterations = max_iterations;
+  }
+  inline int getMaxIterations() const {
+    return params_.max_iterations;
+  }
+
   void setParams(const NdtParams &ndt_params) {
-    this->setTransformationEpsilon(ndt_params.trans_epsilon);
-    this->setStepSize(ndt_params.step_size);
-    this->setResolution(ndt_params.resolution);
-    this->setMaximumIterations(ndt_params.max_iterations);
-    setRegularizationScaleFactor(ndt_params.regularization_scale_factor);
-    setNeighborhoodSearchMethod(ndt_params.search_method);
-    setNumThreads(ndt_params.num_threads);
+    params_ = ndt_params;
   }
 
   NdtParams getParams() const {
-    NdtParams ndt_params;
-    ndt_params.trans_epsilon = transformation_epsilon_;
-    ndt_params.step_size = getStepSize();
-    ndt_params.resolution = getResolution();
-    ndt_params.max_iterations = max_iterations_;
-    ndt_params.regularization_scale_factor = regularization_scale_factor_;
-    ndt_params.search_method = getNeighborhoodSearchMethod();
-    ndt_params.num_threads = num_threads_;
-    return ndt_params;
+    return params_;
   }
 
 protected:
@@ -277,11 +287,9 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
   using pcl::Registration<PointSource, PointTarget>::indices_;
   using pcl::Registration<PointSource, PointTarget>::target_;
   using pcl::Registration<PointSource, PointTarget>::nr_iterations_;
-  using pcl::Registration<PointSource, PointTarget>::max_iterations_;
   using pcl::Registration<PointSource, PointTarget>::previous_transformation_;
   using pcl::Registration<PointSource, PointTarget>::final_transformation_;
   using pcl::Registration<PointSource, PointTarget>::transformation_;
-  using pcl::Registration<PointSource, PointTarget>::transformation_epsilon_;
   using pcl::Registration<PointSource, PointTarget>::converged_;
   using pcl::Registration<PointSource, PointTarget>::corr_dist_threshold_;
   using pcl::Registration<PointSource, PointTarget>::inlier_threshold_;
@@ -303,7 +311,7 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
 
   /** \brief Initiate covariance voxel structure. */
   void inline init() {
-    target_cells_.setLeafSize(resolution_, resolution_, resolution_);
+    target_cells_.setLeafSize(params_.resolution, params_.resolution, params_.resolution);
     target_cells_.setInputCloud(target_);
     // Initiate voxel structure.
     target_cells_.filter(true);
@@ -439,12 +447,6 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
 
   // double fitness_epsilon_;
 
-  /** \brief The side length of voxels. */
-  float resolution_;
-
-  /** \brief The maximum step length. */
-  double step_size_;
-
   /** \brief The ratio of outliers of points w.r.t. a normal distribution, Equation 6.7 [Magnusson 2009]. */
   double outlier_ratio_;
 
@@ -476,19 +478,16 @@ class NormalDistributionsTransform : public pcl::Registration<PointSource, Point
   /** \brief The second order derivative of the transformation of a point w.r.t. the transform vector, \f$ H_E \f$ in Equation 6.20 [Magnusson 2009]. */
   //      Eigen::Matrix<double, 18, 6> point_hessian_;
 
-  int num_threads_;
-
   Eigen::Matrix<double, 6, 6> hessian_;
   std::vector<Eigen::Matrix4f, Eigen::aligned_allocator<Eigen::Matrix4f>> transformation_array_;
   double nearest_voxel_transformation_likelihood_;
 
-  float regularization_scale_factor_;
   boost::optional<Eigen::Matrix4f> regularization_pose_;
   Eigen::Vector3f regularization_pose_translation_;
 
-public:
-  NeighborSearchMethod search_method;
+  NdtParams params_;
 
+public:
   EIGEN_MAKE_ALIGNED_OPERATOR_NEW
 };
 
diff --git a/include/pclomp/ndt_omp_impl.hpp b/include/pclomp/ndt_omp_impl.hpp
index 9a7e6b3b..b8286014 100644
--- a/include/pclomp/ndt_omp_impl.hpp
+++ b/include/pclomp/ndt_omp_impl.hpp
@@ -45,23 +45,25 @@
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 template<typename PointSource, typename PointTarget>
 pclomp::NormalDistributionsTransform<PointSource, PointTarget>::NormalDistributionsTransform()
-    : target_cells_(), resolution_(1.0f), step_size_(0.1), outlier_ratio_(0.55), gauss_d1_(), gauss_d2_(), gauss_d3_(), trans_probability_(), regularization_pose_(boost::none), j_ang_a_(), j_ang_b_(), j_ang_c_(), j_ang_d_(), j_ang_e_(), j_ang_f_(), j_ang_g_(), j_ang_h_(), h_ang_a2_(), h_ang_a3_(), h_ang_b2_(), h_ang_b3_(), h_ang_c2_(), h_ang_c3_(), h_ang_d1_(), h_ang_d2_(), h_ang_d3_(), h_ang_e1_(), h_ang_e2_(), h_ang_e3_(), h_ang_f1_(), h_ang_f2_(), h_ang_f3_() {
+    : target_cells_(), outlier_ratio_(0.55), gauss_d1_(), gauss_d2_(), gauss_d3_(), trans_probability_(), regularization_pose_(boost::none), j_ang_a_(), j_ang_b_(), j_ang_c_(), j_ang_d_(), j_ang_e_(), j_ang_f_(), j_ang_g_(), j_ang_h_(), h_ang_a2_(), h_ang_a3_(), h_ang_b2_(), h_ang_b3_(), h_ang_c2_(), h_ang_c3_(), h_ang_d1_(), h_ang_d2_(), h_ang_d3_(), h_ang_e1_(), h_ang_e2_(), h_ang_e3_(), h_ang_f1_(), h_ang_f2_(), h_ang_f3_() {
   reg_name_ = "NormalDistributionsTransform";
 
+  params_.trans_epsilon = 0.1;
+  params_.step_size = 0.1;
+  params_.resolution = 1.0f;
+  params_.max_iterations = 35;
+  params_.search_method = DIRECT7;
+  params_.num_threads = omp_get_max_threads();
+  params_.regularization_scale_factor = 0.0f;
+
   double gauss_c1, gauss_c2;
 
   // Initializes the gaussian fitting parameters (eq. 6.8) [Magnusson 2009]
   gauss_c1 = 10.0 * (1 - outlier_ratio_);
-  gauss_c2 = outlier_ratio_ / pow(resolution_, 3);
+  gauss_c2 = outlier_ratio_ / pow(params_.resolution, 3);
   gauss_d3_ = -log(gauss_c2);
   gauss_d1_ = -log(gauss_c1 + gauss_c2) - gauss_d3_;
   gauss_d2_ = -2 * log((-log(gauss_c1 * exp(-0.5) + gauss_c2) - gauss_d3_) / gauss_d1_);
-
-  transformation_epsilon_ = 0.1;
-  max_iterations_ = 35;
-
-  search_method = DIRECT7;
-  num_threads_ = omp_get_max_threads();
 }
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -74,7 +76,7 @@ void pclomp::NormalDistributionsTransform<PointSource, PointTarget>::computeTran
 
   // Initializes the gaussian fitting parameters (eq. 6.8) [Magnusson 2009]
   gauss_c1 = 10 * (1 - outlier_ratio_);
-  gauss_c2 = outlier_ratio_ / pow(resolution_, 3);
+  gauss_c2 = outlier_ratio_ / pow(params_.resolution, 3);
   gauss_d3_ = -log(gauss_c2);
   gauss_d1_ = -log(gauss_c1 + gauss_c2) - gauss_d3_;
   gauss_d2_ = -2 * log((-log(gauss_c1 * exp(-0.5) + gauss_c2) - gauss_d3_) / gauss_d1_);
@@ -135,7 +137,7 @@ void pclomp::NormalDistributionsTransform<PointSource, PointTarget>::computeTran
     }
 
     delta_p.normalize();
-    delta_p_norm = computeStepLengthMT(p, delta_p, delta_p_norm, step_size_, transformation_epsilon_ / 2, score, score_gradient, hessian, output);
+    delta_p_norm = computeStepLengthMT(p, delta_p, delta_p_norm, params_.step_size, params_.trans_epsilon / 2, score, score_gradient, hessian, output);
     delta_p *= delta_p_norm;
 
     transformation_ = (Eigen::Translation<float, 3>(static_cast<float>(delta_p(0)), static_cast<float>(delta_p(1)), static_cast<float>(delta_p(2))) * Eigen::AngleAxis<float>(static_cast<float>(delta_p(3)), Eigen::Vector3f::UnitX()) * Eigen::AngleAxis<float>(static_cast<float>(delta_p(4)), Eigen::Vector3f::UnitY()) * Eigen::AngleAxis<float>(static_cast<float>(delta_p(5)), Eigen::Vector3f::UnitZ()))
@@ -148,7 +150,7 @@ void pclomp::NormalDistributionsTransform<PointSource, PointTarget>::computeTran
     // Update Visualizer (untested)
     if(update_visualizer_ != 0) update_visualizer_(output, std::vector<int>(), *target_, std::vector<int>());
 
-    if(nr_iterations_ > max_iterations_ || (nr_iterations_ && (std::fabs(delta_p_norm) < transformation_epsilon_))) {
+    if(nr_iterations_ > params_.max_iterations || (nr_iterations_ && (std::fabs(delta_p_norm) < params_.trans_epsilon))) {
       converged_ = true;
     }
 
@@ -203,11 +205,11 @@ double pclomp::NormalDistributionsTransform<PointSource, PointTarget>::computeDe
   // Precompute Angular Derivatives (eq. 6.19 and 6.21)[Magnusson 2009]
   computeAngleDerivatives(p);
 
-  std::vector<std::vector<TargetGridLeafConstPtr>> neighborhoods(num_threads_);
-  std::vector<std::vector<float>> distancess(num_threads_);
+  std::vector<std::vector<TargetGridLeafConstPtr>> neighborhoods(params_.num_threads);
+  std::vector<std::vector<float>> distancess(params_.num_threads);
 
   // Update gradient and hessian for each point, line 17 in Algorithm 2 [Magnusson 2009]
-#pragma omp parallel for num_threads(num_threads_) schedule(guided, 8)
+#pragma omp parallel for num_threads(params_.num_threads) schedule(guided, 8)
   for(std::size_t idx = 0; idx < input_->points.size(); idx++) {
     int thread_n = omp_get_thread_num();
 
@@ -233,9 +235,9 @@ double pclomp::NormalDistributionsTransform<PointSource, PointTarget>::computeDe
     auto &distances = distancess[thread_n];
 
     // Find neighbors (Radius search has been experimentally faster than direct neighbor checking.
-    switch(search_method) {
+    switch(params_.search_method) {
       case KDTREE:
-        target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+        target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
         break;
       case DIRECT26:
         target_cells_.getNeighborhoodAtPoint(x_trans_pt, neighborhood);
@@ -311,14 +313,14 @@ double pclomp::NormalDistributionsTransform<PointSource, PointTarget>::computeDe
     const float longitudinal_distance = dy * sin_yaw + dx * cos_yaw;
     const auto neighborhood_count_weight = static_cast<float>(total_neighborhood_count);
 
-    regularization_score = -regularization_scale_factor_ * neighborhood_count_weight * longitudinal_distance * longitudinal_distance;
+    regularization_score = -params_.regularization_scale_factor * neighborhood_count_weight * longitudinal_distance * longitudinal_distance;
 
-    regularization_gradient(0, 0) = regularization_scale_factor_ * neighborhood_count_weight * 2.0f * cos_yaw * longitudinal_distance;
-    regularization_gradient(1, 0) = regularization_scale_factor_ * neighborhood_count_weight * 2.0f * sin_yaw * longitudinal_distance;
+    regularization_gradient(0, 0) = params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * cos_yaw * longitudinal_distance;
+    regularization_gradient(1, 0) = params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * sin_yaw * longitudinal_distance;
 
-    regularization_hessian(0, 0) = -regularization_scale_factor_ * neighborhood_count_weight * 2.0f * cos_yaw * cos_yaw;
-    regularization_hessian(0, 1) = -regularization_scale_factor_ * neighborhood_count_weight * 2.0f * cos_yaw * sin_yaw;
-    regularization_hessian(1, 1) = -regularization_scale_factor_ * neighborhood_count_weight * 2.0f * sin_yaw * sin_yaw;
+    regularization_hessian(0, 0) = -params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * cos_yaw * cos_yaw;
+    regularization_hessian(0, 1) = -params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * cos_yaw * sin_yaw;
+    regularization_hessian(1, 1) = -params_.regularization_scale_factor * neighborhood_count_weight * 2.0f * sin_yaw * sin_yaw;
     regularization_hessian(1, 0) = regularization_hessian(0, 1);
 
     score += regularization_score;
@@ -592,9 +594,9 @@ void pclomp::NormalDistributionsTransform<PointSource, PointTarget>::computeHess
     // Find neighbors (Radius search has been experimentally faster than direct neighbor checking.
     std::vector<TargetGridLeafConstPtr> neighborhood;
     std::vector<float> distances;
-    switch(search_method) {
+    switch(params_.search_method) {
       case KDTREE:
-        target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+        target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
         break;
       case DIRECT26:
         target_cells_.getNeighborhoodAtPoint(x_trans_pt, neighborhood);
@@ -935,9 +937,9 @@ double pclomp::NormalDistributionsTransform<PointSource, PointTarget>::calculate
     // Find neighbors (Radius search has been experimentally faster than direct neighbor checking.
     std::vector<TargetGridLeafConstPtr> neighborhood;
     std::vector<float> distances;
-    switch(search_method) {
+    switch(params_.search_method) {
       case KDTREE:
-        target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+        target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
         break;
       case DIRECT26:
         target_cells_.getNeighborhoodAtPoint(x_trans_pt, neighborhood);
@@ -987,9 +989,9 @@ double pclomp::NormalDistributionsTransform<PointSource, PointTarget>::calculate
     // Find neighbors (Radius search has been experimentally faster than direct neighbor checking.
     std::vector<TargetGridLeafConstPtr> neighborhood;
     std::vector<float> distances;
-    switch(search_method) {
+    switch(params_.search_method) {
       case KDTREE:
-        target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+        target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
         break;
       case DIRECT26:
         target_cells_.getNeighborhoodAtPoint(x_trans_pt, neighborhood);
@@ -1041,9 +1043,9 @@ double pclomp::NormalDistributionsTransform<PointSource, PointTarget>::calculate
     // Find neighbors (Radius search has been experimentally faster than direct neighbor checking.
     std::vector<TargetGridLeafConstPtr> neighborhood;
     std::vector<float> distances;
-    switch(search_method) {
+    switch(params_.search_method) {
       case KDTREE:
-        target_cells_.radiusSearch(x_trans_pt, resolution_, neighborhood, distances);
+        target_cells_.radiusSearch(x_trans_pt, params_.resolution, neighborhood, distances);
         break;
       case DIRECT26:
         target_cells_.getNeighborhoodAtPoint(x_trans_pt, neighborhood);