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