1 /**********
2 *
3 * This file is part of the TINA Open Source Image Analysis Environment
4 * henceforth known as TINA
5 *
6 * TINA is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as
8 * published by the Free Software Foundation.
9 *
10 * TINA is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU Lesser General Public License for more details.
14 *
15 * You should have received a copy of the GNU Lesser General Public License
16 * along with TINA; if not, write to the Free Software Foundation, Inc.,
17 * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
18 *
19 **********
20 *
21 * Program : TINA
22 * File : $Source: visPgh_locate.c $
23 * Date : $Date: 2013/05/13 13:50 $
24 * Version : $Revision: 2.1 $
25 *
26 * Author : Legacy TINA; modified NAT/HR
27 */
28 /**
29 * @file Routines for analysis of scene contents using PGH.
30 * @brief Functions for locating individual parts of an object or multiple instances
31 * of specified object in a line based approximation to an edge map.
32 *
33 */
34
35
36 #include <stdio.h>
37 #include <string.h>
38 #include <math.h>
39
40 #include <tina/image/img_GenDef.h>
41 #include <tina/image/img_GenPro.h>
42 #include <tina/image/img_PrcPro.h>
43 #include <tina/math/math_GeomDef.h>
44 #include <tina/math/mathGeom_geom2.h>
45 #include <tina/math/mathTran_rot2.h>
46 #include <tina/math/mathGeom_vec2.h>
47 #include <tina/math/math_TranDef.h>
48 #include <tina/math/mathGeom_mat2.h>
49 #include <tina/file/file_UtilPro.h>
50 #include <tina/geometry/geomGen_free.h>
51 #include <tina/geometry/geomDef.h>
52 #include <tina/geometry/geom_LinePro.h>
53 #include <tina/geometry/geom_PointPro.h>
54 #include <tina/geometry/geom_CurvePro.h>
55 #include <tina/vision/vis_ShapeDef.h>
56 #include <tina/vision/visShape_hough2.h>
57 #include <tina/vision/visShape_hough1.h>
58 #include <tina/vision/visPgh_var.h>
59 #include <tina/vision/visPgh_HT.h>
60 #include <tina/vision/visPgh_hist.h>
61 #include <tina/vision/visPgh_misc.h>
62 #include <tina/vision/visPgh_match.h>
63 #include <tina/vision/visPgh_histfuncs.h>
64 #include <tina/vision/visPgh_locate.h>
65
66
67 double pgh_scale_estimate;
68
69 /* ---------- */
70
71 Point2 *line_to_model(List *geom)
72 {
73 Vec2 translation;
74 Mat2 rotation_matrix;
75 List *ptr;
76 /*List *geom=NULL;*/
77 List *match_list;
78 Line2 *lptr;
79 Imrect *im, *match_hist;
80 double rotation;
81 Hist_ref *scene_hist_ref, *match_hist_ref;
82 Transform2 model_to_scene_trans;
83 Match_ref *match_ref;
84 List *model_pairs_list;
85 List *matched_models_list;
86 Line2 *matched_line;
87 List *picked_geom;
88 Pairs_hist_def *hist_def;
89 Pairs_scale_def *scale_def;
90
91
92 /* (1) Generate a pairwise histogram for the currently
93 selected line. (Or the first line if a number
94 of lines have been selected).
95 (2) Find the best match between this histogram and
96 all of the model histograms.
97 (3) Report this match to the main Tina window.
98 (4) Overlay the model corresponding to the best match
99 onto the scene. Allow the user to transform
100 between the 4 possible locations. */
101
102 model_pairs_list = pairs_model_pairs_list_get();
103 matched_models_list = pairs_matched_models_list_get();
104 matched_line = pairs_matched_line_get();
105 picked_geom = pairs_picked_geom_get();
106 hist_def = pairs_hist_def_get();
107 scale_def = pairs_scale_def_get();
108
109 init_pairs_entry(hist_def->pairs_type,
110 hist_def->dbin_max,
111 hist_def->dbin_size,
112 hist_def->num_abin,
113 hist_def->angle_sigma,
114 hist_def->dist_ramp);
115
116 /*geom = mono_geom_get();*/
117 if (geom==NULL||picked_geom==NULL||model_pairs_list==NULL) return(NULL);
118
119 lptr = picked_geom->to;
120
121 /* Build the pairwise histogram */
122
123 im = build_normalized_pairwise(lptr, NULL, geom, hist_def->dbin_max);
124
125 /* Now search the model database for the best match */
126
127 match_scene_pairs_to_model_list(im, model_pairs_list, scale_def);
128
129 /* Fetch the best (first) match from the histogram matches list */
130
131 scene_hist_ref = hist_ref_get(im);
132 match_list = scene_hist_ref->matches;
133 if (match_list==NULL) return(NULL);
134 match_ref = match_list->to;
135 match_hist = match_ref->hist;
136 match_hist_ref = hist_ref_get(match_hist);
137
138 /* Report the best match */
139
140 printf("Scale:%f Score:%f\n", match_ref->scale_factor, match_ref->dp);
141
142 /* If requested output the scores for matches to this line at different scales */
143
144 /* Now overlay the model over the scene */
145 /* set the p-field of the scene line to its mid-point */
146
147 lptr->p = vec2_midpoint(lptr->p1, lptr->p2);
148
149 /* First transform the model line to fit the scene line */
150
151 /* Initially determine the rotation required to map the model line to the
152 same orientation as the scene line */
153
154 rotation = vec2_angle( lptr->v, match_hist_ref->line_seg->v);
155
156 /* Construct the rotation matrix with the appropiate scale factor */
157
158 rotation_matrix = rot2_with_scale(rotation, match_ref->scale_factor);
159
160 /* Now determine the translation required to map the midpoint of the model
161 line onto the midpoint of the scene line. This mapping has 2 components:
162
163 (1) Mapping the centre point of the model line back to its position before
164 the rotation.
165 (2) Mapping the centre point of the models line from its original position
166 to the midpoint of the scene line. */
167
168 translation = vec2_diff(match_hist_ref->line_seg->p,
169 mat2_vprod(rotation_matrix, match_hist_ref->line_seg->p));
170
171 translation = vec2_sum(translation,
172 vec2_diff(lptr->p, match_hist_ref->line_seg->p));
173
174
175 model_to_scene_trans.R = rotation_matrix;
176 model_to_scene_trans.t = translation;
177
178 /* Now copy the model to the matched model data */
179
180 list_rm(matched_models_list, geom_free);
181 matched_models_list=NULL;
182
183 for(ptr=match_hist_ref->model->model_poly_data;ptr!=NULL;ptr=ptr->next)
184 {
185 lptr = line2_copy(ptr->to);
186 matched_models_list =
187 list_addtoend(matched_models_list, link_alloc(lptr, (void *)NULL/*type*/));
188 if (ptr->to==match_hist_ref->line_seg)
189 matched_line = lptr;
190 }
191
192 pairs_matched_models_list_set(matched_models_list);
193 pairs_matched_line_set(matched_line);
194
195 /* Now transform the copied model data */
196
197 for(ptr=matched_models_list;ptr!=NULL;ptr=ptr->next)
198 line2_transform(ptr->to, model_to_scene_trans);
199
200 if(match_ref->orientation == -1)
201 rotate_matched_model();
202 else
203 {
204 /*draw_matched_models();*/
205 }
206
207 return(point2_make(translation, double_v));
208 }
209
210 /* ---------- */
211
212 List * segment_model(char *model_name)
213 {
214 List *ptr, *match_list, *cut_match_list;
215 List *segmented_lines_list = NULL;
216 Imrect *scene_hist, *match_hist;
217 Hist_ref *scene_hist_ref, *match_hist_ref;
218 int found_match;
219 Match_ref *match_ref;
220 List *scene_pairs_list;
221 Match_cut_def *match_cut_def;
222
223 /* (1) Descend the scene pairwise list adding each
224 scene line to the segmented_lines_list if it
225 has a match corresponding to the specified model.
226 (2) Display the segmented_lines_list. */
227
228 scene_pairs_list = pairs_scene_pairs_list_get();
229 match_cut_def = pairs_match_cut_def_get();
230
231 if (scene_pairs_list==NULL)
232 return(NULL);
233
234 strip_spaces(model_name); /* prepare the model name for comparison later */
235
236 for(ptr=scene_pairs_list;ptr!=NULL;ptr=ptr->next)
237 {
238 scene_hist = ptr->to;
239 scene_hist_ref = prop_get(scene_hist->props, HIST_REF_TYPE);
240 found_match=FALSE;
241
242 if (scene_hist_ref!=NULL)
243 {
244 /* Make a cut copy of the matches list */
245
246 cut_match_list = copy_and_cut_matches_list(scene_hist_ref,
247 match_cut_def);
248
249 for(match_list=cut_match_list;match_list!=NULL;
250 match_list=match_list->next)
251 {
252 match_ref = match_list->to;
253 match_hist = match_ref->hist;
254 if (match_hist!=NULL)
255 {
256 match_hist_ref = prop_get(match_hist->props,
257 HIST_REF_TYPE);
258 if (match_hist_ref!=NULL &&
259 strcmp(model_name, match_hist_ref->model->name)==0)
260 found_match=TRUE;
261 }
262 }
263
264 if (found_match)
265 segmented_lines_list = ref_addtostart(segmented_lines_list,
266 scene_hist_ref->line_seg, LINE2);
267
268 /* Free the generated match list */
269
270 list_rm(cut_match_list, rfree);
271 }
272
273 }/*endfor(ptr)*/
274 return(segmented_lines_list);
275 }
276
277 /* ---------- */
278 /* This function identifies any known models in the mono scene - the hough
279 transforms constructed for 'model_name' are placed onto the image stack */
280
281 Point2 *locate_models_old(char *model_name, List **new_peak_list, Vec2 *shift)
282 {
283 Imrect *im_trans(Imrect *im, double angle, Vec2 *im_center, Vec2 *offset);
284 List *model_ptr;
285 List *peak_list, *peak_ptr;
286 double thres;
287 Imrect *location_hough_space;
288 Imrect *orient_hough;
289 #ifdef _PGHDEBUG
290 Imrect *trans_mono_im; /* HR: 10/05/2013: added ifdef */
291 #endif
292 Model_poly_header *polyh;
293 int peak_count;
294 List *ptr;
295 List *scene_pairs_list;
296 List *matched_models_list;
297 List *model_poly_list;
298 Pairs_hough_def *hough_def;
299 Match_cut_def *match_cut_def;
300 int max_x,max_y;
301 Hough1_peak orient_peak;
302 List *geom;
303 Vec2 origin;
304 Vec2 offset;
305 Point2 *peak_pos;
306 char mess[64];
307 Transform2 poly_trans;
308 double scale_factor;
309 Mat2 *C, rotation_matrix;
310 double theta, lam1, lam2, major, minor;
311
312
313 vec2_x(origin)=(float)0.0;
314 vec2_y(origin)=(float)0.0;
315
316
317 /* (1) For the chosen model in the database, construct a Hough space,
318 and vote for the position using the matched scene histograms.
319 (2) Search for significant peaks which
320 hypothesize the presence of a model in the scene.ref
321 (3) For the chosen peak peak determine the orientation of the model
322 at that location by voting for the angle between each
323 scene line and the model line it has matched to.
324 */
325
326 scene_pairs_list = pairs_scene_pairs_list_get();
327 matched_models_list = pairs_matched_models_list_get();
328 model_poly_list = pairs_model_poly_list_get();
329 hough_def = pairs_hough_def_get();
330 match_cut_def = pairs_match_cut_def_get();
331
332 if (scene_pairs_list==NULL) return(NULL);
333
334 /* Free the matched_models list */
335
336 if(matched_models_list!=NULL) list_rm(matched_models_list, geom_free);
337 matched_models_list=NULL;
338 pairs_matched_models_list_set(matched_models_list);
339
340 for(model_ptr=model_poly_list;model_ptr!=NULL;model_ptr=model_ptr->next)
341 {
342 polyh = model_ptr->to;
343 if (strcmp(model_name, polyh->name)==0) break;
344 }
345 if (model_ptr==NULL) return(NULL);
346
347
348 location_hough_space = hough2_alloc(hough_def->loc_hough2_min_x,
349 hough_def->loc_hough2_min_y,
350 hough_def->loc_hough2_max_x,
351 hough_def->loc_hough2_max_y,
352 hough_def->loc_hough2_binsize_x,
353 hough_def->loc_hough2_binsize_y,
354 double_v);
355
356 /*TEMP*/
357 hough2_fill(location_hough_space, 0.0);
358
359 if (shift !=NULL)
360 polyh->landmark = *shift;
361 else
362 polyh->landmark = polyh->centroid;
363
364 if (hough_def->plot_type==LONG_LINES)
365 {
366 locate_using_long_lines(polyh, location_hough_space,
367 scene_pairs_list, match_cut_def);
368 }
369
370 if (hough_def->plot_type==SHORT_LINES)
371 {
372 locate_using_short_lines(polyh, location_hough_space,
373 scene_pairs_list, match_cut_def);
374 }
375
376 if (hough_def->plot_type==POINTS)
377 {
378 locate_using_points(polyh,
379 location_hough_space, scene_pairs_list, match_cut_def);
380 }
381
382 /* Sqrt the hough transform to approximate a chi_2 */
383
384 sqrt_hough2(location_hough_space);
385
386 /* Search the loc hough2 for peaks */
387
388 thres = imf_max(location_hough_space, &max_y, &max_x)*
389 (hough_def->peak_percent/100.0);
390 peak_list = hough2_locate_peaks(location_hough_space, thres);
391
392 /* file must be compiled as a project with tina-tools library in order to debug */
393 #ifdef _PGHDEBUG
394 stack_push(location_hough_space, IMRECT, im_free);
395 #else
396 im_free(location_hough_space);
397 #endif
398
399 /* Determine model orientation at hough2 peak_no */
400
401 peak_list = sort_list(peak_list, peak_func, NULL);
402
403 for(peak_ptr=peak_list,peak_count=0;peak_ptr!=NULL;peak_count++,peak_ptr=peak_ptr->next)
404 {
405 if(peak_count==hough_def->peak_no) break;
406 }
407
408 if (peak_ptr==NULL) return(NULL);
409
410 sprintf(mess,"hough peak = %4.2f / %4.2f \n",-peak_func(peak_ptr->to),thres);
411 format(mess);
412
413
414 /* Display the error of this peak */
415 C = prop_get(((Point2*)peak_ptr->to)->props, COVAR2);
416 conic_eigen(mat2_xx(*C), mat2_xy(*C), mat2_yy(*C),
417 &theta, &lam1, &lam2);
418 theta = 180.0*theta/PI;
419 /* Assumes binsizex = binsizey */
420 major = sqrt(lam1)*hough_def->loc_hough2_binsize_x;
421 minor = sqrt(lam2)*hough_def->loc_hough2_binsize_x;
422 format("major: %f\nminor: = %f\nangle = %f\n\n", major, minor, theta);
423
424 /* Determine the scale of the model responsible for this peak */
425
426 scale_factor = scale_from_loc(peak_ptr->to, model_ptr->to);
427 printf("Scale:%f\n", scale_factor);
428
429 pgh_scale_estimate = scale_factor; /* Stored for global access */
430
431 orient_hough = hough1_alloc(-PI-(PI/hough_def->orient_binsize),
432 PI+(PI/hough_def->orient_binsize),
433 (PI/hough_def->orient_binsize), float_v);
434
435 angle_from_loc(orient_hough, peak_ptr->to, model_ptr->to);
436 thres = imf_max(orient_hough, &max_y, &max_x)*
437 (hough_def->peak_percent/100.0);
438 orient_peak = hough1_locate_max_peak(orient_hough, thres ,float_v);
439
440 /* Make a copy of the model data */
441
442 polyh=model_ptr->to;
443 geom = poly_copy(polyh->model_poly_data);
444
445 /* Transform to landmark or centroid */
446 rotation_matrix = rot2(0.0);
447 poly_trans.R = rotation_matrix;
448 poly_trans.t = polyh->landmark;
449
450 for(ptr=geom;ptr!=NULL;ptr=ptr->next)
451 line2_transform(ptr->to, poly_trans);
452
453 /* Transform the copied data to the determined angle, position and scale*/
454
455 poly_rotate_and_scale(geom, orient_peak.x, scale_factor, origin);
456 peak_pos = peak_ptr->to;
457 poly_translate(geom, peak_pos->p );
458
459 /* Add the transformed copy to the matched_models list */
460
461 matched_models_list = list_append(matched_models_list, geom);
462
463 /* file must be compiled as a project with tina-tools library in order to debug */
464
465 origin = peak_pos->p;
466 offset = vec2_sum(poly_trans.t,origin);
467
468 #ifdef _PGHDEBUG
469 stack_push(orient_hough, IMRECT, im_free);
470
471 trans_mono_im = im_trans(mono_image_get(), orient_peak.x, &origin, &offset);
472 stack_push(trans_mono_im, IMRECT, im_free);
473 #else
474 im_free(orient_hough);
475 #endif
476
477 pairs_matched_models_list_set(matched_models_list);
478
479 if(new_peak_list != NULL) *new_peak_list = peak_list;
480 return(peak_pos);
481 }
482
483 /* ---------- */
484
485 Point2 *locate_models(char *model_name, List **new_peak_list, Vec2 *shift, double *rot_scale, Vec2 *im_origin, Vec2 *im_offset) /* HR: extra output parameters added */
486 {
487 List *model_ptr;
488 List *peak_list, *peak_ptr;
489 double thres;
490 Imrect *location_hough_space;
491 Imrect *orient_hough;
492 Model_poly_header *polyh;
493 int peak_count;
494 List *ptr;
495 List *scene_pairs_list;
496 List *matched_models_list;
497 List *model_poly_list;
498 Pairs_hough_def *hough_def;
499 Match_cut_def *match_cut_def;
500 int max_x,max_y;
501 Hough1_peak orient_peak;
502 List *geom;
503 Vec2 origin;
504 Vec2 offset;
505 Point2 *peak_pos;
506 char mess[64];
507 Transform2 poly_trans;
508 double scale_factor;
509 Mat2 *C, rotation_matrix;
510 double theta, lam1, lam2, major, minor;
511
512
513 vec2_x(origin)=(float)0.0;
514 vec2_y(origin)=(float)0.0;
515
516
517 /* (1) For the chosen model in the database, construct a Hough space,
518 and vote for the position using the matched scene histograms.
519 (2) Search for significant peaks which
520 hypothesize the presence of a model in the scene.ref
521 (3) For the chosen peak peak determine the orientation of the model
522 at that location by voting for the angle between each
523 scene line and the model line it has matched to.
524 */
525
526 scene_pairs_list = pairs_scene_pairs_list_get();
527 matched_models_list = pairs_matched_models_list_get();
528 model_poly_list = pairs_model_poly_list_get();
529 hough_def = pairs_hough_def_get();
530 match_cut_def = pairs_match_cut_def_get();
531
532 if (scene_pairs_list==NULL) return(NULL);
533
534 /* Free the matched_models list */
535
536 if(matched_models_list!=NULL) list_rm(matched_models_list, geom_free);
537 matched_models_list=NULL;
538 pairs_matched_models_list_set(matched_models_list);
539
540 for(model_ptr=model_poly_list;model_ptr!=NULL;model_ptr=model_ptr->next)
541 {
542 polyh = model_ptr->to;
543 if (strcmp(model_name, polyh->name)==0) break;
544 }
545 if (model_ptr==NULL) return(NULL);
546
547
548 location_hough_space = hough2_alloc(hough_def->loc_hough2_min_x,
549 hough_def->loc_hough2_min_y,
550 hough_def->loc_hough2_max_x,
551 hough_def->loc_hough2_max_y,
552 hough_def->loc_hough2_binsize_x,
553 hough_def->loc_hough2_binsize_y,
554 double_v);
555
556 /*TEMP*/
557 hough2_fill(location_hough_space, 0.0);
558
559 if (shift !=NULL)
560 polyh->landmark = *shift;
561 else
562 polyh->landmark = polyh->centroid;
563
564 if (hough_def->plot_type==LONG_LINES)
565 {
566 locate_using_long_lines(polyh, location_hough_space,
567 scene_pairs_list, match_cut_def);
568 }
569
570 if (hough_def->plot_type==SHORT_LINES)
571 {
572 locate_using_short_lines(polyh, location_hough_space,
573 scene_pairs_list, match_cut_def);
574 }
575
576 if (hough_def->plot_type==POINTS)
577 {
578 locate_using_points(polyh, location_hough_space,
579 scene_pairs_list, match_cut_def);
580 }
581
582 /* Sqrt the hough transform to approximate a chi_2 */
583
584 sqrt_hough2(location_hough_space);
585
586 /* Search the loc hough2 for peaks */
587
588 thres = imf_max(location_hough_space, &max_y, &max_x)*
589 (hough_def->peak_percent/100.0);
590 peak_list = hough2_locate_peaks(location_hough_space, thres);
591
592 /* file must be compiled as a project with tina-tools library in order to debug */
593
594 #ifdef _PGHDEBUG
595 stack_push(location_hough_space, IMRECT, im_free);
596 #else
597 im_free(location_hough_space);
598 #endif
599
600 /* Determine model orientation at hough2 peak_no */
601
602 peak_list = sort_list(peak_list, peak_func, NULL);
603
604 for(peak_ptr=peak_list,peak_count=0;peak_ptr!=NULL;peak_count++,peak_ptr=peak_ptr->next)
605 {
606 if(peak_count==hough_def->peak_no) break;
607 }
608
609 if (peak_ptr==NULL) return(NULL);
610
611 sprintf(mess,"hough peak = %4.2f / %4.2f \n",-peak_func(peak_ptr->to),thres);
612 format(mess);
613
614
615 /* Display the error of this peak */
616 C = prop_get(((Point2*)peak_ptr->to)->props, COVAR2);
617 conic_eigen(mat2_xx(*C), mat2_xy(*C), mat2_yy(*C),
618 &theta, &lam1, &lam2);
619 theta = 180.0*theta/PI;
620 /* Assumes binsizex = binsizey */
621 major = sqrt(lam1)*hough_def->loc_hough2_binsize_x;
622 minor = sqrt(lam2)*hough_def->loc_hough2_binsize_x;
623 format("major: %f\nminor: = %f\nangle = %f\n\n", major, minor, theta);
624
625 /* Determine the scale of the model responsible for this peak */
626
627 scale_factor = scale_from_loc(peak_ptr->to, model_ptr->to); /* HR: check! may return 0.0 as scale! */
628 printf("Scale:%f\n", scale_factor);
629
630 pgh_scale_estimate = scale_factor; /* Stored for global access */
631
632 orient_hough = hough1_alloc(-PI-(PI/hough_def->orient_binsize),
633 PI+(PI/hough_def->orient_binsize),
634 (PI/hough_def->orient_binsize), float_v);
635
636 angle_from_loc(orient_hough, peak_ptr->to, model_ptr->to);
637 thres = imf_max(orient_hough, &max_y, &max_x)*
638 (hough_def->peak_percent/100.0);
639 orient_peak = hough1_locate_max_peak(orient_hough, thres ,float_v);
640
641 /* Make a copy of the model data */
642
643 polyh=model_ptr->to;
644 geom = poly_copy(polyh->model_poly_data);
645
646 /* Transform to landmark or centroid */
647 rotation_matrix = rot2(0.0);
648 poly_trans.R = rotation_matrix;
649 poly_trans.t = polyh->landmark;
650
651 for(ptr=geom;ptr!=NULL;ptr=ptr->next)
652 line2_transform(ptr->to, poly_trans);
653
654 /* Transform the copied data to the determined angle, position and scale*/
655
656 poly_rotate_and_scale(geom, orient_peak.x, scale_factor, origin);
657 peak_pos = peak_ptr->to;
658 poly_translate(geom, peak_pos->p );
659
660 /* Add the transformed copy to the matched_models list */
661
662 matched_models_list = list_append(matched_models_list, geom);
663
664 /* file must be compiled as a project with tina-tools library in order to debug */
665
666 #ifdef _PGHDEBUG
667 stack_push(orient_hough, IMRECT, im_free);
668 #else
669 im_free(orient_hough);
670 #endif
671 /* HR: transformation moved from here */
672 origin = peak_pos->p;
673 offset = vec2_sum(poly_trans.t,origin);
674 /* HR: extra parameters added to pass as arguments */
675 rot_scale[0] = orient_peak.x;
676 rot_scale[1] = scale_factor;
677
678 vec2_x(*im_offset) = vec2_x(offset);
679 vec2_y(*im_offset) = vec2_y(offset);
680
681 vec2_x(*im_origin) = vec2_x(origin);
682 vec2_y(*im_origin) = vec2_y(origin);
683
684 pairs_matched_models_list_set(matched_models_list);
685
686 if(new_peak_list != NULL) *new_peak_list = peak_list;
687 return(peak_pos);
688 }
689
690 static Vec2 array_get_quadmaxf(double **corr_arr , float x, float y, float *px, float *py)
691 {
692 double pixval[3][3], temp;
693 float a, b, c, d, e, f, xs, ys;
694 int i, j, g, h, p, q;
695 Vec2 pxy;
696
697 i = tina_int(x - 1);
698 j = tina_int(y - 1);
699
700 for (g=0; g<3; g++)
701 {
702 for (h=0; h<3; h++)
703 {
704 p = i + h;
705 q = j + g;
706 pixval[g][h] = corr_arr[q][p];
707 }
708 }
709
710 a = (float) pixval[1][1];
711 b = (float) ((pixval[0][2] - pixval[0][0]
712 + pixval[1][2] - pixval[1][0] + pixval[2][2] - pixval[2][0]) / 6.0);
713 c = (float) ((pixval[2][0] - pixval[0][0]
714 + pixval[2][1] - pixval[0][1] + pixval[2][2] - pixval[0][2]) / 6.0);
715 d = (float) ((pixval[0][0] - 2.0 * pixval[0][1] + pixval[0][2]
716 + 3.0 * pixval[1][0] - 6.0 * pixval[1][1] + 3.0 * pixval[1][2]
717 + pixval[2][0] - 2.0 * pixval[2][1] + pixval[2][2]) / 10.0);
718 e = (float) ((pixval[0][0] - pixval[2][0] + pixval[2][2] - pixval[0][2]) / 4.0);
719 f = (float) ((pixval[0][0] + 3.0 * pixval[0][1] + pixval[0][2]
720 - 2.0 * pixval[1][0] - 6.0 * pixval[1][1] - 2.0 * pixval[1][2]
721 + pixval[2][0] + 3.0 * pixval[2][1] + pixval[2][2]) / 10.0);
722
723 temp = 4.0 * d * f - e * e;
724 *px = 1.5 + (double) i + (e * c - 2 * f * b) / temp;
725 *py = 1.5 + (double) j + (e * b - 2 * d * c) / temp;
726 /*xs = (float) ((e * c - 2 * f * b) / temp);
727 ys = (float) ((e * b - 2 * d * c) / temp);
728 *px = x + xs;
729 *py = y + ys;*/
730
731 vec2_x(pxy)=*px;
732 vec2_y(pxy)=*py;
733
734 return (pxy);
735 }
736
737
738 Vec2 corr_code(Imrect *im1, Imrect *im2, Vec2 *peak_p, int window, double range, double sigma, double thresh)
739 {
740 Imrect *imf_gauss_diff(Imrect *im, double sigma, double precision, double factor);
741 Vec2 array_get_quadmaxf(double **corr_arr, float x, float y, float *px, float *py);
742 double **corr_arr, *corr, precision=exp(-5.0);
743 float px, py, data, *row1, *row2, max_corr=0.0;
744 int lx, ux, ly, uy, lx1, ux1, ly1, uy1, pos_x, pos_y;
745 int i, j, m, n, m_max, n_max;
746 Vec2 cv;
747 Imrect *sig_im1 = imf_gauss_diff(im1, sigma, precision, 2.0); /* HR: amended */
748 Imrect *sig_im2 = imf_gauss_diff(im2, sigma, precision, 2.0);
749
750 pos_x = vec2_x(*peak_p);
751 pos_y = vec2_y(*peak_p);
752
753 lx = pos_x - range - window;
754 ly = pos_y - range - window;
755 ux = pos_x + range + window;
756 uy = pos_y + range + window;
757
758 row1 = fvector_alloc(lx, ux);
759 row2 = fvector_alloc(lx, ux);
760
761 corr_arr = darray_alloc((pos_y-range),(pos_x-range),
762 (pos_y+range+1),(pos_x+range+1));
763
764 for (m = - range; m < + range; ++m)
765 {
766 corr= corr_arr[pos_y+m];
767 for (n = - range; n < + range; ++n)
768 {
769 corr[pos_x+n] = 0.0;
770 }
771
772 for (i = pos_y - window; i < pos_y + window; ++i)
773 {
774 im_get_rowf(row1, sig_im1, i+m, lx, ux);
775 im_get_rowf(row2, sig_im2, i, lx, ux);
776
777 for (j = pos_x - window; j < pos_x + window; ++j)
778 {
779 data = row2[j];
780 for (n = - range; n < + range-5; ++n)
781 {
782 corr[pos_x+n] += row1[j+n]*data;
783 n++;
784 corr[pos_x+n] += row1[j+n]*data;
785 n++;
786 corr[pos_x+n] += row1[j+n]*data;
787 n++;
788 corr[pos_x+n] += row1[j+n]*data;
789 n++;
790 corr[pos_x+n] += row1[j+n]*data;
791 }
792 for (; n < range; n++)
793 {
794 corr[pos_x+n] += row1[j+n]*data;
795 }
796 }
797 }
798 }
799
800 for (m = -range+1; m < range-1; ++m)
801 {
802 for (n = -range+1; n < range-1; ++n)
803 {
804 data = corr_arr[pos_y+m][pos_x+n];
805 if (data >= max_corr)
806 {
807 max_corr = data;
808 m_max = m;
809 n_max = n;
810 }
811 }
812 }
813
814 array_get_quadmaxf(corr_arr, pos_x+n_max, pos_y+m_max, &px, &py);
815
816 /* HR: check whether the refined point is inside the range and flag it */
817
818 lx1 = pos_x - range;
819 ly1 = pos_y - range;
820 ux1 = pos_x + range;
821 uy1 = pos_y + range;
822
823 if( (px< lx1) || (px> ux1) || (py< ly1) || (py> uy1) || (fabs(px-pos_x)> thresh) || (fabs(py-pos_y)> thresh) )
824 {
825 format("\n location refinement failed! \n");
826
827 vec2_y(cv) = -1.0;
828 vec2_x(cv) = -1.0;
829
830 if( fabs(px-pos_x)> thresh )
831 {
832 format("shift %5.2f needed in x direction > max shift %5.2f allowed!\n", fabs(px-pos_x), thresh);
833 }
834 if( fabs(py-pos_y)> thresh )
835 {
836 format("shift %5.2f needed in y direction > max shift %5.2f allowed!\n", fabs(py-pos_y), thresh);
837 }
838 }
839 else
840 {
841 vec2_y(cv) = py;
842 vec2_x(cv) = px;
843 }
844
845 fvector_free(row1,lx);
846 fvector_free(row2,lx);
847 darray_free(corr_arr, (pos_y-range), (pos_x-range),
848 (pos_y+range+1), (pos_x+range+1) );
849 return(cv);
850 }
851
852
This page was automatically generated by the
LXR engine.
Visit the LXR main site for more
information.