1 /*
2 * vp_warp.c
3 *
4 * Support routines for 1-pass image warping.
5 *
6 * Copyright (c) 1994 The Board of Trustees of The Leland Stanford
7 * Junior University. All rights reserved.
8 *
9 * Permission to use, copy, modify and distribute this software and its
10 * documentation for any purpose is hereby granted without fee, provided
11 * that the above copyright notice and this permission notice appear in
12 * all copies of this software and that you do not sell the software.
13 * Commercial licensing is available by contacting the author.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS" AND WITHOUT WARRANTY OF ANY KIND,
16 * EXPRESS, IMPLIED OR OTHERWISE, INCLUDING WITHOUT LIMITATION, ANY
17 * WARRANTY OF MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
18 *
19 * Author:
20 * Phil Lacroute
21 * Computer Systems Laboratory
22 * Electrical Engineering Dept.
23 * Stanford University
24 */
25
26 /*
27 * $Date: 1994/12/30 23:52:38 $
28 * $Revision: 1.26 $
29 */
30
31 #include "vp_global.h"
32
33 /* weights for bilirp filter; index with (yfrac, xfrac, tap number) */
34 float VPBilirpWeight[WARP_WEIGHT_ENTRIES][WARP_WEIGHT_ENTRIES][4];
35 static int BilirpWeightsReady = 0; /* true if weight table initialized */
36
37 static void OrderCoords ANSI_ARGS((double coords[4][2], double lft[3][2],
38 double rgt[3][2]));
39
40 /*
41 * VPComputeWarpTables
42 *
43 * Precompute lookup tables for fast filter convolution during warping.
44 */
45
46 void
VPComputeWarpTables()47 VPComputeWarpTables()
48 {
49 float *wptr; /* pointer into weight table */
50
51 int x, y;
52 double in_x, in_y;
53
54 if (BilirpWeightsReady)
55 return;
56
57 #ifdef MEMSPY
58 bin_init(BinNumber(__LINE__, __FILE__, "VPBilirpWeight"), -1, -1,
59 VPBilirpWeight, sizeof(VPBilirpWeight), "VPBilirpWeight");
60 #endif
61
62 wptr = &VPBilirpWeight[0][0][0];
63 for (y = 0; y < WARP_WEIGHT_ENTRIES; y++) {
64 in_y = (double)y / (WARP_WEIGHT_ENTRIES-1);
65 for (x = 0; x < WARP_WEIGHT_ENTRIES; x++) {
66 in_x = (double)x / (WARP_WEIGHT_ENTRIES-1);
67 *wptr++ = (1. - in_x)*(1. - in_y);
68 *wptr++ = in_x * (1. - in_y);
69 *wptr++ = (1. - in_x) * in_y;
70 *wptr++ = 1. - wptr[-1] - wptr[-2] - wptr[-3];
71 }
72 }
73
74 BilirpWeightsReady = 1;
75 }
76
77 /*
78 * VPAffineComputeImageOverlap
79 *
80 * Compute the intersection of the intermediate image and the final image.
81 * This intersection is broken up into a series of trapezoids whose
82 * bases are parallel to the scanlines of the final image (for ease
83 * of scan conversion). Two sets of trapezoids are computed: one set
84 * gives the region of the final image for which the filter kernel is
85 * completely contained within the intermedaite image. The second set of
86 * trapezoids gives the region of the final image for which the filter
87 * kernel overlaps at least some of the intermedaite image, but may or may not
88 * be completely contained. Any region of the final image which is not
89 * within the second set of trapezoids is not affected by the intermediate
90 * image at all and should be set to the background color.
91 */
92
93 void
VPAffineImageOverlap(in_width,in_height,out_width,out_height,warp_matrix,filter_width,full_overlap,part_overlap)94 VPAffineImageOverlap(in_width, in_height, out_width, out_height,
95 warp_matrix, filter_width, full_overlap, part_overlap)
96 int in_width, in_height; /* input (intermediate) image size */
97 int out_width, out_height; /* output (final) image size */
98 vpMatrix3 warp_matrix; /* affine transformation from input image
99 to output image */
100 double filter_width; /* filter kernel width in input image
101 coordinates */
102 Trapezoid full_overlap[9]; /* portion of the intersection with full
103 filter kernel overlap */
104 Trapezoid part_overlap[9]; /* portion of the intersection with partial
105 filter kernel overlap */
106 {
107 double int_lft[3][2]; /* coordinates bounding the left side of the
108 full_overlap region in output coordinates,
109 sorted from top (lowest y) to bottom */
110 double int_rgt[3][2]; /* right side of full_overlap region */
111 double ext_lft[3][2]; /* left side of part_overlap region */
112 double ext_rgt[3][2]; /* right side of part_overlap region */
113 double coords[4][2];
114 double inset;
115 int ilft, irgt, elft, ergt;
116 int region;
117 int lasty, nexty, y;
118
119 /* compute the bounding box of the full_overlap region and store it
120 in int_lft and int_rgt */
121 inset = -1.0 + filter_width / 2.0 + 1.0e-5;
122 coords[0][0] = warp_matrix[0][0] * inset +
123 warp_matrix[0][1] * inset +
124 warp_matrix[0][2];
125 coords[0][1] = warp_matrix[1][0] * inset +
126 warp_matrix[1][1] * inset +
127 warp_matrix[1][2];
128 coords[1][0] = warp_matrix[0][0] * (in_width - 1 - inset) +
129 warp_matrix[0][1] * inset +
130 warp_matrix[0][2];
131 coords[1][1] = warp_matrix[1][0] * (in_width - 1 - inset) +
132 warp_matrix[1][1] * inset +
133 warp_matrix[1][2];
134 coords[2][0] = warp_matrix[0][0] * (in_width - 1 - inset) +
135 warp_matrix[0][1] * (in_height - 1 - inset) +
136 warp_matrix[0][2];
137 coords[2][1] = warp_matrix[1][0] * (in_width - 1 - inset) +
138 warp_matrix[1][1] * (in_height - 1 - inset) +
139 warp_matrix[1][2];
140 coords[3][0] = warp_matrix[0][0] * inset +
141 warp_matrix[0][1] * (in_height - 1 - inset) +
142 warp_matrix[0][2];
143 coords[3][1] = warp_matrix[1][0] * inset +
144 warp_matrix[1][1] * (in_height - 1 - inset) +
145 warp_matrix[1][2];
146 OrderCoords(coords, int_lft, int_rgt);
147
148 /* compute the bounding box of the part_overlap region and store it
149 in int_lft and int_rgt */
150 inset = -filter_width / 2.0;
151 coords[0][0] = warp_matrix[0][0] * inset +
152 warp_matrix[0][1] * inset +
153 warp_matrix[0][2];
154 coords[0][1] = warp_matrix[1][0] * inset +
155 warp_matrix[1][1] * inset +
156 warp_matrix[1][2];
157 coords[1][0] = warp_matrix[0][0] * (in_width - 1 - inset) +
158 warp_matrix[0][1] * inset +
159 warp_matrix[0][2];
160 coords[1][1] = warp_matrix[1][0] * (in_width - 1 - inset) +
161 warp_matrix[1][1] * inset +
162 warp_matrix[1][2];
163 coords[2][0] = warp_matrix[0][0] * (in_width - 1 - inset) +
164 warp_matrix[0][1] * (in_height - 1 - inset) +
165 warp_matrix[0][2];
166 coords[2][1] = warp_matrix[1][0] * (in_width - 1 - inset) +
167 warp_matrix[1][1] * (in_height - 1 - inset) +
168 warp_matrix[1][2];
169 coords[3][0] = warp_matrix[0][0] * inset +
170 warp_matrix[0][1] * (in_height - 1 - inset) +
171 warp_matrix[0][2];
172 coords[3][1] = warp_matrix[1][0] * inset +
173 warp_matrix[1][1] * (in_height - 1 - inset) +
174 warp_matrix[1][2];
175 OrderCoords(coords, ext_lft, ext_rgt);
176
177 for (ilft = 0; ilft < 3 && int_lft[ilft][1] <= 0.; ilft++);
178 for (irgt = 0; irgt < 3 && int_rgt[irgt][1] <= 0.; irgt++);
179 for (elft = 0; elft < 3 && ext_lft[elft][1] <= 0.; elft++);
180 for (ergt = 0; ergt < 3 && ext_rgt[ergt][1] <= 0.; ergt++);
181 region = 0;
182 lasty = -1;
183 while (lasty < out_height-1) {
184 ASSERT(region < 9);
185
186 /* find nexty */
187 nexty = out_height - 1;
188 if (ilft < 3) {
189 y = (int)floor(int_lft[ilft][1]);
190 if (nexty > y)
191 nexty = y;
192 }
193 if (irgt < 3) {
194 y = (int)floor(int_rgt[irgt][1]);
195 if (nexty > y)
196 nexty = y;
197 }
198 if (elft < 3) {
199 y = (int)floor(ext_lft[elft][1]);
200 if (nexty > y)
201 nexty = y;
202 }
203 if (ergt < 3) {
204 y = (int)floor(ext_rgt[ergt][1]);
205 if (nexty > y)
206 nexty = y;
207 }
208 ASSERT((ilft == 0 && (int)floor(int_lft[0][1]) >= nexty) ||
209 (ilft == 3 && (int)floor(int_lft[2][1]) <= lasty) ||
210 (((int)floor(int_lft[ilft-1][1]) <= lasty || lasty == -1)
211 && (int)floor(int_lft[ilft][1]) >= nexty));
212 ASSERT((irgt == 0 && (int)floor(int_rgt[0][1]) >= nexty) ||
213 (irgt == 3 && (int)floor(int_rgt[2][1]) <= lasty) ||
214 (((int)floor(int_rgt[irgt-1][1]) <= lasty || lasty == -1)
215 && (int)floor(int_rgt[irgt][1]) >= nexty));
216 ASSERT((elft == 0 && (int)floor(ext_lft[0][1]) >= nexty) ||
217 (elft == 3 && (int)floor(ext_lft[2][1]) <= lasty) ||
218 (((int)floor(ext_lft[elft-1][1]) <= lasty || lasty == -1)
219 && (int)floor(ext_lft[elft][1]) >= nexty));
220 ASSERT((ergt == 0 && (int)floor(ext_rgt[0][1]) >= nexty) ||
221 (ergt == 3 && (int)floor(ext_rgt[2][1]) <= lasty) ||
222 (((int)floor(ext_rgt[ergt-1][1]) <= lasty || lasty == -1)
223 && (int)floor(ext_rgt[ergt][1]) >= nexty));
224 full_overlap[region].miny = lasty + 1;
225 full_overlap[region].maxy = nexty;
226 part_overlap[region].miny = lasty + 1;
227 part_overlap[region].maxy = nexty;
228 if (ilft == 0 || ilft == 3) {
229 /* this trapezoid does not intersect full_overlap */
230 full_overlap[region].x_top_lft = 0;
231 full_overlap[region].x_top_rgt = -1;
232 full_overlap[region].x_incr_lft = 0;
233 full_overlap[region].x_incr_rgt = 0;
234 } else {
235 full_overlap[region].x_incr_lft =
236 (int_lft[ilft][0] - int_lft[ilft-1][0]) /
237 (int_lft[ilft][1] - int_lft[ilft-1][1]);
238 full_overlap[region].x_top_lft =
239 int_lft[ilft-1][0] + (lasty+1 - int_lft[ilft-1][1]) *
240 full_overlap[region].x_incr_lft;
241 full_overlap[region].x_incr_rgt =
242 (int_rgt[irgt][0] - int_rgt[irgt-1][0]) /
243 (int_rgt[irgt][1] - int_rgt[irgt-1][1]);
244 full_overlap[region].x_top_rgt =
245 int_rgt[irgt-1][0] + (lasty+1 - int_rgt[irgt-1][1]) *
246 full_overlap[region].x_incr_rgt;
247 }
248 if (elft == 0 || elft == 3) {
249 /* this trapezoid does not intersect part_overlap */
250 part_overlap[region].x_top_lft = 0;
251 part_overlap[region].x_top_rgt = -1;
252 part_overlap[region].x_incr_lft = 0;
253 part_overlap[region].x_incr_rgt = 0;
254 } else {
255 part_overlap[region].x_incr_lft =
256 (ext_lft[elft][0] - ext_lft[elft-1][0]) /
257 (ext_lft[elft][1] - ext_lft[elft-1][1]);
258 part_overlap[region].x_top_lft =
259 ext_lft[elft-1][0] + (lasty+1 - ext_lft[elft-1][1]) *
260 part_overlap[region].x_incr_lft;
261 part_overlap[region].x_incr_rgt =
262 (ext_rgt[ergt][0] - ext_rgt[ergt-1][0]) /
263 (ext_rgt[ergt][1] - ext_rgt[ergt-1][1]);
264 part_overlap[region].x_top_rgt =
265 ext_rgt[ergt-1][0] + (lasty+1 - ext_rgt[ergt-1][1]) *
266 part_overlap[region].x_incr_rgt;
267 }
268 ASSERT(!(full_overlap[region].x_top_lft <=
269 full_overlap[region].x_top_rgt &&
270 part_overlap[region].x_top_lft >
271 part_overlap[region].x_top_rgt));
272 for (; ilft < 3 && (int)floor(int_lft[ilft][1]) <= nexty; ilft++);
273 for (; irgt < 3 && (int)floor(int_rgt[irgt][1]) <= nexty; irgt++);
274 for (; elft < 3 && (int)floor(ext_lft[elft][1]) <= nexty; elft++);
275 for (; ergt < 3 && (int)floor(ext_rgt[ergt][1]) <= nexty; ergt++);
276 region++;
277 lasty = nexty;
278 }
279 for (; region < 9; region++) {
280 full_overlap[region].miny = out_height;
281 full_overlap[region].maxy = out_height;
282 part_overlap[region].miny = out_height;
283 part_overlap[region].maxy = out_height;
284 full_overlap[region].x_top_lft = 0;
285 full_overlap[region].x_top_rgt = -1;
286 full_overlap[region].x_incr_lft = 0;
287 full_overlap[region].x_incr_rgt = 0;
288 part_overlap[region].x_top_lft = 0;
289 part_overlap[region].x_top_rgt = -1;
290 part_overlap[region].x_incr_lft = 0;
291 part_overlap[region].x_incr_rgt = 0;
292 }
293
294 #ifdef DEBUG_OVERLAP
295 for (region = 0; region < 9; region++) {
296 printf("region %d: y = %d to %d, [%d+%g [%d+%g %d+%g] %d+%g]\n",
297 region,
298 full_overlap[region].miny,
299 full_overlap[region].maxy,
300 (int)part_overlap[region].x_top_lft,
301 part_overlap[region].x_incr_lft,
302 (int)full_overlap[region].x_top_lft,
303 full_overlap[region].x_incr_lft,
304 (int)full_overlap[region].x_top_rgt,
305 full_overlap[region].x_incr_rgt,
306 (int)part_overlap[region].x_top_rgt,
307 part_overlap[region].x_incr_rgt);
308 }
309 #endif
310 }
311
312 /*
313 * OrderCoords
314 *
315 * Sort an array of coordinates.
316 */
317
318 static void
OrderCoords(coords,lft,rgt)319 OrderCoords(coords, lft, rgt)
320 double coords[4][2]; /* inputs */
321 double lft[3][2]; /* outputs */
322 double rgt[3][2];
323 {
324 int index;
325 double swap_buf;
326 double sorted_coords[4][2];
327 double xmid;
328
329 /* sort the coordinates as follows:
330 coords[0]: smallest y value; if there is a tie, then smallest x value
331 to break the tie; this is the top or top left corner
332 coords[1]: next coordinate in CCW order; this is the left or bottom
333 left corner
334 coords[2]: next coordinate in CCW order; this is the bottom or
335 bottom right corner
336 coords[3]: next coordinate in CCW order; this is the right or
337 top right corner */
338
339 /* search for coords[0] */
340 index = 0;
341 if (coords[1][1] < coords[index][1] ||
342 (coords[1][1] == coords[index][1] && coords[1][0] < coords[index][0]))
343 index = 1;
344 if (coords[2][1] < coords[index][1] ||
345 (coords[2][1] == coords[index][1] && coords[2][0] < coords[index][0]))
346 index = 2;
347 if (coords[3][1] < coords[index][1] ||
348 (coords[3][1] == coords[index][1] && coords[3][0] < coords[index][0]))
349 index = 3;
350 sorted_coords[0][0] = coords[index][0];
351 sorted_coords[0][1] = coords[index][1];
352
353 /* coords[2] is the opposite corner */
354 sorted_coords[2][0] = coords[(index+2)%4][0];
355 sorted_coords[2][1] = coords[(index+2)%4][1];
356
357 /* determine which of the remaining two coordinates is to the left
358 of the line joining coords[0] and coords[2]; this coordinate
359 is coords[1], and the final remaining coordinate is coords[3] */
360 index = (index + 1) % 4;
361 if (fabs(sorted_coords[0][1] - sorted_coords[2][1]) < VP_EPS) {
362 /* input image is degenerate (transforms to a horizontal line) */
363 lft[0][0] = sorted_coords[0][0];
364 lft[0][1] = sorted_coords[0][1];
365 lft[1][0] = sorted_coords[0][0];
366 lft[1][1] = sorted_coords[0][1];
367 lft[2][0] = sorted_coords[0][0];
368 lft[2][1] = sorted_coords[0][1];
369 rgt[0][0] = sorted_coords[2][0];
370 rgt[0][1] = sorted_coords[2][1];
371 rgt[1][0] = sorted_coords[2][0];
372 rgt[1][1] = sorted_coords[2][1];
373 rgt[2][0] = sorted_coords[2][0];
374 rgt[2][1] = sorted_coords[2][1];
375 return;
376 }
377 xmid = sorted_coords[0][0] + (coords[index][1] - sorted_coords[0][1]) *
378 (sorted_coords[2][0] - sorted_coords[0][0]) /
379 (sorted_coords[2][1] - sorted_coords[0][1]);
380 if (coords[index][0] < xmid) {
381 sorted_coords[1][0] = coords[index][0];
382 sorted_coords[1][1] = coords[index][1];
383 sorted_coords[3][0] = coords[(index+2)%4][0];
384 sorted_coords[3][1] = coords[(index+2)%4][1];
385 } else {
386 sorted_coords[1][0] = coords[(index+2)%4][0];
387 sorted_coords[1][1] = coords[(index+2)%4][1];
388 sorted_coords[3][0] = coords[index][0];
389 sorted_coords[3][1] = coords[index][1];
390 }
391
392 /* store the results in the output array */
393 lft[0][0] = sorted_coords[0][0];
394 lft[0][1] = sorted_coords[0][1];
395 lft[1][0] = sorted_coords[1][0];
396 lft[1][1] = sorted_coords[1][1];
397 if (sorted_coords[1][1] == sorted_coords[2][1]) {
398 lft[2][0] = sorted_coords[1][0];
399 lft[2][1] = sorted_coords[1][1];
400 } else {
401 lft[2][0] = sorted_coords[2][0];
402 lft[2][1] = sorted_coords[2][1];
403 }
404 if (sorted_coords[0][1] == sorted_coords[3][1]) {
405 rgt[0][0] = sorted_coords[3][0];
406 rgt[0][1] = sorted_coords[3][1];
407 rgt[1][0] = sorted_coords[2][0];
408 rgt[1][1] = sorted_coords[2][1];
409 rgt[2][0] = sorted_coords[2][0];
410 rgt[2][1] = sorted_coords[2][1];
411 } else {
412 rgt[0][0] = sorted_coords[0][0];
413 rgt[0][1] = sorted_coords[0][1];
414 rgt[1][0] = sorted_coords[3][0];
415 rgt[1][1] = sorted_coords[3][1];
416 rgt[2][0] = sorted_coords[2][0];
417 rgt[2][1] = sorted_coords[2][1];
418 }
419 }
420