1 #include <iostream>
2 #include <iomanip>
3 #include <bvrml/bvrml_write.h>
4 #include "boxm2_vecf_orbit_param_stats.h"
5 #include "boxm2_vecf_plot_orbit.h"
6 #include <vnl/algo/vnl_svd.h>
7 #ifdef _MSC_VER
8 #  include "vcl_msvc_warnings.h"
9 #endif
average_params()10 void boxm2_vecf_orbit_param_stats::average_params(){
11   for(auto & pit : param_map_){
12     const std::string& pid = pit.first;
13     boxm2_vecf_orbit_params avg;
14     boxm2_vecf_orbit_params& left_params = (pit.second).first;
15     boxm2_vecf_orbit_params& right_params = (pit.second).second;
16     avg.eye_radius_= 0.5*(left_params.eye_radius_ + right_params.eye_radius_);
17     devs_[pid]["eye_radius"]=std::fabs(left_params.eye_radius_-avg.eye_radius_);
18     avg.iris_radius_= 0.5*(left_params.iris_radius_ + right_params.iris_radius_);
19     devs_[pid]["iris_radius"]=std::fabs(left_params.iris_radius_-avg.iris_radius_);
20     avg.medial_socket_radius_coef_ = 0.5*(left_params.medial_socket_radius_coef_ + right_params.medial_socket_radius_coef_);
21     devs_[pid]["medial_socket_radius_coef"]=std::fabs(left_params.medial_socket_radius_coef_-avg.medial_socket_radius_coef_);
22     avg.lateral_socket_radius_coef_ = 0.5*(left_params.lateral_socket_radius_coef_ + right_params.lateral_socket_radius_coef_);
23     devs_[pid]["lateral_socket_radius_coef"]=std::fabs(left_params.lateral_socket_radius_coef_-avg.lateral_socket_radius_coef_);
24     avg.inferior_lid_radius_offset_ = 0.5*(left_params.inferior_lid_radius_offset_ + right_params.inferior_lid_radius_offset_);
25     devs_[pid]["inferior_lid_radius_offset"]=std::fabs(left_params.inferior_lid_radius_offset_-avg.inferior_lid_radius_offset_);
26     avg.eyelid_crease_scale_y_coef_ = 0.5*(left_params.eyelid_crease_scale_y_coef_ + right_params.eyelid_crease_scale_y_coef_);
27     devs_[pid]["eyelid_crease_scale_y_coef"]=std::fabs(left_params.eyelid_crease_scale_y_coef_-avg.eyelid_crease_scale_y_coef_);
28     double h_to_w = left_params.height_to_width_ratio();
29     h_to_w += right_params.height_to_width_ratio();
30     avg.height_to_width_ratio_ = 0.5*h_to_w;
31     devs_[pid]["height_to_width_ratio"]=std::fabs(left_params.height_to_width_ratio_-avg.height_to_width_ratio_);
32     avg.scale_x_coef_ =  0.5*(left_params.scale_x_coef_ + right_params.scale_x_coef_);
33     devs_[pid]["scale_x_coef"]=std::fabs(left_params.scale_x_coef_-avg.scale_x_coef_);
34     avg.scale_y_coef_ =  0.5*(left_params.scale_y_coef_ + right_params.scale_y_coef_);
35     devs_[pid]["scale_y_coef"]=std::fabs(left_params.scale_y_coef_-avg.scale_y_coef_);
36     avg.mid_eyelid_crease_z_ = 0.5*(left_params.mid_eyelid_crease_z_ + right_params.mid_eyelid_crease_z_);
37     devs_[pid]["mid_eyelid_crease_z"]=std::fabs(left_params.mid_eyelid_crease_z_-avg.mid_eyelid_crease_z_);
38     param_avg_[pit.first] = avg;
39   }
40 }
generate_stats()41 void boxm2_vecf_orbit_param_stats::generate_stats(){
42   //average left and right vals
43   this->average_params();
44 }
print_stats()45 void boxm2_vecf_orbit_param_stats::print_stats() {
46   std::cout << "patient_id  eye_radius iris_radius lid_radius_off height/width cr_scale_y crease_z_rat\n";
47   for(auto & pit : param_avg_){
48     const boxm2_vecf_orbit_params& pr = pit.second;
49     const std::string& pid = pit.first;
50     std::cout << pid << ' ' << pr.eye_radius_ << ' '<< devs_[pid]["eye_radius"]<< ' '
51              << pr.iris_radius_ << ' '<< devs_[pid]["iris_radius"] << ' '
52              << pr.inferior_lid_radius_offset_ << ' ' << devs_[pid]["inferior_lid_radius_offset"];
53     double rat = pr.height_to_width_ratio_;
54     std::cout << ' '<< rat << ' '<< devs_[pid]["height_to_width_ratio"]
55              << ' '<< pr.eyelid_crease_scale_y() << ' ' << devs_[pid]["eyelid_crease_scale_y_coef"]*pr.eye_radius_
56              << ' ' << pr.mid_eyelid_crease_z_/pr.lid_sph_.radius() << ' ' << devs_[pid]["mid_eyelid_crease_z"]/pr.lid_sph_.radius()<< '\n';
57   }
58 }
print_xy_fitting_error()59 void boxm2_vecf_orbit_param_stats::print_xy_fitting_error(){
60   std::cout << "patient_id  left_inf_xy left_sup_xy left_sup_crease right_inf_xy right_sup_xy right_sup_crease_xy \n";
61   for(auto & pit : param_map_){
62     const std::string& pid = pit.first;
63     boxm2_vecf_orbit_params& lp = (pit.second).first;
64     boxm2_vecf_orbit_params& rp = (pit.second).second;
65 
66     std::cout << pid << ' ' << lp.inferior_margin_xy_error_ << ' '
67              << lp.superior_margin_xy_error_ << ' '
68                    << lp.superior_crease_xy_error_ << ' '
69              << rp.inferior_margin_xy_error_ << ' '
70                    << rp.superior_margin_xy_error_ << ' '
71                    << rp.superior_crease_xy_error_ << '\n';
72  }
73 }
74 
print_xyz_fitting_error()75 void boxm2_vecf_orbit_param_stats::print_xyz_fitting_error(){
76   std::cout << "patient_id  left_inf_xyz left_sup_xyz left_sup_crease right_inf_xyz right_sup_xyz right_sup_crease_xyz \n";
77   for(auto & pit : param_map_){
78     const std::string& pid = pit.first;
79     boxm2_vecf_orbit_params& lp = (pit.second).first;
80     boxm2_vecf_orbit_params& rp = (pit.second).second;
81 
82     std::cout << pid << ' ' << lp.inferior_margin_xyz_error_ << ' '
83              << lp.superior_margin_xyz_error_ << ' '
84                    << lp.superior_crease_xyz_error_ << ' '
85              << rp.inferior_margin_xyz_error_ << ' '
86                    << rp.superior_margin_xyz_error_ << ' '
87                    << rp.superior_crease_xyz_error_ << '\n';
88  }
89 }
merge_margins_and_crease()90 bool boxm2_vecf_orbit_param_stats::merge_margins_and_crease(){
91   for(auto & pit : param_map_){
92     const std::string& pid = pit.first;
93     boxm2_vecf_orbit_params lp = (pit.second).first;
94     boxm2_vecf_orbit_params rp = (pit.second).second;
95 
96     // plot left eye
97     double xm_min_left = lp.x_min()-10.0;
98     double xm_max_left = lp.x_max()+10.0;
99     std::vector<vgl_point_3d<double> > left_inf_pts, left_sup_pts, left_crease_pts;
100     boxm2_vecf_plot_orbit::plot_inferior_margin(lp, false, xm_min_left, xm_max_left, left_inf_pts);
101     boxm2_vecf_plot_orbit::plot_superior_margin(lp, false, xm_min_left, xm_max_left, left_sup_pts);
102     boxm2_vecf_plot_orbit::plot_crease(lp, false, xm_min_left, xm_max_left, left_crease_pts);
103     int left_imin = -1, left_imax = -1;
104     bool success = boxm2_vecf_plot_orbit::plot_limits(left_inf_pts, left_sup_pts, left_imin, left_imax);
105     if(!success){
106       std::cout << "Find left plot limits failed for "<< pid << "\n";
107       continue;
108     }
109     // get left medial canthus (as origin)
110     vgl_point_3d<double> plc_inf = left_inf_pts[left_imin];
111     vgl_point_3d<double> plc_sup = left_sup_pts[left_imin];
112 #if 1// medial canthus origin (leads to excessive scatter)
113     vgl_vector_3d<double> vlc(0.5*(plc_inf.x() + plc_sup.x()),
114                               0.5*(plc_inf.y() + plc_sup.y()),
115                               0.5*(plc_inf.z() + plc_sup.z()));
116 #endif
117     // use sclera-derived z origin and palpebral x-y origin instead
118     //vgl_vector_3d<double> vlc(0.0, 0.0, lp.eyelid_radius());
119     // plot right eye
120     double xm_min_right = rp.x_min()-10.0;
121     double xm_max_right = rp.x_max()+10.0;
122     std::vector<vgl_point_3d<double> > right_inf_pts, right_sup_pts, right_crease_pts;
123     boxm2_vecf_plot_orbit::plot_inferior_margin(rp, false, xm_min_right, xm_max_right, right_inf_pts);
124     boxm2_vecf_plot_orbit::plot_superior_margin(rp, false, xm_min_right, xm_max_right, right_sup_pts);
125     boxm2_vecf_plot_orbit::plot_crease(rp, false, xm_min_right, xm_max_right, right_crease_pts);
126     int right_imin = -1, right_imax = -1;
127     success = boxm2_vecf_plot_orbit::plot_limits(right_inf_pts, right_sup_pts, right_imin, right_imax);
128     if(!success){
129       std::cout << "Find right plot limits failed\n";
130       return false;
131     }
132     // get right medial canthus (as origin)
133     vgl_point_3d<double> prc_inf = right_inf_pts[right_imin];
134     vgl_point_3d<double> prc_sup = right_sup_pts[right_imin];
135 #if 1// medial canthus origin (leads to excessive scatter)
136     vgl_vector_3d<double> vrc(0.5*(prc_inf.x() + prc_sup.x()),
137                               0.5*(prc_inf.y() + prc_sup.y()),
138                               0.5*(prc_inf.z() + prc_sup.z()));
139 #endif
140     // use sclera-derived z origin and palpebral x-y origin instead
141     //vgl_vector_3d<double> vrc(0.0, 0.0, rp.eyelid_radius());
142     // merge margins and crease
143     std::vector<vgl_point_3d<double> > inf_pts, sup_pts, crease_pts;
144     for(int i = left_imin; i<=left_imax; ++i){
145       inf_pts.push_back(left_inf_pts[i]-vlc);
146       sup_pts.push_back(left_sup_pts[i]-vlc);
147       crease_pts.push_back(left_crease_pts[i]-vlc);
148     }
149     for(int i = right_imin; i<=right_imax; ++i){
150       inf_pts.push_back(right_inf_pts[i]-vrc);
151       sup_pts.push_back(right_sup_pts[i]-vrc);
152       crease_pts.push_back(right_crease_pts[i]-vrc);
153     }
154     merged_inf_margin_[pid]=inf_pts;
155     merged_sup_margin_[pid]=sup_pts;
156     merged_crease_[pid]=crease_pts;
157 
158   }
159   return true;
160 }
plot_merged_margins(std::ofstream & os,unsigned sample_skip)161 bool boxm2_vecf_orbit_param_stats::plot_merged_margins(std::ofstream& os, unsigned sample_skip){
162   if(!os){
163     std::cout << "bad ostream in boxm2_vecf_orbit_param_stats\n";
164     return false;
165   }
166   bvrml_write::write_vrml_header(os);
167   float r = 0.5f; // error range 1.0 mm
168   // write inferior margins
169   for(auto & pit : merged_inf_margin_){
170      std::vector<vgl_point_3d<double> > inf_pts = pit.second;
171      auto n = static_cast<unsigned>(inf_pts.size());
172      for(unsigned i = 0; i<n; i+=sample_skip){
173        vgl_point_3d<double> pd = inf_pts[i];
174        vgl_point_3d<float> pf(static_cast<float>(pd.x()),
175                               static_cast<float>(pd.y()),
176                               static_cast<float>(pd.z()));
177        vgl_sphere_3d<float> sp(pf, r);
178        bvrml_write::write_vrml_sphere(os, sp, 1.0f, 1.0f, 0.0f);
179      }
180   }
181   // write superior margins
182   for(auto & pit : merged_sup_margin_){
183     std::vector<vgl_point_3d<double> > sup_pts = pit.second;
184     auto n = static_cast<unsigned>(sup_pts.size());
185     for(unsigned i = 0; i<n; i+=sample_skip){
186       vgl_point_3d<double> pd = sup_pts[i];
187       vgl_point_3d<float> pf(static_cast<float>(pd.x()),
188                              static_cast<float>(pd.y()),
189                              static_cast<float>(pd.z()));
190       vgl_sphere_3d<float> sp(pf, r);
191       bvrml_write::write_vrml_sphere(os, sp, 1.0f, 0.0f, 0.0f);
192     }
193   }
194   // write creases
195   for(auto & pit : merged_crease_){
196     std::vector<vgl_point_3d<double> > crease_pts = pit.second;
197     auto n = static_cast<unsigned>(crease_pts.size());
198     for(unsigned i = 0; i<n; i+=sample_skip){
199       vgl_point_3d<double> pd = crease_pts[i];
200       vgl_point_3d<float> pf(static_cast<float>(pd.x()),
201                              static_cast<float>(pd.y()),
202                              static_cast<float>(pd.z()));
203       vgl_sphere_3d<float> sp(pf, r);
204       bvrml_write::write_vrml_sphere(os, sp, 0.0f, 1.0f, 1.0f);
205     }
206   }
207   os.close();
208   return true;
209 }
compute_feature_vectors()210 void boxm2_vecf_orbit_param_stats::compute_feature_vectors(){
211   unsigned dim = 12;
212   for(auto & mit : merged_inf_margin_){
213     std::string pid = mit.first;
214     std::vector<vgl_point_3d<double> > inf_pts = mit.second;
215     std::vector<vgl_point_3d<double> > sup_pts = merged_sup_margin_[pid];
216     std::vector<vgl_point_3d<double> > crease_pts = merged_crease_[pid];
217     vgl_point_3d<double> pl = inf_pts[inf_pts.size()-1];
218     vgl_point_3d<double> pmini, pmaxs, pmaxc;
219     unsigned n = 0;
220 #if 1
221     n = static_cast<unsigned>(inf_pts.size());
222     double yi = 1000.0;
223     for(unsigned i = 0; i<n; ++i){
224       vgl_point_3d<double> pi = inf_pts[i];
225       if(pi.y()<yi){
226         yi = pi.y();
227         pmini = pi;
228       }
229     }
230 #endif
231     n = static_cast<unsigned>(sup_pts.size());
232     double ys = -1000.0;
233     for(unsigned i = 0; i<n; ++i){
234       vgl_point_3d<double> ps = sup_pts[i];
235       if(ps.y()>ys){
236         ys = ps.y();
237         pmaxs = ps;
238       }
239     }
240     n = static_cast<unsigned>(crease_pts.size());
241     double yc = -1000.0;
242     for(unsigned i = 0; i<n; ++i){
243       vgl_point_3d<double> pc = crease_pts[i];
244       if(pc.y()>yc){
245         yc = pc.y();
246         pmaxc = pc;
247       }
248     }
249     vnl_matrix<double> fv(dim,1,0.0);
250 #if 0
251     double alpha = (ys-yi)/pl.x();
252     double zeta = 0.5*(pmaxc.z()-pl.z())/8.7;
253     double kappa = (yc-ys)/ys;
254     double eta = (pmaxs.z()-pmaxc.z())/pmaxs.z();
255     fv[0][0]= alpha;     fv[1][0]= zeta;
256     fv[2][0]= kappa;     fv[3][0]= eta;
257 #endif
258 #if 1
259     fv[0][0]=pl.x(); fv[1][0]=pl.y(); fv[2][0]= pl.z();
260     fv[3][0]=pmaxs.x(); fv[4][0]=pmaxs.y(); fv[5][0]= pmaxs.z();
261     fv[6][0]=pmaxc.x(); fv[7][0]=pmaxc.y(); fv[8][0]= pmaxc.z();
262     fv[9][0]=pmini.x(); fv[10][0]=pmini.y(); fv[11][0]= pmini.z();
263 #endif
264     feature_vectors_[pid]=fv;
265   }
266 }
267 
compute_covariance_matrix()268 void boxm2_vecf_orbit_param_stats::compute_covariance_matrix(){
269   unsigned dim = 12; //dimension of the feature vector
270   // compute mean vector
271   mean_ = vnl_matrix<double>(dim, 1, 0.0);
272   double nv = 0.0;
273   for(auto fit = feature_vectors_.begin();
274       fit != feature_vectors_.end(); ++fit, nv+=1.0)
275     mean_ += fit->second;
276   mean_/=nv;
277   std::cout << "mean feature vector\n";
278   for(unsigned i = 0; i<dim; ++i)
279     std::cout << std::setprecision(2) << mean_[i][0] << ' ';
280   std::cout << '\n';
281   cov_ = vnl_matrix<double>(dim, dim, 0.0);
282   auto m = static_cast<unsigned>(feature_vectors_.size());
283   // the data matrix - rows denote the elements of the feature vector, columns each sample
284   vnl_matrix<double> X(dim, m);
285   unsigned k = 0;
286   for(auto fit = feature_vectors_.begin();
287       fit != feature_vectors_.end(); ++fit, ++k){
288     vnl_matrix<double> x = (fit->second - mean_);
289     for(unsigned r = 0; r<static_cast<unsigned>(x.rows()); ++r)
290       X[r][k]=x[r][0];
291   }
292   vnl_svd<double> svd(X);
293   vnl_matrix<double> W = svd.W();
294   vnl_matrix<double> Wsq = W*W/m;
295   std::cout << Wsq << '\n';
296   vnl_matrix<double> U = svd.U();
297   vnl_matrix<double> Ut = U.transpose();
298   // the primary eigenvector
299   std::cout << "primary eigenvector\n";
300   for(unsigned i = 0; i<dim; ++i)
301     std::cout << std::setprecision(2) << Ut[0][i] << ' ';
302   std::cout << '\n';
303   std::cout << "secondary eigenvector\n";
304   for(unsigned i = 0; i<dim; ++i)
305     std::cout << std::setprecision(2) << Ut[1][i] << ' ';
306   std::cout << '\n';
307   cov_ = U*Wsq*Ut;
308   // a m x m matrix only the first dim rows are meaningful.
309   // the columns are the m data samples transformed so as
310   // to produce a diagonal covariance matrix
311   vnl_matrix<double> Xw = Ut*X;
312   k = 0;
313   for(auto fit = feature_vectors_.begin();
314       fit != feature_vectors_.end(); ++fit, ++k){
315     std::string pid = fit->first;
316     std::cout << pid << ' ' << Xw[0][k] << ' ' << Xw[1][k] << '\n';
317   }
318 }
separation_stats()319 void boxm2_vecf_orbit_param_stats::separation_stats(){
320   double avg_ratio = 0.0;
321   double ns = 0.0;
322   for(auto pit =  param_map_.begin();
323       pit!=param_map_.end(); ++pit, ns += 1.0){
324     const std::string& pid = pit->first;
325     boxm2_vecf_orbit_params& left_params = (pit->second).first;
326     boxm2_vecf_orbit_params& right_params = (pit->second).second;
327     double avg_eye_radius= 0.5*(left_params.eye_radius_ + right_params.eye_radius_);
328     vgl_vector_3d<double> org_l(left_params.x_trans(), left_params.y_trans(), left_params.z_trans());
329     vgl_vector_3d<double> org_r(right_params.x_trans(), right_params.y_trans(), right_params.z_trans());
330     vgl_vector_3d<double> dif = org_l - org_r;
331     double sep = dif.length();
332     double ratio = sep/avg_eye_radius;
333     avg_ratio += ratio;
334     std::cout << pid << ' ' << avg_eye_radius << ' ' << sep << ' ' << ratio << '\n';
335   }
336   std::cout << "Average separation ratio " << avg_ratio/ns << '\n';
337 }
338