transform.c

Go to the documentation of this file.
00001 /* @(#)transform.c      2.1  6/26/87 */
00002 /****************************************************************
00003 This  file  contains    routines    which    perform    (affine?)
00004 transformations  from  one  coordinate  system  into another. The
00005 second system may be translated, stretched, and rotated  relative
00006 to  the  first.   The  input  system is system "a" and the output
00007 system is "b"
00008 
00009 note: uses sqrt() from math library
00010 *****************************************************************
00011 compute_transformation_coef (ax,ay,bx,by,use,n)
00012 
00013     double ax[], ay[];       coordinate from system a
00014     double bx[], by[];       coordinate from system b
00015     char   use[];            use point flags
00016     int n;                   number of points in ax,ay,bx,by
00017 
00018 The first step is to compute coefficients for a set of  equations
00019 which  are then used to convert from the one system to the other.
00020 A set of x,y points from both systems is input into the  equation
00021 generator  which  determines the equation coefficients which most
00022 nearly represent the original points. These coefficients are kept
00023 in a static variables internal to this file.
00024 
00025 note: use[i] must be true for ax[i],ay[i],bx[i],by[i] to be used
00026       in the equation
00027 
00028       also, the total number of used points must be 4 or larger
00029 
00030 returns:
00031    -2  less than 4 used points.
00032    -1  couldn't solve the equation. points probably colinear
00033        and need to be spread out more.
00034     1  ok
00035 *****************************************************************
00036 Then points from one system may be converted into the  second  by
00037 use of one of the two equation routines.
00038 
00039 transform_a_into_b (ax,ay,bx,by)
00040 
00041     double ax,ay;            input point from system a
00042     double *bx,*by;          resultant point in system b
00043 
00044 transform_b_into_a (bx,by,ax,ay)
00045 
00046     double bx,by;            input point from system b
00047     double *ax,*ay;          resultant point in system a
00048 *****************************************************************
00049 Residual analysis on the equation can be run to test how well
00050 the equations work.  Either test how well b is predicted by a
00051 or vice versa.
00052 
00053 residuals_a_predicts_b (ax,ay,bx,by,use,n,residuals,rms)
00054 residuals_b_predicts_a (ax,ay,bx,by,use,n,residuals,rms)
00055 
00056     double ax[], ay[];       coordinate from system a
00057     double bx[], by[];       coordinate from system b
00058     char   use[];            use point flags
00059     int n;                   number of points in ax,ay,bx,by
00060     double residual[]        residual error for each point
00061     double *rms;             overall root mean square error
00062 ****************************************************************/
00063 
00064 #include <math.h>
00065 #include <grass/libtrans.h>
00066 
00067 /* the coefficients */
00068 static double A0,A1,A2,A3,A4,A5;
00069 static double B0,B1,B2,B3,B4,B5;
00070 
00071 static int resid(
00072     double *,double *,double *,double *,int *,int,double *,double *,int);
00073 
00074 int compute_transformation_coef(
00075     double ax[],double ay[],double bx[],double by[],
00076     int use[],int n)
00077 {
00078     int i;
00079     int j;
00080     int count;
00081     double aa[3];
00082     double aar[3];
00083     double bb[3];
00084     double bbr[3];
00085 
00086     double cc[3][3];
00087     double x;
00088 
00089     count = 0;
00090     for (i = 0; i < n; i++)
00091         if (use[i])
00092             count++;
00093     if (count < 4)
00094         return -2; /* must have at least 4 points */
00095 
00096     for (i = 0; i < 3; i++)
00097     {
00098         aa[i] = bb[i] = 0.0;
00099 
00100         for (j = 0; j < 3; j++)
00101             cc[i][j] = 0.0;
00102     }
00103 
00104     for (i = 0; i < n; i++)
00105     {
00106         if (!use[i]) continue;  /* skip this point */
00107         cc[0][0] += 1;
00108         cc[0][1] += bx[i];
00109         cc[0][2] += by[i];
00110 
00111         cc[1][1] += bx[i] * bx[i];
00112         cc[1][2] += bx[i] * by[i];
00113         cc[2][2] += by[i] * by[i];
00114 
00115         aa[0] += ay[i];
00116         aa[1] += ay[i] * bx[i];
00117         aa[2] += ay[i] * by[i];
00118 
00119         bb[0] += ax[i];
00120         bb[1] += ax[i] * bx[i];
00121         bb[2] += ax[i] * by[i];
00122     }
00123 
00124     cc[1][0] = cc[0][1];
00125     cc[2][0] = cc[0][2];
00126     cc[2][1] = cc[1][2];
00127 
00128 /* aa and bb are solved */
00129         
00130         if ( inverse (cc) < 0)
00131                 return (-1) ;
00132         if ( m_mult ( cc, aa, aar) < 0  ||  m_mult ( cc, bb, bbr) < 0)
00133                 return (-1) ;
00134 
00135 
00136 
00137 /* the equation coefficients */
00138 
00139     B0 = aar[0];
00140     B1 = aar[1];
00141     B2 = aar[2];
00142 
00143     B3 = bbr[0];
00144     B4 = bbr[1];
00145     B5 = bbr[2];
00146 
00147 /* the inverse equation */
00148 
00149     x = B2 * B4 - B1 * B5 ;
00150 
00151     if( ! x)
00152         return (-1) ;
00153 
00154     A0 = (B1 * B3 - B0 * B4) / x ;
00155     A1 = -B1 / x ;
00156     A2 =  B4 / x ;
00157     A3 = (B0 * B5 - B2 * B3) / x ;
00158     A4 = B2 / x ;
00159     A5 = -B5 / x ;
00160 
00161     return 1;
00162 }
00163 
00164 int transform_a_into_b(
00165     double ax,double ay,
00166     double *bx,double *by)
00167 {
00168     *by = A0 + A1 * ax + A2 * ay ;
00169     *bx = A3 + A4 * ax + A5 * ay ;
00170 
00171     return 0;
00172 }
00173 
00174 int transform_b_into_a(
00175     double bx,double by,
00176     double *ax,double *ay)
00177 {
00178     *ay = B0  + B1 * bx + B2 * by ;
00179     *ax = B3  + B4 * bx + B5 * by ;
00180 
00181     return 0;
00182 }
00183 
00184 /**************************************************************
00185 These routines are internal to this source code
00186 
00187 solve (a, b)
00188     double a[3][3]
00189     double b[3]
00190 
00191     equation solver used by compute_transformation_coef()
00192 **************************************************************/
00193 
00194 /*  #define abs(xx) (xx >= 0 ? xx : -xx)  */
00195 /*      #define N 3  */
00196 
00197 
00198 int residuals_a_predicts_b (
00199     double ax[],double ay[],double bx[],double by[],
00200     int use[],
00201     int n,
00202     double residuals[],
00203     double *rms)
00204 {
00205     resid (ax,ay,bx,by,use,n,residuals,rms,1);
00206 
00207     return 0;
00208 }
00209 
00210 int residuals_b_predicts_a(
00211     double ax[],double ay[],double bx[],double by[],
00212     int use[],
00213     int n,
00214     double residuals[],
00215     double *rms)
00216 {
00217     resid (ax,ay,bx,by,use,n,residuals,rms,0);
00218 
00219     return 0;
00220 }
00221 
00222 static int resid(
00223     double ax[],double ay[],double bx[],double by[],
00224     int use[],int n,
00225     double residuals[],
00226     double *rms,
00227     int atob)
00228 {
00229     double x,y;
00230     int i;
00231     int count;
00232     double sum;
00233     double delta;
00234     double dx,dy;
00235 
00236     count = 0;
00237     sum = 0.0;
00238     for (i=0; i < n; i++)
00239     {
00240         if (!use[i]) continue;
00241 
00242         count++;
00243         if (atob)
00244         {
00245             transform_a_into_b (ax[i], ay[i], &x, &y);
00246             dx = x - bx[i];
00247             dy = y - by[i];
00248         }
00249         else
00250         {
00251             transform_b_into_a (bx[i], by[i], &x, &y);
00252             dx = x - ax[i];
00253             dy = y - ay[i];
00254         }
00255 
00256         delta = dx * dx + dy * dy;
00257         residuals[i] = sqrt (delta);
00258         sum += delta;
00259     }
00260     *rms = sqrt (sum / count);
00261 
00262     return 0;
00263 }

Generated on Sun Apr 6 17:32:44 2008 for GRASS by  doxygen 1.5.5