GCC Code Coverage Report

Directory: ./
File: libfprint/nbis/mindtct/maps.c
Date: 2024-09-16 14:36:32
Exec Total Coverage
Lines: 581 640 90.8%
Functions: 20 20 100.0%
Branches: 273 320 85.3%
Line Branch Exec Source
1 /*******************************************************************************
2
3 License:
4 This software and/or related materials was developed at the National Institute
5 of Standards and Technology (NIST) by employees of the Federal Government
6 in the course of their official duties. Pursuant to title 17 Section 105
7 of the United States Code, this software is not subject to copyright
8 protection and is in the public domain.
9
10 This software and/or related materials have been determined to be not subject
11 to the EAR (see Part 734.3 of the EAR for exact details) because it is
12 a publicly available technology and software, and is freely distributed
13 to any interested party with no licensing requirements. Therefore, it is
14 permissible to distribute this software as a free download from the internet.
15
16 Disclaimer:
17 This software and/or related materials was developed to promote biometric
18 standards and biometric technology testing for the Federal Government
19 in accordance with the USA PATRIOT Act and the Enhanced Border Security
20 and Visa Entry Reform Act. Specific hardware and software products identified
21 in this software were used in order to perform the software development.
22 In no case does such identification imply recommendation or endorsement
23 by the National Institute of Standards and Technology, nor does it imply that
24 the products and equipment identified are necessarily the best available
25 for the purpose.
26
27 This software and/or related materials are provided "AS-IS" without warranty
28 of any kind including NO WARRANTY OF PERFORMANCE, MERCHANTABILITY,
29 NO WARRANTY OF NON-INFRINGEMENT OF ANY 3RD PARTY INTELLECTUAL PROPERTY
30 or FITNESS FOR A PARTICULAR PURPOSE or for any purpose whatsoever, for the
31 licensed product, however used. In no event shall NIST be liable for any
32 damages and/or costs, including but not limited to incidental or consequential
33 damages of any kind, including economic damage or injury to property and lost
34 profits, regardless of whether NIST shall be advised, have reason to know,
35 or in fact shall know of the possibility.
36
37 By using this software, you agree to bear all risk relating to quality,
38 use and performance of the software and/or related materials. You agree
39 to hold the Government harmless from any claim arising from your use
40 of the software.
41
42 *******************************************************************************/
43
44
45 /***********************************************************************
46 LIBRARY: LFS - NIST Latent Fingerprint System
47
48 FILE: MAPS.C
49 AUTHOR: Michael D. Garris
50 DATE: 03/16/1999
51 UPDATED: 10/04/1999 Version 2 by MDG
52 UPDATED: 10/26/1999 by MDG
53 To permit margin blocks to be flagged in
54 low contrast and low flow maps.
55 UPDATED: 03/16/2005 by MDG
56
57 Contains routines responsible for computing various block-based
58 maps (including directional ridge flow maps) from an arbitrarily-
59 sized image as part of the NIST Latent Fingerprint System (LFS).
60
61 ***********************************************************************
62 ROUTINES:
63 gen_image_maps()
64 gen_initial_maps()
65 interpolate_direction_map()
66 morph_TF_map()
67 pixelize_map()
68 smooth_direction_map()
69 gen_high_curve_map()
70 gen_imap()
71 gen_initial_imap()
72 primary_dir_test()
73 secondary_fork_test()
74 remove_incon_dirs()
75 test_top_edge()
76 test_right_edge()
77 test_bottom_edge()
78 test_left_edge()
79 remove_dir()
80 average_8nbr_dir()
81 num_valid_8nbrs()
82 smooth_imap()
83 gen_nmap()
84 vorticity()
85 accum_nbr_vorticity()
86 curvature()
87
88 ***********************************************************************/
89
90 #include <stdio.h>
91 #include <lfs.h>
92 #include <morph.h>
93 #include <log.h>
94
95 /*************************************************************************
96 **************************************************************************
97 #cat: gen_image_maps - Computes a set of image maps based on Version 2
98 #cat: of the NIST LFS System. The first map is a Direction Map
99 #cat: which is a 2D vector of integer directions, where each
100 #cat: direction represents the dominant ridge flow in a block of
101 #cat: the input grayscale image. The Low Contrast Map flags
102 #cat: blocks with insufficient contrast. The Low Flow Map flags
103 #cat: blocks with insufficient ridge flow. The High Curve Map
104 #cat: flags blocks containing high curvature. This routine will
105 #cat: generate maps for an arbitrarily sized, non-square, image.
106
107 Input:
108 pdata - padded input image data (8 bits [0..256) grayscale)
109 pw - padded width (in pixels) of the input image
110 ph - padded height (in pixels) of the input image
111 dir2rad - lookup table for converting integer directions
112 dftwaves - structure containing the DFT wave forms
113 dftgrids - structure containing the rotated pixel grid offsets
114 lfsparms - parameters and thresholds for controlling LFS
115 Output:
116 odmap - points to the created Direction Map
117 olcmap - points to the created Low Contrast Map
118 olfmap - points to the Low Ridge Flow Map
119 ohcmap - points to the High Curvature Map
120 omw - width (in blocks) of the maps
121 omh - height (in blocks) of the maps
122 Return Code:
123 Zero - successful completion
124 Negative - system error
125 **************************************************************************/
126 48 int gen_image_maps(int **odmap, int **olcmap, int **olfmap, int **ohcmap,
127 int *omw, int *omh,
128 unsigned char *pdata, const int pw, const int ph,
129 const DIR2RAD *dir2rad, const DFTWAVES *dftwaves,
130 const ROTGRIDS *dftgrids, const LFSPARMS *lfsparms)
131 {
132 48 int *direction_map, *low_contrast_map, *low_flow_map, *high_curve_map;
133 48 int mw, mh, iw, ih;
134 48 int *blkoffs;
135 48 int ret; /* return code */
136
137 /* 1. Compute block offsets for the entire image, accounting for pad */
138 /* Block_offsets() assumes square block (grid), so ERROR otherwise. */
139
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
 48 if(dftgrids->grid_w != dftgrids->grid_h){
140 ✗ fprintf(stderr,
141 "ERROR : gen_image_maps : DFT grids must be square\n");
142 ✗ return(-540);
143 }
144 /* Compute unpadded image dimensions. */
145 48 iw = pw - (dftgrids->pad<<1);
146 48 ih = ph - (dftgrids->pad<<1);
147
1/2
✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
 48 if((ret = block_offsets(&blkoffs, &mw, &mh, iw, ih,
148 48 dftgrids->pad, lfsparms->blocksize))){
149 return(ret);
150 }
151
152 /* 2. Generate initial Direction Map and Low Contrast Map*/
153
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
 48 if((ret = gen_initial_maps(&direction_map, &low_contrast_map,
154 &low_flow_map, blkoffs, mw, mh,
155 pdata, pw, ph, dftwaves, dftgrids, lfsparms))){
156 /* Free memory allocated to this point. */
157 ✗ g_free(blkoffs);
158 ✗ return(ret);
159 }
160
161
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
 48 if((ret = morph_TF_map(low_flow_map, mw, mh, lfsparms))){
162 ✗ g_free(direction_map);
163 ✗ g_free(low_contrast_map);
164 ✗ g_free(low_flow_map);
165 ✗ return(ret);
166 }
167
168 /* 3. Remove directions that are inconsistent with neighbors */
169 48 remove_incon_dirs(direction_map, mw, mh, dir2rad, lfsparms);
170
171
172 /* 4. Smooth Direction Map values with their neighbors */
173 48 smooth_direction_map(direction_map, low_contrast_map, mw, mh,
174 dir2rad, lfsparms);
175
176 /* 5. Interpolate INVALID direction blocks with their valid neighbors. */
177
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
 48 if((ret = interpolate_direction_map(direction_map, low_contrast_map,
178 mw, mh, lfsparms))){
179 ✗ g_free(direction_map);
180 ✗ g_free(low_contrast_map);
181 ✗ g_free(low_flow_map);
182 ✗ return(ret);
183 }
184
185 /* May be able to skip steps 6 and/or 7 if computation time */
186 /* is a critical factor. */
187
188 /* 6. Remove directions that are inconsistent with neighbors */
189 48 remove_incon_dirs(direction_map, mw, mh, dir2rad, lfsparms);
190
191 /* 7. Smooth Direction Map values with their neighbors. */
192 48 smooth_direction_map(direction_map, low_contrast_map, mw, mh,
193 dir2rad, lfsparms);
194
195 /* 8. Set the Direction Map values in the image margin to INVALID. */
196 48 set_margin_blocks(direction_map, mw, mh, INVALID_DIR);
197
198 /* 9. Generate High Curvature Map from interpolated Direction Map. */
199
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
 48 if((ret = gen_high_curve_map(&high_curve_map, direction_map, mw, mh,
200 lfsparms))){
201 ✗ g_free(direction_map);
202 ✗ g_free(low_contrast_map);
203 ✗ g_free(low_flow_map);
204 ✗ return(ret);
205 }
206
207 /* Deallocate working memory. */
208 48 g_free(blkoffs);
209
210 48 *odmap = direction_map;
211 48 *olcmap = low_contrast_map;
212 48 *olfmap = low_flow_map;
213 48 *ohcmap = high_curve_map;
214 48 *omw = mw;
215 48 *omh = mh;
216 48 return(0);
217 }
218
219 /*************************************************************************
220 **************************************************************************
221 #cat: gen_initial_maps - Creates an initial Direction Map from the given
222 #cat: input image. It very important that the image be properly
223 #cat: padded so that rotated grids along the boundary of the image
224 #cat: do not access unkown memory. The rotated grids are used by a
225 #cat: DFT-based analysis to determine the integer directions
226 #cat: in the map. Typically this initial vector of directions will
227 #cat: subsequently have weak or inconsistent directions removed
228 #cat: followed by a smoothing process. The resulting Direction
229 #cat: Map contains valid directions >= 0 and INVALID values = -1.
230 #cat: This routine also computes and returns 2 other image maps.
231 #cat: The Low Contrast Map flags blocks in the image with
232 #cat: insufficient contrast. Blocks with low contrast have a
233 #cat: corresponding direction of INVALID in the Direction Map.
234 #cat: The Low Flow Map flags blocks in which the DFT analyses
235 #cat: could not determine a significant ridge flow. Blocks with
236 #cat: low ridge flow also have a corresponding direction of
237 #cat: INVALID in the Direction Map.
238
239 Input:
240 blkoffs - offsets to the pixel origin of each block in the padded image
241 mw - number of blocks horizontally in the padded input image
242 mh - number of blocks vertically in the padded input image
243 pdata - padded input image data (8 bits [0..256) grayscale)
244 pw - width (in pixels) of the padded input image
245 ph - height (in pixels) of the padded input image
246 dftwaves - structure containing the DFT wave forms
247 dftgrids - structure containing the rotated pixel grid offsets
248 lfsparms - parameters and thresholds for controlling LFS
249 Output:
250 odmap - points to the newly created Direction Map
251 olcmap - points to the newly created Low Contrast Map
252 Return Code:
253 Zero - successful completion
254 Negative - system error
255 **************************************************************************/
256 48 int gen_initial_maps(int **odmap, int **olcmap, int **olfmap,
257 int *blkoffs, const int mw, const int mh,
258 unsigned char *pdata, const int pw, const int ph,
259 const DFTWAVES *dftwaves, const ROTGRIDS *dftgrids,
260 const LFSPARMS *lfsparms)
261 {
262 48 int *direction_map, *low_contrast_map, *low_flow_map;
263 48 int bi, bsize, blkdir;
264 48 int *wis, *powmax_dirs;
265 48 double **powers, *powmaxs, *pownorms;
266 48 int nstats;
267 48 int ret; /* return code */
268 48 int dft_offset;
269 48 int xminlimit, xmaxlimit, yminlimit, ymaxlimit;
270 48 int win_x, win_y, low_contrast_offset;
271
272 48 print2log("INITIAL MAP\n");
273
274 /* Compute total number of blocks in map */
275
2/4
✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 48 times.
 48 ASSERT_INT_MUL(mw, mh);
276 48 bsize = mw * mh;
277
278 /* Allocate Direction Map memory */
279 48 direction_map = (int *)g_malloc(bsize * sizeof(int));
280 /* Initialize the Direction Map to INVALID (-1). */
281 48 memset(direction_map, INVALID_DIR, bsize * sizeof(int));
282
283 /* Allocate Low Contrast Map memory */
284 48 low_contrast_map = (int *)g_malloc(bsize * sizeof(int));
285 /* Initialize the Low Contrast Map to FALSE (0). */
286 48 memset(low_contrast_map, 0, bsize * sizeof(int));
287
288 /* Allocate Low Ridge Flow Map memory */
289 48 low_flow_map = (int *)g_malloc(bsize * sizeof(int));
290 /* Initialize the Low Flow Map to FALSE (0). */
291 48 memset(low_flow_map, 0, bsize * sizeof(int));
292
293 /* Allocate DFT directional power vectors */
294
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
 48 if((ret = alloc_dir_powers(&powers, dftwaves->nwaves, dftgrids->ngrids))){
295 /* Free memory allocated to this point. */
296 ✗ g_free(direction_map);
297 ✗ g_free(low_contrast_map);
298 ✗ g_free(low_flow_map);
299 ✗ return(ret);
300 }
301
302 /* Allocate DFT power statistic arrays */
303 /* Compute length of statistics arrays. Statistics not needed */
304 /* for the first DFT wave, so the length is number of waves - 1. */
305 48 nstats = dftwaves->nwaves - 1;
306
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 48 times.
 48 if((ret = alloc_power_stats(&wis, &powmaxs, &powmax_dirs,
307 &pownorms, nstats))){
308 /* Free memory allocated to this point. */
309 ✗ g_free(direction_map);
310 ✗ g_free(low_contrast_map);
311 ✗ g_free(low_flow_map);
312 ✗ free_dir_powers(powers, dftwaves->nwaves);
313 ✗ return(ret);
314 }
315
316 /* Compute special window origin limits for determining low contrast. */
317 /* These pixel limits avoid analyzing the padded borders of the image. */
318 48 xminlimit = dftgrids->pad;
319 48 yminlimit = dftgrids->pad;
320 48 xmaxlimit = pw - dftgrids->pad - lfsparms->windowsize - 1;
321 48 ymaxlimit = ph - dftgrids->pad - lfsparms->windowsize - 1;
322
323 /* Foreach block in image ... */
324
2/2
✓ Branch 0 taken 50730 times.
✓ Branch 1 taken 48 times.
 50778 for(bi = 0; bi < bsize; bi++){
325 /* Adjust block offset from pointing to block origin to pointing */
326 /* to surrounding window origin. */
327 50730 dft_offset = blkoffs[bi] - (lfsparms->windowoffset * pw) -
328 lfsparms->windowoffset;
329
330 /* Compute pixel coords of window origin. */
331 50730 win_x = dft_offset % pw;
332 50730 win_y = (int)(dft_offset / pw);
333
334 /* Make sure the current window does not access padded image pixels */
335 /* for analyzing low contrast. */
336 50730 win_x = max(xminlimit, win_x);
337 50730 win_x = min(xmaxlimit, win_x);
338 50730 win_y = max(yminlimit, win_y);
339 50730 win_y = min(ymaxlimit, win_y);
340 50730 low_contrast_offset = (win_y * pw) + win_x;
341
342 50730 print2log(" BLOCK %2d (%2d, %2d) ", bi, bi%mw, bi/mw);
343
344 /* If block is low contrast ... */
345
2/2
✓ Branch 1 taken 4591 times.
✓ Branch 2 taken 46139 times.
 50730 if((ret = low_contrast_block(low_contrast_offset, lfsparms->windowsize,
346 pdata, pw, ph, lfsparms))){
347 /* If system error ... */
348
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4591 times.
 4591 if(ret < 0){
349 ✗ g_free(direction_map);
350 ✗ g_free(low_contrast_map);
351 ✗ g_free(low_flow_map);
352 ✗ free_dir_powers(powers, dftwaves->nwaves);
353 ✗ g_free(wis);
354 ✗ g_free(powmaxs);
355 ✗ g_free(powmax_dirs);
356 ✗ g_free(pownorms);
357 ✗ return(ret);
358 }
359
360 /* Otherwise, block is low contrast ... */
361 4591 print2log("LOW CONTRAST\n");
362 4591 low_contrast_map[bi] = TRUE;
363 /* Direction Map's block is already set to INVALID. */
364 }
365 /* Otherwise, sufficient contrast for DFT processing ... */
366 else {
367 46139 print2log("\n");
368
369 /* Compute DFT powers */
370
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 46139 times.
 46139 if((ret = dft_dir_powers(powers, pdata, low_contrast_offset, pw, ph,
371 dftwaves, dftgrids))){
372 /* Free memory allocated to this point. */
373 ✗ g_free(direction_map);
374 ✗ g_free(low_contrast_map);
375 ✗ g_free(low_flow_map);
376 ✗ free_dir_powers(powers, dftwaves->nwaves);
377 ✗ g_free(wis);
378 ✗ g_free(powmaxs);
379 ✗ g_free(powmax_dirs);
380 ✗ g_free(pownorms);
381 ✗ return(ret);
382 }
383
384 /* Compute DFT power statistics, skipping first applied DFT */
385 /* wave. This is dependent on how the primary and secondary */
386 /* direction tests work below. */
387
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 46139 times.
 46139 if((ret = dft_power_stats(wis, powmaxs, powmax_dirs, pownorms, powers,
388 46139 1, dftwaves->nwaves, dftgrids->ngrids))){
389 /* Free memory allocated to this point. */
390 ✗ g_free(direction_map);
391 ✗ g_free(low_contrast_map);
392 ✗ g_free(low_flow_map);
393 ✗ free_dir_powers(powers, dftwaves->nwaves);
394 ✗ g_free(wis);
395 ✗ g_free(powmaxs);
396 ✗ g_free(powmax_dirs);
397 ✗ g_free(pownorms);
398 ✗ return(ret);
399 }
400
401 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
402 { int _w;
403 fprintf(logfp, " Power\n");
404 for(_w = 0; _w < nstats; _w++){
405 /* Add 1 to wis[w] to create index to original g_dft_coefs[] */
406 fprintf(logfp, " wis[%d] %d %12.3f %2d %9.3f %12.3f\n",
407 _w, wis[_w]+1,
408 powmaxs[wis[_w]], powmax_dirs[wis[_w]], pownorms[wis[_w]],
409 powers[0][powmax_dirs[wis[_w]]]);
410 }
411 }
412 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
413
414 /* Conduct primary direction test */
415 46139 blkdir = primary_dir_test(powers, wis, powmaxs, powmax_dirs,
416 pownorms, nstats, lfsparms);
417
418
2/2
✓ Branch 0 taken 35558 times.
✓ Branch 1 taken 10581 times.
 46139 if(blkdir != INVALID_DIR)
419 35558 direction_map[bi] = blkdir;
420 else{
421 /* Conduct secondary (fork) direction test */
422 10581 blkdir = secondary_fork_test(powers, wis, powmaxs, powmax_dirs,
423 pownorms, nstats, lfsparms);
424
2/2
✓ Branch 0 taken 473 times.
✓ Branch 1 taken 10108 times.
 10581 if(blkdir != INVALID_DIR)
425 473 direction_map[bi] = blkdir;
426 /* Otherwise current direction in Direction Map remains INVALID */
427 else
428 /* Flag the block as having LOW RIDGE FLOW. */
429 10108 low_flow_map[bi] = TRUE;
430 }
431
432 } /* End DFT */
433 } /* bi */
434
435 /* Deallocate working memory */
436 48 free_dir_powers(powers, dftwaves->nwaves);
437 48 g_free(wis);
438 48 g_free(powmaxs);
439 48 g_free(powmax_dirs);
440 48 g_free(pownorms);
441
442 48 *odmap = direction_map;
443 48 *olcmap = low_contrast_map;
444 48 *olfmap = low_flow_map;
445 48 return(0);
446 }
447
448 /*************************************************************************
449 **************************************************************************
450 #cat: interpolate_direction_map - Take a Direction Map and Low Contrast
451 #cat: Map and attempts to fill in INVALID directions in the
452 #cat: Direction Map based on a blocks valid neighbors. The
453 #cat: valid neighboring directions are combined in a weighted
454 #cat: average inversely proportional to their distance from
455 #cat: the block being interpolated. Low Contrast blocks are
456 #cat: used to prempt the search for a valid neighbor in a
457 #cat: specific direction, which keeps the process from
458 #cat: interpolating directions for blocks in the background and
459 #cat: and perimeter of the fingerprint in the image.
460
461 Input:
462 direction_map - map of blocks containing directional ridge flow
463 low_contrast_map - map of blocks flagged as LOW CONTRAST
464 mw - number of blocks horizontally in the maps
465 mh - number of blocks vertically in the maps
466 lfsparms - parameters and thresholds for controlling LFS
467 Output:
468 direction_map - contains the newly interpolated results
469 Return Code:
470 Zero - successful completion
471 Negative - system error
472 **************************************************************************/
473 48 int interpolate_direction_map(int *direction_map, int *low_contrast_map,
474 const int mw, const int mh, const LFSPARMS *lfsparms)
475 {
476 48 int x, y, new_dir;
477 48 int n_dir, e_dir, s_dir, w_dir;
478 48 int n_dist = 0, e_dist = 0, s_dist = 0, w_dist = 0, total_dist;
479 48 int n_found, e_found, s_found, w_found, total_found;
480 48 int n_delta = 0, e_delta = 0, s_delta = 0, w_delta = 0, total_delta;
481 48 int nbr_x, nbr_y;
482 48 int *omap, *dptr, *cptr, *optr;
483 48 double avr_dir;
484
485 48 print2log("INTERPOLATE DIRECTION MAP\n");
486
487 /* Allocate output (interpolated) Direction Map. */
488
1/2
✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
 48 ASSERT_SIZE_MUL(mw, mh);
489
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
 48 ASSERT_SIZE_MUL(mw * mh, sizeof(int));
490 48 omap = (int *)g_malloc(mw * mh * sizeof(int));
491
492 /* Set pointers to the first block in the maps. */
493 48 dptr = direction_map;
494 48 cptr = low_contrast_map;
495 48 optr = omap;
496
497 /* Foreach block in the maps ... */
498
2/2
✓ Branch 1 taken 1673 times.
✓ Branch 2 taken 48 times.
 1769 for(y = 0; y < mh; y++){
499
2/2
✓ Branch 0 taken 50730 times.
✓ Branch 1 taken 1673 times.
 52403 for(x = 0; x < mw; x++){
500
501 /* If image block is NOT LOW CONTRAST and has INVALID direction ... */
502
4/4
✓ Branch 0 taken 46139 times.
✓ Branch 1 taken 4591 times.
✓ Branch 2 taken 12738 times.
✓ Branch 3 taken 33401 times.
 50730 if((!*cptr) && (*dptr == INVALID_DIR)){
503
504 /* Set neighbor accumulators to 0. */
505 12738 total_found = 0;
506 12738 total_dist = 0;
507
508 /* Find north neighbor. */
509
2/2
✓ Branch 1 taken 7763 times.
✓ Branch 2 taken 4975 times.
 12738 if((n_found = find_valid_block(&n_dir, &nbr_x, &nbr_y,
510 direction_map, low_contrast_map,
511 x, y, mw, mh, 0, -1)) == FOUND){
512 /* Compute north distance. */
513 7763 n_dist = y - nbr_y;
514 /* Accumulate neighbor distance. */
515 7763 total_dist += n_dist;
516 /* Bump number of neighbors found. */
517 7763 total_found++;
518 }
519
520 /* Find east neighbor. */
521
2/2
✓ Branch 1 taken 8641 times.
✓ Branch 2 taken 4097 times.
 12738 if((e_found = find_valid_block(&e_dir, &nbr_x, &nbr_y,
522 direction_map, low_contrast_map,
523 x, y, mw, mh, 1, 0)) == FOUND){
524 /* Compute east distance. */
525 8641 e_dist = nbr_x - x;
526 /* Accumulate neighbor distance. */
527 8641 total_dist += e_dist;
528 /* Bump number of neighbors found. */
529 8641 total_found++;
530 }
531
532 /* Find south neighbor. */
533
2/2
✓ Branch 1 taken 9466 times.
✓ Branch 2 taken 3272 times.
 12738 if((s_found = find_valid_block(&s_dir, &nbr_x, &nbr_y,
534 direction_map, low_contrast_map,
535 x, y, mw, mh, 0, 1)) == FOUND){
536 /* Compute south distance. */
537 9466 s_dist = nbr_y - y;
538 /* Accumulate neighbor distance. */
539 9466 total_dist += s_dist;
540 /* Bump number of neighbors found. */
541 9466 total_found++;
542 }
543
544 /* Find west neighbor. */
545
2/2
✓ Branch 1 taken 9448 times.
✓ Branch 2 taken 3290 times.
 12738 if((w_found = find_valid_block(&w_dir, &nbr_x, &nbr_y,
546 direction_map, low_contrast_map,
547 x, y, mw, mh, -1, 0)) == FOUND){
548 /* Compute west distance. */
549 9448 w_dist = x - nbr_x;
550 /* Accumulate neighbor distance. */
551 9448 total_dist += w_dist;
552 /* Bump number of neighbors found. */
553 9448 total_found++;
554 }
555
556 /* If a sufficient number of neighbors found (Ex. 2) ... */
557
2/2
✓ Branch 0 taken 10891 times.
✓ Branch 1 taken 1847 times.
 12738 if(total_found >= lfsparms->min_interpolate_nbrs){
558
559 /* Accumulate weighted sum of neighboring directions */
560 /* inversely related to the distance from current block. */
561 10891 total_delta = 0.0;
562 /* If neighbor found to the north ... */
563
2/2
✓ Branch 0 taken 7373 times.
✓ Branch 1 taken 3518 times.
 10891 if(n_found){
564 7373 n_delta = total_dist - n_dist;
565 7373 total_delta += n_delta;
566 }
567 /* If neighbor found to the east ... */
568
2/2
✓ Branch 0 taken 8280 times.
✓ Branch 1 taken 2611 times.
 10891 if(e_found){
569 8280 e_delta = total_dist - e_dist;
570 8280 total_delta += e_delta;
571 }
572 /* If neighbor found to the south ... */
573
2/2
✓ Branch 0 taken 9188 times.
✓ Branch 1 taken 1703 times.
 10891 if(s_found){
574 9188 s_delta = total_dist - s_dist;
575 9188 total_delta += s_delta;
576 }
577 /* If neighbor found to the west ... */
578
2/2
✓ Branch 0 taken 8941 times.
✓ Branch 1 taken 1950 times.
 10891 if(w_found){
579 8941 w_delta = total_dist - w_dist;
580 8941 total_delta += w_delta;
581 }
582
583 10891 avr_dir = 0.0;
584
585
2/2
✓ Branch 0 taken 7373 times.
✓ Branch 1 taken 3518 times.
 10891 if(n_found){
586 7373 avr_dir += (n_dir*(n_delta/(double)total_delta));
587 }
588
2/2
✓ Branch 0 taken 8280 times.
✓ Branch 1 taken 2611 times.
 10891 if(e_found){
589 8280 avr_dir += (e_dir*(e_delta/(double)total_delta));
590 }
591
2/2
✓ Branch 0 taken 9188 times.
✓ Branch 1 taken 1703 times.
 10891 if(s_found){
592 9188 avr_dir += (s_dir*(s_delta/(double)total_delta));
593 }
594
2/2
✓ Branch 0 taken 8941 times.
✓ Branch 1 taken 1950 times.
 10891 if(w_found){
595 8941 avr_dir += (w_dir*(w_delta/(double)total_delta));
596 }
597
598 /* Need to truncate precision so that answers are consistent */
599 /* on different computer architectures when rounding doubles. */
600
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10891 times.
 10891 avr_dir = trunc_dbl_precision(avr_dir, TRUNC_SCALE);
601
602 /* Assign interpolated direction to output Direction Map. */
603
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10891 times.
 10891 new_dir = sround(avr_dir);
604
605 10891 print2log(" Block %d,%d INTERP numnbs=%d newdir=%d\n",
606 x, y, total_found, new_dir);
607
608 10891 *optr = new_dir;
609 }
610 else{
611 /* Otherwise, the direction remains INVALID. */
612 1847 *optr = *dptr;
613 }
614 }
615 else{
616 /* Otherwise, assign the current direction to the output block. */
617 37992 *optr = *dptr;
618 }
619
620 /* Bump to the next block in the maps ... */
621 50730 dptr++;
622 50730 cptr++;
623 50730 optr++;
624 }
625 }
626
627 /* Copy the interpolated directions into the input map. */
628 48 memcpy(direction_map, omap, mw*mh*sizeof(int));
629 /* Deallocate the working memory. */
630 48 g_free(omap);
631
632 /* Return normally. */
633 48 return(0);
634 }
635
636 /*************************************************************************
637 **************************************************************************
638 #cat: morph_tf_map - Takes a 2D vector of TRUE and FALSE values integers
639 #cat: and dialates and erodes the map in an attempt to fill
640 #cat: in voids in the map.
641
642 Input:
643 tfmap - vector of integer block values
644 mw - width (in blocks) of the map
645 mh - height (in blocks) of the map
646 lfsparms - parameters and thresholds for controlling LFS
647 Output:
648 tfmap - resulting morphed map
649 **************************************************************************/
650 48 int morph_TF_map(int *tfmap, const int mw, const int mh,
651 const LFSPARMS *lfsparms)
652 {
653 48 unsigned char *cimage, *mimage, *cptr;
654 48 int *mptr;
655 48 int i;
656
657
2/4
✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 48 times.
 48 ASSERT_INT_MUL(mw, mh);
658
659 /* Convert TRUE/FALSE map into a binary byte image. */
660 48 cimage = (unsigned char *)g_malloc(mw * mh);
661
662 48 mimage = (unsigned char *)g_malloc(mw * mh);
663
664 48 cptr = cimage;
665 48 mptr = tfmap;
666
2/2
✓ Branch 1 taken 50730 times.
✓ Branch 2 taken 48 times.
 50826 for(i = 0; i < mw*mh; i++){
667 50730 *cptr++ = *mptr++;
668 }
669
670 48 dilate_charimage_2(cimage, mimage, mw, mh);
671 48 dilate_charimage_2(mimage, cimage, mw, mh);
672 48 erode_charimage_2(cimage, mimage, mw, mh);
673 48 erode_charimage_2(mimage, cimage, mw, mh);
674
675 48 cptr = cimage;
676 48 mptr = tfmap;
677
2/2
✓ Branch 1 taken 50730 times.
✓ Branch 2 taken 48 times.
 50826 for(i = 0; i < mw*mh; i++){
678 50730 *mptr++ = *cptr++;
679 }
680
681 48 g_free(cimage);
682 48 g_free(mimage);
683
684 48 return(0);
685 }
686
687 /*************************************************************************
688 **************************************************************************
689 #cat: pixelize_map - Takes a block image map and assigns each pixel in the
690 #cat: image its corresponding block value. This allows block
691 #cat: values in maps to be directly accessed via pixel addresses.
692
693 Input:
694 iw - the width (in pixels) of the corresponding image
695 ih - the height (in pixels) of the corresponding image
696 imap - input block image map
697 mw - the width (in blocks) of the map
698 mh - the height (in blocks) of the map
699 blocksize - the dimension (in pixels) of each block
700 Output:
701 omap - points to the resulting pixelized map
702 Return Code:
703 Zero - successful completion
704 Negative - system error
705 **************************************************************************/
706 192 int pixelize_map(int **omap, const int iw, const int ih,
707 int *imap, const int mw, const int mh, const int blocksize)
708 {
709 192 int *pmap;
710 192 int ret, x, y;
711 192 int *blkoffs, bw, bh, bi;
712 192 int *spptr, *pptr;
713
714
1/2
✓ Branch 0 taken 192 times.
✗ Branch 1 not taken.
 192 ASSERT_SIZE_MUL(iw, ih);
715
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 192 times.
 192 ASSERT_SIZE_MUL(iw * ih, sizeof(int));
716
717 192 pmap = (int *)g_malloc(iw * ih * sizeof(int));
718
719
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 192 times.
 192 if((ret = block_offsets(&blkoffs, &bw, &bh, iw, ih, 0, blocksize))){
720 ✗ g_free(pmap);
721 ✗ return(ret);
722 }
723
724
2/4
✓ Branch 0 taken 192 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 192 times.
 192 if((bw != mw) || (bh != mh)){
725 ✗ g_free(blkoffs);
726 ✗ g_free(pmap);
727 ✗ fprintf(stderr,
728 "ERROR : pixelize_map : block dimensions do not match\n");
729 ✗ return(-591);
730 }
731
732
2/2
✓ Branch 0 taken 202920 times.
✓ Branch 1 taken 192 times.
 203112 for(bi = 0; bi < mw*mh; bi++){
733 202920 spptr = pmap + blkoffs[bi];
734
2/2
✓ Branch 0 taken 1623360 times.
✓ Branch 1 taken 202920 times.
 1826280 for(y = 0; y < blocksize; y++){
735 pptr = spptr;
736
2/2
✓ Branch 0 taken 12986880 times.
✓ Branch 1 taken 1623360 times.
 14610240 for(x = 0; x < blocksize; x++){
737 12986880 *pptr++ = imap[bi];
738 }
739 1623360 spptr += iw;
740 }
741 }
742
743 /* Deallocate working memory. */
744 192 g_free(blkoffs);
745 /* Assign pixelized map to output pointer. */
746 192 *omap = pmap;
747
748 /* Return normally. */
749 192 return(0);
750 }
751
752 /*************************************************************************
753 **************************************************************************
754 #cat: smooth_direction_map - Takes a vector of integer directions and smooths
755 #cat: them by analyzing the direction of adjacent neighbors.
756
757 Input:
758 direction_map - vector of integer block values
759 mw - width (in blocks) of the map
760 mh - height (in blocks) of the map
761 dir2rad - lookup table for converting integer directions
762 lfsparms - parameters and thresholds for controlling LFS
763 Output:
764 imap - vector of smoothed input values
765 **************************************************************************/
766 96 void smooth_direction_map(int *direction_map, int *low_contrast_map,
767 const int mw, const int mh,
768 const DIR2RAD *dir2rad, const LFSPARMS *lfsparms)
769 {
770 96 int mx, my;
771 96 int *dptr, *cptr;
772 96 int avrdir, nvalid;
773 96 double dir_strength;
774
775 96 print2log("SMOOTH DIRECTION MAP\n");
776
777 /* Assign pointers to beginning of both maps. */
778 96 dptr = direction_map;
779 96 cptr = low_contrast_map;
780
781 /* Foreach block in maps ... */
782
2/2
✓ Branch 1 taken 3346 times.
✓ Branch 2 taken 96 times.
 3538 for(my = 0; my < mh; my++){
783
2/2
✓ Branch 0 taken 101460 times.
✓ Branch 1 taken 3346 times.
 104806 for(mx = 0; mx < mw; mx++){
784 /* If the current block does NOT have LOW CONTRAST ... */
785
2/2
✓ Branch 0 taken 92278 times.
✓ Branch 1 taken 9182 times.
 101460 if(!*cptr){
786
787 /* Compute average direction from neighbors, returning the */
788 /* number of valid neighbors used in the computation, and */
789 /* the "strength" of the average direction. */
790 92278 average_8nbr_dir(&avrdir, &dir_strength, &nvalid,
791 direction_map, mx, my, mw, mh, dir2rad);
792
793 /* If average direction strength is strong enough */
794 /* (Ex. thresh==0.2)... */
795
2/2
✓ Branch 0 taken 83572 times.
✓ Branch 1 taken 8706 times.
 92278 if(dir_strength >= lfsparms->dir_strength_min){
796 /* If Direction Map direction is valid ... */
797
2/2
✓ Branch 0 taken 74131 times.
✓ Branch 1 taken 9441 times.
 83572 if(*dptr != INVALID_DIR){
798 /* Conduct valid neighbor test (Ex. thresh==3)... */
799
1/2
✓ Branch 0 taken 74131 times.
✗ Branch 1 not taken.
 74131 if(nvalid >= lfsparms->rmv_valid_nbr_min){
800
801 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
802 fprintf(logfp, " BLOCK %2d (%2d, %2d)\n",
803 mx+(my*mw), mx, my);
804 fprintf(logfp, " Average NBR : %2d %6.3f %d\n",
805 avrdir, dir_strength, nvalid);
806 fprintf(logfp, " 1. Valid NBR (%d >= %d)\n",
807 nvalid, lfsparms->rmv_valid_nbr_min);
808 fprintf(logfp, " Valid Direction = %d\n", *dptr);
809 fprintf(logfp, " Smoothed Direction = %d\n", avrdir);
810 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
811
812 /* Reassign valid direction with average direction. */
813 74131 *dptr = avrdir;
814 }
815 }
816 /* Otherwise direction is invalid ... */
817 else{
818 /* Even if DIRECTION_MAP value is invalid, if number of */
819 /* valid neighbors is big enough (Ex. thresh==7)... */
820
2/2
✓ Branch 0 taken 1309 times.
✓ Branch 1 taken 8132 times.
 9441 if(nvalid >= lfsparms->smth_valid_nbr_min){
821
822 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
823 fprintf(logfp, " BLOCK %2d (%2d, %2d)\n",
824 mx+(my*mw), mx, my);
825 fprintf(logfp, " Average NBR : %2d %6.3f %d\n",
826 avrdir, dir_strength, nvalid);
827 fprintf(logfp, " 2. Invalid NBR (%d >= %d)\n",
828 nvalid, lfsparms->smth_valid_nbr_min);
829 fprintf(logfp, " Invalid Direction = %d\n", *dptr);
830 fprintf(logfp, " Smoothed Direction = %d\n", avrdir);
831 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
832
833 /* Assign invalid direction with average direction. */
834 1309 *dptr = avrdir;
835 }
836 }
837 }
838 }
839 /* Otherwise, block has LOW CONTRAST, so keep INVALID direction. */
840
841 /* Bump to next block in maps. */
842 101460 dptr++;
843 101460 cptr++;
844 }
845 }
846 96 }
847
848 /*************************************************************************
849 **************************************************************************
850 #cat: gen_high_curve_map - Takes a Direction Map and generates a new map
851 #cat: that flags blocks with HIGH CURVATURE.
852
853 Input:
854 direction_map - map of blocks containing directional ridge flow
855 mw - the width (in blocks) of the map
856 mh - the height (in blocks) of the map
857 lfsparms - parameters and thresholds for controlling LFS
858 Output:
859 ohcmap - points to the created High Curvature Map
860 Return Code:
861 Zero - successful completion
862 Negative - system error
863 **************************************************************************/
864 48 int gen_high_curve_map(int **ohcmap, int *direction_map,
865 const int mw, const int mh, const LFSPARMS *lfsparms)
866 {
867 48 int *high_curve_map, mapsize;
868 48 int *hptr, *dptr;
869 48 int bx, by;
870 48 int nvalid, cmeasure, vmeasure;
871
872
2/4
✓ Branch 0 taken 48 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 48 times.
 48 ASSERT_INT_MUL(mw, mh);
873 48 mapsize = mw*mh;
874
875 /* Allocate High Curvature Map. */
876
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
 48 ASSERT_SIZE_MUL(mapsize, sizeof(int));
877 48 high_curve_map = (int *)g_malloc(mapsize * sizeof(int));
878 /* Initialize High Curvature Map to FALSE (0). */
879 48 memset(high_curve_map, 0, mapsize*sizeof(int));
880
881 48 hptr = high_curve_map;
882 48 dptr = direction_map;
883
884 /* Foreach row in maps ... */
885
2/2
✓ Branch 0 taken 1673 times.
✓ Branch 1 taken 48 times.
 1721 for(by = 0; by < mh; by++){
886 /* Foreach column in maps ... */
887
2/2
✓ Branch 0 taken 50730 times.
✓ Branch 1 taken 1673 times.
 52403 for(bx = 0; bx < mw; bx++){
888
889 /* Count number of valid neighbors around current block ... */
890 50730 nvalid = num_valid_8nbrs(direction_map, bx, by, mw, mh);
891
892 /* If valid neighbors exist ... */
893
2/2
✓ Branch 0 taken 45618 times.
✓ Branch 1 taken 5112 times.
 50730 if(nvalid > 0){
894 /* If current block's direction is INVALID ... */
895
2/2
✓ Branch 0 taken 6667 times.
✓ Branch 1 taken 38951 times.
 45618 if(*dptr == INVALID_DIR){
896 /* If a sufficient number of VALID neighbors exists ... */
897
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 6651 times.
 6667 if(nvalid >= lfsparms->vort_valid_nbr_min){
898 /* Measure vorticity of neighbors. */
899 32 vmeasure = vorticity(direction_map, bx, by, mw, mh,
900 16 lfsparms->num_directions);
901 /* If vorticity is sufficiently high ... */
902
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 16 times.
 16 if(vmeasure >= lfsparms->highcurv_vorticity_min)
903 /* Flag block as HIGH CURVATURE. */
904 ✗ *hptr = TRUE;
905 }
906 }
907 /* Otherwise block has valid direction ... */
908 else{
909 /* Measure curvature around the valid block. */
910 77902 cmeasure = curvature(direction_map, bx, by, mw, mh,
911 38951 lfsparms->num_directions);
912 /* If curvature is sufficiently high ... */
913
2/2
✓ Branch 0 taken 1866 times.
✓ Branch 1 taken 37085 times.
 38951 if(cmeasure >= lfsparms->highcurv_curvature_min)
914 1866 *hptr = TRUE;
915 }
916 } /* Else (nvalid <= 0) */
917
918 /* Bump pointers to next block in maps. */
919 50730 dptr++;
920 50730 hptr++;
921
922 } /* bx */
923 } /* by */
924
925 /* Assign High Curvature Map to output pointer. */
926 48 *ohcmap = high_curve_map;
927
928 /* Return normally. */
929 48 return(0);
930 }
931
932 /*************************************************************************
933 **************************************************************************
934 #cat: gen_imap - Computes an IMAP, which is a 2D vector of integer directions,
935 #cat: where each direction represents the dominant ridge flow in
936 #cat: a block of the input grayscale image. This routine will
937 #cat: generate an IMAP for arbitrarily sized, non-square, images.
938
939 Input:
940 pdata - padded input image data (8 bits [0..256) grayscale)
941 pw - padded width (in pixels) of the input image
942 ph - padded height (in pixels) of the input image
943 dir2rad - lookup table for converting integer directions
944 dftwaves - structure containing the DFT wave forms
945 dftgrids - structure containing the rotated pixel grid offsets
946 lfsparms - parameters and thresholds for controlling LFS
947 Output:
948 optr - points to the created IMAP
949 ow - width (in blocks) of the IMAP
950 oh - height (in blocks) of the IMAP
951 Return Code:
952 Zero - successful completion
953 Negative - system error
954 **************************************************************************/
955
956 /*************************************************************************
957 **************************************************************************
958 #cat: gen_initial_imap - Creates an initial IMAP from the given input image.
959 #cat: It very important that the image be properly padded so
960 #cat: that rotated grids along the boudary of the image do not
961 #cat: access unkown memory. The rotated grids are used by a
962 #cat: DFT-based analysis to determine the integer directions
963 #cat: in the IMAP. Typically this initial vector of directions will
964 #cat: subsequently have weak or inconsistent directions removed
965 #cat: followed by a smoothing process.
966
967 Input:
968 blkoffs - offsets to the pixel origin of each block in the padded image
969 mw - number of blocks horizontally in the padded input image
970 mh - number of blocks vertically in the padded input image
971 pdata - padded input image data (8 bits [0..256) grayscale)
972 pw - width (in pixels) of the padded input image
973 ph - height (in pixels) of the padded input image
974 dftwaves - structure containing the DFT wave forms
975 dftgrids - structure containing the rotated pixel grid offsets
976 lfsparms - parameters and thresholds for controlling LFS
977 Output:
978 optr - points to the newly created IMAP
979 Return Code:
980 Zero - successful completion
981 Negative - system error
982 **************************************************************************/
983
984 /*************************************************************************
985 **************************************************************************
986 #cat: primary_dir_test - Applies the primary set of criteria for selecting
987 #cat: an IMAP integer direction from a set of DFT results
988 #cat: computed from a block of image data
989
990 Input:
991 powers - DFT power computed from each (N) wave frequencies at each
992 rotation direction in the current image block
993 wis - sorted order of the highest N-1 frequency power statistics
994 powmaxs - maximum power for each of the highest N-1 frequencies
995 powmax_dirs - directions associated with each of the N-1 maximum powers
996 pownorms - normalized power for each of the highest N-1 frequencies
997 nstats - N-1 wave frequencies (where N is the length of g_dft_coefs)
998 lfsparms - parameters and thresholds for controlling LFS
999 Return Code:
1000 Zero or Positive - The selected IMAP integer direction
1001 INVALID_DIR - IMAP Integer direction could not be determined
1002 **************************************************************************/
1003 46139 int primary_dir_test(double **powers, const int *wis,
1004 const double *powmaxs, const int *powmax_dirs,
1005 const double *pownorms, const int nstats,
1006 const LFSPARMS *lfsparms)
1007 {
1008 46139 int w;
1009
1010 46139 print2log(" Primary\n");
1011
1012 /* Look at max power statistics in decreasing order ... */
1013
2/2
✓ Branch 1 taken 70845 times.
✓ Branch 2 taken 10581 times.
 127565 for(w = 0; w < nstats; w++){
1014 /* 1. Test magnitude of current max power (Ex. Thresh==100000) */
1015
2/2
✓ Branch 0 taken 64207 times.
✓ Branch 1 taken 6638 times.
 70845 if((powmaxs[wis[w]] > lfsparms->powmax_min) &&
1016 /* 2. Test magnitude of normalized max power (Ex. Thresh==3.8) */
1017
2/2
✓ Branch 0 taken 36391 times.
✓ Branch 1 taken 27816 times.
 64207 (pownorms[wis[w]] > lfsparms->pownorm_min) &&
1018 /* 3. Test magnitude of power of lowest DFT frequency at current */
1019 /* max power direction and make sure it is not too big. */
1020 /* (Ex. Thresh==50000000) */
1021
2/2
✓ Branch 0 taken 35558 times.
✓ Branch 1 taken 833 times.
 36391 (powers[0][powmax_dirs[wis[w]]] <= lfsparms->powmax_max)){
1022
1023 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
1024 /* Add 1 to wis[w] to create index to original g_dft_coefs[] */
1025 fprintf(logfp,
1026 " Selected Wave = %d\n", wis[w]+1);
1027 fprintf(logfp,
1028 " 1. Power Magnitude (%12.3f > %12.3f)\n",
1029 powmaxs[wis[w]], lfsparms->powmax_min);
1030 fprintf(logfp,
1031 " 2. Norm Power Magnitude (%9.3f > %9.3f)\n",
1032 pownorms[wis[w]], lfsparms->pownorm_min);
1033 fprintf(logfp,
1034 " 3. Low Freq Wave Magnitude (%12.3f <= %12.3f)\n",
1035 powers[0][powmax_dirs[wis[w]]], lfsparms->powmax_max);
1036 fprintf(logfp,
1037 " PASSED\n");
1038 fprintf(logfp,
1039 " Selected Direction = %d\n",
1040 powmax_dirs[wis[w]]);
1041 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
1042
1043 /* If ALL 3 criteria met, return current max power direction. */
1044 35558 return(powmax_dirs[wis[w]]);
1045
1046 }
1047 }
1048
1049 /* Otherwise test failed. */
1050
1051 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
1052 fprintf(logfp, " 1. Power Magnitude ( > %12.3f)\n",
1053 lfsparms->powmax_min);
1054 fprintf(logfp, " 2. Norm Power Magnitude ( > %9.3f)\n",
1055 lfsparms->pownorm_min);
1056 fprintf(logfp, " 3. Low Freq Wave Magnitude ( <= %12.3f)\n",
1057 lfsparms->powmax_max);
1058 fprintf(logfp, " FAILED\n");
1059 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
1060
1061 return(INVALID_DIR);
1062 }
1063
1064 /*************************************************************************
1065 **************************************************************************
1066 #cat: secondary_fork_test - Applies a secondary set of criteria for selecting
1067 #cat: an IMAP integer direction from a set of DFT results
1068 #cat: computed from a block of image data. This test
1069 #cat: analyzes the strongest power statistics associated
1070 #cat: with a given frequency and direction and analyses
1071 #cat: small changes in direction to the left and right to
1072 #cat: determine if the block contains a "fork".
1073
1074 Input:
1075 powers - DFT power computed from each (N) wave frequencies at each
1076 rotation direction in the current image block
1077 wis - sorted order of the highest N-1 frequency power statistics
1078 powmaxs - maximum power for each of the highest N-1 frequencies
1079 powmax_dirs - directions associated with each of the N-1 maximum powers
1080 pownorms - normalized power for each of the highest N-1 frequencies
1081 nstats - N-1 wave frequencies (where N is the length of g_dft_coefs)
1082 lfsparms - parameters and thresholds for controlling LFS
1083 Return Code:
1084 Zero or Positive - The selected IMAP integer direction
1085 INVALID_DIR - IMAP Integer direction could not be determined
1086 **************************************************************************/
1087 10581 int secondary_fork_test(double **powers, const int *wis,
1088 const double *powmaxs, const int *powmax_dirs,
1089 const double *pownorms, const int nstats,
1090 const LFSPARMS *lfsparms)
1091 {
1092 10581 int ldir, rdir;
1093 10581 double fork_pownorm_min, fork_pow_thresh;
1094
1095 #ifdef LOG_REPORT
1096 { int firstpart = 0; /* Flag to determine if passed 1st part ... */
1097 fprintf(logfp, " Secondary\n");
1098 #endif
1099
1100 /* Relax the normalized power threshold under fork conditions. */
1101 10581 fork_pownorm_min = lfsparms->fork_pct_pownorm * lfsparms->pownorm_min;
1102
1103 /* 1. Test magnitude of largest max power (Ex. Thresh==100000) */
1104
2/2
✓ Branch 0 taken 9500 times.
✓ Branch 1 taken 1081 times.
 10581 if((powmaxs[wis[0]] > lfsparms->powmax_min) &&
1105 /* 2. Test magnitude of corresponding normalized power */
1106 /* (Ex. Thresh==2.85) */
1107
2/2
✓ Branch 0 taken 3919 times.
✓ Branch 1 taken 5581 times.
 9500 (pownorms[wis[0]] >= fork_pownorm_min) &&
1108 /* 3. Test magnitude of power of lowest DFT frequency at largest */
1109 /* max power direction and make sure it is not too big. */
1110 /* (Ex. Thresh==50000000) */
1111
2/2
✓ Branch 0 taken 376 times.
✓ Branch 1 taken 5205 times.
 5581 (powers[0][powmax_dirs[wis[0]]] <= lfsparms->powmax_max)){
1112
1113 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
1114 /* First part passed ... */
1115 firstpart = 1;
1116 fprintf(logfp,
1117 " Selected Wave = %d\n", wis[0]+1);
1118 fprintf(logfp,
1119 " 1. Power Magnitude (%12.3f > %12.3f)\n",
1120 powmaxs[wis[0]], lfsparms->powmax_min);
1121 fprintf(logfp,
1122 " 2. Norm Power Magnitude (%9.3f >= %9.3f)\n",
1123 pownorms[wis[0]], fork_pownorm_min);
1124 fprintf(logfp,
1125 " 3. Low Freq Wave Magnitude (%12.3f <= %12.3f)\n",
1126 powers[0][powmax_dirs[wis[0]]], lfsparms->powmax_max);
1127 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
1128
1129 /* Add FORK_INTERVALs to current direction modulo NDIRS */
1130 5205 rdir = (powmax_dirs[wis[0]] + lfsparms->fork_interval) %
1131 5205 lfsparms->num_directions;
1132
1133 /* Subtract FORK_INTERVALs from direction modulo NDIRS */
1134 /* For example, FORK_INTERVAL==2 & NDIRS==16, then */
1135 /* ldir = (dir - (16-2)) % 16 */
1136 /* which keeps result in proper modulo range. */
1137 5205 ldir = (powmax_dirs[wis[0]] + lfsparms->num_directions -
1138 lfsparms->fork_interval) % lfsparms->num_directions;
1139
1140 5205 print2log(" Left = %d, Current = %d, Right = %d\n",
1141 ldir, powmax_dirs[wis[0]], rdir);
1142
1143 /* Set forked angle threshold to be a % of the max directional */
1144 /* power. (Ex. thresh==0.7*powmax) */
1145 5205 fork_pow_thresh = powmaxs[wis[0]] * lfsparms->fork_pct_powmax;
1146
1147 /* Look up and test the computed power for the left and right */
1148 /* fork directions.s */
1149 /* The power stats (and thus wis) are on the range [0..nwaves-1) */
1150 /* as the statistics for the first DFT wave are not included. */
1151 /* The original power vectors exist for ALL DFT waves, therefore */
1152 /* wis indices must be added by 1 before addressing the original */
1153 /* powers vector. */
1154 /* LFS permits one and only one of the fork angles to exceed */
1155 /* the relative power threshold. */
1156
2/2
✓ Branch 0 taken 290 times.
✓ Branch 1 taken 4915 times.
 5205 if(((powers[wis[0]+1][ldir] <= fork_pow_thresh) ||
1157
4/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 289 times.
✓ Branch 2 taken 4915 times.
✓ Branch 3 taken 289 times.
 5205 (powers[wis[0]+1][rdir] <= fork_pow_thresh)) &&
1158 4915 ((powers[wis[0]+1][ldir] > fork_pow_thresh) ||
1159
2/2
✓ Branch 0 taken 4731 times.
✓ Branch 1 taken 184 times.
 4915 (powers[wis[0]+1][rdir] > fork_pow_thresh))){
1160
1161 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
1162 fprintf(logfp,
1163 " 4. Left Power Magnitude (%12.3f > %12.3f)\n",
1164 powers[wis[0]+1][ldir], fork_pow_thresh);
1165 fprintf(logfp,
1166 " 5. Right Power Magnitude (%12.3f > %12.3f)\n",
1167 powers[wis[0]+1][rdir], fork_pow_thresh);
1168 fprintf(logfp, " PASSED\n");
1169 fprintf(logfp,
1170 " Selected Direction = %d\n", powmax_dirs[wis[0]]);
1171 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
1172
1173 /* If ALL the above criteria hold, then return the direction */
1174 /* of the largest max power. */
1175 473 return(powmax_dirs[wis[0]]);
1176 }
1177 }
1178
1179 /* Otherwise test failed. */
1180
1181 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
1182 if(!firstpart){
1183 fprintf(logfp,
1184 " 1. Power Magnitude ( > %12.3f)\n",
1185 lfsparms->powmax_min);
1186 fprintf(logfp,
1187 " 2. Norm Power Magnitude ( > %9.3f)\n",
1188 fork_pownorm_min);
1189 fprintf(logfp,
1190 " 3. Low Freq Wave Magnitude ( <= %12.3f)\n",
1191 lfsparms->powmax_max);
1192 }
1193 else{
1194 fprintf(logfp, " 4. Left Power Magnitude (%12.3f > %12.3f)\n",
1195 powers[wis[0]+1][ldir], fork_pow_thresh);
1196 fprintf(logfp, " 5. Right Power Magnitude (%12.3f > %12.3f)\n",
1197 powers[wis[0]+1][rdir], fork_pow_thresh);
1198 }
1199 fprintf(logfp, " FAILED\n");
1200 } /* Close scope of firstpart */
1201 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
1202
1203 return(INVALID_DIR);
1204 }
1205
1206 /*************************************************************************
1207 **************************************************************************
1208 #cat: remove_incon_dirs - Takes a vector of integer directions and removes
1209 #cat: individual directions that are too weak or inconsistent.
1210 #cat: Directions are tested from the center of the IMAP working
1211 #cat: outward in concentric squares, and the process resets to
1212 #cat: the center and continues until no changes take place during
1213 #cat: a complete pass.
1214
1215 Input:
1216 imap - vector of IMAP integer directions
1217 mw - width (in blocks) of the IMAP
1218 mh - height (in blocks) of the IMAP
1219 dir2rad - lookup table for converting integer directions
1220 lfsparms - parameters and thresholds for controlling LFS
1221 Output:
1222 imap - vector of pruned input values
1223 **************************************************************************/
1224 96 void remove_incon_dirs(int *imap, const int mw, const int mh,
1225 const DIR2RAD *dir2rad, const LFSPARMS *lfsparms)
1226 {
1227 96 int cx, cy;
1228 96 int *iptr;
1229 96 int nremoved;
1230 96 int lbox, rbox, tbox, bbox;
1231
1232 #ifdef LOG_REPORT
1233 { int numpass = 0;
1234 fprintf(logfp, "REMOVE MAP\n");
1235 #endif
1236
1237 /* Compute center coords of IMAP */
1238 96 cx = mw>>1;
1239 96 cy = mh>>1;
1240
1241 /* Do pass, while directions have been removed in a pass ... */
1242 328 do{
1243
1244 #ifdef LOG_REPORT
1245 /* Count number of complete passes through IMAP */
1246 ++numpass;
1247 fprintf(logfp, " PASS = %d\n", numpass);
1248 #endif
1249
1250 /* Reinitialize number of removed directions to 0 */
1251 328 nremoved = 0;
1252
1253 /* Start at center */
1254 328 iptr = imap + (cy * mw) + cx;
1255 /* If valid IMAP direction and test for removal is true ... */
1256
4/4
✓ Branch 0 taken 196 times.
✓ Branch 1 taken 132 times.
✓ Branch 2 taken 127 times.
✓ Branch 3 taken 5 times.
 460 if((*iptr != INVALID_DIR)&&
1257 132 (remove_dir(imap, cx, cy, mw, mh, dir2rad, lfsparms))){
1258
1259 /* Set to INVALID */
1260 5 *iptr = INVALID_DIR;
1261 /* Bump number of removed IMAP directions */
1262 5 nremoved++;
1263 }
1264
1265 /* Initialize side indices of concentric boxes */
1266 328 lbox = cx-1;
1267 328 tbox = cy-1;
1268 328 rbox = cx+1;
1269 328 bbox = cy+1;
1270
1271 /* Grow concentric boxes, until ALL edges of imap are exceeded */
1272
4/4
✓ Branch 0 taken 5195 times.
✓ Branch 1 taken 1271 times.
✓ Branch 2 taken 943 times.
✓ Branch 3 taken 328 times.
 6466 while((lbox >= 0) || (rbox < mw) || (tbox >= 0) || (bbox < mh)){
1273
1274 /* test top edge of box */
1275
2/2
✓ Branch 0 taken 5718 times.
✓ Branch 1 taken 420 times.
 6138 if(tbox >= 0)
1276 5718 nremoved += test_top_edge(lbox, tbox, rbox, bbox, imap, mw, mh,
1277 dir2rad, lfsparms);
1278
1279 /* test right edge of box */
1280
2/2
✓ Branch 0 taken 4899 times.
✓ Branch 1 taken 1239 times.
 6138 if(rbox < mw)
1281 4899 nremoved += test_right_edge(lbox, tbox, rbox, bbox, imap, mw, mh,
1282 dir2rad, lfsparms);
1283
1284 /* test bottom edge of box */
1285
2/2
✓ Branch 0 taken 5443 times.
✓ Branch 1 taken 695 times.
 6138 if(bbox < mh)
1286 5443 nremoved += test_bottom_edge(lbox, tbox, rbox, bbox, imap, mw, mh,
1287 dir2rad, lfsparms);
1288
1289 /* test left edge of box */
1290
2/2
✓ Branch 0 taken 5195 times.
✓ Branch 1 taken 943 times.
 6138 if(lbox >=0)
1291 5195 nremoved += test_left_edge(lbox, tbox, rbox, bbox, imap, mw, mh,
1292 dir2rad, lfsparms);
1293
1294 /* Resize current box */
1295 6138 lbox--;
1296 6138 tbox--;
1297 6138 rbox++;
1298 6138 bbox++;
1299 }
1300
2/2
✓ Branch 0 taken 232 times.
✓ Branch 1 taken 96 times.
 328 }while(nremoved);
1301
1302 #ifdef LOG_REPORT
1303 } /* Close scope of numpass */
1304 #endif
1305
1306 96 }
1307
1308 /*************************************************************************
1309 **************************************************************************
1310 #cat: test_top_edge - Walks the top edge of a concentric square in the IMAP,
1311 #cat: testing directions along the way to see if they should
1312 #cat: be removed due to being too weak or inconsistent with
1313 #cat: respect to their adjacent neighbors.
1314
1315 Input:
1316 lbox - left edge of current concentric square
1317 tbox - top edge of current concentric square
1318 rbox - right edge of current concentric square
1319 bbox - bottom edge of current concentric square
1320 imap - vector of IMAP integer directions
1321 mw - width (in blocks) of the IMAP
1322 mh - height (in blocks) of the IMAP
1323 dir2rad - lookup table for converting integer directions
1324 lfsparms - parameters and thresholds for controlling LFS
1325 Return Code:
1326 Positive - direction should be removed from IMAP
1327 Zero - direction should NOT be remove from IMAP
1328 **************************************************************************/
1329 5718 int test_top_edge(const int lbox, const int tbox, const int rbox,
1330 const int bbox, int *imap, const int mw, const int mh,
1331 const DIR2RAD *dir2rad, const LFSPARMS *lfsparms)
1332 {
1333 5718 int bx, by, sx, ex;
1334 5718 int *iptr, *sptr, *eptr;
1335 5718 int nremoved;
1336
1337 /* Initialize number of directions removed on edge to 0 */
1338 5718 nremoved = 0;
1339
1340 /* Set start pointer to top-leftmost point of box, or set it to */
1341 /* the leftmost point in the IMAP row (0), whichever is larger. */
1342 5718 sx = max(lbox, 0);
1343 5718 sptr = imap + (tbox*mw) + sx;
1344
1345 /* Set end pointer to either 1 point short of the top-rightmost */
1346 /* point of box, or set it to the rightmost point in the IMAP */
1347 /* row (lastx=mw-1), whichever is smaller. */
1348
2/2
✓ Branch 0 taken 4730 times.
✓ Branch 1 taken 988 times.
 5718 ex = min(rbox-1, mw-1);
1349 5718 eptr = imap + (tbox*mw) + ex;
1350
1351 /* For each point on box's edge ... */
1352 5718 for(iptr = sptr, bx = sx, by = tbox;
1353
2/2
✓ Branch 0 taken 97987 times.
✓ Branch 1 taken 5718 times.
 103705 iptr <= eptr;
1354 97987 iptr++, bx++){
1355 /* If valid IMAP direction and test for removal is true ... */
1356
4/4
✓ Branch 0 taken 34527 times.
✓ Branch 1 taken 63460 times.
✓ Branch 2 taken 62080 times.
✓ Branch 3 taken 1380 times.
 161447 if((*iptr != INVALID_DIR)&&
1357 63460 (remove_dir(imap, bx, by, mw, mh, dir2rad, lfsparms))){
1358 /* Set to INVALID */
1359 1380 *iptr = INVALID_DIR;
1360 /* Bump number of removed IMAP directions */
1361 1380 nremoved++;
1362 }
1363 }
1364
1365 /* Return the number of directions removed on edge */
1366 5718 return(nremoved);
1367 }
1368
1369 /*************************************************************************
1370 **************************************************************************
1371 #cat: test_right_edge - Walks the right edge of a concentric square in the
1372 #cat: IMAP, testing directions along the way to see if they
1373 #cat: should be removed due to being too weak or inconsistent
1374 #cat: with respect to their adjacent neighbors.
1375
1376 Input:
1377 lbox - left edge of current concentric square
1378 tbox - top edge of current concentric square
1379 rbox - right edge of current concentric square
1380 bbox - bottom edge of current concentric square
1381 imap - vector of IMAP integer directions
1382 mw - width (in blocks) of the IMAP
1383 mh - height (in blocks) of the IMAP
1384 dir2rad - lookup table for converting integer directions
1385 lfsparms - parameters and thresholds for controlling LFS
1386 Return Code:
1387 Positive - direction should be removed from IMAP
1388 Zero - direction should NOT be remove from IMAP
1389 **************************************************************************/
1390 4899 int test_right_edge(const int lbox, const int tbox, const int rbox,
1391 const int bbox, int *imap, const int mw, const int mh,
1392 const DIR2RAD *dir2rad, const LFSPARMS *lfsparms)
1393 {
1394 4899 int bx, by, sy, ey;
1395 4899 int *iptr, *sptr, *eptr;
1396 4899 int nremoved;
1397
1398 /* Initialize number of directions removed on edge to 0 */
1399 4899 nremoved = 0;
1400
1401 /* Set start pointer to top-rightmost point of box, or set it to */
1402 /* the topmost point in IMAP column (0), whichever is larger. */
1403 4899 sy = max(tbox, 0);
1404 4899 sptr = imap + (sy*mw) + rbox;
1405
1406 /* Set end pointer to either 1 point short of the bottom- */
1407 /* rightmost point of box, or set it to the bottommost point */
1408 /* in the IMAP column (lasty=mh-1), whichever is smaller. */
1409
2/2
✓ Branch 0 taken 4508 times.
✓ Branch 1 taken 391 times.
 4899 ey = min(bbox-1,mh-1);
1410 4899 eptr = imap + (ey*mw) + rbox;
1411
1412 /* For each point on box's edge ... */
1413 4899 for(iptr = sptr, bx = rbox, by = sy;
1414
2/2
✓ Branch 0 taken 79343 times.
✓ Branch 1 taken 4899 times.
 84242 iptr <= eptr;
1415 79343 iptr+=mw, by++){
1416 /* If valid IMAP direction and test for removal is true ... */
1417
4/4
✓ Branch 0 taken 25784 times.
✓ Branch 1 taken 53559 times.
✓ Branch 2 taken 52622 times.
✓ Branch 3 taken 937 times.
 132902 if((*iptr != INVALID_DIR)&&
1418 53559 (remove_dir(imap, bx, by, mw, mh, dir2rad, lfsparms))){
1419 /* Set to INVALID */
1420 937 *iptr = INVALID_DIR;
1421 /* Bump number of removed IMAP directions */
1422 937 nremoved++;
1423 }
1424 }
1425
1426 /* Return the number of directions removed on edge */
1427 4899 return(nremoved);
1428 }
1429
1430 /*************************************************************************
1431 **************************************************************************
1432 #cat: test_bottom_edge - Walks the bottom edge of a concentric square in the
1433 #cat: IMAP, testing directions along the way to see if they
1434 #cat: should be removed due to being too weak or inconsistent
1435 #cat: with respect to their adjacent neighbors.
1436 Input:
1437 lbox - left edge of current concentric square
1438 tbox - top edge of current concentric square
1439 rbox - right edge of current concentric square
1440 bbox - bottom edge of current concentric square
1441 imap - vector of IMAP integer directions
1442 mw - width (in blocks) of the IMAP
1443 mh - height (in blocks) of the IMAP
1444 dir2rad - lookup table for converting integer directions
1445 lfsparms - parameters and thresholds for controlling LFS
1446 Return Code:
1447 Positive - direction should be removed from IMAP
1448 Zero - direction should NOT be remove from IMAP
1449 **************************************************************************/
1450 5443 int test_bottom_edge(const int lbox, const int tbox, const int rbox,
1451 const int bbox, int *imap, const int mw, const int mh,
1452 const DIR2RAD *dir2rad, const LFSPARMS *lfsparms)
1453 {
1454 5443 int bx, by, sx, ex;
1455 5443 int *iptr, *sptr, *eptr;
1456 5443 int nremoved;
1457
1458 /* Initialize number of directions removed on edge to 0 */
1459 5443 nremoved = 0;
1460
1461 /* Set start pointer to bottom-rightmost point of box, or set it to the */
1462 /* rightmost point in the IMAP ROW (lastx=mw-1), whichever is smaller. */
1463 5443 sx = min(rbox, mw-1);
1464 5443 sptr = imap + (bbox*mw) + sx;
1465
1466 /* Set end pointer to either 1 point short of the bottom- */
1467 /* lefttmost point of box, or set it to the leftmost point */
1468 /* in the IMAP row (x=0), whichever is larger. */
1469 5443 ex = max(lbox-1, 0);
1470 5443 eptr = imap + (bbox*mw) + ex;
1471
1472 /* For each point on box's edge ... */
1473 5443 for(iptr = sptr, bx = sx, by = bbox;
1474
2/2
✓ Branch 0 taken 98932 times.
✓ Branch 1 taken 5443 times.
 104375 iptr >= eptr;
1475 98932 iptr--, bx--){
1476 /* If valid IMAP direction and test for removal is true ... */
1477
4/4
✓ Branch 0 taken 26222 times.
✓ Branch 1 taken 72710 times.
✓ Branch 2 taken 71958 times.
✓ Branch 3 taken 752 times.
 171642 if((*iptr != INVALID_DIR)&&
1478 72710 (remove_dir(imap, bx, by, mw, mh, dir2rad, lfsparms))){
1479 /* Set to INVALID */
1480 752 *iptr = INVALID_DIR;
1481 /* Bump number of removed IMAP directions */
1482 752 nremoved++;
1483 }
1484 }
1485
1486 /* Return the number of directions removed on edge */
1487 5443 return(nremoved);
1488 }
1489
1490 /*************************************************************************
1491 **************************************************************************
1492 #cat: test_left_edge - Walks the left edge of a concentric square in the IMAP,
1493 #cat: testing directions along the way to see if they should
1494 #cat: be removed due to being too weak or inconsistent with
1495 #cat: respect to their adjacent neighbors.
1496
1497 Input:
1498 lbox - left edge of current concentric square
1499 tbox - top edge of current concentric square
1500 rbox - right edge of current concentric square
1501 bbox - bottom edge of current concentric square
1502 imap - vector of IMAP integer directions
1503 mw - width (in blocks) of the IMAP
1504 mh - height (in blocks) of the IMAP
1505 dir2rad - lookup table for converting integer directions
1506 lfsparms - parameters and thresholds for controlling LFS
1507 Return Code:
1508 Positive - direction should be removed from IMAP
1509 Zero - direction should NOT be remove from IMAP
1510 **************************************************************************/
1511 5195 int test_left_edge(const int lbox, const int tbox, const int rbox,
1512 const int bbox, int *imap, const int mw, const int mh,
1513 const DIR2RAD *dir2rad, const LFSPARMS *lfsparms)
1514 {
1515 5195 int bx, by, sy, ey;
1516 5195 int *iptr, *sptr, *eptr;
1517 5195 int nremoved;
1518
1519 /* Initialize number of directions removed on edge to 0 */
1520 5195 nremoved = 0;
1521
1522 /* Set start pointer to bottom-leftmost point of box, or set it to */
1523 /* the bottommost point in IMAP column (lasty=mh-1), whichever */
1524 /* is smaller. */
1525 5195 sy = min(bbox, mh-1);
1526 5195 sptr = imap + (sy*mw) + lbox;
1527
1528 /* Set end pointer to either 1 point short of the top-leftmost */
1529 /* point of box, or set it to the topmost point in the IMAP */
1530 /* column (y=0), whichever is larger. */
1531 5195 ey = max(tbox-1, 0);
1532 5195 eptr = imap + (ey*mw) + lbox;
1533
1534 /* For each point on box's edge ... */
1535 5195 for(iptr = sptr, bx = lbox, by = sy;
1536
2/2
✓ Branch 0 taken 97241 times.
✓ Branch 1 taken 5195 times.
 102436 iptr >= eptr;
1537 97241 iptr-=mw, by--){
1538 /* If valid IMAP direction and test for removal is true ... */
1539
4/4
✓ Branch 0 taken 26919 times.
✓ Branch 1 taken 70322 times.
✓ Branch 2 taken 69174 times.
✓ Branch 3 taken 1148 times.
 167563 if((*iptr != INVALID_DIR)&&
1540 70322 (remove_dir(imap, bx, by, mw, mh, dir2rad, lfsparms))){
1541 /* Set to INVALID */
1542 1148 *iptr = INVALID_DIR;
1543 /* Bump number of removed IMAP directions */
1544 1148 nremoved++;
1545 }
1546 }
1547
1548 /* Return the number of directions removed on edge */
1549 5195 return(nremoved);
1550 }
1551
1552 /*************************************************************************
1553 **************************************************************************
1554 #cat: remove_dir - Determines if an IMAP direction should be removed based
1555 #cat: on analyzing its adjacent neighbors
1556
1557 Input:
1558 imap - vector of IMAP integer directions
1559 mx - IMAP X-coord of the current direction being tested
1560 my - IMPA Y-coord of the current direction being tested
1561 mw - width (in blocks) of the IMAP
1562 mh - height (in blocks) of the IMAP
1563 dir2rad - lookup table for converting integer directions
1564 lfsparms - parameters and thresholds for controlling LFS
1565 Return Code:
1566 Positive - direction should be removed from IMAP
1567 Zero - direction should NOT be remove from IMAP
1568 **************************************************************************/
1569 260183 int remove_dir(int *imap, const int mx, const int my,
1570 const int mw, const int mh, const DIR2RAD *dir2rad,
1571 const LFSPARMS *lfsparms)
1572 {
1573 260183 int avrdir, nvalid, dist;
1574 260183 double dir_strength;
1575
1576 /* Compute average direction from neighbors, returning the */
1577 /* number of valid neighbors used in the computation, and */
1578 /* the "strength" of the average direction. */
1579 260183 average_8nbr_dir(&avrdir, &dir_strength, &nvalid, imap, mx, my, mw, mh,
1580 dir2rad);
1581
1582 /* Conduct valid neighbor test (Ex. thresh==3) */
1583
2/2
✓ Branch 0 taken 258033 times.
✓ Branch 1 taken 2150 times.
 260183 if(nvalid < lfsparms->rmv_valid_nbr_min){
1584
1585 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
1586 fprintf(logfp, " BLOCK %2d (%2d, %2d)\n",
1587 mx+(my*mw), mx, my);
1588 fprintf(logfp, " Average NBR : %2d %6.3f %d\n",
1589 avrdir, dir_strength, nvalid);
1590 fprintf(logfp, " 1. Valid NBR (%d < %d)\n",
1591 nvalid, lfsparms->rmv_valid_nbr_min);
1592 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
1593
1594 return(1);
1595 }
1596
1597 /* If stregnth of average neighbor direction is large enough to */
1598 /* put credence in ... (Ex. thresh==0.2) */
1599
2/2
✓ Branch 0 taken 248735 times.
✓ Branch 1 taken 9298 times.
 258033 if(dir_strength >= lfsparms->dir_strength_min){
1600
1601 /* Conduct direction distance test (Ex. thresh==3) */
1602 /* Compute minimum absolute distance between current and */
1603 /* average directions accounting for wrapping from 0 to NDIRS. */
1604 248735 dist = abs(avrdir - *(imap+(my*mw)+mx));
1605 248735 dist = min(dist, dir2rad->ndirs-dist);
1606
2/2
✓ Branch 0 taken 2072 times.
✓ Branch 1 taken 246663 times.
 248735 if(dist > lfsparms->dir_distance_max){
1607
1608 #ifdef LOG_REPORT /*vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv*/
1609 fprintf(logfp, " BLOCK %2d (%2d, %2d)\n",
1610 mx+(my*mw), mx, my);
1611 fprintf(logfp, " Average NBR : %2d %6.3f %d\n",
1612 avrdir, dir_strength, nvalid);
1613 fprintf(logfp, " 1. Valid NBR (%d < %d)\n",
1614 nvalid, lfsparms->rmv_valid_nbr_min);
1615 fprintf(logfp, " 2. Direction Strength (%6.3f >= %6.3f)\n",
1616 dir_strength, lfsparms->dir_strength_min);
1617 fprintf(logfp, " Current Dir = %d, Average Dir = %d\n",
1618 *(imap+(my*mw)+mx), avrdir);
1619 fprintf(logfp, " 3. Direction Distance (%d > %d)\n",
1620 dist, lfsparms->dir_distance_max);
1621 #endif /*^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*/
1622
1623 2072 return(2);
1624 }
1625 }
1626
1627 /* Otherwise, the strength of the average direciton is not strong enough */
1628 /* to put credence in, so leave the current block's directon alone. */
1629
1630 return(0);
1631 }
1632
1633 /*************************************************************************
1634 **************************************************************************
1635 #cat: average_8nbr_dir - Given an IMAP direction, computes an average
1636 #cat: direction from its adjacent 8 neighbors returning
1637 #cat: the average direction, its strength, and the
1638 #cat: number of valid direction in the neighborhood.
1639
1640 Input:
1641 imap - vector of IMAP integer directions
1642 mx - IMAP X-coord of the current direction
1643 my - IMPA Y-coord of the current direction
1644 mw - width (in blocks) of the IMAP
1645 mh - height (in blocks) of the IMAP
1646 dir2rad - lookup table for converting integer directions
1647 Output:
1648 avrdir - the average direction computed from neighbors
1649 dir_strenght - the strength of the average direction
1650 nvalid - the number of valid directions used to compute the
1651 average
1652 **************************************************************************/
1653 352461 void average_8nbr_dir(int *avrdir, double *dir_strength, int *nvalid,
1654 int *imap, const int mx, const int my,
1655 const int mw, const int mh,
1656 const DIR2RAD *dir2rad)
1657 {
1658 352461 int *iptr;
1659 352461 int e,w,n,s;
1660 352461 double cospart, sinpart;
1661 352461 double pi2, pi_factor, theta;
1662 352461 double avr;
1663
1664 /* Compute neighbor coordinates to current IMAP direction */
1665 352461 e = mx+1; /* East */
1666 352461 w = mx-1; /* West */
1667 352461 n = my-1; /* North */
1668 352461 s = my+1; /* South */
1669
1670 /* Intialize accumulators */
1671 352461 *nvalid = 0;
1672 352461 cospart = 0.0;
1673 352461 sinpart = 0.0;
1674
1675 /* 1. Test NW */
1676 /* If NW point within IMAP boudaries ... */
1677
2/2
✓ Branch 0 taken 332827 times.
✓ Branch 1 taken 19634 times.
 352461 if((w >= 0) && (n >= 0)){
1678 332827 iptr = imap + (n*mw) + w;
1679 /* If valid direction ... */
1680
2/2
✓ Branch 0 taken 291383 times.
✓ Branch 1 taken 41444 times.
 332827 if(*iptr != INVALID_DIR){
1681 /* Accumulate cosine and sine components of the direction */
1682 291383 cospart += dir2rad->cos[*iptr];
1683 291383 sinpart += dir2rad->sin[*iptr];
1684 /* Bump number of accumulated directions */
1685 291383 (*nvalid)++;
1686 }
1687 }
1688
1689 /* 2. Test N */
1690 /* If N point within IMAP boudaries ... */
1691
2/2
✓ Branch 0 taken 343874 times.
✓ Branch 1 taken 8587 times.
 352461 if(n >= 0){
1692 343874 iptr = imap + (n*mw) + mx;
1693 /* If valid direction ... */
1694
2/2
✓ Branch 0 taken 305668 times.
✓ Branch 1 taken 38206 times.
 343874 if(*iptr != INVALID_DIR){
1695 /* Accumulate cosine and sine components of the direction */
1696 305668 cospart += dir2rad->cos[*iptr];
1697 305668 sinpart += dir2rad->sin[*iptr];
1698 /* Bump number of accumulated directions */
1699 305668 (*nvalid)++;
1700 }
1701 }
1702
1703 /* 3. Test NE */
1704 /* If NE point within IMAP boudaries ... */
1705
2/2
✓ Branch 0 taken 334128 times.
✓ Branch 1 taken 18333 times.
 352461 if((e < mw) && (n >= 0)){
1706 334128 iptr = imap + (n*mw) + e;
1707 /* If valid direction ... */
1708
2/2
✓ Branch 0 taken 291706 times.
✓ Branch 1 taken 42422 times.
 334128 if(*iptr != INVALID_DIR){
1709 /* Accumulate cosine and sine components of the direction */
1710 291706 cospart += dir2rad->cos[*iptr];
1711 291706 sinpart += dir2rad->sin[*iptr];
1712 /* Bump number of accumulated directions */
1713 291706 (*nvalid)++;
1714 }
1715 }
1716
1717 /* 4. Test E */
1718 /* If E point within IMAP boudaries ... */
1719
2/2
✓ Branch 0 taken 342501 times.
✓ Branch 1 taken 9960 times.
 352461 if(e < mw){
1720 342501 iptr = imap + (my*mw) + e;
1721 /* If valid direction ... */
1722
2/2
✓ Branch 0 taken 303000 times.
✓ Branch 1 taken 39501 times.
 342501 if(*iptr != INVALID_DIR){
1723 /* Accumulate cosine and sine components of the direction */
1724 303000 cospart += dir2rad->cos[*iptr];
1725 303000 sinpart += dir2rad->sin[*iptr];
1726 /* Bump number of accumulated directions */
1727 303000 (*nvalid)++;
1728 }
1729 }
1730
1731 /* 5. Test SE */
1732 /* If SE point within IMAP boudaries ... */
1733
2/2
✓ Branch 0 taken 331886 times.
✓ Branch 1 taken 20575 times.
 352461 if((e < mw) && (s < mh)){
1734 331886 iptr = imap + (s*mw) + e;
1735 /* If valid direction ... */
1736
2/2
✓ Branch 0 taken 290021 times.
✓ Branch 1 taken 41865 times.
 331886 if(*iptr != INVALID_DIR){
1737 /* Accumulate cosine and sine components of the direction */
1738 290021 cospart += dir2rad->cos[*iptr];
1739 290021 sinpart += dir2rad->sin[*iptr];
1740 /* Bump number of accumulated directions */
1741 290021 (*nvalid)++;
1742 }
1743 }
1744
1745 /* 6. Test S */
1746 /* If S point within IMAP boudaries ... */
1747
2/2
✓ Branch 0 taken 341524 times.
✓ Branch 1 taken 10937 times.
 352461 if(s < mh){
1748 341524 iptr = imap + (s*mw) + mx;
1749 /* If valid direction ... */
1750
2/2
✓ Branch 0 taken 303870 times.
✓ Branch 1 taken 37654 times.
 341524 if(*iptr != INVALID_DIR){
1751 /* Accumulate cosine and sine components of the direction */
1752 303870 cospart += dir2rad->cos[*iptr];
1753 303870 sinpart += dir2rad->sin[*iptr];
1754 /* Bump number of accumulated directions */
1755 303870 (*nvalid)++;
1756 }
1757 }
1758
1759 /* 7. Test SW */
1760 /* If SW point within IMAP boudaries ... */
1761
2/2
✓ Branch 0 taken 330415 times.
✓ Branch 1 taken 22046 times.
 352461 if((w >= 0) && (s < mh)){
1762 330415 iptr = imap + (s*mw) + w;
1763 /* If valid direction ... */
1764
2/2
✓ Branch 0 taken 290298 times.
✓ Branch 1 taken 40117 times.
 330415 if(*iptr != INVALID_DIR){
1765 /* Accumulate cosine and sine components of the direction */
1766 290298 cospart += dir2rad->cos[*iptr];
1767 290298 sinpart += dir2rad->sin[*iptr];
1768 /* Bump number of accumulated directions */
1769 290298 (*nvalid)++;
1770 }
1771 }
1772
1773 /* 8. Test W */
1774 /* If W point within IMAP boudaries ... */
1775
2/2
✓ Branch 0 taken 341061 times.
✓ Branch 1 taken 11400 times.
 352461 if(w >= 0){
1776 341061 iptr = imap + (my*mw) + w;
1777 /* If valid direction ... */
1778
2/2
✓ Branch 0 taken 304074 times.
✓ Branch 1 taken 36987 times.
 341061 if(*iptr != INVALID_DIR){
1779 /* Accumulate cosine and sine components of the direction */
1780 304074 cospart += dir2rad->cos[*iptr];
1781 304074 sinpart += dir2rad->sin[*iptr];
1782 /* Bump number of accumulated directions */
1783 304074 (*nvalid)++;
1784 }
1785 }
1786
1787 /* If there were no neighbors found with valid direction ... */
1788
2/2
✓ Branch 0 taken 6620 times.
✓ Branch 1 taken 345841 times.
 352461 if(*nvalid == 0){
1789 /* Return INVALID direction. */
1790 6620 *dir_strength = 0;
1791 6620 *avrdir = INVALID_DIR;
1792 6620 return;
1793 }
1794
1795 /* Compute averages of accumulated cosine and sine direction components */
1796 345841 cospart /= (double)(*nvalid);
1797 345841 sinpart /= (double)(*nvalid);
1798
1799 /* Compute directional strength as hypotenuse (without sqrt) of average */
1800 /* cosine and sine direction components. Believe this value will be on */
1801 /* the range of [0 .. 1]. */
1802 345841 *dir_strength = (cospart * cospart) + (sinpart * sinpart);
1803 /* Need to truncate precision so that answers are consistent */
1804 /* on different computer architectures when comparing doubles. */
1805 345841 *dir_strength = trunc_dbl_precision(*dir_strength, TRUNC_SCALE);
1806
1807 /* If the direction strength is not sufficiently high ... */
1808
2/2
✓ Branch 0 taken 11731 times.
✓ Branch 1 taken 334110 times.
 345841 if(*dir_strength < DIR_STRENGTH_MIN){
1809 /* Return INVALID direction. */
1810 11731 *dir_strength = 0;
1811 11731 *avrdir = INVALID_DIR;
1812 11731 return;
1813 }
1814
1815 /* Compute angle (in radians) from Arctan of avarage */
1816 /* cosine and sine direction components. I think this order */
1817 /* is necessary because 0 direction is vertical and positive */
1818 /* direction is clockwise. */
1819 334110 theta = atan2(sinpart, cospart);
1820
1821 /* Atan2 returns theta on range [-PI..PI]. Adjust theta so that */
1822 /* it is on the range [0..2PI]. */
1823 334110 pi2 = 2*M_PI;
1824 334110 theta += pi2;
1825 334110 theta = fmod(theta, pi2);
1826
1827 /* Pi_factor sets the period of the trig functions to NDIRS units in x. */
1828 /* For example, if NDIRS==16, then pi_factor = 2(PI/16) = .3926... */
1829 /* Dividing theta (in radians) by this factor ((1/pi_factor)==2.546...) */
1830 /* will produce directions on the range [0..NDIRS]. */
1831 334110 pi_factor = pi2/(double)dir2rad->ndirs; /* 2(M_PI/ndirs) */
1832
1833 /* Round off the direction and return it as an average direction */
1834 /* for the neighborhood. */
1835 334110 avr = theta / pi_factor;
1836 /* Need to truncate precision so that answers are consistent */
1837 /* on different computer architectures when rounding doubles. */
1838
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 334110 times.
 334110 avr = trunc_dbl_precision(avr, TRUNC_SCALE);
1839
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 334110 times.
 334110 *avrdir = sround(avr);
1840
1841 /* Really do need to map values > NDIRS back onto [0..NDIRS) range. */
1842 334110 *avrdir %= dir2rad->ndirs;
1843 }
1844
1845 /*************************************************************************
1846 **************************************************************************
1847 #cat: num_valid_8nbrs - Given a block in an IMAP, counts the number of
1848 #cat: immediate neighbors that have a valid IMAP direction.
1849
1850 Input:
1851 imap - 2-D vector of directional ridge flows
1852 mx - horizontal coord of current block in IMAP
1853 my - vertical coord of current block in IMAP
1854 mw - width (in blocks) of the IMAP
1855 mh - height (in blocks) of the IMAP
1856 Return Code:
1857 Non-negative - the number of valid IMAP neighbors
1858 **************************************************************************/
1859 51995 int num_valid_8nbrs(int *imap, const int mx, const int my,
1860 const int mw, const int mh)
1861 {
1862 51995 int e_ind, w_ind, n_ind, s_ind;
1863 51995 int nvalid;
1864
1865 /* Initialize VALID IMAP counter to zero. */
1866 51995 nvalid = 0;
1867
1868 /* Compute neighbor coordinates to current IMAP direction */
1869 51995 e_ind = mx+1; /* East index */
1870 51995 w_ind = mx-1; /* West index */
1871 51995 n_ind = my-1; /* North index */
1872 51995 s_ind = my+1; /* South index */
1873
1874 /* 1. Test NW IMAP value. */
1875 /* If neighbor indices are within IMAP boundaries and it is VALID ... */
1876
4/4
✓ Branch 0 taken 3623 times.
✓ Branch 1 taken 48372 times.
✓ Branch 2 taken 8837 times.
✓ Branch 3 taken 39535 times.
 51995 if((w_ind >= 0) && (n_ind >= 0) && (*(imap + (n_ind*mw) + w_ind) >= 0))
1877 /* Bump VALID counter. */
1878 51995 nvalid++;
1879
1880 /* 2. Test N IMAP value. */
1881
4/4
✓ Branch 0 taken 50174 times.
✓ Branch 1 taken 1821 times.
✓ Branch 2 taken 39404 times.
✓ Branch 3 taken 10770 times.
 51995 if((n_ind >= 0) && (*(imap + (n_ind*mw) + mx) >= 0))
1882 39404 nvalid++;
1883
1884 /* 3. Test NE IMAP value. */
1885
4/4
✓ Branch 0 taken 48519 times.
✓ Branch 1 taken 3476 times.
✓ Branch 2 taken 39703 times.
✓ Branch 3 taken 8816 times.
 51995 if((n_ind >= 0) && (e_ind < mw) && (*(imap + (n_ind*mw) + e_ind) >= 0))
1886 39703 nvalid++;
1887
1888 /* 4. Test E IMAP value. */
1889
4/4
✓ Branch 0 taken 50292 times.
✓ Branch 1 taken 1703 times.
✓ Branch 2 taken 39460 times.
✓ Branch 3 taken 10832 times.
 51995 if((e_ind < mw) && (*(imap + (my*mw) + e_ind) >= 0))
1890 39460 nvalid++;
1891
1892 /* 5. Test SE IMAP value. */
1893
4/4
✓ Branch 0 taken 48645 times.
✓ Branch 1 taken 3350 times.
✓ Branch 2 taken 39825 times.
✓ Branch 3 taken 8820 times.
 51995 if((e_ind < mw) && (s_ind < mh) && (*(imap + (s_ind*mw) + e_ind) >= 0))
1894 39825 nvalid++;
1895
1896 /* 6. Test S IMAP value. */
1897
4/4
✓ Branch 0 taken 50300 times.
✓ Branch 1 taken 1695 times.
✓ Branch 2 taken 39564 times.
✓ Branch 3 taken 10736 times.
 51995 if((s_ind < mh) && (*(imap + (s_ind*mw) + mx) >= 0))
1898 39564 nvalid++;
1899
1900 /* 7. Test SW IMAP value. */
1901
4/4
✓ Branch 0 taken 48498 times.
✓ Branch 1 taken 3497 times.
✓ Branch 2 taken 39749 times.
✓ Branch 3 taken 8749 times.
 51995 if((w_ind >= 0) && (s_ind < mh) && (*(imap + (s_ind*mw) + w_ind) >= 0))
1902 39749 nvalid++;
1903
1904 /* 8. Test W IMAP value. */
1905
4/4
✓ Branch 0 taken 50145 times.
✓ Branch 1 taken 1850 times.
✓ Branch 2 taken 39314 times.
✓ Branch 3 taken 10831 times.
 51995 if((w_ind >= 0) && (*(imap + (my*mw) + w_ind) >= 0))
1906 39314 nvalid++;
1907
1908 /* Return number of neighbors with VALID IMAP values. */
1909 51995 return(nvalid);
1910 }
1911
1912 /*************************************************************************
1913 **************************************************************************
1914 #cat: smooth_imap - Takes a vector of integer directions and smooths them
1915 #cat: by analyzing the direction of adjacent neighbors.
1916
1917 Input:
1918 imap - vector of IMAP integer directions
1919 mw - width (in blocks) of the IMAP
1920 mh - height (in blocks) of the IMAP
1921 dir2rad - lookup table for converting integer directions
1922 lfsparms - parameters and thresholds for controlling LFS
1923 Output:
1924 imap - vector of smoothed input values
1925 **************************************************************************/
1926
1927 /*************************************************************************
1928 **************************************************************************
1929 #cat: gen_nmap - Computes an NMAP from its associated 2D vector of integer
1930 #cat: directions (IMAP). Each value in the NMAP either represents
1931 #cat: a direction of dominant ridge flow in a block of the input
1932 #cat: grayscale image, or it contains a codes describing why such
1933 #cat: a direction was not procuded.
1934 #cat: For example, blocks near areas of high-curvature (such as
1935 #cat: with cores and deltas) will not produce reliable IMAP
1936 #cat: directions.
1937
1938 Input:
1939 imap - associated input vector of IMAP directions
1940 mw - the width (in blocks) of the IMAP
1941 mh - the height (in blocks) of the IMAP
1942 lfsparms - parameters and thresholds for controlling LFS
1943 Output:
1944 optr - points to the created NMAP
1945 Return Code:
1946 Zero - successful completion
1947 Negative - system error
1948 **************************************************************************/
1949
1950 /*************************************************************************
1951 **************************************************************************
1952 #cat: vorticity - Measures the amount of cummulative curvature incurred
1953 #cat: among the IMAP neighbors of the given block.
1954
1955 Input:
1956 imap - 2D vector of ridge flow directions
1957 mx - horizontal coord of current IMAP block
1958 my - vertical coord of current IMAP block
1959 mw - width (in blocks) of the IMAP
1960 mh - height (in blocks) of the IMAP
1961 ndirs - number of possible directions in the IMAP
1962 Return Code:
1963 Non-negative - the measured vorticity among the neighbors
1964 **************************************************************************/
1965 16 int vorticity(int *imap, const int mx, const int my,
1966 const int mw, const int mh, const int ndirs)
1967 {
1968 16 int e_ind, w_ind, n_ind, s_ind;
1969 16 int nw_val, n_val, ne_val, e_val, se_val, s_val, sw_val, w_val;
1970 16 int vmeasure;
1971
1972 /* Compute neighbor coordinates to current IMAP direction */
1973 16 e_ind = mx+1; /* East index */
1974 16 w_ind = mx-1; /* West index */
1975 16 n_ind = my-1; /* North index */
1976 16 s_ind = my+1; /* South index */
1977
1978 /* 1. Get NW IMAP value. */
1979 /* If neighbor indices are within IMAP boundaries ... */
1980
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if((w_ind >= 0) && (n_ind >= 0))
1981 /* Set neighbor value to IMAP value. */
1982 16 nw_val = *(imap + (n_ind*mw) + w_ind);
1983 else
1984 /* Otherwise, set the neighbor value to INVALID. */
1985 nw_val = INVALID_DIR;
1986
1987 /* 2. Get N IMAP value. */
1988
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if(n_ind >= 0)
1989 16 n_val = *(imap + (n_ind*mw) + mx);
1990 else
1991 n_val = INVALID_DIR;
1992
1993 /* 3. Get NE IMAP value. */
1994
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if((n_ind >= 0) && (e_ind < mw))
1995 16 ne_val = *(imap + (n_ind*mw) + e_ind);
1996 else
1997 ne_val = INVALID_DIR;
1998
1999 /* 4. Get E IMAP value. */
2000
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if(e_ind < mw)
2001 16 e_val = *(imap + (my*mw) + e_ind);
2002 else
2003 e_val = INVALID_DIR;
2004
2005 /* 5. Get SE IMAP value. */
2006
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if((e_ind < mw) && (s_ind < mh))
2007 16 se_val = *(imap + (s_ind*mw) + e_ind);
2008 else
2009 se_val = INVALID_DIR;
2010
2011 /* 6. Get S IMAP value. */
2012
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if(s_ind < mh)
2013 16 s_val = *(imap + (s_ind*mw) + mx);
2014 else
2015 s_val = INVALID_DIR;
2016
2017 /* 7. Get SW IMAP value. */
2018
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if((w_ind >= 0) && (s_ind < mh))
2019 16 sw_val = *(imap + (s_ind*mw) + w_ind);
2020 else
2021 sw_val = INVALID_DIR;
2022
2023 /* 8. Get W IMAP value. */
2024
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 if(w_ind >= 0)
2025 16 w_val = *(imap + (my*mw) + w_ind);
2026 else
2027 w_val = INVALID_DIR;
2028
2029 /* Now that we have all IMAP neighbors, accumulate vorticity between */
2030 /* the neighboring directions. */
2031
2032 /* Initialize vorticity accumulator to zero. */
2033 16 vmeasure = 0;
2034
2035 /* 1. NW & N */
2036 16 accum_nbr_vorticity(&vmeasure, nw_val, n_val, ndirs);
2037
2038 /* 2. N & NE */
2039 16 accum_nbr_vorticity(&vmeasure, n_val, ne_val, ndirs);
2040
2041 /* 3. NE & E */
2042 16 accum_nbr_vorticity(&vmeasure, ne_val, e_val, ndirs);
2043
2044 /* 4. E & SE */
2045 16 accum_nbr_vorticity(&vmeasure, e_val, se_val, ndirs);
2046
2047 /* 5. SE & S */
2048 16 accum_nbr_vorticity(&vmeasure, se_val, s_val, ndirs);
2049
2050 /* 6. S & SW */
2051 16 accum_nbr_vorticity(&vmeasure, s_val, sw_val, ndirs);
2052
2053 /* 7. SW & W */
2054 16 accum_nbr_vorticity(&vmeasure, sw_val, w_val, ndirs);
2055
2056 /* 8. W & NW */
2057 16 accum_nbr_vorticity(&vmeasure, w_val, nw_val, ndirs);
2058
2059 /* Return the accumulated vorticity measure. */
2060 16 return(vmeasure);
2061 }
2062
2063 /*************************************************************************
2064 **************************************************************************
2065 #cat: accum_nbor_vorticity - Accumlates the amount of curvature measures
2066 #cat: between neighboring IMAP blocks.
2067
2068 Input:
2069 dir1 - first neighbor's integer IMAP direction
2070 dir2 - second neighbor's integer IMAP direction
2071 ndirs - number of possible IMAP directions
2072 Output:
2073 vmeasure - accumulated vorticity among neighbors measured so far
2074 **************************************************************************/
2075 128 void accum_nbr_vorticity(int *vmeasure, const int dir1, const int dir2,
2076 const int ndirs)
2077 {
2078 128 int dist;
2079
2080 /* Measure difference in direction between a pair of neighboring */
2081 /* directions. */
2082 /* If both neighbors are not equal and both are VALID ... */
2083
4/4
✓ Branch 0 taken 62 times.
✓ Branch 1 taken 66 times.
✓ Branch 2 taken 51 times.
✓ Branch 3 taken 11 times.
 128 if((dir1 != dir2) && (dir1 >= 0)&&(dir2 >= 0)){
2084 /* Measure the clockwise distance from the first to the second */
2085 /* directions. */
2086 51 dist = dir2 - dir1;
2087 /* If dist is negative, then clockwise distance must wrap around */
2088 /* the high end of the direction range. For example: */
2089 /* dir1 = 8 */
2090 /* dir2 = 3 */
2091 /* and ndirs = 16 */
2092 /* 3 - 8 = -5 */
2093 /* so 16 - 5 = 11 (the clockwise distance from 8 to 3) */
2094
2/2
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 23 times.
 51 if(dist < 0)
2095 28 dist += ndirs;
2096 /* If the change in clockwise direction is larger than 90 degrees as */
2097 /* in total the total number of directions covers 180 degrees. */
2098
2/2
✓ Branch 0 taken 28 times.
✓ Branch 1 taken 23 times.
 51 if(dist > (ndirs>>1))
2099 /* Decrement the vorticity measure. */
2100 28 (*vmeasure)--;
2101 else
2102 /* Otherwise, bump the vorticity measure. */
2103 23 (*vmeasure)++;
2104 }
2105 /* Otherwise both directions are either equal or */
2106 /* one or both directions are INVALID, so ignore. */
2107 128 }
2108
2109 /*************************************************************************
2110 **************************************************************************
2111 #cat: curvature - Measures the largest change in direction between the
2112 #cat: current IMAP direction and its immediate neighbors.
2113
2114 Input:
2115 imap - 2D vector of ridge flow directions
2116 mx - horizontal coord of current IMAP block
2117 my - vertical coord of current IMAP block
2118 mw - width (in blocks) of the IMAP
2119 mh - height (in blocks) of the IMAP
2120 ndirs - number of possible directions in the IMAP
2121 Return Code:
2122 Non-negative - maximum change in direction found (curvature)
2123 Negative - No valid neighbor found to measure change in direction
2124 **************************************************************************/
2125 38951 int curvature(int *imap, const int mx, const int my,
2126 const int mw, const int mh, const int ndirs)
2127 {
2128 38951 int *iptr;
2129 38951 int e_ind, w_ind, n_ind, s_ind;
2130 38951 int nw_val, n_val, ne_val, e_val, se_val, s_val, sw_val, w_val;
2131 38951 int cmeasure, dist;
2132
2133 /* Compute neighbor coordinates to current IMAP direction */
2134 38951 e_ind = mx+1; /* East index */
2135 38951 w_ind = mx-1; /* West index */
2136 38951 n_ind = my-1; /* North index */
2137 38951 s_ind = my+1; /* South index */
2138
2139 /* 1. Get NW IMAP value. */
2140 /* If neighbor indices are within IMAP boundaries ... */
2141
1/2
✓ Branch 0 taken 38951 times.
✗ Branch 1 not taken.
 38951 if((w_ind >= 0) && (n_ind >= 0))
2142 /* Set neighbor value to IMAP value. */
2143 38951 nw_val = *(imap + (n_ind*mw) + w_ind);
2144 else
2145 /* Otherwise, set the neighbor value to INVALID. */
2146 nw_val = INVALID_DIR;
2147
2148 /* 2. Get N IMAP value. */
2149
1/2
✓ Branch 0 taken 38951 times.
✗ Branch 1 not taken.
 38951 if(n_ind >= 0)
2150 38951 n_val = *(imap + (n_ind*mw) + mx);
2151 else
2152 n_val = INVALID_DIR;
2153
2154 /* 3. Get NE IMAP value. */
2155
1/2
✓ Branch 0 taken 38951 times.
✗ Branch 1 not taken.
 38951 if((n_ind >= 0) && (e_ind < mw))
2156 38951 ne_val = *(imap + (n_ind*mw) + e_ind);
2157 else
2158 ne_val = INVALID_DIR;
2159
2160 /* 4. Get E IMAP value. */
2161
1/2
✓ Branch 0 taken 38951 times.
✗ Branch 1 not taken.
 38951 if(e_ind < mw)
2162 38951 e_val = *(imap + (my*mw) + e_ind);
2163 else
2164 e_val = INVALID_DIR;
2165
2166 /* 5. Get SE IMAP value. */
2167
1/2
✓ Branch 0 taken 38951 times.
✗ Branch 1 not taken.
 38951 if((e_ind < mw) && (s_ind < mh))
2168 38951 se_val = *(imap + (s_ind*mw) + e_ind);
2169 else
2170 se_val = INVALID_DIR;
2171
2172 /* 6. Get S IMAP value. */
2173
1/2
✓ Branch 0 taken 38951 times.
✗ Branch 1 not taken.
 38951 if(s_ind < mh)
2174 38951 s_val = *(imap + (s_ind*mw) + mx);
2175 else
2176 s_val = INVALID_DIR;
2177
2178 /* 7. Get SW IMAP value. */
2179
1/2
✓ Branch 0 taken 38951 times.
✗ Branch 1 not taken.
 38951 if((w_ind >= 0) && (s_ind < mh))
2180 38951 sw_val = *(imap + (s_ind*mw) + w_ind);
2181 else
2182 sw_val = INVALID_DIR;
2183
2184 /* 8. Get W IMAP value. */
2185
1/2
✓ Branch 0 taken 38951 times.
✗ Branch 1 not taken.
 38951 if(w_ind >= 0)
2186 38951 w_val = *(imap + (my*mw) + w_ind);
2187 else
2188 w_val = INVALID_DIR;
2189
2190 /* Now that we have all IMAP neighbors, determine largest change in */
2191 /* direction from current block to each of its 8 VALID neighbors. */
2192
2193 /* Initialize pointer to current IMAP value. */
2194 38951 iptr = imap + (my*mw) + mx;
2195
2196 /* Initialize curvature measure to negative as closest_dir_dist() */
2197 /* always returns -1=INVALID or a positive value. */
2198 38951 cmeasure = -1;
2199
2200 /* 1. With NW */
2201 /* Compute closest distance between neighboring directions. */
2202 38951 dist = closest_dir_dist(*iptr, nw_val, ndirs);
2203 /* Keep track of maximum. */
2204 38951 if(dist > cmeasure)
2205 cmeasure = dist;
2206
2207 /* 2. With N */
2208 38951 dist = closest_dir_dist(*iptr, n_val, ndirs);
2209 38951 if(dist > cmeasure)
2210 cmeasure = dist;
2211
2212 /* 3. With NE */
2213 38951 dist = closest_dir_dist(*iptr, ne_val, ndirs);
2214 38951 if(dist > cmeasure)
2215 cmeasure = dist;
2216
2217 /* 4. With E */
2218 38951 dist = closest_dir_dist(*iptr, e_val, ndirs);
2219 38951 if(dist > cmeasure)
2220 cmeasure = dist;
2221
2222 /* 5. With SE */
2223 38951 dist = closest_dir_dist(*iptr, se_val, ndirs);
2224 38951 if(dist > cmeasure)
2225 cmeasure = dist;
2226
2227 /* 6. With S */
2228 38951 dist = closest_dir_dist(*iptr, s_val, ndirs);
2229 38951 if(dist > cmeasure)
2230 cmeasure = dist;
2231
2232 /* 7. With SW */
2233 38951 dist = closest_dir_dist(*iptr, sw_val, ndirs);
2234 38951 if(dist > cmeasure)
2235 cmeasure = dist;
2236
2237 /* 8. With W */
2238 38951 dist = closest_dir_dist(*iptr, w_val, ndirs);
2239 38951 if(dist > cmeasure)
2240 cmeasure = dist;
2241
2242 /* Return maximum difference between current block's IMAP direction */
2243 /* and the rest of its VALID neighbors. */
2244 38951 return(cmeasure);
2245 }
2246