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  * VPWarpA111N
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 VPWarpA111N (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         float gray_acc; int gray_acc_int;
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 < ((1) + (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 * ((1) + (1)));
00323                 out_ptr += lft_zero_cnt * ((1) + (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         gray_acc = 0;
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         gray_acc += (wptr[0]) * (in_ptr[0].clrflt);
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         gray_acc += (wptr[2]) * (in_ptr[in_width].clrflt);
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         gray_acc += (wptr[1]) * (in_ptr[1].clrflt);
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         gray_acc += (wptr[3]) * (in_ptr[in_width + 1].clrflt);
00354         opc_acc += (wptr[3]) * (in_ptr[in_width + 1].opcflt);;
00355                     }
00356                 }
00357                 
00358             
00359         gray_acc_int = gray_acc;
00360         if (gray_acc_int > 255)
00361             gray_acc_int = 255;
00362         ((out_ptr)[0]) = gray_acc_int;
00363             
00364         opc_acc_int = opc_acc * (float)255.;
00365         if (opc_acc_int > 255)
00366             opc_acc_int = 255;
00367         ((out_ptr)[1]) = opc_acc_int;;
00368                 out_ptr += ((1) + (1));
00369                 xfrac += xfrac_incr;
00370                 yfrac += yfrac_incr;
00371                 if (xfrac < 0) {
00372                     xfrac &= 0x7fffffff;
00373                     in_x_int += in_x_incr_int + in_x_incr_dlt;
00374                     in_ptr += in_x_incr_int + in_x_incr_dlt;
00375                 } else {
00376                     in_x_int += in_x_incr_int;
00377                     in_ptr += in_x_incr_int;
00378                 }
00379                 if (yfrac < 0) {
00380                     yfrac &= 0x7fffffff;
00381                     in_y_int += in_y_incr_int + in_y_incr_dlt;
00382                     in_ptr += in_width * (in_y_incr_int + in_y_incr_dlt);
00383                 } else {
00384                     in_y_int += in_y_incr_int;
00385                     in_ptr += in_width * in_y_incr_int;
00386                 }
00387             }
00388 
00389             /* process output pixels affected by four input pixels */
00390             for (x = lft_zero_cnt + lft_edge_cnt;
00391                  x < lft_zero_cnt + lft_edge_cnt + full_cnt; x++) {
00392                 ASSERT(in_x_int >= 0 && in_x_int < in_width-1);
00393                 ASSERT(in_y_int >= 0 && in_y_int < in_height-1);
00394                 ASSERT((GrayIntPixel *)(((char *)in_image + in_y_int *
00395                                 in_bytes_per_scan)) + in_x_int == in_ptr);
00396                 wptr = VPBilirpWeight[FIX31_TO_WGTIND(yfrac)]
00397                                      [FIX31_TO_WGTIND(xfrac)];
00398                 
00399         
00400         gray_acc = (wptr[0]) * (in_ptr[0].clrflt) +
00401                    (wptr[2]) * (in_ptr[in_width].clrflt) +      
00402                    (wptr[1]) * (in_ptr[1].clrflt) +
00403                    (wptr[3]) * (in_ptr[in_width+1].clrflt);
00404         
00405         opc_acc = (wptr[0]) * (in_ptr[0].opcflt) +
00406                   (wptr[2]) * (in_ptr[in_width].opcflt) +       
00407                   (wptr[1]) * (in_ptr[1].opcflt) +
00408                   (wptr[3]) * (in_ptr[in_width+1].opcflt);;
00409                 
00410             
00411         gray_acc_int = gray_acc;
00412         if (gray_acc_int > 255)
00413             gray_acc_int = 255;
00414         ((out_ptr)[0]) = gray_acc_int;
00415             
00416         opc_acc_int = opc_acc * (float)255.;
00417         if (opc_acc_int > 255)
00418             opc_acc_int = 255;
00419         ((out_ptr)[1]) = opc_acc_int;;
00420                 out_ptr += ((1) + (1));
00421                 xfrac += xfrac_incr;
00422                 yfrac += yfrac_incr;
00423                 if (xfrac < 0) {
00424                     xfrac &= 0x7fffffff;
00425                     in_x_int += in_x_incr_int + in_x_incr_dlt;
00426                     in_ptr += in_x_incr_int + in_x_incr_dlt;
00427                 } else {
00428                     in_x_int += in_x_incr_int;
00429                     in_ptr += in_x_incr_int;
00430                 }
00431                 if (yfrac < 0) {
00432                     yfrac &= 0x7fffffff;
00433                     in_y_int += in_y_incr_int + in_y_incr_dlt;
00434                     in_ptr += in_width * (in_y_incr_int + in_y_incr_dlt);
00435                 } else {
00436                     in_y_int += in_y_incr_int;
00437                     in_ptr += in_width * in_y_incr_int;
00438                 }
00439             }
00440 
00441             /* process right edge case pixels */
00442             for (x = lft_zero_cnt + lft_edge_cnt + full_cnt;
00443                  x < lft_zero_cnt + lft_edge_cnt + full_cnt + rgt_edge_cnt;
00444                  x++) {
00445                 wptr = VPBilirpWeight[FIX31_TO_WGTIND(yfrac)]
00446                                      [FIX31_TO_WGTIND(xfrac)];
00447                 
00448         gray_acc = 0;
00449         opc_acc = 0;;
00450                 if (in_x_int >= 0 && in_x_int < in_width) {
00451                     if (in_y_int >= 0 && in_y_int < in_height) {
00452                         
00453         gray_acc += (wptr[0]) * (in_ptr[0].clrflt);
00454         opc_acc += (wptr[0]) * (in_ptr[0].opcflt);;
00455                     }
00456                     if (in_y_int+1 >= 0 && in_y_int+1 < in_height) {
00457                         
00458         gray_acc += (wptr[2]) * (in_ptr[in_width].clrflt);
00459         opc_acc += (wptr[2]) * (in_ptr[in_width].opcflt);;
00460                     }
00461                 }
00462                 if (in_x_int+1 >= 0 && in_x_int+1 < in_width) {
00463                     if (in_y_int >= 0 && in_y_int < in_height) {
00464                         
00465         gray_acc += (wptr[1]) * (in_ptr[1].clrflt);
00466         opc_acc += (wptr[1]) * (in_ptr[1].opcflt);;
00467                     }
00468                     if (in_y_int+1 >= 0 && in_y_int+1 < in_height) {
00469                         
00470         gray_acc += (wptr[3]) * (in_ptr[in_width + 1].clrflt);
00471         opc_acc += (wptr[3]) * (in_ptr[in_width + 1].opcflt);;
00472                     }
00473                 }
00474                 
00475             
00476         gray_acc_int = gray_acc;
00477         if (gray_acc_int > 255)
00478             gray_acc_int = 255;
00479         ((out_ptr)[0]) = gray_acc_int;
00480             
00481         opc_acc_int = opc_acc * (float)255.;
00482         if (opc_acc_int > 255)
00483             opc_acc_int = 255;
00484         ((out_ptr)[1]) = opc_acc_int;;
00485                 out_ptr += ((1) + (1));
00486                 xfrac += xfrac_incr;
00487                 yfrac += yfrac_incr;
00488                 if (xfrac < 0) {
00489                     xfrac &= 0x7fffffff;
00490                     in_x_int += in_x_incr_int + in_x_incr_dlt;
00491                     in_ptr += in_x_incr_int + in_x_incr_dlt;
00492                 } else {
00493                     in_x_int += in_x_incr_int;
00494                     in_ptr += in_x_incr_int;
00495                 }
00496                 if (yfrac < 0) {
00497                     yfrac &= 0x7fffffff;
00498                     in_y_int += in_y_incr_int + in_y_incr_dlt;
00499                     in_ptr += in_width * (in_y_incr_int + in_y_incr_dlt);
00500                 } else {
00501                     in_y_int += in_y_incr_int;
00502                     in_ptr += in_width * in_y_incr_int;
00503                 }
00504             }
00505 
00506             /* zero out unaffected pixels on right edge of scan */
00507             if (rgt_zero_cnt > 0) {
00508                 bzero(out_ptr, rgt_zero_cnt * ((1) + (1)));
00509                 out_ptr += rgt_zero_cnt * ((1) + (1));
00510             }
00511 
00512             /* go on to next scan */
00513             scans_to_next_vertex--;
00514             out_scan_y++;
00515             out_ptr += out_bytes_per_scan - out_width * ((1) + (1));
00516             x_lft_full += full_overlap[region].x_incr_lft;
00517             x_rgt_full += full_overlap[region].x_incr_rgt;
00518             x_lft_part += part_overlap[region].x_incr_lft;
00519             x_rgt_part += part_overlap[region].x_incr_rgt;
00520         } /* next scanline in region */
00521     } /* next region */
00522     ASSERT(out_scan_y == out_height);
00523 }