~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/vision/visPgh_HTcovar.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  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 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 General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU 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  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  *
 25  * Program :    TINA
 26  * File    :  $Source: /home/tina/cvs/tina-libs/tina/vision/visPgh_HTcovar.c,v $
 27  * Date    :  $Date: 2007/02/15 01:52:29 $
 28  * Version :  $Revision: 1.2 $
 29  * CVS Id  :  $Id: visPgh_HTcovar.c,v 1.2 2007/02/15 01:52:29 paul Exp $
 30  *
 31  * Author  :  A.Ashbrook, NAT.
 32  */
 33 /**
 34  * @file Routines for computing the covariance on the location of a rigid object from
 35  *       correspondences between line fragments established using PGH (or other methods).
 36  * @brief Assumes that the quadratic shape of the peak in the location hough transform 
 37  *        is a good description of the true likelihood.
 38  *
 39 */
 40 
 41 #include "visPgh_HTcovar.h"
 42 
 43 #include <math.h>
 44 #include <float.h>
 45 
 46 
 47 #undef MAXFLOAT
 48 #define MAXFLOAT FLT_MAX
 49 #define MINFLOAT FLT_MIN 
 50 
 51 
 52 
 53 /* ---------- Code from hough1_float.c APA 4/2/94 ---------- */
 54 
 55 /* ---------- Code from covariance.c APA 4/2/94 ---------- */
 56 
 57 void  get_covariance(Line2 *ue_lin1_s, Line2 *ue_lin2_s,
 58                      double *a, double *b, double *c, double *d)
 59 {
 60      double  matrix[8][2], matrix_T[2][8], inter_covar_matrix[2][2];
 61      double  det;
 62      float   ax1, ay1, ax2, ay2, bx1, by1, bx2, by2;
 63 
 64 /* Does a chi squared test by propagating the measurement errors through 
 65    if the chi squared value is high then the centroid is probably near the
 66    previously calculated centroid */
 67 
 68 /* Get data about the position of the two lines */
 69 
 70      ax1 = vec2_x(ue_lin1_s->p1);
 71      ay1 = vec2_y(ue_lin1_s->p1);
 72      ax2 = vec2_x(ue_lin1_s->p2);
 73      ay2 = vec2_y(ue_lin1_s->p2);
 74      bx1 = vec2_x(ue_lin2_s->p1);
 75      by1 = vec2_y(ue_lin2_s->p1);
 76      bx2 = vec2_x(ue_lin2_s->p2);
 77      by2 = vec2_y(ue_lin2_s->p2);
 78 
 79 /* First find the covariance matrix for the position of the centroid, then
 80    find the transform */
 81 
 82    partial_derv_covariance(ax1, ay1, ax2, ay2, bx1, by1, bx2, by2, &matrix[0][0], &matrix[0][1]);
 83    partial_derv_covariance(ax2, ay2, ax1, ay1, bx1, by1, bx2, by2, &matrix[2][0], &matrix[2][1]);
 84    partial_derv_covariance(ay1, ax1, ay2, ax2, by1, bx1, by2, bx2, &matrix[1][1], &matrix[1][0]);
 85    partial_derv_covariance(ay2, ax2, ay1, ax1, by1, bx1, by2, bx2, &matrix[3][1], &matrix[3][0]);
 86    partial_derv_covariance(bx1, by1, bx2, by2, ax1, ay1, ax2, ay2, &matrix[4][0], &matrix[4][1]);
 87    partial_derv_covariance(bx2, by2, bx1, by1, ax1, ay1, ax2, ay2, &matrix[6][0], &matrix[6][1]);
 88    partial_derv_covariance(by1, bx1, by2, bx2, ay1, ax1, ay2, ax2, &matrix[5][1], &matrix[5][0]);
 89    partial_derv_covariance(by2, bx2, by1, bx1, ay1, ax1, ay2, ax2, &matrix[7][1], &matrix[7][0]);
 90 
 91 /* Transpose the matrix */
 92 
 93      matrix_T[0][0] = matrix[0][0];
 94      matrix_T[0][1] = matrix[1][0];
 95      matrix_T[0][1] = matrix[1][0];
 96      matrix_T[0][2] = matrix[2][0];
 97      matrix_T[0][3] = matrix[3][0];
 98      matrix_T[0][4] = matrix[4][0];
 99      matrix_T[0][5] = matrix[5][0];
100      matrix_T[0][6] = matrix[6][0];
101      matrix_T[0][7] = matrix[7][0];
102    
103      matrix_T[1][0] = matrix[0][1];
104      matrix_T[1][1] = matrix[1][1];
105      matrix_T[1][2] = matrix[2][1];
106      matrix_T[1][3] = matrix[3][1];
107      matrix_T[1][4] = matrix[4][1];
108      matrix_T[1][5] = matrix[5][1];
109      matrix_T[1][6] = matrix[6][1];
110      matrix_T[1][7] = matrix[7][1];
111    
112 /* Find the propagated covariance matrix for the intersection point */
113 
114      inter_covar_matrix[0][0] = (matrix_T[0][0]*matrix[0][0]) + (matrix_T[0][1]*matrix[1][0]);
115      inter_covar_matrix[0][0] = inter_covar_matrix[0][0] + (matrix_T[0][2]*matrix[2][0]) 
116                                 + (matrix_T[0][3]*matrix[3][0]);
117      inter_covar_matrix[0][0] = inter_covar_matrix[0][0] + (matrix_T[0][4]*matrix[4][0]) 
118                                 + (matrix_T[0][5]*matrix[5][0]);
119      inter_covar_matrix[0][0] = inter_covar_matrix[0][0] + (matrix_T[0][6]*matrix[6][0]) 
120                                 + (matrix_T[0][7]*matrix[7][0]);
121 
122      inter_covar_matrix[0][1] = (matrix_T[0][0]*matrix[0][1]) + (matrix_T[0][1]*matrix[1][1]);
123      inter_covar_matrix[0][1] = inter_covar_matrix[0][1] + (matrix_T[0][2]*matrix[2][1]) 
124                                 + (matrix_T[0][3]*matrix[3][1]);
125      inter_covar_matrix[0][1] = inter_covar_matrix[0][1] + (matrix_T[0][4]*matrix[4][1]) 
126                                 + (matrix_T[0][5]*matrix[5][1]);
127      inter_covar_matrix[0][1] = inter_covar_matrix[0][1] + (matrix_T[0][6]*matrix[6][1]) 
128                                 + (matrix_T[0][7]*matrix[7][1]);
129 
130      inter_covar_matrix[1][0] = (matrix_T[1][0]*matrix[0][0]) + (matrix_T[1][1]*matrix[1][0]);
131      inter_covar_matrix[1][0] = inter_covar_matrix[1][0] + (matrix_T[1][2]*matrix[2][0]) 
132                                 + (matrix_T[1][3]*matrix[3][0]);
133      inter_covar_matrix[1][0] = inter_covar_matrix[1][0] + (matrix_T[1][4]*matrix[4][0])
134                                 + (matrix_T[1][5]*matrix[5][0]);
135      inter_covar_matrix[1][0] = inter_covar_matrix[1][0] + (matrix_T[1][6]*matrix[6][0]) 
136                                 + (matrix_T[1][7]*matrix[7][0]);
137 
138      inter_covar_matrix[1][1] = (matrix_T[1][0]*matrix[0][1]) + (matrix_T[1][1]*matrix[1][1]);
139      inter_covar_matrix[1][1] = inter_covar_matrix[1][1]+ (matrix_T[1][2]*matrix[2][1]) 
140                                 + (matrix_T[1][3]*matrix[3][1]);     inter_covar_matrix[1][1] = inter_covar_matrix[1][1] +
141 (matrix_T[1][4]*matrix[4][1]) 
142                                 + (matrix_T[1][5]*matrix[5][1]);
143      inter_covar_matrix[1][1] = inter_covar_matrix[1][1] + (matrix_T[1][6]*matrix[6][1])
144                                 + (matrix_T[1][7]*matrix[7][1]);
145 
146     *a = inter_covar_matrix[0][0];
147     *b = inter_covar_matrix[0][1];
148     *c = inter_covar_matrix[1][0];
149     *d = inter_covar_matrix[1][1];
150 
151      det = ((*a) * (*d)) - ((*c) * (*b));
152 
153      inter_covar_matrix[0][0] =  (*d) / det;
154      inter_covar_matrix[0][1] = -((*b) / det);
155      inter_covar_matrix[1][0] = -((*c) / det);
156      inter_covar_matrix[1][1] =  (*a) / det;
157 
158     *a = inter_covar_matrix[0][0];
159     *b = inter_covar_matrix[0][1];
160     *c = inter_covar_matrix[1][0];
161     *d = inter_covar_matrix[1][1];
162 }
163 
164 /* ---------- */
165 
166 void partial_derv_covariance(float ax1, float ay1, float ax2, float ay2,
167                              float bx1, float by1, float bx2, float by2,
168                              double *matrix_1, double *matrix_2)
169 {
170      double v, u, derv;
171 
172 /* Compute the values of u and v for differentiation of a quotient */
173 
174      v = (((ay2 - ay1)*(bx2 - bx1)) - ((by2 - by1)*(ax2 - ax1)));
175      u = (((by1*bx2) - (bx1 * by2))*(ax2 - ax1)) - (((ay1*ax2)-(ax1*ay2))*(bx2-bx1));
176 
177 /* Computes the parital derviative with respect to ax1 */
178 
179      derv = (((-(by1 * bx2) + (bx1*by2) + (ay2*bx2) - (ay2*bx1))*v)-(u*(by2-by1)))/(v*v);
180     
181     *matrix_1 = derv;
182      if (fabs(bx2 - bx1)<MINFLOAT ) *matrix_2 = MAXFLOAT;
183      else
184          *matrix_2 = ((by2 - by1) / (bx2 - bx1)) * derv; 
185 }
186 
187 /* ---------- */
188 
189 double  chi_squared_test(Vec2 *int_pt, Point2 *loc, double *a, double *b,
190                          double *c, double *d)
191 {
192      double   delta_x, delta_y, chi_sq;
193 
194      delta_x = vec2_x(loc->p) - vec2_x(*int_pt);
195      delta_y = vec2_y(loc->p) - vec2_y(*int_pt);
196      chi_sq = ((*a) * delta_x * delta_x) + ((*c) * delta_x * delta_y) + 
197               ((*b) * delta_x * delta_y) + ((*d) * delta_y * delta_y);
198 
199 return(chi_sq);
200 }
201 
202 /* ---------- */
203 
204 
205 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.