00001 /*
00002  * DO NOT EDIT THIS FILE! It was created automatically by m4.
00003  */
00004 
00005 /*
00006  * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
00007  * Junior University.  All rights reserved.
00008  * 
00009  * Permission to use, copy, modify and distribute this software and its
00010  * documentation for any purpose is hereby granted without fee, provided
00011  * that the above copyright notice and this permission notice appear in
00012  * all copies of this software and that you do not sell the software.
00013  * Commercial licensing is available by contacting the author.
00014  * 
00015  * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
00016  * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
00017  * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
00018  * 
00019  * Author:
00020  *    Phil Lacroute
00021  *    Computer Systems Laboratory
00022  *    Electrical Engineering Dept.
00023  *    Stanford University
00024  */
00025 /*
00026  * vp_warpA.m4
00027  *
00028  * One-pass image warping routine for affine transformations.
00029  *
00030  * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
00031  * Junior University.  All rights reserved.
00032  *
00033  * Permission to use, copy, modify and distribute this software and its
00034  * documentation for any purpose is hereby granted without fee, provided
00035  * that the above copyright notice and this permission notice appear in
00036  * all copies of this software and that you do not sell the software.
00037  * Commercial licensing is available by contacting the author.
00038  * 
00039  * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
00040  * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
00041  * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
00042  *
00043  * Author:
00044  *    Phil Lacroute
00045  *    Computer Systems Laboratory
00046  *    Electrical Engineering Dept.
00047  *    Stanford University
00048  */
00049 
00050 /*
00051  * $Date: 2001/12/17 16:16:23 $
00052  * $Revision: 1.1 $
00053  */
00054 
00055 #include "vp_global.h"
00056 
00057 
00058 
00059     
00060     
00061     
00062     
00063     
00064 
00065     
00066     
00067     
00068 
00069 
00070 
00071 
00072 
00073     
00074     
00075     
00076     
00077     
00078 
00079 
00080 
00081 
00082 
00083 
00084 
00085 
00086 
00087     
00088     
00089     
00090     
00091     
00092 
00093 
00094 
00095 
00096 
00097 
00098 
00099 
00100         
00101 
00102 /* convert a float in the interval [0-1) to a 31-bit fixed point */
00103 #define FLTFRAC_TO_FIX31(f)     ((int)((f) * 2147483648.))
00104 
00105 /* convert a 31-bit fixed point to a weight table index */
00106 #define FIX31_TO_WGTIND(f)      ((f) >> (31 - WARP_WEIGHT_INDEX_BITS))
00107 
00108 extern float VPBilirpWeight[WARP_WEIGHT_ENTRIES][WARP_WEIGHT_ENTRIES][4];
00109 
00110 /*
00111  * VPWarpA101N
00112  *
00113  * One-pass warper.  Transforms in_image to out_image according to
00114  * the affine warp specified by warp_matrix.  The resampling filter
00115  * is a bilirp (suitable for upsampling only).
00116  */
00117 
00118 void
00119 VPWarpA101N (in_image, in_width, in_height, in_bytes_per_scan,
00120           out_image, out_width, out_height, out_bytes_per_scan,
00121           warp_matrix)
00122 GrayIntPixel *in_image;         /* input image data */
00123 int in_width;                   /* size of input image */
00124 int in_height;
00125 int in_bytes_per_scan;          /* bytes per scanline in input image */
00126 char *out_image;                /* output image data */
00127 int out_width;                  /* size of output image */
00128 int out_height;
00129 int out_bytes_per_scan;         /* bytes per scanline in output image */
00130 vpMatrix3 warp_matrix;          /* [ outx ]                 [ inx ] */
00131                                 /* [ outy ] = warp_matrix * [ iny ] */
00132                                 /* [   1  ]                 [  1  ] */
00133 {
00134     Trapezoid full_overlap[9];  /* description of the area of overlap
00135                                    of output image (shrunk by the size
00136                                    of the filter kernel) with input image */
00137     Trapezoid part_overlap[9];  /* description of the area of overlap
00138                                    of output image (unlarged by the size
00139                                    of the filter kernel) with input image */
00140     int region;                 /* index into full/part_overlap */
00141     char *out_ptr;              /* pointer to current pixel of output image */
00142     int out_scan_y;             /* coordinate of current output scanline */
00143     int scans_to_next_vertex;   /* number of scans left to process before
00144                                    the next vertex is reached */
00145     GrayIntPixel *in_ptr;       /* pointer to current pixel of input image */
00146     double x_lft_full, x_rgt_full; /* intersection of scan with full_overlap */
00147     double x_lft_part, x_rgt_part; /* intersection of scan with part_overlap */
00148     int no_full_pixels;         /* true if full_overlap is empty for scan */
00149     double in_x, in_y;          /* exact coordinates in the input image of
00150                                    the current output image pixel */
00151     int in_x_int, in_y_int;     /* coordinates of the nearest input image
00152                                    pixel to the upper-left of the current
00153                                    output image pixel */
00154     int xfrac, yfrac;           /* in_x - in_x_int and in_y - in_y_int,
00155                                    stored as a fixed-point number with 31 bits
00156                                    of fraction */
00157     int xfrac_incr, yfrac_incr; /* increments to xfrac and yfrac to give
00158                                    the fractions for the next output image
00159                                    pixel in the current scan */
00160     double in_x_incr, in_y_incr;/* increments to in_x and in_y to give the
00161                                    input image coordinates of the next
00162                                    output image pixel in the current scan 
00163                                    (equal to dx_in/dx_out and dy_in/dx_out) */
00164     int in_x_incr_int, in_y_incr_int; /* integer part of in_x/y_incr */
00165     int in_x_incr_dlt, in_y_incr_dlt; /* sign of in_x/y_incr */
00166     float *wptr;                /* pointer into weight table */
00167     int lft_zero_cnt;           /* # zero pixels on left edge of scan */
00168     int lft_edge_cnt;           /* # pixels on left w/ part filter overlap */
00169     int full_cnt;               /* # pixels w/ full filter overlap */
00170     int rgt_edge_cnt;           /* # pixels on rgt w/ part filter overlap */
00171     int rgt_zero_cnt;           /* # zero pixels on right edge of scan */
00172     int x;                      /* pixel index */
00173     
00174         
00175         float opc_acc; int opc_acc_int;         /* pixel accumulator */
00176     double denom;
00177     int c;
00178 
00179 #ifdef DEBUG
00180     {
00181         int y;
00182 
00183         for (y = 0; y < out_height; y++) {
00184             out_ptr = out_image + y*out_bytes_per_scan;
00185             for (x = 0; x < out_width; x++) {
00186                 for (c = 0; c < ((0) + (1)); c++)
00187                     *out_ptr++ = 255;
00188             }
00189         }
00190     }
00191 #endif
00192 
00193     /* initialize tables */
00194     VPComputeWarpTables();
00195 
00196     /* compute the intersection of the input image and the output image */
00197     /* filter width = 2.0 in input image space (triangle filter) */
00198     VPAffineImageOverlap(in_width, in_height, out_width, out_height,
00199                          warp_matrix, 2., full_overlap, part_overlap);
00200 
00201     /* compute the output image */
00202     out_ptr = out_image;
00203     out_scan_y = 0;
00204     denom = 1. / (warp_matrix[0][0] * warp_matrix[1][1] -
00205                   warp_matrix[0][1] * warp_matrix[1][0]);
00206     in_x_incr = warp_matrix[1][1]*denom;
00207     in_y_incr = -warp_matrix[1][0]*denom;
00208     if (in_x_incr < 0) {
00209         in_x_incr_int = (int)ceil(in_x_incr);
00210         in_x_incr_dlt = -1;
00211     } else {
00212         in_x_incr_int = (int)floor(in_x_incr);
00213         in_x_incr_dlt = 1;
00214     }
00215     if (in_y_incr < 0) {
00216         in_y_incr_int = (int)ceil(in_y_incr);
00217         in_y_incr_dlt = -1;
00218     } else {
00219         in_y_incr_int = (int)floor(in_y_incr);
00220         in_y_incr_dlt = 1;
00221     }
00222     xfrac_incr = FLTFRAC_TO_FIX31(in_x_incr - in_x_incr_int);
00223     yfrac_incr = FLTFRAC_TO_FIX31(in_y_incr - in_y_incr_int);
00224     for (region = 0; region < 9; region++) {
00225         /* check for empty region */
00226         if (part_overlap[region].miny >= out_height) {
00227             break;
00228         }
00229 
00230         /* check if this region of the output image is unaffected by
00231            the input image */
00232         if (part_overlap[region].x_top_lft >
00233             part_overlap[region].x_top_rgt) {
00234             c = (part_overlap[region].maxy - part_overlap[region].miny + 1) *
00235                 out_bytes_per_scan;
00236             bzero(out_ptr, c);
00237             out_ptr += c;
00238             out_scan_y += part_overlap[region].maxy -
00239                           part_overlap[region].miny + 1;
00240             continue;
00241         }
00242 
00243         /* process scanlines of this region */
00244         scans_to_next_vertex = part_overlap[region].maxy -
00245                                part_overlap[region].miny + 1;
00246         x_lft_full = full_overlap[region].x_top_lft;
00247         x_rgt_full = full_overlap[region].x_top_rgt;
00248         x_lft_part = part_overlap[region].x_top_lft;
00249         x_rgt_part = part_overlap[region].x_top_rgt;
00250         if (x_lft_full > x_rgt_full)
00251             no_full_pixels = 1;
00252         else
00253             no_full_pixels = 0;
00254         ASSERT(scans_to_next_vertex > 0);
00255         ASSERT(out_scan_y == part_overlap[region].miny);
00256         while (scans_to_next_vertex > 0) {
00257             /* compute the portions of the scanline which are zero
00258                and which intersect the full and partially-full regions */
00259             lft_zero_cnt = (int)floor(x_lft_part);
00260             if (lft_zero_cnt < 0)
00261                 lft_zero_cnt = 0;
00262             else if (lft_zero_cnt > out_width)
00263                 lft_zero_cnt = out_width;
00264             if (no_full_pixels) {
00265                 lft_edge_cnt = (int)ceil(x_rgt_part);
00266                 if (lft_edge_cnt < 0)
00267                     lft_edge_cnt = 0;
00268                 else if (lft_edge_cnt > out_width)
00269                     lft_edge_cnt = out_width;
00270                 lft_edge_cnt -= lft_zero_cnt;
00271                 if (lft_edge_cnt < 0)
00272                     lft_edge_cnt = 0;
00273                 full_cnt = 0;
00274                 rgt_edge_cnt = 0;
00275                 rgt_zero_cnt = out_width - lft_zero_cnt - lft_edge_cnt;
00276             } else {
00277                 lft_edge_cnt = (int)ceil(x_lft_full);
00278                 if (lft_edge_cnt < 0)
00279                     lft_edge_cnt = 0;
00280                 else if (lft_edge_cnt > out_width)
00281                     lft_edge_cnt = out_width;
00282                 lft_edge_cnt -= lft_zero_cnt;
00283                 if (lft_edge_cnt < 0)
00284                     lft_edge_cnt = 0;
00285                 full_cnt = (int)floor(x_rgt_full);
00286                 if (full_cnt < 0)
00287                     full_cnt = 0;
00288                 else if (full_cnt > out_width)
00289                     full_cnt = out_width;
00290                 full_cnt -= lft_edge_cnt + lft_zero_cnt;
00291                 if (full_cnt < 0)
00292                     full_cnt = 0;
00293                 rgt_edge_cnt = (int)ceil(x_rgt_part);
00294                 if (rgt_edge_cnt < 0)
00295                     rgt_edge_cnt = 0;
00296                 else if (rgt_edge_cnt > out_width)
00297                     rgt_edge_cnt = out_width;
00298                 rgt_edge_cnt -= full_cnt + lft_edge_cnt + lft_zero_cnt;
00299                 if (rgt_edge_cnt < 0)
00300                     rgt_edge_cnt = 0;
00301                 rgt_zero_cnt = out_width - lft_zero_cnt - lft_edge_cnt - 
00302                                full_cnt - rgt_edge_cnt;
00303             }
00304 
00305             /* reverse map the first left-edge output pixel coordinate into
00306                the input image coordinate system */
00307             in_x = ((lft_zero_cnt - warp_matrix[0][2]) * warp_matrix[1][1] -
00308                     (out_scan_y - warp_matrix[1][2])*warp_matrix[0][1])*denom;
00309             in_y = (-(lft_zero_cnt - warp_matrix[0][2]) * warp_matrix[1][0] +
00310                     (out_scan_y - warp_matrix[1][2])*warp_matrix[0][0])*denom;
00311             in_x_int = (int)floor(in_x);
00312             in_y_int = (int)floor(in_y);
00313             in_ptr = (GrayIntPixel *)(((char *)in_image + in_y_int *
00314                                        in_bytes_per_scan)) + in_x_int;
00315 
00316             /* compute the weight lookup table indices and increments */
00317             xfrac = FLTFRAC_TO_FIX31(in_x - in_x_int);
00318             yfrac = FLTFRAC_TO_FIX31(in_y - in_y_int);
00319 
00320             /* zero out unaffected pixels on left edge of scan */
00321             if (lft_zero_cnt > 0) {
00322                 bzero(out_ptr, lft_zero_cnt * ((0) + (1)));
00323                 out_ptr += lft_zero_cnt * ((0) + (1));
00324             }
00325 
00326             /* process left edge case pixels */
00327             for (x = lft_zero_cnt; x < lft_zero_cnt + lft_edge_cnt; x++) {
00328                 wptr = VPBilirpWeight[FIX31_TO_WGTIND(yfrac)]
00329                                      [FIX31_TO_WGTIND(xfrac)];
00330                 
00331         
00332         opc_acc = 0;;
00333                 if (in_x_int >= 0 && in_x_int < in_width) {
00334                     if (in_y_int >= 0 && in_y_int < in_height) {
00335                         
00336         
00337         opc_acc += (wptr[0]) * (in_ptr[0].opcflt);;
00338                     }
00339                     if (in_y_int+1 >= 0 && in_y_int+1 < in_height) {
00340                         
00341         
00342         opc_acc += (wptr[2]) * (in_ptr[in_width].opcflt);;
00343                     }
00344                 }
00345                 if (in_x_int+1 >= 0 && in_x_int+1 < in_width) {
00346                     if (in_y_int >= 0 && in_y_int < in_height) {
00347                         
00348         
00349         opc_acc += (wptr[1]) * (in_ptr[1].opcflt);;
00350                     }
00351                     if (in_y_int+1 >= 0 && in_y_int+1 < in_height) {
00352                         
00353         
00354         opc_acc += (wptr[3]) * (in_ptr[in_width + 1].opcflt);;
00355                     }
00356                 }
00357                 
00358             
00359             
00360         opc_acc_int = opc_acc * (float)255.;
00361         if (opc_acc_int > 255)
00362             opc_acc_int = 255;
00363         ((out_ptr)[0]) = opc_acc_int;;
00364                 out_ptr += ((0) + (1));
00365                 xfrac += xfrac_incr;
00366                 yfrac += yfrac_incr;
00367                 if (xfrac < 0) {
00368                     xfrac &= 0x7fffffff;
00369                     in_x_int += in_x_incr_int + in_x_incr_dlt;
00370                     in_ptr += in_x_incr_int + in_x_incr_dlt;
00371                 } else {
00372                     in_x_int += in_x_incr_int;
00373                     in_ptr += in_x_incr_int;
00374                 }
00375                 if (yfrac < 0) {
00376                     yfrac &= 0x7fffffff;
00377                     in_y_int += in_y_incr_int + in_y_incr_dlt;
00378                     in_ptr += in_width * (in_y_incr_int + in_y_incr_dlt);
00379                 } else {
00380                     in_y_int += in_y_incr_int;
00381                     in_ptr += in_width * in_y_incr_int;
00382                 }
00383             }
00384 
00385             /* process output pixels affected by four input pixels */
00386             for (x = lft_zero_cnt + lft_edge_cnt;
00387                  x < lft_zero_cnt + lft_edge_cnt + full_cnt; x++) {
00388                 ASSERT(in_x_int >= 0 && in_x_int < in_width-1);
00389                 ASSERT(in_y_int >= 0 && in_y_int < in_height-1);
00390                 ASSERT((GrayIntPixel *)(((char *)in_image + in_y_int *
00391                                 in_bytes_per_scan)) + in_x_int == in_ptr);
00392                 wptr = VPBilirpWeight[FIX31_TO_WGTIND(yfrac)]
00393                                      [FIX31_TO_WGTIND(xfrac)];
00394                 
00395         
00396         
00397         opc_acc = (wptr[0]) * (in_ptr[0].opcflt) +
00398                   (wptr[2]) * (in_ptr[in_width].opcflt) +       
00399                   (wptr[1]) * (in_ptr[1].opcflt) +
00400                   (wptr[3]) * (in_ptr[in_width+1].opcflt);;
00401                 
00402             
00403             
00404         opc_acc_int = opc_acc * (float)255.;
00405         if (opc_acc_int > 255)
00406             opc_acc_int = 255;
00407         ((out_ptr)[0]) = opc_acc_int;;
00408                 out_ptr += ((0) + (1));
00409                 xfrac += xfrac_incr;
00410                 yfrac += yfrac_incr;
00411                 if (xfrac < 0) {
00412                     xfrac &= 0x7fffffff;
00413                     in_x_int += in_x_incr_int + in_x_incr_dlt;
00414                     in_ptr += in_x_incr_int + in_x_incr_dlt;
00415                 } else {
00416                     in_x_int += in_x_incr_int;
00417                     in_ptr += in_x_incr_int;
00418                 }
00419                 if (yfrac < 0) {
00420                     yfrac &= 0x7fffffff;
00421                     in_y_int += in_y_incr_int + in_y_incr_dlt;
00422                     in_ptr += in_width * (in_y_incr_int + in_y_incr_dlt);
00423                 } else {
00424                     in_y_int += in_y_incr_int;
00425                     in_ptr += in_width * in_y_incr_int;
00426                 }
00427             }
00428 
00429             /* process right edge case pixels */
00430             for (x = lft_zero_cnt + lft_edge_cnt + full_cnt;
00431                  x < lft_zero_cnt + lft_edge_cnt + full_cnt + rgt_edge_cnt;
00432                  x++) {
00433                 wptr = VPBilirpWeight[FIX31_TO_WGTIND(yfrac)]
00434                                      [FIX31_TO_WGTIND(xfrac)];
00435                 
00436         
00437         opc_acc = 0;;
00438                 if (in_x_int >= 0 && in_x_int < in_width) {
00439                     if (in_y_int >= 0 && in_y_int < in_height) {
00440                         
00441         
00442         opc_acc += (wptr[0]) * (in_ptr[0].opcflt);;
00443                     }
00444                     if (in_y_int+1 >= 0 && in_y_int+1 < in_height) {
00445                         
00446         
00447         opc_acc += (wptr[2]) * (in_ptr[in_width].opcflt);;
00448                     }
00449                 }
00450                 if (in_x_int+1 >= 0 && in_x_int+1 < in_width) {
00451                     if (in_y_int >= 0 && in_y_int < in_height) {
00452                         
00453         
00454         opc_acc += (wptr[1]) * (in_ptr[1].opcflt);;
00455                     }
00456                     if (in_y_int+1 >= 0 && in_y_int+1 < in_height) {
00457                         
00458         
00459         opc_acc += (wptr[3]) * (in_ptr[in_width + 1].opcflt);;
00460                     }
00461                 }
00462                 
00463             
00464             
00465         opc_acc_int = opc_acc * (float)255.;
00466         if (opc_acc_int > 255)
00467             opc_acc_int = 255;
00468         ((out_ptr)[0]) = opc_acc_int;;
00469                 out_ptr += ((0) + (1));
00470                 xfrac += xfrac_incr;
00471                 yfrac += yfrac_incr;
00472                 if (xfrac < 0) {
00473                     xfrac &= 0x7fffffff;
00474                     in_x_int += in_x_incr_int + in_x_incr_dlt;
00475                     in_ptr += in_x_incr_int + in_x_incr_dlt;
00476                 } else {
00477                     in_x_int += in_x_incr_int;
00478                     in_ptr += in_x_incr_int;
00479                 }
00480                 if (yfrac < 0) {
00481                     yfrac &= 0x7fffffff;
00482                     in_y_int += in_y_incr_int + in_y_incr_dlt;
00483                     in_ptr += in_width * (in_y_incr_int + in_y_incr_dlt);
00484                 } else {
00485                     in_y_int += in_y_incr_int;
00486                     in_ptr += in_width * in_y_incr_int;
00487                 }
00488             }
00489 
00490             /* zero out unaffected pixels on right edge of scan */
00491             if (rgt_zero_cnt > 0) {
00492                 bzero(out_ptr, rgt_zero_cnt * ((0) + (1)));
00493                 out_ptr += rgt_zero_cnt * ((0) + (1));
00494             }
00495 
00496             /* go on to next scan */
00497             scans_to_next_vertex--;
00498             out_scan_y++;
00499             out_ptr += out_bytes_per_scan - out_width * ((0) + (1));
00500             x_lft_full += full_overlap[region].x_incr_lft;
00501             x_rgt_full += full_overlap[region].x_incr_rgt;
00502             x_lft_part += part_overlap[region].x_incr_lft;
00503             x_rgt_part += part_overlap[region].x_incr_rgt;
00504         } /* next scanline in region */
00505     } /* next region */
00506     ASSERT(out_scan_y == out_height);
00507 }