summaryrefslogtreecommitdiffstats
path: root/SD-VBS/benchmarks/texture_synthesis/src/matlab/MEX/edges.c
diff options
context:
space:
mode:
Diffstat (limited to 'SD-VBS/benchmarks/texture_synthesis/src/matlab/MEX/edges.c')
-rwxr-xr-xSD-VBS/benchmarks/texture_synthesis/src/matlab/MEX/edges.c647
1 files changed, 647 insertions, 0 deletions
diff --git a/SD-VBS/benchmarks/texture_synthesis/src/matlab/MEX/edges.c b/SD-VBS/benchmarks/texture_synthesis/src/matlab/MEX/edges.c
new file mode 100755
index 0000000..98b377d
--- /dev/null
+++ b/SD-VBS/benchmarks/texture_synthesis/src/matlab/MEX/edges.c
@@ -0,0 +1,647 @@
1/*
2;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
3;;; File: edges.c
4;;; Author: Eero Simoncelli
5;;; Description: Boundary handling routines for use with convolve.c
6;;; Creation Date: Spring 1987.
7;;; MODIFIED, 6/96, to operate on double float arrays.
8;;; MODIFIED by dgp, 4/1/97, to support THINK C.
9;;; MODIFIED, 8/97: reflect1, reflect2, repeat, extend upgraded to
10;;; work properly for non-symmetric filters. Added qreflect2 to handle
11;;; even-length QMF's which broke under the reflect2 modification.
12;;; ----------------------------------------------------------------
13;;; Object-Based Vision and Image Understanding System (OBVIUS),
14;;; Copyright 1988, Vision Science Group, Media Laboratory,
15;;; Massachusetts Institute of Technology.
16;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
17*/
18
19/* This file contains functions which determine how edges are to be
20handled when performing convolutions of images with linear filters.
21Any edge handling function which is local and linear may be defined,
22except (unfortunately) constants cannot be added. So to treat the
23edges as if the image is surrounded by a gray field, you must paste it
24into a gray image, convolve, and crop it out... The main convolution
25functions are called internal_reduce and internal_expand and are
26defined in the file convolve.c. The idea is that the convolution
27function calls the edge handling function which computes a new filter
28based on the old filter and the distance to the edge of the image.
29For example, reflection is done by reflecting the filter through the
30appropriate axis and summing. Currently defined functions are listed
31below.
32*/
33
34
35#include <stdio.h>
36#include <math.h>
37#include <string.h>
38#include "convolve.h"
39
40#define sgn(a) ( ((a)>0)?1:(((a)<0)?-1:0) )
41#define clip(a,mn,mx) ( ((a)<(mn))?(mn):(((a)>=(mx))?(mx-1):(a)) )
42
43int reflect1(), reflect2(), qreflect2(), repeat(), zero(), Extend(), nocompute();
44int ereflect(), predict();
45
46/* Lookup table matching a descriptive string to the edge-handling function */
47#if !THINK_C
48 static EDGE_HANDLER edge_foos[] =
49 {
50 { "dont-compute", nocompute }, /* zero output for filter touching edge */
51 { "zero", zero }, /* zero outside of image */
52 { "repeat", repeat }, /* repeat edge pixel */
53 { "reflect1", reflect1 }, /* reflect about edge pixels */
54 { "reflect2", reflect2 }, /* reflect image, including edge pixels */
55 { "qreflect2", qreflect2 }, /* reflect image, including edge pixels
56 for even-length QMF decompositions */
57 { "extend", Extend }, /* extend (reflect & invert) */
58 { "ereflect", ereflect }, /* orthogonal QMF reflection */
59 };
60#else
61 /*
62 This is really stupid, but THINK C won't allow initialization of static variables in
63 a code resource with string addresses. So we do it this way.
64 The 68K code for a MATLAB 4 MEX file can only be created by THINK C.
65 However, for MATLAB 5, we'll be able to use Metrowerks CodeWarrior for both 68K and PPC, so this
66 cludge can be dropped when we drop support for MATLAB 4.
67 Denis Pelli, 4/1/97.
68 */
69 static EDGE_HANDLER edge_foos[8];
70
71 void InitializeTable(EDGE_HANDLER edge_foos[])
72 {
73 static int i=0;
74
75 if(i>0) return;
76 edge_foos[i].name="dont-compute";
77 edge_foos[i++].func=nocompute;
78 edge_foos[i].name="zero";
79 edge_foos[i++].func=zero;
80 edge_foos[i].name="repeat";
81 edge_foos[i++].func=repeat;
82 edge_foos[i].name="reflect1";
83 edge_foos[i++].func=reflect1;
84 edge_foos[i].name="reflect2";
85 edge_foos[i++].func=reflect2;
86 edge_foos[i].name="qreflect2";
87 edge_foos[i++].func=qreflect2;
88 edge_foos[i].name="extend";
89 edge_foos[i++].func=Extend;
90 edge_foos[i].name="ereflect";
91 edge_foos[i++].func=ereflect;
92 }
93#endif
94
95/*
96Function looks up an edge handler id string in the structure above, and
97returns the associated function
98*/
99fptr edge_function(char *edges)
100 {
101 int i;
102
103#if THINK_C
104 InitializeTable(edge_foos);
105#endif
106 for (i = 0; i<sizeof(edge_foos)/sizeof(EDGE_HANDLER); i++)
107 if (strcmp(edges,edge_foos[i].name) IS 0)
108 return(edge_foos[i].func);
109 printf("Error: '%s' is not the name of a valid edge-handler!\n",edges);
110 for (i=0; i<sizeof(edge_foos)/sizeof(EDGE_HANDLER); i++)
111 {
112 if (i IS 0) printf(" Options are: ");
113 else printf(", ");
114 printf("%s",edge_foos[i].name);
115 }
116 printf("\n");
117 return(0);
118 }
119
120/*
121---------------- EDGE HANDLER ARGUMENTS ------------------------
122filt - array of filter taps.
123x_dim, y_dim - x and y dimensions of filt.
124x_pos - position of filter relative to the horizontal image edges. Negative
125 values indicate left edge, positive indicate right edge. Zero
126 indicates that the filter is not touching either edge. An absolute
127 value of 1 indicates that the edge tap of the filter is over the
128 edge pixel of the image.
129y_pos - analogous to x_pos.
130result - array where the resulting filter will go. The edge
131 of this filter will be aligned with the image for application...
132r_or_e - equal to one of the two constants EXPAND or REDUCE.
133--------------------------------------------------------------------
134*/
135
136
137/* --------------------------------------------------------------------
138nocompute() - Return zero for values where filter hangs over the edge.
139*/
140
141int nocompute(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
142 register double *filt, *result;
143 register int x_dim;
144 int y_dim, x_pos, y_pos, r_or_e;
145 {
146 register int i;
147 register int size = x_dim*y_dim;
148
149 if ( (x_pos>1) OR (x_pos<-1) OR (y_pos>1) OR (y_pos<-1) )
150 for (i=0; i<size; i++) result[i] = 0.0;
151 else
152 for (i=0; i<size; i++) result[i] = filt[i];
153
154 return(0);
155 }
156
157/* --------------------------------------------------------------------
158zero() - Zero outside of image. Discontinuous, but adds zero energy. */
159
160int zero(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
161 register double *filt, *result;
162 register int x_dim;
163 int y_dim, x_pos, y_pos, r_or_e;
164 {
165 register int y_filt,x_filt, y_res,x_res;
166 int filt_sz = x_dim*y_dim;
167 int x_start = ((x_pos>0)?(x_pos-1):((x_pos<0)?(x_pos+1):0));
168 int y_start = x_dim * ((y_pos>0)?(y_pos-1):((y_pos<0)?(y_pos+1):0));
169 int i;
170
171 for (i=0; i<filt_sz; i++) result[i] = 0.0;
172
173 for (y_filt=0, y_res=y_start;
174 y_filt<filt_sz;
175 y_filt+=x_dim, y_res+=x_dim)
176 if ((y_res >= 0) AND (y_res < filt_sz))
177 for (x_filt=y_filt, x_res=x_start;
178 x_filt<y_filt+x_dim;
179 x_filt++, x_res++)
180 if ((x_res >= 0) AND (x_res < x_dim))
181 result[y_res+x_res] = filt[x_filt];
182 return(0);
183 }
184
185
186/* --------------------------------------------------------------------
187reflect1() - Reflection through the edge pixels. Continuous, but
188discontinuous first derivative. This is the right thing to do if you
189are subsampling by 2, since it maintains parity (even pixels positions
190remain even, odd ones remain odd).
191*/
192
193int reflect1(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
194 register double *filt, *result;
195 register int x_dim;
196 int y_dim, x_pos, y_pos, r_or_e;
197 {
198 int filt_sz = x_dim*y_dim;
199 register int y_filt,x_filt, y_res, x_res;
200 register int x_base = (x_pos>0)?(x_dim-1):0;
201 register int y_base = x_dim * ((y_pos>0)?(y_dim-1):0);
202 int x_overhang = (x_pos>0)?(x_pos-1):((x_pos<0)?(x_pos+1):0);
203 int y_overhang = x_dim * ((y_pos>0)?(y_pos-1):((y_pos<0)?(y_pos+1):0));
204 int i;
205 int mx_pos = (x_pos<0)?(x_dim/2):((x_dim-1)/2);
206 int my_pos = x_dim * ((y_pos<0)?(y_dim/2):((y_dim-1)/2));
207
208 for (i=0; i<filt_sz; i++) result[i] = 0.0;
209
210 if (r_or_e IS REDUCE)
211 for (y_filt=0, y_res=y_overhang-y_base;
212 y_filt<filt_sz;
213 y_filt+=x_dim, y_res+=x_dim)
214 for (x_filt=y_filt, x_res=x_overhang-x_base;
215 x_filt<y_filt+x_dim;
216 x_filt++, x_res++)
217 result[ABS(y_base-ABS(y_res)) + ABS(x_base-ABS(x_res))]
218 += filt[x_filt];
219 else {
220 y_overhang = ABS(y_overhang);
221 x_overhang = ABS(x_overhang);
222 for (y_res=y_base, y_filt = y_base-y_overhang;
223 y_filt > y_base-filt_sz;
224 y_filt-=x_dim, y_res-=x_dim)
225 {
226 for (x_res=x_base, x_filt=x_base-x_overhang;
227 x_filt > x_base-x_dim;
228 x_res--, x_filt--)
229 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
230 if ((x_overhang ISNT mx_pos) AND (x_pos ISNT 0))
231 for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang;
232 x_filt > x_base-x_dim;
233 x_res--, x_filt--)
234 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
235 }
236 if ((y_overhang ISNT my_pos) AND (y_pos ISNT 0))
237 for (y_res=y_base, y_filt = y_base-2*my_pos+y_overhang;
238 y_filt > y_base-filt_sz;
239 y_filt-=x_dim, y_res-=x_dim)
240 {
241 for (x_res=x_base, x_filt=x_base-x_overhang;
242 x_filt > x_base-x_dim;
243 x_res--, x_filt--)
244 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
245 if ((x_overhang ISNT mx_pos) AND (x_pos ISNT 0))
246 for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang;
247 x_filt > x_base-x_dim;
248 x_res--, x_filt--)
249 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
250 }
251 }
252
253 return(0);
254 }
255
256/* --------------------------------------------------------------------
257reflect2() - Reflect image at boundary. The edge pixel is repeated,
258then the next pixel, etc. Continuous, but discontinuous first
259derivative.
260*/
261
262int reflect2(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
263 register double *filt, *result;
264 register int x_dim;
265 int y_dim, x_pos, y_pos, r_or_e;
266 {
267 int filt_sz = x_dim*y_dim;
268 register int y_filt,x_filt, y_res, x_res;
269 register int x_base = (x_pos>0)?(x_dim-1):0;
270 register int y_base = x_dim * ((y_pos>0)?(y_dim-1):0);
271 int x_overhang = (x_pos>0)?(x_pos-1):((x_pos<0)?(x_pos+1):0);
272 int y_overhang = x_dim * ((y_pos>0)?(y_pos-1):((y_pos<0)?(y_pos+1):0));
273 int i;
274 int mx_pos = (x_pos<0)?(x_dim/2):((x_dim-1)/2);
275 int my_pos = x_dim * ((y_pos<0)?(y_dim/2):((y_dim-1)/2));
276
277 for (i=0; i<filt_sz; i++) result[i] = 0.0;
278
279 if (r_or_e IS REDUCE)
280 for (y_filt=0, y_res = (y_overhang-y_base) - ((y_pos>0)?x_dim:0);
281 y_filt<filt_sz;
282 y_filt+=x_dim, y_res+=x_dim)
283 {
284 if (y_res IS 0) y_res+=x_dim;
285 for (x_filt=y_filt, x_res = (x_overhang-x_base) - ((x_pos>0)?1:0);
286 x_filt<y_filt+x_dim;
287 x_filt++, x_res++)
288 {
289 if (x_res IS 0) x_res++;
290 result[ABS(y_base-ABS(y_res)+x_dim) + ABS(x_base-ABS(x_res)+1)]
291 += filt[x_filt];
292 }
293 }
294 else {
295 y_overhang = ABS(y_overhang);
296 x_overhang = ABS(x_overhang);
297 for (y_res=y_base, y_filt = y_base-y_overhang;
298 y_filt > y_base-filt_sz;
299 y_filt-=x_dim, y_res-=x_dim)
300 {
301 for (x_res=x_base, x_filt=x_base-x_overhang;
302 x_filt > x_base-x_dim;
303 x_res--, x_filt--)
304 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
305 if (x_pos ISNT 0)
306 for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang-1;
307 x_filt > x_base-x_dim;
308 x_res--, x_filt--)
309 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
310 }
311 if (y_pos ISNT 0)
312 for (y_res=y_base, y_filt = y_base-2*my_pos+y_overhang-x_dim;
313 y_filt > y_base-filt_sz;
314 y_filt-=x_dim, y_res-=x_dim)
315 {
316 for (x_res=x_base, x_filt=x_base-x_overhang;
317 x_filt > x_base-x_dim;
318 x_res--, x_filt--)
319 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
320 if (x_pos ISNT 0)
321 for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang-1;
322 x_filt > x_base-x_dim;
323 x_res--, x_filt--)
324 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
325 }
326 }
327
328 return(0);
329 }
330
331
332/* --------------------------------------------------------------------
333qreflect2() - Modified version of reflect2 that works properly for
334even-length QMF filters.
335*/
336
337int qreflect2(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
338 double *filt, *result;
339 int x_dim, y_dim, x_pos, y_pos, r_or_e;
340 {
341 reflect2(filt,x_dim,y_dim,x_pos,y_pos,result,0);
342 return(0);
343 }
344
345/* --------------------------------------------------------------------
346repeat() - repeat edge pixel. Continuous, with discontinuous first
347derivative.
348*/
349
350int repeat(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
351 register double *filt, *result;
352 register int x_dim;
353 int y_dim, x_pos, y_pos, r_or_e;
354 {
355 register int y_filt,x_filt, y_res,x_res, y_tmp, x_tmp;
356 register int x_base = (x_pos>0)?(x_dim-1):0;
357 register int y_base = x_dim * ((y_pos>0)?(y_dim-1):0);
358 int x_overhang = ((x_pos>0)?(x_pos-1):((x_pos<0)?(x_pos+1):0));
359 int y_overhang = x_dim * ((y_pos>0)?(y_pos-1):((y_pos<0)?(y_pos+1):0));
360 int filt_sz = x_dim*y_dim;
361 int mx_pos = (x_dim/2);
362 int my_pos = x_dim * (y_dim/2);
363 int i;
364
365 for (i=0; i<filt_sz; i++) result[i] = 0.0;
366
367 if (r_or_e IS REDUCE)
368 for (y_filt=0, y_res=y_overhang;
369 y_filt<filt_sz;
370 y_filt+=x_dim, y_res+=x_dim)
371 for (x_filt=y_filt, x_res=x_overhang;
372 x_filt<y_filt+x_dim;
373 x_filt++, x_res++)
374 /* Clip the index: */
375 result[((y_res>=0)?((y_res<filt_sz)?y_res:(filt_sz-x_dim)):0)
376 + ((x_res>=0)?((x_res<x_dim)?x_res:(x_dim-1)):0)]
377 += filt[x_filt];
378 else {
379 if ((y_base-y_overhang) ISNT my_pos)
380 for (y_res=y_base, y_filt=y_base-ABS(y_overhang);
381 y_filt > y_base-filt_sz;
382 y_filt-=x_dim, y_res-=x_dim)
383 if ((x_base-x_overhang) ISNT mx_pos)
384 for (x_res=x_base, x_filt=x_base-ABS(x_overhang);
385 x_filt > x_base-x_dim;
386 x_res--, x_filt--)
387 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
388 else /* ((x_base-x_overhang) IS mx_pos) */
389 for (x_res=x_base, x_filt=x_base-ABS(x_overhang);
390 x_filt > x_base-x_dim;
391 x_filt--, x_res--)
392 for(x_tmp=x_filt; x_tmp > x_base-x_dim; x_tmp--)
393 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_tmp)];
394 else /* ((y_base-y_overhang) IS my_pos) */
395 for (y_res=y_base, y_filt=y_base-ABS(y_overhang);
396 y_filt > y_base-filt_sz;
397 y_filt-=x_dim, y_res-=x_dim)
398 for (y_tmp=y_filt; y_tmp > y_base-filt_sz; y_tmp-=x_dim)
399 if ((x_base-x_overhang) ISNT mx_pos)
400 for (x_res=x_base, x_filt=x_base-ABS(x_overhang);
401 x_filt > x_base-x_dim;
402 x_filt--, x_res--)
403 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_tmp)+ABS(x_filt)];
404 else /* ((x_base-x_overhang) IS mx_pos) */
405 for (x_res=x_base, x_filt=x_base-ABS(x_overhang);
406 x_filt > x_base-x_dim;
407 x_filt--, x_res--)
408 for (x_tmp=x_filt; x_tmp > x_base-x_dim; x_tmp--)
409 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_tmp)+ABS(x_tmp)];
410 } /* End, if r_or_e IS EXPAND */
411
412 return(0);
413 }
414
415/* --------------------------------------------------------------------
416extend() - Extend image by reflecting and inverting about edge pixel
417value. Maintains continuity in intensity AND first derivative (but
418not higher derivs).
419*/
420
421int Extend(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
422 register double *filt, *result;
423 register int x_dim;
424 int y_dim, x_pos, y_pos, r_or_e;
425 {
426 int filt_sz = x_dim*y_dim;
427 register int y_filt,x_filt, y_res,x_res, y_tmp, x_tmp;
428 register int x_base = (x_pos>0)?(x_dim-1):0;
429 register int y_base = x_dim * ((y_pos>0)?(y_dim-1):0);
430 int x_overhang = (x_pos>0)?(x_pos-1):((x_pos<0)?(x_pos+1):0);
431 int y_overhang = x_dim * ((y_pos>0)?(y_pos-1):((y_pos<0)?(y_pos+1):0));
432 int mx_pos = (x_pos<0)?(x_dim/2):((x_dim-1)/2);
433 int my_pos = x_dim * ((y_pos<0)?(y_dim/2):((y_dim-1)/2));
434 int i;
435
436 for (i=0; i<filt_sz; i++) result[i] = 0.0;
437
438 /* Modeled on reflect1() */
439 if (r_or_e IS REDUCE)
440 for (y_filt=0, y_res=y_overhang;
441 y_filt<filt_sz;
442 y_filt+=x_dim, y_res+=x_dim)
443 if ((y_res>=0) AND (y_res<filt_sz))
444 for (x_filt=y_filt, x_res=x_overhang;
445 x_filt<y_filt+x_dim;
446 x_filt++, x_res++)
447 if ((x_res>=0) AND (x_res<x_dim))
448 result[y_res+x_res] += filt[x_filt];
449 else
450 {
451 result[y_res+ABS(x_base-ABS(x_res-x_base))] -= filt[x_filt];
452 result[y_res+x_base] += 2*filt[x_filt];
453 }
454 else
455 for (x_filt=y_filt, x_res=x_overhang;
456 x_filt<y_filt+x_dim;
457 x_filt++, x_res++)
458 if ((x_res>=0) AND (x_res<x_dim))
459 {
460 result[ABS(y_base-ABS(y_res-y_base))+x_res] -= filt[x_filt];
461 result[y_base+x_res] += 2*filt[x_filt];
462 }
463 else
464 {
465 result[ABS(y_base-ABS(y_res-y_base))+ABS(x_base-ABS(x_res-x_base))]
466 -= filt[x_filt];
467 result[y_base+x_base] += 2*filt[x_filt];
468 }
469 else { /* r_or_e ISNT REDUCE */
470 y_overhang = ABS(y_overhang);
471 x_overhang = ABS(x_overhang);
472 for (y_res=y_base, y_filt = y_base-y_overhang;
473 y_filt > y_base-filt_sz;
474 y_filt-=x_dim, y_res-=x_dim)
475 {
476 for (x_res=x_base, x_filt=x_base-x_overhang;
477 x_filt > x_base-x_dim;
478 x_res--, x_filt--)
479 result[ABS(y_res)+ABS(x_res)] += filt[ABS(y_filt)+ABS(x_filt)];
480 if (x_pos ISNT 0)
481 if (x_overhang ISNT mx_pos)
482 for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang;
483 x_filt > x_base-x_dim;
484 x_res--, x_filt--)
485 result[ABS(y_res)+ABS(x_res)] -= filt[ABS(y_filt)+ABS(x_filt)];
486 else /* x_overhang IS mx_pos */
487 for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang-1;
488 x_filt > x_base-x_dim;
489 x_res--, x_filt--)
490 for (x_tmp=x_filt; x_tmp > x_base-x_dim; x_tmp--)
491 result[ABS(y_res)+ABS(x_res)] += 2*filt[ABS(y_filt)+ABS(x_tmp)];
492 }
493 if (y_pos ISNT 0)
494 if (y_overhang ISNT my_pos)
495 for (y_res=y_base, y_filt = y_base-2*my_pos+y_overhang;
496 y_filt > y_base-filt_sz;
497 y_filt-=x_dim, y_res-=x_dim)
498 {
499 for (x_res=x_base, x_filt=x_base-x_overhang;
500 x_filt > x_base-x_dim;
501 x_res--, x_filt--)
502 result[ABS(y_res)+ABS(x_res)] -= filt[ABS(y_filt)+ABS(x_filt)];
503 if ((x_pos ISNT 0) AND (x_overhang ISNT mx_pos))
504 for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang;
505 x_filt > x_base-x_dim;
506 x_res--, x_filt--)
507 result[ABS(y_res)+ABS(x_res)] -= filt[ABS(y_filt)+ABS(x_filt)];
508 }
509 else /* y_overhang IS my_pos */
510 for (y_res=y_base, y_filt = y_base-2*my_pos+y_overhang-x_dim;
511 y_filt > y_base-filt_sz;
512 y_res-=x_dim, y_filt-=x_dim)
513 for (y_tmp=y_filt; y_tmp > y_base-filt_sz; y_tmp-=x_dim)
514 {
515 for (x_res=x_base, x_filt=x_base-x_overhang;
516 x_filt > x_base-x_dim;
517 x_res--, x_filt--)
518 result[ABS(y_res)+ABS(x_res)] += 2*filt[ABS(y_tmp)+ABS(x_filt)];
519 if ((x_pos ISNT 0) AND (x_overhang IS mx_pos))
520 for (x_res=x_base, x_filt=x_base-2*mx_pos+x_overhang-1;
521 x_filt > x_base-x_dim;
522 x_res--, x_filt--)
523 for (x_tmp=x_filt; x_tmp > x_base-x_dim; x_tmp--)
524 result[ABS(y_res)+ABS(x_res)] += 2*filt[ABS(y_tmp)+ABS(x_tmp)];
525 }
526 } /* r_or_e ISNT REDUCE */
527
528 return(0);
529 }
530
531/* --------------------------------------------------------------------
532predict() - Simple prediction. Like zero, but multiplies the result
533by the reciprocal of the percentage of filter being used. (i.e. if
53450% of the filter is hanging over the edge of the image, multiply the
535taps being used by 2). */
536
537int predict(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
538 register double *filt, *result;
539 register int x_dim;
540 int y_dim, x_pos, y_pos, r_or_e;
541 {
542 register int y_filt,x_filt, y_res,x_res;
543 register double taps_used = 0.0; /* int *** */
544 register double fraction = 0.0;
545 int filt_sz = x_dim*y_dim;
546 int x_start = ((x_pos>0)?(x_pos-1):((x_pos<0)?(x_pos+1):0));
547 int y_start = x_dim * ((y_pos>0)?(y_pos-1):((y_pos<0)?(y_pos+1):0));
548 int i;
549
550 for (i=0; i<filt_sz; i++) result[i] = 0.0;
551
552 for (y_filt=0, y_res=y_start;
553 y_filt<filt_sz;
554 y_filt+=x_dim, y_res+=x_dim)
555 if ((y_res >= 0) AND (y_res < filt_sz))
556 for (x_filt=y_filt, x_res=x_start;
557 x_filt<y_filt+x_dim;
558 x_filt++, x_res++)
559 if ((x_res >= 0) AND (x_res < x_dim))
560 {
561 result[y_res+x_res] = filt[x_filt];
562 taps_used += ABS(filt[x_filt]);
563 }
564
565 if (r_or_e IS REDUCE)
566 {
567 /* fraction = ( (double) filt_sz ) / ( (double) taps_used ); */
568 for (i=0; i<filt_sz; i++) fraction += ABS(filt[i]);
569 fraction = ( fraction / taps_used );
570 for (i=0; i<filt_sz; i++) result[i] *= fraction;
571 }
572 return(0);
573 }
574
575
576/* --------------------------------------------------------------------
577Reflect, multiplying tap of filter which is at the edge of the image
578by root 2. This maintains orthogonality of odd-length linear-phase
579QMF filters, but it is not useful for most applications, since it
580alters the DC level. */
581
582int ereflect(filt,x_dim,y_dim,x_pos,y_pos,result,r_or_e)
583 register double *filt, *result;
584 register int x_dim;
585 int y_dim, x_pos, y_pos, r_or_e;
586 {
587 register int y_filt,x_filt, y_res,x_res;
588 register int x_base = (x_pos>0)?(x_dim-1):0;
589 register int y_base = x_dim * ((y_pos>0)?(y_dim-1):0);
590 int filt_sz = x_dim*y_dim;
591 int x_overhang = (x_pos>1)?(x_pos-x_dim):((x_pos<-1)?(x_pos+1):0);
592 int y_overhang = x_dim * ( (y_pos>1)?(y_pos-y_dim):((y_pos<-1)?(y_pos+1):0) );
593 int i;
594 double norm,onorm;
595
596 for (i=0; i<filt_sz; i++) result[i] = 0.0;
597
598 /* reflect at boundary */
599 for (y_filt=0, y_res=y_overhang;
600 y_filt<filt_sz;
601 y_filt+=x_dim, y_res+=x_dim)
602 for (x_filt=y_filt, x_res=x_overhang;
603 x_filt<y_filt+x_dim;
604 x_filt++, x_res++)
605 result[ABS(y_base-ABS(y_res)) + ABS(x_base-ABS(x_res))]
606 += filt[x_filt];
607
608 /* now multiply edge by root 2 */
609 if (x_pos ISNT 0)
610 for (y_filt=x_base; y_filt<filt_sz; y_filt+=x_dim)
611 result[y_filt] *= ROOT2;
612 if (y_pos ISNT 0)
613 for (x_filt=y_base; x_filt<y_base+x_dim; x_filt++)
614 result[x_filt] *= ROOT2;
615
616 /* now normalize to norm of original filter */
617 for (norm=0.0,i=0; i<filt_sz; i++)
618 norm += (result[i]*result[i]);
619 norm=sqrt(norm);
620
621 for (onorm=0.0,i=0; i<filt_sz; i++)
622 onorm += (filt[i]*filt[i]);
623 onorm = sqrt(onorm);
624
625 norm = norm/onorm;
626 for (i=0; i<filt_sz; i++)
627 result[i] /= norm;
628 return(0);
629 }
630
631
632/* ------- printout stuff for testing ------------------------------
633 printf("Xpos: %d, Ypos: %d", x_pos, y_pos);
634 for (y_filt=0; y_filt<y_dim; y_filt++)
635 {
636 printf("\n");
637 for (x_filt=0; x_filt<x_dim; x_filt++)
638 printf("%6.1f", result[y_filt*x_dim+x_filt]);
639 }
640 printf("\n");
641*/
642
643
644
645/* Local Variables: */
646/* buffer-read-only: t */
647/* End: */
generated by cgit v1.2.2 (git 2.25.0) at 2025-04-30 18:26:18 -0400