00001 /*
00002  * vp_warp.c
00003  *
00004  * Support routines for 1-pass image warping.
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 /*
00027  * $Date: 2001/12/17 16:16:24 $
00028  * $Revision: 1.1 $
00029  */
00030 
00031 #include "vp_global.h"
00032 
00033 /* weights for bilirp filter; index with (yfrac, xfrac, tap number) */
00034 float VPBilirpWeight[WARP_WEIGHT_ENTRIES][WARP_WEIGHT_ENTRIES][4];
00035 static int BilirpWeightsReady = 0;      /* true if weight table initialized */
00036 
00037 static void OrderCoords ANSI_ARGS((double coords[4][2], double lft[3][2],
00038     double rgt[3][2]));
00039 
00040 /*
00041  * VPComputeWarpTables
00042  *
00043  * Precompute lookup tables for fast filter convolution during warping.
00044  */
00045 
00046 void
00047 VPComputeWarpTables()
00048 {
00049     float *wptr;        /* pointer into weight table */
00050 
00051     int x, y;
00052     double in_x, in_y;
00053 
00054     if (BilirpWeightsReady)
00055         return;
00056 
00057 #ifdef MEMSPY
00058     bin_init(BinNumber(__LINE__, __FILE__, "VPBilirpWeight"), -1, -1,
00059              VPBilirpWeight, sizeof(VPBilirpWeight), "VPBilirpWeight");
00060 #endif
00061 
00062     wptr = &VPBilirpWeight[0][0][0];
00063     for (y = 0; y < WARP_WEIGHT_ENTRIES; y++) {
00064         in_y = (double)y / (WARP_WEIGHT_ENTRIES-1);
00065         for (x = 0; x < WARP_WEIGHT_ENTRIES; x++) {
00066             in_x = (double)x / (WARP_WEIGHT_ENTRIES-1);
00067             *wptr++ = (1. - in_x)*(1. - in_y);
00068             *wptr++ = in_x * (1. - in_y);
00069             *wptr++ = (1. - in_x) * in_y;
00070             *wptr++ = 1. - wptr[-1] - wptr[-2] - wptr[-3];
00071         }
00072     }
00073 
00074     BilirpWeightsReady = 1;
00075 }
00076 
00077 /*
00078  * VPAffineComputeImageOverlap
00079  *
00080  * Compute the intersection of the intermediate image and the final image.
00081  * This intersection is broken up into a series of trapezoids whose
00082  * bases are parallel to the scanlines of the final image (for ease
00083  * of scan conversion).  Two sets of trapezoids are computed: one set
00084  * gives the region of the final image for which the filter kernel is
00085  * completely contained within the intermedaite image.  The second set of
00086  * trapezoids gives the region of the final image for which the filter
00087  * kernel overlaps at least some of the intermedaite image, but may or may not
00088  * be completely contained.  Any region of the final image which is not
00089  * within the second set of trapezoids is not affected by the intermediate
00090  * image at all and should be set to the background color.
00091  */
00092 
00093 void
00094 VPAffineImageOverlap(in_width, in_height, out_width, out_height,
00095                      warp_matrix, filter_width, full_overlap, part_overlap)
00096 int in_width, in_height;        /* input (intermediate) image size */
00097 int out_width, out_height;      /* output (final) image size */
00098 vpMatrix3 warp_matrix;          /* affine transformation from input image
00099                                    to output image */
00100 double filter_width;            /* filter kernel width in input image
00101                                    coordinates */
00102 Trapezoid full_overlap[9]; /* portion of the intersection with full
00103                                       filter kernel overlap */
00104 Trapezoid part_overlap[9]; /* portion of the intersection with partial
00105                                       filter kernel overlap */
00106 {
00107     double int_lft[3][2];       /* coordinates bounding the left side of the
00108                                    full_overlap region in output coordinates,
00109                                    sorted from top (lowest y) to bottom */
00110     double int_rgt[3][2];       /* right side of full_overlap region */
00111     double ext_lft[3][2];       /* left side of part_overlap region */
00112     double ext_rgt[3][2];       /* right side of part_overlap region */
00113     double coords[4][2];
00114     double inset;
00115     int ilft, irgt, elft, ergt;
00116     int region;
00117     int lasty, nexty, y;
00118 
00119     /* compute the bounding box of the full_overlap region and store it
00120        in int_lft and int_rgt */
00121     inset = -1.0 + filter_width / 2.0 + 1.0e-5;
00122     coords[0][0] = warp_matrix[0][0] * inset +
00123                    warp_matrix[0][1] * inset + 
00124                    warp_matrix[0][2];
00125     coords[0][1] = warp_matrix[1][0] * inset +
00126                    warp_matrix[1][1] * inset + 
00127                    warp_matrix[1][2];
00128     coords[1][0] = warp_matrix[0][0] * (in_width - 1 - inset) +
00129                    warp_matrix[0][1] * inset + 
00130                    warp_matrix[0][2];
00131     coords[1][1] = warp_matrix[1][0] * (in_width - 1 - inset) +
00132                    warp_matrix[1][1] * inset + 
00133                    warp_matrix[1][2];
00134     coords[2][0] = warp_matrix[0][0] * (in_width - 1 - inset) +
00135                    warp_matrix[0][1] * (in_height - 1 - inset) + 
00136                    warp_matrix[0][2];
00137     coords[2][1] = warp_matrix[1][0] * (in_width - 1 - inset) +
00138                    warp_matrix[1][1] * (in_height - 1 - inset) + 
00139                    warp_matrix[1][2];
00140     coords[3][0] = warp_matrix[0][0] * inset +
00141                    warp_matrix[0][1] * (in_height - 1 - inset) + 
00142                    warp_matrix[0][2];
00143     coords[3][1] = warp_matrix[1][0] * inset +
00144                    warp_matrix[1][1] * (in_height - 1 - inset) + 
00145                    warp_matrix[1][2];
00146     OrderCoords(coords, int_lft, int_rgt);
00147 
00148     /* compute the bounding box of the part_overlap region and store it
00149        in int_lft and int_rgt */
00150     inset = -filter_width / 2.0;
00151     coords[0][0] = warp_matrix[0][0] * inset +
00152                    warp_matrix[0][1] * inset + 
00153                    warp_matrix[0][2];
00154     coords[0][1] = warp_matrix[1][0] * inset +
00155                    warp_matrix[1][1] * inset + 
00156                    warp_matrix[1][2];
00157     coords[1][0] = warp_matrix[0][0] * (in_width - 1 - inset) +
00158                    warp_matrix[0][1] * inset + 
00159                    warp_matrix[0][2];
00160     coords[1][1] = warp_matrix[1][0] * (in_width - 1 - inset) +
00161                    warp_matrix[1][1] * inset + 
00162                    warp_matrix[1][2];
00163     coords[2][0] = warp_matrix[0][0] * (in_width - 1 - inset) +
00164                    warp_matrix[0][1] * (in_height - 1 - inset) + 
00165                    warp_matrix[0][2];
00166     coords[2][1] = warp_matrix[1][0] * (in_width - 1 - inset) +
00167                    warp_matrix[1][1] * (in_height - 1 - inset) + 
00168                    warp_matrix[1][2];
00169     coords[3][0] = warp_matrix[0][0] * inset +
00170                    warp_matrix[0][1] * (in_height - 1 - inset) + 
00171                    warp_matrix[0][2];
00172     coords[3][1] = warp_matrix[1][0] * inset +
00173                    warp_matrix[1][1] * (in_height - 1 - inset) + 
00174                    warp_matrix[1][2];
00175     OrderCoords(coords, ext_lft, ext_rgt);
00176 
00177     for (ilft = 0; ilft < 3 && int_lft[ilft][1] <= 0.; ilft++);
00178     for (irgt = 0; irgt < 3 && int_rgt[irgt][1] <= 0.; irgt++);
00179     for (elft = 0; elft < 3 && ext_lft[elft][1] <= 0.; elft++);
00180     for (ergt = 0; ergt < 3 && ext_rgt[ergt][1] <= 0.; ergt++);
00181     region = 0;
00182     lasty = -1;
00183     while (lasty < out_height-1) {
00184         ASSERT(region < 9);
00185 
00186         /* find nexty */
00187         nexty = out_height - 1;
00188         if (ilft < 3) {
00189             y = (int)floor(int_lft[ilft][1]);
00190             if (nexty > y)
00191                 nexty = y;
00192         }
00193         if (irgt < 3) {
00194             y = (int)floor(int_rgt[irgt][1]);
00195             if (nexty > y)
00196                 nexty = y;
00197         }
00198         if (elft < 3) {
00199             y = (int)floor(ext_lft[elft][1]);
00200             if (nexty > y)
00201                 nexty = y;
00202         }
00203         if (ergt < 3) {
00204             y = (int)floor(ext_rgt[ergt][1]);
00205             if (nexty > y)
00206                 nexty = y;
00207         }
00208         ASSERT((ilft == 0 && (int)floor(int_lft[0][1]) >= nexty) ||
00209                (ilft == 3 && (int)floor(int_lft[2][1]) <= lasty) ||
00210                (((int)floor(int_lft[ilft-1][1]) <= lasty || lasty == -1)
00211                 && (int)floor(int_lft[ilft][1]) >= nexty));
00212         ASSERT((irgt == 0 && (int)floor(int_rgt[0][1]) >= nexty) ||
00213                (irgt == 3 && (int)floor(int_rgt[2][1]) <= lasty) ||
00214                (((int)floor(int_rgt[irgt-1][1]) <= lasty || lasty == -1)
00215                 && (int)floor(int_rgt[irgt][1]) >= nexty));
00216         ASSERT((elft == 0 && (int)floor(ext_lft[0][1]) >= nexty) ||
00217                (elft == 3 && (int)floor(ext_lft[2][1]) <= lasty) ||
00218                (((int)floor(ext_lft[elft-1][1]) <= lasty || lasty == -1)
00219                 && (int)floor(ext_lft[elft][1]) >= nexty));
00220         ASSERT((ergt == 0 && (int)floor(ext_rgt[0][1]) >= nexty) ||
00221                (ergt == 3 && (int)floor(ext_rgt[2][1]) <= lasty) ||
00222                (((int)floor(ext_rgt[ergt-1][1]) <= lasty || lasty == -1)
00223                 && (int)floor(ext_rgt[ergt][1]) >= nexty));
00224         full_overlap[region].miny = lasty + 1;
00225         full_overlap[region].maxy = nexty;
00226         part_overlap[region].miny = lasty + 1;
00227         part_overlap[region].maxy = nexty;
00228         if (ilft == 0 || ilft == 3) {
00229             /* this trapezoid does not intersect full_overlap */
00230             full_overlap[region].x_top_lft = 0;
00231             full_overlap[region].x_top_rgt = -1;
00232             full_overlap[region].x_incr_lft = 0;
00233             full_overlap[region].x_incr_rgt = 0;
00234         } else {
00235             full_overlap[region].x_incr_lft = 
00236                 (int_lft[ilft][0] - int_lft[ilft-1][0]) /
00237                 (int_lft[ilft][1] - int_lft[ilft-1][1]);
00238             full_overlap[region].x_top_lft =
00239                 int_lft[ilft-1][0] + (lasty+1 - int_lft[ilft-1][1]) *
00240                 full_overlap[region].x_incr_lft;
00241             full_overlap[region].x_incr_rgt = 
00242                 (int_rgt[irgt][0] - int_rgt[irgt-1][0]) /
00243                 (int_rgt[irgt][1] - int_rgt[irgt-1][1]);
00244             full_overlap[region].x_top_rgt =
00245                 int_rgt[irgt-1][0] + (lasty+1 - int_rgt[irgt-1][1]) *
00246                 full_overlap[region].x_incr_rgt;
00247         }
00248         if (elft == 0 || elft == 3) {
00249             /* this trapezoid does not intersect part_overlap */
00250             part_overlap[region].x_top_lft = 0;
00251             part_overlap[region].x_top_rgt = -1;
00252             part_overlap[region].x_incr_lft = 0;
00253             part_overlap[region].x_incr_rgt = 0;
00254         } else {
00255             part_overlap[region].x_incr_lft = 
00256                 (ext_lft[elft][0] - ext_lft[elft-1][0]) /
00257                 (ext_lft[elft][1] - ext_lft[elft-1][1]);
00258             part_overlap[region].x_top_lft =
00259                 ext_lft[elft-1][0] + (lasty+1 - ext_lft[elft-1][1]) *
00260                 part_overlap[region].x_incr_lft;
00261             part_overlap[region].x_incr_rgt = 
00262                 (ext_rgt[ergt][0] - ext_rgt[ergt-1][0]) /
00263                 (ext_rgt[ergt][1] - ext_rgt[ergt-1][1]);
00264             part_overlap[region].x_top_rgt =
00265                 ext_rgt[ergt-1][0] + (lasty+1 - ext_rgt[ergt-1][1]) *
00266                 part_overlap[region].x_incr_rgt;
00267         }
00268         ASSERT(!(full_overlap[region].x_top_lft <= 
00269                  full_overlap[region].x_top_rgt &&
00270                  part_overlap[region].x_top_lft >
00271                  part_overlap[region].x_top_rgt));
00272         for (; ilft < 3 && (int)floor(int_lft[ilft][1]) <= nexty; ilft++);
00273         for (; irgt < 3 && (int)floor(int_rgt[irgt][1]) <= nexty; irgt++);
00274         for (; elft < 3 && (int)floor(ext_lft[elft][1]) <= nexty; elft++);
00275         for (; ergt < 3 && (int)floor(ext_rgt[ergt][1]) <= nexty; ergt++);
00276         region++;
00277         lasty = nexty;
00278     }
00279     for (; region < 9; region++) {
00280         full_overlap[region].miny = out_height;
00281         full_overlap[region].maxy = out_height;
00282         part_overlap[region].miny = out_height;
00283         part_overlap[region].maxy = out_height;
00284         full_overlap[region].x_top_lft = 0;
00285         full_overlap[region].x_top_rgt = -1;
00286         full_overlap[region].x_incr_lft = 0;
00287         full_overlap[region].x_incr_rgt = 0;
00288         part_overlap[region].x_top_lft = 0;
00289         part_overlap[region].x_top_rgt = -1;
00290         part_overlap[region].x_incr_lft = 0;
00291         part_overlap[region].x_incr_rgt = 0;
00292     }
00293 
00294 #ifdef DEBUG_OVERLAP
00295     for (region = 0; region < 9; region++) {
00296         printf("region %d: y = %d to %d, [%d+%g [%d+%g %d+%g] %d+%g]\n",
00297                region,
00298                full_overlap[region].miny,
00299                full_overlap[region].maxy,
00300                (int)part_overlap[region].x_top_lft,
00301                part_overlap[region].x_incr_lft,
00302                (int)full_overlap[region].x_top_lft,
00303                full_overlap[region].x_incr_lft,
00304                (int)full_overlap[region].x_top_rgt,
00305                full_overlap[region].x_incr_rgt,
00306                (int)part_overlap[region].x_top_rgt,
00307                part_overlap[region].x_incr_rgt);
00308     }
00309 #endif
00310 }
00311 
00312 /*
00313  * OrderCoords
00314  *
00315  * Sort an array of coordinates.
00316  */
00317 
00318 static void
00319 OrderCoords(coords, lft, rgt)
00320 double coords[4][2];    /* inputs */
00321 double lft[3][2];       /* outputs */
00322 double rgt[3][2];
00323 {
00324     int index;
00325     double swap_buf;
00326     double sorted_coords[4][2];
00327     double xmid;
00328 
00329     /* sort the coordinates as follows:
00330        coords[0]: smallest y value; if there is a tie, then smallest x value
00331                   to break the tie; this is the top or top left corner
00332        coords[1]: next coordinate in CCW order; this is the left or bottom
00333                   left corner
00334        coords[2]: next coordinate in CCW order; this is the bottom or
00335                   bottom right corner
00336        coords[3]: next coordinate in CCW order; this is the right or
00337                   top right corner */
00338 
00339     /* search for coords[0] */
00340     index = 0;
00341     if (coords[1][1] < coords[index][1] ||
00342         (coords[1][1] == coords[index][1] && coords[1][0] < coords[index][0]))
00343         index = 1;
00344     if (coords[2][1] < coords[index][1] ||
00345         (coords[2][1] == coords[index][1] && coords[2][0] < coords[index][0]))
00346         index = 2;
00347     if (coords[3][1] < coords[index][1] ||
00348         (coords[3][1] == coords[index][1] && coords[3][0] < coords[index][0]))
00349         index = 3;
00350     sorted_coords[0][0] = coords[index][0];
00351     sorted_coords[0][1] = coords[index][1];
00352 
00353     /* coords[2] is the opposite corner */
00354     sorted_coords[2][0] = coords[(index+2)%4][0];
00355     sorted_coords[2][1] = coords[(index+2)%4][1];
00356 
00357     /* determine which of the remaining two coordinates is to the left
00358        of the line joining coords[0] and coords[2]; this coordinate
00359        is coords[1], and the final remaining coordinate is coords[3] */
00360     index = (index + 1) % 4;
00361     if (fabs(sorted_coords[0][1] - sorted_coords[2][1]) < VP_EPS) {
00362         /* input image is degenerate (transforms to a horizontal line) */
00363         lft[0][0] = sorted_coords[0][0];
00364         lft[0][1] = sorted_coords[0][1];
00365         lft[1][0] = sorted_coords[0][0];
00366         lft[1][1] = sorted_coords[0][1];
00367         lft[2][0] = sorted_coords[0][0];
00368         lft[2][1] = sorted_coords[0][1];
00369         rgt[0][0] = sorted_coords[2][0];
00370         rgt[0][1] = sorted_coords[2][1];
00371         rgt[1][0] = sorted_coords[2][0];
00372         rgt[1][1] = sorted_coords[2][1];
00373         rgt[2][0] = sorted_coords[2][0];
00374         rgt[2][1] = sorted_coords[2][1];
00375         return;
00376     }
00377     xmid = sorted_coords[0][0] + (coords[index][1] - sorted_coords[0][1]) *
00378            (sorted_coords[2][0] - sorted_coords[0][0]) /
00379            (sorted_coords[2][1] - sorted_coords[0][1]);
00380     if (coords[index][0] < xmid) {
00381         sorted_coords[1][0] = coords[index][0];
00382         sorted_coords[1][1] = coords[index][1];
00383         sorted_coords[3][0] = coords[(index+2)%4][0];
00384         sorted_coords[3][1] = coords[(index+2)%4][1];
00385     } else {
00386         sorted_coords[1][0] = coords[(index+2)%4][0];
00387         sorted_coords[1][1] = coords[(index+2)%4][1];
00388         sorted_coords[3][0] = coords[index][0];
00389         sorted_coords[3][1] = coords[index][1];
00390     }
00391 
00392     /* store the results in the output array */
00393     lft[0][0] = sorted_coords[0][0];
00394     lft[0][1] = sorted_coords[0][1];
00395     lft[1][0] = sorted_coords[1][0];
00396     lft[1][1] = sorted_coords[1][1];
00397     if (sorted_coords[1][1] == sorted_coords[2][1]) {
00398         lft[2][0] = sorted_coords[1][0];
00399         lft[2][1] = sorted_coords[1][1];
00400     } else {
00401         lft[2][0] = sorted_coords[2][0];
00402         lft[2][1] = sorted_coords[2][1];
00403     }
00404     if (sorted_coords[0][1] == sorted_coords[3][1]) {
00405         rgt[0][0] = sorted_coords[3][0];
00406         rgt[0][1] = sorted_coords[3][1];
00407         rgt[1][0] = sorted_coords[2][0];
00408         rgt[1][1] = sorted_coords[2][1];
00409         rgt[2][0] = sorted_coords[2][0];
00410         rgt[2][1] = sorted_coords[2][1];
00411     } else {
00412         rgt[0][0] = sorted_coords[0][0];
00413         rgt[0][1] = sorted_coords[0][1];
00414         rgt[1][0] = sorted_coords[3][0];
00415         rgt[1][1] = sorted_coords[3][1];
00416         rgt[2][0] = sorted_coords[2][0];
00417         rgt[2][1] = sorted_coords[2][1];
00418     }
00419 }