Added distortion correction post processor with associated test

src/base/DICe_FieldEnums.cpp

@@ -108,6 +108,10 @@ void create_string_maps()
   field_name_string[field_enums::PROJECTION_AUG_Y]                                   = "PROJECTION_AUG_Y";
   field_name_string[field_enums::SUBSET_DISPLACEMENT_X]                              = "DISPLACEMENT_X"; // the string is different from the enum because of legacy field names
   field_name_string[field_enums::SUBSET_DISPLACEMENT_Y]                              = "DISPLACEMENT_Y"; // the string is different from the enum because of legacy field names
+  field_name_string[field_enums::DIST_CORRECTED_SUBSET_DISP_X]                       = "DIST_CORRECTED_SUBSET_DISP_X";
+  field_name_string[field_enums::DIST_CORRECTED_SUBSET_DISP_Y]                       = "DIST_CORRECTED_SUBSET_DISP_Y";
+  field_name_string[field_enums::DIST_CORRECTED_SUBSET_X]                            = "DIST_CORRECTED_SUBSET_X";
+  field_name_string[field_enums::DIST_CORRECTED_SUBSET_Y]                            = "DIST_CORRECTED_SUBSET_Y";
   field_name_string[field_enums::STEREO_DISPLACEMENT_X]                              = "STEREO_DISPLACEMENT_X";
   field_name_string[field_enums::STEREO_DISPLACEMENT_Y]                              = "STEREO_DISPLACEMENT_Y";
   field_name_string[field_enums::MODEL_DISPLACEMENT_X]                               = "MODEL_DISPLACEMENT_X";

src/base/DICe_FieldEnums.h

@@ -94,6 +94,10 @@ Field_Name
   PROJECTION_AUG_Y,
   SUBSET_DISPLACEMENT_X,
   SUBSET_DISPLACEMENT_Y,
+  DIST_CORRECTED_SUBSET_X,
+  DIST_CORRECTED_SUBSET_Y,
+  DIST_CORRECTED_SUBSET_DISP_X,
+  DIST_CORRECTED_SUBSET_DISP_Y,
   STEREO_DISPLACEMENT_X,
   STEREO_DISPLACEMENT_Y,
   MODEL_DISPLACEMENT_X,
@@ -476,6 +480,14 @@ const Field_Spec SUBSET_DISPLACEMENT_Y_FS(field_enums::SCALAR_FIELD_TYPE,field_e
 /// field spec
 const Field_Spec SUBSET_DISPLACEMENT_Y_NM1_FS(field_enums::SCALAR_FIELD_TYPE,field_enums::SUBSET_DISPLACEMENT_Y,field_enums::NODE_RANK,field_enums::STATE_N_MINUS_ONE,false,true);
 /// field spec
+const Field_Spec DIST_CORRECTED_SUBSET_DISP_X_FS(field_enums::SCALAR_FIELD_TYPE,field_enums::DIST_CORRECTED_SUBSET_DISP_X,field_enums::NODE_RANK,field_enums::NO_FIELD_STATE,true,true);
+/// field spec
+const Field_Spec DIST_CORRECTED_SUBSET_DISP_Y_FS(field_enums::SCALAR_FIELD_TYPE,field_enums::DIST_CORRECTED_SUBSET_DISP_Y,field_enums::NODE_RANK,field_enums::NO_FIELD_STATE,true,true);
+/// field spec
+const Field_Spec DIST_CORRECTED_SUBSET_X_FS(field_enums::SCALAR_FIELD_TYPE,field_enums::DIST_CORRECTED_SUBSET_X,field_enums::NODE_RANK,field_enums::NO_FIELD_STATE,true,true);
+/// field spec
+const Field_Spec DIST_CORRECTED_SUBSET_Y_FS(field_enums::SCALAR_FIELD_TYPE,field_enums::DIST_CORRECTED_SUBSET_Y,field_enums::NODE_RANK,field_enums::NO_FIELD_STATE,true,true);
+/// field spec
 const Field_Spec STEREO_SUBSET_DISPLACEMENT_X_FS(field_enums::SCALAR_FIELD_TYPE,field_enums::STEREO_DISPLACEMENT_X,field_enums::NODE_RANK,field_enums::NO_FIELD_STATE,true,true);
 /// field spec
 const Field_Spec STEREO_SUBSET_DISPLACEMENT_Y_FS(field_enums::SCALAR_FIELD_TYPE,field_enums::STEREO_DISPLACEMENT_Y,field_enums::NODE_RANK,field_enums::NO_FIELD_STATE,true,true);
@@ -709,7 +721,7 @@ const Field_Spec ROT_TRANS_3D_TRANS_Y_FS(field_enums::SCALAR_FIELD_TYPE, field_e
 const Field_Spec ROT_TRANS_3D_TRANS_Z_FS(field_enums::SCALAR_FIELD_TYPE, field_enums::ROT_TRANS_3D_TRANS_Z, field_enums::NODE_RANK, field_enums::NO_FIELD_STATE, true);
 
 /// the number of fields that have been defined (must be set at compile time)
-const int_t num_fields_defined = 142;
+const int_t num_fields_defined = 146;
 
 /// array of all the valid field specs
 const field_enums::Field_Spec fs_spec_vec[num_fields_defined] = {
@@ -736,6 +748,10 @@ const field_enums::Field_Spec fs_spec_vec[num_fields_defined] = {
     SUBSET_DISPLACEMENT_X_NM1_FS,
     SUBSET_DISPLACEMENT_Y_FS,
     SUBSET_DISPLACEMENT_Y_NM1_FS,
+    DIST_CORRECTED_SUBSET_DISP_X_FS,
+    DIST_CORRECTED_SUBSET_DISP_Y_FS,
+    DIST_CORRECTED_SUBSET_X_FS,
+    DIST_CORRECTED_SUBSET_Y_FS,
     STEREO_SUBSET_DISPLACEMENT_X_FS,
     STEREO_SUBSET_DISPLACEMENT_Y_FS,
     MODEL_DISPLACEMENT_X_FS,

src/core/DICe_PostProcessor.cpp

@@ -822,6 +822,57 @@ Plotly_Contour_Post_Processor::execute(Teuchos::RCP<Image> ref_img, Teuchos::RCP
   DEBUG_MSG("Plotly_Contour_Post_Processor execute() end");
 }
 
+Distortion_Correction_Post_Processor::Distortion_Correction_Post_Processor(const Teuchos::RCP<Teuchos::ParameterList> & params) :
+  Post_Processor(post_process_distortion_correction){
+  field_specs_.push_back(DICe::field_enums::DIST_CORRECTED_SUBSET_X_FS);
+  field_specs_.push_back(DICe::field_enums::DIST_CORRECTED_SUBSET_Y_FS);
+  field_specs_.push_back(DICe::field_enums::DIST_CORRECTED_SUBSET_DISP_X_FS);
+  field_specs_.push_back(DICe::field_enums::DIST_CORRECTED_SUBSET_DISP_Y_FS);
+  DEBUG_MSG("Enabling post processor Distortion_Correction_Post_Processor with associated fields:");
+  for(size_t i=0;i<field_specs_.size();++i){
+    DEBUG_MSG(field_specs_[i].get_name_label());
+  }
+  std::string camera_sys_file = params->get<std::string>(camera_system_file);
+  DEBUG_MSG("Distortion_Correction_Post_Processor(): camera system file: " << camera_sys_file);
+  tri_ = Teuchos::rcp(new DICe::Triangulation(camera_sys_file));
+}
+
+void
+Distortion_Correction_Post_Processor::execute(Teuchos::RCP<Image> ref_img, Teuchos::RCP<Image> def_img){
+  DEBUG_MSG("Distortion_Correction_Post_Processor execute() begin");
+  Teuchos::RCP<DICe::MultiField> subset_x = mesh_->get_field(DICe::field_enums::SUBSET_COORDINATES_X_FS);
+  Teuchos::RCP<DICe::MultiField> subset_y = mesh_->get_field(DICe::field_enums::SUBSET_COORDINATES_Y_FS);
+  Teuchos::RCP<DICe::MultiField> subset_u = mesh_->get_field(DICe::field_enums::SUBSET_DISPLACEMENT_X_FS);
+  Teuchos::RCP<DICe::MultiField> subset_v = mesh_->get_field(DICe::field_enums::SUBSET_DISPLACEMENT_Y_FS);
+  Teuchos::RCP<DICe::MultiField> subset_x_prime = mesh_->get_field(DICe::field_enums::DIST_CORRECTED_SUBSET_X_FS);
+  Teuchos::RCP<DICe::MultiField> subset_y_prime = mesh_->get_field(DICe::field_enums::DIST_CORRECTED_SUBSET_Y_FS);
+  Teuchos::RCP<DICe::MultiField> subset_u_prime = mesh_->get_field(DICe::field_enums::DIST_CORRECTED_SUBSET_DISP_X_FS);
+  Teuchos::RCP<DICe::MultiField> subset_v_prime = mesh_->get_field(DICe::field_enums::DIST_CORRECTED_SUBSET_DISP_Y_FS);
+
+  std::vector<cv::Point2f> undist_points;
+  std::vector<cv::Point2f> undist_def_points;
+  const int_t num_subsets = mesh_->get_scalar_node_overlap_map()->get_num_local_elements();
+  std::vector<cv::Point2f> dist_points(num_subsets);
+  std::vector<cv::Point2f> dist_def_points(num_subsets);
+  for(int_t i=0;i<num_subsets;++i){
+    dist_points[i] = cv::Point2f(subset_x->local_value(i),subset_y->local_value(i));
+    dist_def_points[i] = cv::Point2f(subset_x->local_value(i)+subset_u->local_value(i),subset_y->local_value(i)+subset_v->local_value(i));
+  }
+  cv::Mat M = tri_->camera_matrix(0);
+  cv::Mat D = tri_->distortion_matrix(0);
+  cv::Mat M_prime = cv::Mat::zeros(3,4,CV_64F);
+  M.copyTo(M_prime(cv::Rect(0,0,3,3)));
+  cv::undistortPoints(dist_points,undist_points,M,D,cv::Mat(),M_prime);
+  cv::undistortPoints(dist_def_points,undist_def_points,M,D,cv::Mat(),M_prime);
+
+  for(int_t i=0;i<num_subsets;++i){
+    subset_x_prime->local_value(i) = undist_points[i].x;
+    subset_y_prime->local_value(i) = undist_points[i].y;
+    subset_u_prime->local_value(i) = undist_def_points[i].x - undist_points[i].x;
+    subset_v_prime->local_value(i) = undist_def_points[i].y - undist_points[i].y;
+  }
+  DEBUG_MSG("Distortion_Correction_Post_Processor execute() end");
+}
 
 VSG_Strain_Post_Processor::VSG_Strain_Post_Processor(const Teuchos::RCP<Teuchos::ParameterList> & params) :
   Post_Processor(post_process_vsg_strain){

src/core/DICe_PostProcessor.h

@@ -45,6 +45,7 @@
 #include <DICe.h>
 #include <DICe_Mesh.h>
 #include <DICe_PointCloud.h>
+#include <DICe_Triangulation.h>
 
 #include <Teuchos_ParameterList.hpp>
 
@@ -82,17 +83,20 @@ const char * const post_process_live_plots = "post_process_live_plots";
 const char * const post_process_plotly_contour = "post_process_plotly_contour";
 /// String parameter name
 const char * const plotly_contour_grid_step = "plotly_contour_grid_step";
+/// String parameter name
+const char * const post_process_distortion_correction = "post_process_distortion_correction";
 
 
 /// Number of post processor options
-const int_t num_valid_post_processor_params = 5;
+const int_t num_valid_post_processor_params = 6;
 /// Set of all the valid post processors
 const char * const valid_post_processor_params[num_valid_post_processor_params] = {
   post_process_vsg_strain,
   post_process_nlvc_strain,
   post_process_plotly_contour,
   post_process_altitude,
-  post_process_crack_locator
+  post_process_crack_locator,
+  post_process_distortion_correction
 };
 
 /// String field name
@@ -213,7 +217,7 @@ Post_Processor{
 };
 
 /// \class DICe::VSG_Strain_Post_Processor
-/// \brief A specific instance of post processor that computes vurtual strain gauge (VSG) strain
+/// \brief A specific instance of post processor that computes virtual strain gauge (VSG) strain
 ///
 /// The VSG strain is computed by doing a least-squares fit of the data and computing the
 /// strain using the coefficients of the fitted polynomial. It is well know that for large
@@ -259,6 +263,46 @@ VSG_Strain_Post_Processor : public Post_Processor {
   int_t window_size_;
 };
 
+/// \class DICe::Distortion_Correction_Post_Processor
+/// \brief A specific instance of post processor that corrects the subset displacements for image distortions
+///
+/// This post processor corrects the subset displacement values for lens distortions and outputs the new corrected
+/// displacement field in a new field called DIST_CORRECTED_SUBSET_DISP_X and ..._Y
+class DICE_LIB_DLL_EXPORT
+Distortion_Correction_Post_Processor : public Post_Processor {
+
+public:
+
+  /// Default constructor
+  /// \param params Pointer to the set of parameters for this post processor
+  Distortion_Correction_Post_Processor(const Teuchos::RCP<Teuchos::ParameterList> & params);
+
+  /// Virtual destructor
+  virtual ~Distortion_Correction_Post_Processor(){}
+
+  /// See base clase docutmentation
+  virtual int_t strain_window_size(){return -1.0;}
+
+  /// See base clase docutmentation
+  virtual void set_params(const Teuchos::RCP<Teuchos::ParameterList> & params){}
+
+  /// Collect the neighborhoods of each of the points
+  virtual void pre_execution_tasks(){};
+
+  /// Execute the post processor
+  virtual void execute(Teuchos::RCP<Image> ref_img, Teuchos::RCP<Image> def_img);
+
+  /// See base class documentation
+  using Post_Processor::field_specs;
+
+protected:
+
+  /// Triangulation to use for correcting displacements
+  Teuchos::RCP<DICe::Triangulation> tri_;
+
+
+};
+
 /// \class DICe::NLVC_Strain_Post_Processor
 /// \brief A specific instance of post processor that computes Nonlocal Vector Calculus (NLVC) strains
 ///

src/core/DICe_Schema.cpp

@@ -79,22 +79,22 @@ Schema::Schema(const std::string & input_file_name,
   const std::string & params_file_name){
   // create a parameter list from the selected file
   Teuchos::RCP<Teuchos::ParameterList> corr_params = read_correlation_params(params_file_name);
-  default_constructor_tasks(corr_params);
   // create a parameter list from the selected file
   Teuchos::RCP<Teuchos::ParameterList> input_params = read_input_params(input_file_name);
+  default_constructor_tasks(corr_params,input_params);
   initialize(input_params,corr_params);
 }
 
 Schema::Schema(const Teuchos::RCP<Teuchos::ParameterList> & input_params,
   const Teuchos::RCP<Teuchos::ParameterList> & correlation_params){
-  default_constructor_tasks(correlation_params);
+  default_constructor_tasks(correlation_params,input_params);
   initialize(input_params,correlation_params);
 }
 
 Schema::Schema(const Teuchos::RCP<Teuchos::ParameterList> & input_params,
   const Teuchos::RCP<Teuchos::ParameterList> & correlation_params,
   const Teuchos::RCP<Schema> & schema){
-  default_constructor_tasks(correlation_params);
+  default_constructor_tasks(correlation_params,input_params);
   initialize(input_params,schema);
 }
 
@@ -445,7 +445,8 @@ Schema::set_ref_image(Teuchos::RCP<Image> img){
 }
 
 void
-Schema::default_constructor_tasks(const Teuchos::RCP<Teuchos::ParameterList> & params){
+Schema::default_constructor_tasks(const Teuchos::RCP<Teuchos::ParameterList> & corr_params,
+    const Teuchos::RCP<Teuchos::ParameterList> & input_params){
   global_num_subsets_ = 0;
   local_num_subsets_ = 0;
   subset_dim_ = -1;
@@ -462,8 +463,15 @@ Schema::default_constructor_tasks(const Teuchos::RCP<Teuchos::ParameterList> & p
   normalize_gamma_with_active_pixels_ = false;
   gauss_filter_images_ = false;
   gauss_filter_mask_size_ = 7;
-  init_params_ = params==Teuchos::null ? Teuchos::rcp(new Teuchos::ParameterList()):
-    Teuchos::rcp(new Teuchos::ParameterList(*params));
+  init_params_ = corr_params==Teuchos::null ? Teuchos::rcp(new Teuchos::ParameterList()):
+    Teuchos::rcp(new Teuchos::ParameterList(*corr_params));
+  if(input_params!=Teuchos::null){
+    // copy over the camera system file to so both the input and correlation params have it
+    // could be needed by the post processors
+    if(input_params->isParameter(DICe::camera_system_file)){
+      init_params_->set(DICe::camera_system_file,input_params->get<std::string>(DICe::camera_system_file));
+    }
+  }
   comm_ = Teuchos::rcp(new MultiField_Comm());
   path_file_names_ = Teuchos::rcp(new std::map<int_t,std::string>());
   optical_flow_flags_ = Teuchos::rcp(new std::map<int_t,bool>());
@@ -478,7 +486,7 @@ Schema::default_constructor_tasks(const Teuchos::RCP<Teuchos::ParameterList> & p
   read_full_images_ = false;
   sort_txt_output_ = false;
   threshold_block_size_ = -1;
-  set_params(params);
+  set_params(corr_params);
   prev_imgs_.push_back(Teuchos::null);
   def_imgs_.push_back(Teuchos::null);
   has_extents_ = false;
@@ -771,6 +779,13 @@ Schema::set_params(const Teuchos::RCP<Teuchos::ParameterList> & params){
     Teuchos::RCP<VSG_Strain_Post_Processor> vsg_ptr = Teuchos::rcp (new VSG_Strain_Post_Processor(ppParams));
     post_processors_.push_back(vsg_ptr);
   }
+  if(diceParams->isParameter(DICe::post_process_distortion_correction)){
+    Teuchos::RCP<Teuchos::ParameterList> ppParams = Teuchos::rcp( new Teuchos::ParameterList());
+    TEUCHOS_TEST_FOR_EXCEPTION(!init_params_->isParameter(camera_system_file),std::runtime_error,"Error, Distortion Correction Post Processor requires a camera system file");
+    ppParams->set<std::string>(camera_system_file,init_params_->get<std::string>(camera_system_file));
+    Teuchos::RCP<Distortion_Correction_Post_Processor> dc_ptr = Teuchos::rcp (new Distortion_Correction_Post_Processor(ppParams));
+    post_processors_.push_back(dc_ptr);
+  }
   if(diceParams->isParameter(DICe::post_process_plotly_contour)){
     Teuchos::ParameterList sublist = diceParams->sublist(DICe::post_process_plotly_contour);
     Teuchos::RCP<Teuchos::ParameterList> ppParams = Teuchos::rcp( new Teuchos::ParameterList());
@@ -977,10 +992,6 @@ Schema::initialize(const Teuchos::RCP<Teuchos::ParameterList> & input_params,
   const std::string output_prefix = input_params->get<std::string>(DICe::output_prefix,"DICe_solution");
   init_params_->set(DICe::output_prefix,output_prefix);
   init_params_->set(DICe::output_folder,output_folder);
-  if(input_params->isParameter(DICe::camera_system_file)){
-    const std::string camera_sys_file = input_params->get<std::string>(DICe::camera_system_file);
-    init_params_->set(DICe::camera_system_file,camera_sys_file);
-  }
 
   if(analysis_type_==GLOBAL_DIC){
     // create the computational mesh:

src/core/DICe_Schema.h

@@ -1152,7 +1152,8 @@ Schema {
 
   /// \brief Sets the default values for the schema's member data and other initialization tasks
   /// \param params Optional correlation parameters
-  void default_constructor_tasks(const Teuchos::RCP<Teuchos::ParameterList> & params);
+  void default_constructor_tasks(const Teuchos::RCP<Teuchos::ParameterList> & corr_params,
+      const Teuchos::RCP<Teuchos::ParameterList> & input_params=Teuchos::null);
 
   /// \brief Create an exodus mesh for output
   /// \param decomp pointer to a decomposition

src/core/DICe_Triangulation.cpp

@@ -271,6 +271,12 @@ Triangulation::load_calibration_parameters(const std::string & param_file_name){
   DEBUG_MSG("Triangulation::load_calibration_parameters(): Parsing calibration parameters from file: " << param_file_name);
 
   camera_system_ = Teuchos::rcp(new Camera_System(param_file_name));
+
+  // set the camera intrinsic parameters
+  for(size_t i=0;i<camera_system_->num_cameras();++i)
+    for(size_t j=0;j<Camera::MAX_CAM_INTRINSIC_PARAM;++j)
+    cal_intrinsics_[i][j] = (*camera_system_->camera(i)->intrinsics())[j];
+
   if(camera_system_->num_cameras()==1) return;
 
   TEUCHOS_TEST_FOR_EXCEPTION(camera_system_->num_cameras()!=2,std::runtime_error,"");
@@ -296,11 +302,6 @@ Triangulation::load_calibration_parameters(const std::string & param_file_name){
     cam_0_to_cam_1_ = cam_0_to_cam_1_ * cam_0_to_world_;
   }
 
-  // set the camera intrinsic parameters
-  for(size_t i=0;i<camera_system_->num_cameras();++i)
-    for(size_t j=0;j<Camera::MAX_CAM_INTRINSIC_PARAM;++j)
-    cal_intrinsics_[i][j] = (*camera_system_->camera(i)->intrinsics())[j];
-
 #ifdef DICE_DEBUG_MSG
   std::cout << *camera_system_.get() << std::endl;
 #endif

tests/regression/CMakeLists.txt

@@ -15,6 +15,7 @@ SET(REGRESSION_TESTS
   simple_tracking
   reset_lost_subsets
   epipolar_distortion
+  distortion_correction_pp
 )
 IF(NOT UNIX)
   SET(REGRESSION_TESTS

tests/regression/distortion_correction_pp/cal.xml

@@ -0,0 +1,20 @@
+<ParameterList>
+  <Parameter name="xml_file_format" type="string" value="DICe_xml_camera_system_file" />  
+  <Parameter name="system_type_3D" type="string" value="OPENCV" />  
+  <ParameterList name="CAMERA 0">
+    <Parameter name="CAMERA_ID" type="string" value="left" />  
+    <Parameter name="CX" type="double" value="959.5" />  
+    <Parameter name="CY" type="double" value="539.5" />  
+    <Parameter name="FX" type="double" value="2666.6666666666667" />  
+    <Parameter name="FY" type="double" value="2666.6666666666667" />  
+    <Parameter name="FS" type="double" value="0.0" />  
+    <Parameter name="TX" type="double" value="0.0" />  
+    <Parameter name="TY" type="double" value="0.0" />  
+    <Parameter name="TZ" type="double" value="15.0" />  
+    <Parameter name="K1" type="double" value="-2.361956054191e-01" />
+    <Parameter name="K2" type="double" value="-9.858048851165e-02" />
+    <Parameter name="K3" type="double" value="-3.652319258538e-01" />
+    <Parameter name="LENS_DISTORTION_MODEL" type="string" value="OPENCV_LENS_DISTORTION" />
+    <Parameter name="IMAGE_HEIGHT_WIDTH" type="string" value="{ 1080, 1920 }" />  
+  </ParameterList>
+</ParameterList>