GCC Code Coverage Report

Directory: ./
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 90.9% 577 / 0 / 635
Functions: 100.0% 20 / 0 / 20
Branches: 85.6% 274 / 0 / 320
libfprint/nbis/mindtct/maps.c
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 50 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 50 int *direction_map, *low_contrast_map, *low_flow_map, *high_curve_map;
133 50 int mw, mh, iw, ih;
134 50 int *blkoffs;
135 50 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 2 → 3 not taken.
✓ Branch 2 → 5 taken 50 times.
 50 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 50 iw = pw - (dftgrids->pad<<1);
146 50 ih = ph - (dftgrids->pad<<1);
147
1/2
✓ Branch 6 → 7 taken 50 times.
✗ Branch 6 → 36 not taken.
 50 if((ret = block_offsets(&blkoffs, &mw, &mh, iw, ih,
148 dftgrids->pad, lfsparms->blocksize))){
149 return(ret);
150 }
151
152 /* 2. Generate initial Direction Map and Low Contrast Map*/
153
1/2
✗ Branch 8 → 9 not taken.
✓ Branch 8 → 11 taken 50 times.
 50 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 12 → 13 not taken.
✓ Branch 12 → 17 taken 50 times.
 50 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 50 remove_incon_dirs(direction_map, mw, mh, dir2rad, lfsparms);
170
171
172 /* 4. Smooth Direction Map values with their neighbors */
173 50 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 20 → 21 not taken.
✓ Branch 20 → 25 taken 50 times.
 50 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 50 remove_incon_dirs(direction_map, mw, mh, dir2rad, lfsparms);
190
191 /* 7. Smooth Direction Map values with their neighbors. */
192 50 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 50 set_margin_blocks(direction_map, mw, mh, INVALID_DIR);
197
198 /* 9. Generate High Curvature Map from interpolated Direction Map. */
199
1/2
✗ Branch 29 → 30 not taken.
✓ Branch 29 → 34 taken 50 times.
 50 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 50 g_free(blkoffs);
209
210 50 *odmap = direction_map;
211 50 *olcmap = low_contrast_map;
212 50 *olfmap = low_flow_map;
213 50 *ohcmap = high_curve_map;
214 50 *omw = mw;
215 50 *omh = mh;
216 50 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 50 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 50 int *direction_map, *low_contrast_map, *low_flow_map;
263 50 int bi, bsize, blkdir;
264 50 int *wis, *powmax_dirs;
265 50 double **powers, *powmaxs, *pownorms;
266 50 int nstats;
267 50 int ret; /* return code */
268 50 int dft_offset;
269 50 int xminlimit, xmaxlimit, yminlimit, ymaxlimit;
270 50 int win_x, win_y, low_contrast_offset;
271
272 50 print2log("INITIAL MAP\n");
273
274 /* Compute total number of blocks in map */
275
2/4
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 50 times.
✓ Branch 5 → 6 taken 50 times.
✗ Branch 5 → 11 not taken.
 50 ASSERT_INT_MUL(mw, mh);
276 50 bsize = mw * mh;
277
278 /* Allocate Direction Map memory */
279 50 direction_map = (int *)g_malloc(bsize * sizeof(int));
280 /* Initialize the Direction Map to INVALID (-1). */
281 50 memset(direction_map, INVALID_DIR, bsize * sizeof(int));
282
283 /* Allocate Low Contrast Map memory */
284 50 low_contrast_map = (int *)g_malloc(bsize * sizeof(int));
285 /* Initialize the Low Contrast Map to FALSE (0). */
286 50 memset(low_contrast_map, 0, bsize * sizeof(int));
287
288 /* Allocate Low Ridge Flow Map memory */
289 50 low_flow_map = (int *)g_malloc(bsize * sizeof(int));
290 /* Initialize the Low Flow Map to FALSE (0). */
291 50 memset(low_flow_map, 0, bsize * sizeof(int));
292
293 /* Allocate DFT directional power vectors */
294
1/2
✗ Branch 10 → 12 not taken.
✓ Branch 10 → 16 taken 50 times.
 50 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 50 nstats = dftwaves->nwaves - 1;
306
1/2
✗ Branch 17 → 18 not taken.
✓ Branch 17 → 23 taken 50 times.
 50 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 50 xminlimit = dftgrids->pad;
319 50 yminlimit = dftgrids->pad;
320 50 xmaxlimit = pw - dftgrids->pad - lfsparms->windowsize - 1;
321 50 ymaxlimit = ph - dftgrids->pad - lfsparms->windowsize - 1;
322
323 /* Foreach block in image ... */
324
2/2
✓ Branch 70 → 24 taken 55834 times.
✓ Branch 70 → 71 taken 50 times.
 55884 for(bi = 0; bi < bsize; bi++){
325 /* Adjust block offset from pointing to block origin to pointing */
326 /* to surrounding window origin. */
327 55834 dft_offset = blkoffs[bi] - (lfsparms->windowoffset * pw) -
328 lfsparms->windowoffset;
329
330 /* Compute pixel coords of window origin. */
331 55834 win_x = dft_offset % pw;
332 55834 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 55834 win_x = max(xminlimit, win_x);
337 55834 win_x = min(xmaxlimit, win_x);
338 55834 win_y = max(yminlimit, win_y);
339 55834 win_y = min(ymaxlimit, win_y);
340 55834 low_contrast_offset = (win_y * pw) + win_x;
341
342 55834 print2log(" BLOCK %2d (%2d, %2d) ", bi, bi%mw, bi/mw);
343
344 /* If block is low contrast ... */
345
2/2
✓ Branch 26 → 27 taken 7286 times.
✓ Branch 26 → 39 taken 48548 times.
 55834 if((ret = low_contrast_block(low_contrast_offset, lfsparms->windowsize,
346 pdata, pw, ph, lfsparms))){
347 /* If system error ... */
348
1/2
✗ Branch 27 → 28 not taken.
✓ Branch 27 → 37 taken 7286 times.
 7286 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 7286 print2log("LOW CONTRAST\n");
362 7286 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 48548 print2log("\n");
368
369 /* Compute DFT powers */
370
1/2
✗ Branch 41 → 42 not taken.
✓ Branch 41 → 51 taken 48548 times.
 48548 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 52 → 53 not taken.
✓ Branch 52 → 62 taken 48548 times.
 48548 if((ret = dft_power_stats(wis, powmaxs, powmax_dirs, pownorms, powers,
388 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 48548 blkdir = primary_dir_test(powers, wis, powmaxs, powmax_dirs,
416 pownorms, nstats, lfsparms);
417
418
2/2
✓ Branch 63 → 64 taken 36693 times.
✓ Branch 63 → 65 taken 11855 times.
 48548 if(blkdir != INVALID_DIR)
419 36693 direction_map[bi] = blkdir;
420 else{
421 /* Conduct secondary (fork) direction test */
422 11855 blkdir = secondary_fork_test(powers, wis, powmaxs, powmax_dirs,
423 pownorms, nstats, lfsparms);
424
2/2
✓ Branch 66 → 67 taken 477 times.
✓ Branch 66 → 68 taken 11378 times.
 11855 if(blkdir != INVALID_DIR)
425 477 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 11378 low_flow_map[bi] = TRUE;
430 }
431
432 } /* End DFT */
433 } /* bi */
434
435 /* Deallocate working memory */
436 50 free_dir_powers(powers, dftwaves->nwaves);
437 50 g_free(wis);
438 50 g_free(powmaxs);
439 50 g_free(powmax_dirs);
440 50 g_free(pownorms);
441
442 50 *odmap = direction_map;
443 50 *olcmap = low_contrast_map;
444 50 *olfmap = low_flow_map;
445 50 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 50 int interpolate_direction_map(int *direction_map, int *low_contrast_map,
474 const int mw, const int mh, const LFSPARMS *lfsparms)
475 {
476 50 int x, y, new_dir;
477 50 int n_dir, e_dir, s_dir, w_dir;
478 50 int n_dist = 0, e_dist = 0, s_dist = 0, w_dist = 0, total_dist;
479 50 int n_found, e_found, s_found, w_found, total_found;
480 50 int n_delta = 0, e_delta = 0, s_delta = 0, w_delta = 0, total_delta;
481 50 int nbr_x, nbr_y;
482 50 int *omap, *dptr, *cptr, *optr;
483 50 double avr_dir;
484
485 50 print2log("INTERPOLATE DIRECTION MAP\n");
486
487 /* Allocate output (interpolated) Direction Map. */
488
1/2
✗ Branch 3 → 4 not taken.
✓ Branch 3 → 5 taken 50 times.
 50 ASSERT_SIZE_MUL(mw, mh);
489
1/2
✓ Branch 5 → 6 taken 50 times.
✗ Branch 5 → 8 not taken.
 50 ASSERT_SIZE_MUL(mw * mh, sizeof(int));
490 50 omap = (int *)g_malloc(mw * mh * sizeof(int));
491
492 /* Set pointers to the first block in the maps. */
493 50 dptr = direction_map;
494 50 cptr = low_contrast_map;
495 50 optr = omap;
496
497 /* Foreach block in the maps ... */
498
2/2
✓ Branch 53 → 54 taken 1812 times.
✓ Branch 53 → 55 taken 50 times.
 1912 for(y = 0; y < mh; y++){
499
2/2
✓ Branch 51 → 9 taken 55834 times.
✓ Branch 51 → 52 taken 1812 times.
 57646 for(x = 0; x < mw; x++){
500
501 /* If image block is NOT LOW CONTRAST and has INVALID direction ... */
502
4/4
✓ Branch 9 → 10 taken 48548 times.
✓ Branch 9 → 49 taken 7286 times.
✓ Branch 10 → 11 taken 14164 times.
✓ Branch 10 → 49 taken 34384 times.
 55834 if((!*cptr) && (*dptr == INVALID_DIR)){
503
504 /* Set neighbor accumulators to 0. */
505 14164 total_found = 0;
506 14164 total_dist = 0;
507
508 /* Find north neighbor. */
509
2/2
✓ Branch 12 → 13 taken 8765 times.
✓ Branch 12 → 14 taken 5399 times.
 14164 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 8765 n_dist = y - nbr_y;
514 /* Accumulate neighbor distance. */
515 8765 total_dist += n_dist;
516 /* Bump number of neighbors found. */
517 8765 total_found++;
518 }
519
520 /* Find east neighbor. */
521
2/2
✓ Branch 15 → 16 taken 9688 times.
✓ Branch 15 → 17 taken 4476 times.
 14164 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 9688 e_dist = nbr_x - x;
526 /* Accumulate neighbor distance. */
527 9688 total_dist += e_dist;
528 /* Bump number of neighbors found. */
529 9688 total_found++;
530 }
531
532 /* Find south neighbor. */
533
2/2
✓ Branch 18 → 19 taken 10511 times.
✓ Branch 18 → 20 taken 3653 times.
 14164 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 10511 s_dist = nbr_y - y;
538 /* Accumulate neighbor distance. */
539 10511 total_dist += s_dist;
540 /* Bump number of neighbors found. */
541 10511 total_found++;
542 }
543
544 /* Find west neighbor. */
545
2/2
✓ Branch 21 → 22 taken 10281 times.
✓ Branch 21 → 23 taken 3883 times.
 14164 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 10281 w_dist = x - nbr_x;
550 /* Accumulate neighbor distance. */
551 10281 total_dist += w_dist;
552 /* Bump number of neighbors found. */
553 10281 total_found++;
554 }
555
556 /* If a sufficient number of neighbors found (Ex. 2) ... */
557
2/2
✓ Branch 23 → 24 taken 12173 times.
✓ Branch 23 → 48 taken 1991 times.
 14164 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 12173 total_delta = 0.0;
562 /* If neighbor found to the north ... */
563
2/2
✓ Branch 24 → 25 taken 8345 times.
✓ Branch 24 → 26 taken 3828 times.
 12173 if(n_found){
564 8345 n_delta = total_dist - n_dist;
565 8345 total_delta += n_delta;
566 }
567 /* If neighbor found to the east ... */
568
2/2
✓ Branch 26 → 27 taken 9260 times.
✓ Branch 26 → 28 taken 2913 times.
 12173 if(e_found){
569 9260 e_delta = total_dist - e_dist;
570 9260 total_delta += e_delta;
571 }
572 /* If neighbor found to the south ... */
573
2/2
✓ Branch 28 → 29 taken 10205 times.
✓ Branch 28 → 30 taken 1968 times.
 12173 if(s_found){
574 10205 s_delta = total_dist - s_dist;
575 10205 total_delta += s_delta;
576 }
577 /* If neighbor found to the west ... */
578
2/2
✓ Branch 30 → 31 taken 9765 times.
✓ Branch 30 → 32 taken 2408 times.
 12173 if(w_found){
579 9765 w_delta = total_dist - w_dist;
580 9765 total_delta += w_delta;
581 }
582
583 12173 avr_dir = 0.0;
584
585
2/2
✓ Branch 32 → 33 taken 8345 times.
✓ Branch 32 → 34 taken 3828 times.
 12173 if(n_found){
586 8345 avr_dir += (n_dir*(n_delta/(double)total_delta));
587 }
588
2/2
✓ Branch 34 → 35 taken 9260 times.
✓ Branch 34 → 36 taken 2913 times.
 12173 if(e_found){
589 9260 avr_dir += (e_dir*(e_delta/(double)total_delta));
590 }
591
2/2
✓ Branch 36 → 37 taken 10205 times.
✓ Branch 36 → 38 taken 1968 times.
 12173 if(s_found){
592 10205 avr_dir += (s_dir*(s_delta/(double)total_delta));
593 }
594
2/2
✓ Branch 38 → 39 taken 9765 times.
✓ Branch 38 → 40 taken 2408 times.
 12173 if(w_found){
595 9765 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 40 → 41 not taken.
✓ Branch 40 → 42 taken 12173 times.
 12173 avr_dir = trunc_dbl_precision(avr_dir, TRUNC_SCALE);
601
602 /* Assign interpolated direction to output Direction Map. */
603
1/2
✗ Branch 43 → 44 not taken.
✓ Branch 43 → 45 taken 12173 times.
 12173 new_dir = sround(avr_dir);
604
605 12173 print2log(" Block %d,%d INTERP numnbs=%d newdir=%d\n",
606 x, y, total_found, new_dir);
607
608 12173 *optr = new_dir;
609 }
610 else{
611 /* Otherwise, the direction remains INVALID. */
612 1991 *optr = *dptr;
613 }
614 }
615 else{
616 /* Otherwise, assign the current direction to the output block. */
617 41670 *optr = *dptr;
618 }
619
620 /* Bump to the next block in the maps ... */
621 55834 dptr++;
622 55834 cptr++;
623 55834 optr++;
624 }
625 }
626
627 /* Copy the interpolated directions into the input map. */
628 50 memcpy(direction_map, omap, mw*mh*sizeof(int));
629 /* Deallocate the working memory. */
630 50 g_free(omap);
631
632 /* Return normally. */
633 50 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 50 int morph_TF_map(int *tfmap, const int mw, const int mh,
651 const LFSPARMS *lfsparms)
652 {
653 50 unsigned char *cimage, *mimage, *cptr;
654 50 int *mptr;
655 50 int i;
656
657
2/4
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 50 times.
✓ Branch 4 → 5 taken 50 times.
✗ Branch 4 → 8 not taken.
 50 ASSERT_INT_MUL(mw, mh);
658
659 /* Convert TRUE/FALSE map into a binary byte image. */
660 50 cimage = (unsigned char *)g_malloc(mw * mh);
661
662 50 mimage = (unsigned char *)g_malloc(mw * mh);
663
664 50 cptr = cimage;
665 50 mptr = tfmap;
666
2/2
✓ Branch 10 → 9 taken 55834 times.
✓ Branch 10 → 11 taken 50 times.
 55934 for(i = 0; i < mw*mh; i++){
667 55834 *cptr++ = *mptr++;
668 }
669
670 50 dilate_charimage_2(cimage, mimage, mw, mh);
671 50 dilate_charimage_2(mimage, cimage, mw, mh);
672 50 erode_charimage_2(cimage, mimage, mw, mh);
673 50 erode_charimage_2(mimage, cimage, mw, mh);
674
675 50 cptr = cimage;
676 50 mptr = tfmap;
677
2/2
✓ Branch 17 → 16 taken 55834 times.
✓ Branch 17 → 18 taken 50 times.
 55934 for(i = 0; i < mw*mh; i++){
678 55834 *mptr++ = *cptr++;
679 }
680
681 50 g_free(cimage);
682 50 g_free(mimage);
683
684 50 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 200 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 200 int *pmap;
710 200 int ret, x, y;
711 200 int *blkoffs, bw, bh, bi;
712 200 int *spptr, *pptr;
713
714
1/2
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 200 times.
 200 ASSERT_SIZE_MUL(iw, ih);
715
1/2
✓ Branch 4 → 5 taken 200 times.
✗ Branch 4 → 8 not taken.
 200 ASSERT_SIZE_MUL(iw * ih, sizeof(int));
716
717 200 pmap = (int *)g_malloc(iw * ih * sizeof(int));
718
719
1/2
✗ Branch 7 → 9 not taken.
✓ Branch 7 → 11 taken 200 times.
 200 if((ret = block_offsets(&blkoffs, &bw, &bh, iw, ih, 0, blocksize))){
720 ✗ g_free(pmap);
721 ✗ return(ret);
722 }
723
724
2/4
✓ Branch 11 → 12 taken 200 times.
✗ Branch 11 → 13 not taken.
✗ Branch 12 → 13 not taken.
✓ Branch 12 → 23 taken 200 times.
 200 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 23 → 17 taken 223336 times.
✓ Branch 23 → 24 taken 200 times.
 223536 for(bi = 0; bi < mw*mh; bi++){
733 223336 spptr = pmap + blkoffs[bi];
734
2/2
✓ Branch 21 → 19 taken 1786688 times.
✓ Branch 21 → 22 taken 223336 times.
 2010024 for(y = 0; y < blocksize; y++){
735 pptr = spptr;
736
2/2
✓ Branch 19 → 18 taken 14293504 times.
✓ Branch 19 → 20 taken 1786688 times.
 16080192 for(x = 0; x < blocksize; x++){
737 14293504 *pptr++ = imap[bi];
738 }
739 1786688 spptr += iw;
740 }
741 }
742
743 /* Deallocate working memory. */
744 200 g_free(blkoffs);
745 /* Assign pixelized map to output pointer. */
746 200 *omap = pmap;
747
748 /* Return normally. */
749 200 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 100 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 100 int mx, my;
771 100 int *dptr, *cptr;
772 100 int avrdir, nvalid;
773 100 double dir_strength;
774
775 100 print2log("SMOOTH DIRECTION MAP\n");
776
777 /* Assign pointers to beginning of both maps. */
778 100 dptr = direction_map;
779 100 cptr = low_contrast_map;
780
781 /* Foreach block in maps ... */
782
2/2
✓ Branch 15 → 16 taken 3624 times.
✓ Branch 15 → 17 taken 100 times.
 3824 for(my = 0; my < mh; my++){
783
2/2
✓ Branch 13 → 4 taken 111668 times.
✓ Branch 13 → 14 taken 3624 times.
 115292 for(mx = 0; mx < mw; mx++){
784 /* If the current block does NOT have LOW CONTRAST ... */
785
2/2
✓ Branch 4 → 5 taken 97096 times.
✓ Branch 4 → 12 taken 14572 times.
 111668 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 97096 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 6 → 7 taken 87102 times.
✓ Branch 6 → 12 taken 9994 times.
 97096 if(dir_strength >= lfsparms->dir_strength_min){
796 /* If Direction Map direction is valid ... */
797
2/2
✓ Branch 7 → 8 taken 77113 times.
✓ Branch 7 → 10 taken 9989 times.
 87102 if(*dptr != INVALID_DIR){
798 /* Conduct valid neighbor test (Ex. thresh==3)... */
799
1/2
✓ Branch 8 → 9 taken 77113 times.
✗ Branch 8 → 12 not taken.
 77113 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 77113 *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 10 → 11 taken 1492 times.
✓ Branch 10 → 12 taken 8497 times.
 9989 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 1492 *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 111668 dptr++;
843 111668 cptr++;
844 }
845 }
846 100 }
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 50 int gen_high_curve_map(int **ohcmap, int *direction_map,
865 const int mw, const int mh, const LFSPARMS *lfsparms)
866 {
867 50 int *high_curve_map, mapsize;
868 50 int *hptr, *dptr;
869 50 int bx, by;
870 50 int nvalid, cmeasure, vmeasure;
871
872
2/4
✗ Branch 2 → 3 not taken.
✓ Branch 2 → 4 taken 50 times.
✓ Branch 4 → 5 taken 50 times.
✗ Branch 4 → 6 not taken.
 50 ASSERT_INT_MUL(mw, mh);
873 50 mapsize = mw*mh;
874
875 /* Allocate High Curvature Map. */
876
1/2
✓ Branch 5 → 7 taken 50 times.
✗ Branch 5 → 9 not taken.
 50 ASSERT_SIZE_MUL(mapsize, sizeof(int));
877 50 high_curve_map = (int *)g_malloc(mapsize * sizeof(int));
878 /* Initialize High Curvature Map to FALSE (0). */
879 50 memset(high_curve_map, 0, mapsize*sizeof(int));
880
881 50 hptr = high_curve_map;
882 50 dptr = direction_map;
883
884 /* Foreach row in maps ... */
885
2/2
✓ Branch 23 → 24 taken 1812 times.
✓ Branch 23 → 25 taken 50 times.
 1862 for(by = 0; by < mh; by++){
886 /* Foreach column in maps ... */
887
2/2
✓ Branch 21 → 10 taken 55834 times.
✓ Branch 21 → 22 taken 1812 times.
 57646 for(bx = 0; bx < mw; bx++){
888
889 /* Count number of valid neighbors around current block ... */
890 55834 nvalid = num_valid_8nbrs(direction_map, bx, by, mw, mh);
891
892 /* If valid neighbors exist ... */
893
2/2
✓ Branch 11 → 12 taken 48110 times.
✓ Branch 11 → 20 taken 7724 times.
 55834 if(nvalid > 0){
894 /* If current block's direction is INVALID ... */
895
2/2
✓ Branch 12 → 13 taken 6892 times.
✓ Branch 12 → 17 taken 41218 times.
 48110 if(*dptr == INVALID_DIR){
896 /* If a sufficient number of VALID neighbors exists ... */
897
2/2
✓ Branch 13 → 14 taken 24 times.
✓ Branch 13 → 20 taken 6868 times.
 6892 if(nvalid >= lfsparms->vort_valid_nbr_min){
898 /* Measure vorticity of neighbors. */
899 24 vmeasure = vorticity(direction_map, bx, by, mw, mh,
900 lfsparms->num_directions);
901 /* If vorticity is sufficiently high ... */
902
2/2
✓ Branch 15 → 16 taken 3 times.
✓ Branch 15 → 20 taken 21 times.
 24 if(vmeasure >= lfsparms->highcurv_vorticity_min)
903 /* Flag block as HIGH CURVATURE. */
904 3 *hptr = TRUE;
905 }
906 }
907 /* Otherwise block has valid direction ... */
908 else{
909 /* Measure curvature around the valid block. */
910 41218 cmeasure = curvature(direction_map, bx, by, mw, mh,
911 lfsparms->num_directions);
912 /* If curvature is sufficiently high ... */
913
2/2
✓ Branch 18 → 19 taken 1887 times.
✓ Branch 18 → 20 taken 39331 times.
 41218 if(cmeasure >= lfsparms->highcurv_curvature_min)
914 1887 *hptr = TRUE;
915 }
916 } /* Else (nvalid <= 0) */
917
918 /* Bump pointers to next block in maps. */
919 55834 dptr++;
920 55834 hptr++;
921
922 } /* bx */
923 } /* by */
924
925 /* Assign High Curvature Map to output pointer. */
926 50 *ohcmap = high_curve_map;
927
928 /* Return normally. */
929 50 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 48548 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 48548 int w;
1009
1010 48548 print2log(" Primary\n");
1011
1012 /* Look at max power statistics in decreasing order ... */
1013
2/2
✓ Branch 8 → 4 taken 76333 times.
✓ Branch 8 → 9 taken 11855 times.
 136736 for(w = 0; w < nstats; w++){
1014 /* 1. Test magnitude of current max power (Ex. Thresh==100000) */
1015
2/2
✓ Branch 4 → 5 taken 68654 times.
✓ Branch 4 → 7 taken 7679 times.
 76333 if((powmaxs[wis[w]] > lfsparms->powmax_min) &&
1016 /* 2. Test magnitude of normalized max power (Ex. Thresh==3.8) */
1017
2/2
✓ Branch 5 → 6 taken 37564 times.
✓ Branch 5 → 7 taken 31090 times.
 68654 (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 6 → 7 taken 871 times.
✓ Branch 6 → 9 taken 36693 times.
 37564 (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 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 11855 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 11855 int ldir, rdir;
1093 11855 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 11855 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 2 → 3 taken 10643 times.
✓ Branch 2 → 11 taken 1212 times.
 11855 if((powmaxs[wis[0]] > lfsparms->powmax_min) &&
1105 /* 2. Test magnitude of corresponding normalized power */
1106 /* (Ex. Thresh==2.85) */
1107
2/2
✓ Branch 3 → 4 taken 6183 times.
✓ Branch 3 → 11 taken 4460 times.
 10643 (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 4 → 5 taken 5787 times.
✓ Branch 4 → 11 taken 396 times.
 6183 (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 5787 rdir = (powmax_dirs[wis[0]] + lfsparms->fork_interval) %
1131 5787 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 5787 ldir = (powmax_dirs[wis[0]] + lfsparms->num_directions -
1138 lfsparms->fork_interval) % lfsparms->num_directions;
1139
1140 5787 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 5787 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 6 → 7 taken 268 times.
✓ Branch 6 → 8 taken 5519 times.
 5787 if(((powers[wis[0]+1][ldir] <= fork_pow_thresh) ||
1157
4/4
✓ Branch 7 → 8 taken 260 times.
✓ Branch 7 → 11 taken 8 times.
✓ Branch 8 → 9 taken 5519 times.
✓ Branch 8 → 10 taken 260 times.
 5787 (powers[wis[0]+1][rdir] <= fork_pow_thresh)) &&
1158 5519 ((powers[wis[0]+1][ldir] > fork_pow_thresh) ||
1159
2/2
✓ Branch 9 → 10 taken 217 times.
✓ Branch 9 → 11 taken 5302 times.
 5519 (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 477 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 100 void remove_incon_dirs(int *imap, const int mw, const int mh,
1225 const DIR2RAD *dir2rad, const LFSPARMS *lfsparms)
1226 {
1227 100 int cx, cy;
1228 100 int *iptr;
1229 100 int nremoved;
1230 100 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 100 cx = mw>>1;
1239 100 cy = mh>>1;
1240
1241 /* Do pass, while directions have been removed in a pass ... */
1242 355 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 355 nremoved = 0;
1252
1253 /* Start at center */
1254 355 iptr = imap + (cy * mw) + cx;
1255 /* If valid IMAP direction and test for removal is true ... */
1256
4/4
✓ Branch 3 → 4 taken 227 times.
✓ Branch 3 → 7 taken 128 times.
✓ Branch 5 → 6 taken 12 times.
✓ Branch 5 → 7 taken 215 times.
 582 if((*iptr != INVALID_DIR)&&
1257 227 (remove_dir(imap, cx, cy, mw, mh, dir2rad, lfsparms))){
1258
1259 /* Set to INVALID */
1260 12 *iptr = INVALID_DIR;
1261 /* Bump number of removed IMAP directions */
1262 12 nremoved++;
1263 }
1264
1265 /* Initialize side indices of concentric boxes */
1266 355 lbox = cx-1;
1267 355 tbox = cy-1;
1268 355 rbox = cx+1;
1269 355 bbox = cy+1;
1270
1271 /* Grow concentric boxes, until ALL edges of imap are exceeded */
1272
4/4
✓ Branch 21 → 8 taken 5676 times.
✓ Branch 21 → 22 taken 1572 times.
✓ Branch 22 → 8 taken 1217 times.
✓ Branch 22 → 23 taken 355 times.
 7248 while((lbox >= 0) || (rbox < mw) || (tbox >= 0) || (bbox < mh)){
1273
1274 /* test top edge of box */
1275
2/2
✓ Branch 8 → 9 taken 6465 times.
✓ Branch 8 → 11 taken 428 times.
 6893 if(tbox >= 0)
1276 6465 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 11 → 12 taken 5353 times.
✓ Branch 11 → 14 taken 1540 times.
 6893 if(rbox < mw)
1281 5353 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 14 → 15 taken 6171 times.
✓ Branch 14 → 17 taken 722 times.
 6893 if(bbox < mh)
1286 6171 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 17 → 18 taken 5676 times.
✓ Branch 17 → 20 taken 1217 times.
 6893 if(lbox >=0)
1291 5676 nremoved += test_left_edge(lbox, tbox, rbox, bbox, imap, mw, mh,
1292 dir2rad, lfsparms);
1293
1294 /* Resize current box */
1295 6893 lbox--;
1296 6893 tbox--;
1297 6893 rbox++;
1298 6893 bbox++;
1299 }
1300
2/2
✓ Branch 23 → 24 taken 255 times.
✓ Branch 23 → 25 taken 100 times.
 355 }while(nremoved);
1301
1302 #ifdef LOG_REPORT
1303 } /* Close scope of numpass */
1304 #endif
1305
1306 100 }
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 6465 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 6465 int bx, by, sx, ex;
1334 6465 int *iptr, *sptr, *eptr;
1335 6465 int nremoved;
1336
1337 /* Initialize number of directions removed on edge to 0 */
1338 6465 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 6465 sx = max(lbox, 0);
1343 6465 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 2 → 3 taken 5184 times.
✓ Branch 2 → 4 taken 1281 times.
 6465 ex = min(rbox-1, mw-1);
1349 6465 eptr = imap + (tbox*mw) + ex;
1350
1351 /* For each point on box's edge ... */
1352 6465 for(iptr = sptr, bx = sx, by = tbox;
1353
2/2
✓ Branch 11 → 6 taken 116819 times.
✓ Branch 11 → 12 taken 6465 times.
 123284 iptr <= eptr;
1354 116819 iptr++, bx++){
1355 /* If valid IMAP direction and test for removal is true ... */
1356
4/4
✓ Branch 6 → 7 taken 70301 times.
✓ Branch 6 → 10 taken 46518 times.
✓ Branch 8 → 9 taken 1445 times.
✓ Branch 8 → 10 taken 68856 times.
 187120 if((*iptr != INVALID_DIR)&&
1357 70301 (remove_dir(imap, bx, by, mw, mh, dir2rad, lfsparms))){
1358 /* Set to INVALID */
1359 1445 *iptr = INVALID_DIR;
1360 /* Bump number of removed IMAP directions */
1361 1445 nremoved++;
1362 }
1363 }
1364
1365 /* Return the number of directions removed on edge */
1366 6465 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 5353 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 5353 int bx, by, sy, ey;
1395 5353 int *iptr, *sptr, *eptr;
1396 5353 int nremoved;
1397
1398 /* Initialize number of directions removed on edge to 0 */
1399 5353 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 5353 sy = max(tbox, 0);
1404 5353 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 2 → 3 taken 4954 times.
✓ Branch 2 → 4 taken 399 times.
 5353 ey = min(bbox-1,mh-1);
1410 5353 eptr = imap + (ey*mw) + rbox;
1411
1412 /* For each point on box's edge ... */
1413 5353 for(iptr = sptr, bx = rbox, by = sy;
1414
2/2
✓ Branch 11 → 6 taken 87473 times.
✓ Branch 11 → 12 taken 5353 times.
 92826 iptr <= eptr;
1415 87473 iptr+=mw, by++){
1416 /* If valid IMAP direction and test for removal is true ... */
1417
4/4
✓ Branch 6 → 7 taken 57064 times.
✓ Branch 6 → 10 taken 30409 times.
✓ Branch 8 → 9 taken 996 times.
✓ Branch 8 → 10 taken 56068 times.
 144537 if((*iptr != INVALID_DIR)&&
1418 57064 (remove_dir(imap, bx, by, mw, mh, dir2rad, lfsparms))){
1419 /* Set to INVALID */
1420 996 *iptr = INVALID_DIR;
1421 /* Bump number of removed IMAP directions */
1422 996 nremoved++;
1423 }
1424 }
1425
1426 /* Return the number of directions removed on edge */
1427 5353 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 6171 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 6171 int bx, by, sx, ex;
1455 6171 int *iptr, *sptr, *eptr;
1456 6171 int nremoved;
1457
1458 /* Initialize number of directions removed on edge to 0 */
1459 6171 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 6171 sx = min(rbox, mw-1);
1464 6171 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 6171 ex = max(lbox-1, 0);
1470 6171 eptr = imap + (bbox*mw) + ex;
1471
1472 /* For each point on box's edge ... */
1473 6171 for(iptr = sptr, bx = sx, by = bbox;
1474
2/2
✓ Branch 8 → 3 taken 117998 times.
✓ Branch 8 → 9 taken 6171 times.
 124169 iptr >= eptr;
1475 117998 iptr--, bx--){
1476 /* If valid IMAP direction and test for removal is true ... */
1477
4/4
✓ Branch 3 → 4 taken 79720 times.
✓ Branch 3 → 7 taken 38278 times.
✓ Branch 5 → 6 taken 975 times.
✓ Branch 5 → 7 taken 78745 times.
 197718 if((*iptr != INVALID_DIR)&&
1478 79720 (remove_dir(imap, bx, by, mw, mh, dir2rad, lfsparms))){
1479 /* Set to INVALID */
1480 975 *iptr = INVALID_DIR;
1481 /* Bump number of removed IMAP directions */
1482 975 nremoved++;
1483 }
1484 }
1485
1486 /* Return the number of directions removed on edge */
1487 6171 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 5676 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 5676 int bx, by, sy, ey;
1516 5676 int *iptr, *sptr, *eptr;
1517 5676 int nremoved;
1518
1519 /* Initialize number of directions removed on edge to 0 */
1520 5676 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 5676 sy = min(bbox, mh-1);
1526 5676 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 5676 ey = max(tbox-1, 0);
1532 5676 eptr = imap + (ey*mw) + lbox;
1533
1534 /* For each point on box's edge ... */
1535 5676 for(iptr = sptr, bx = lbox, by = sy;
1536
2/2
✓ Branch 8 → 3 taken 107247 times.
✓ Branch 8 → 9 taken 5676 times.
 112923 iptr >= eptr;
1537 107247 iptr-=mw, by--){
1538 /* If valid IMAP direction and test for removal is true ... */
1539
4/4
✓ Branch 3 → 4 taken 75402 times.
✓ Branch 3 → 7 taken 31845 times.
✓ Branch 5 → 6 taken 1071 times.
✓ Branch 5 → 7 taken 74331 times.
 182649 if((*iptr != INVALID_DIR)&&
1540 75402 (remove_dir(imap, bx, by, mw, mh, dir2rad, lfsparms))){
1541 /* Set to INVALID */
1542 1071 *iptr = INVALID_DIR;
1543 /* Bump number of removed IMAP directions */
1544 1071 nremoved++;
1545 }
1546 }
1547
1548 /* Return the number of directions removed on edge */
1549 5676 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 282714 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 282714 int avrdir, nvalid, dist;
1574 282714 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 282714 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 3 → 4 taken 280197 times.
✓ Branch 3 → 7 taken 2517 times.
 282714 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 4 → 5 taken 269137 times.
✓ Branch 4 → 7 taken 11060 times.
 280197 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 269137 dist = abs(avrdir - *(imap+(my*mw)+mx));
1605 269137 dist = min(dist, dir2rad->ndirs-dist);
1606
2/2
✓ Branch 5 → 6 taken 1982 times.
✓ Branch 5 → 7 taken 267155 times.
 269137 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 1982 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 379810 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 379810 int *iptr;
1659 379810 int e,w,n,s;
1660 379810 double cospart, sinpart;
1661 379810 double pi2, pi_factor, theta;
1662 379810 double avr;
1663
1664 /* Compute neighbor coordinates to current IMAP direction */
1665 379810 e = mx+1; /* East */
1666 379810 w = mx-1; /* West */
1667 379810 n = my-1; /* North */
1668 379810 s = my+1; /* South */
1669
1670 /* Intialize accumulators */
1671 379810 *nvalid = 0;
1672 379810 cospart = 0.0;
1673 379810 sinpart = 0.0;
1674
1675 /* 1. Test NW */
1676 /* If NW point within IMAP boudaries ... */
1677
2/2
✓ Branch 2 → 3 taken 358743 times.
✓ Branch 2 → 5 taken 21067 times.
 379810 if((w >= 0) && (n >= 0)){
1678 358743 iptr = imap + (n*mw) + w;
1679 /* If valid direction ... */
1680
2/2
✓ Branch 3 → 4 taken 314122 times.
✓ Branch 3 → 5 taken 44621 times.
 358743 if(*iptr != INVALID_DIR){
1681 /* Accumulate cosine and sine components of the direction */
1682 314122 cospart += dir2rad->cos[*iptr];
1683 314122 sinpart += dir2rad->sin[*iptr];
1684 /* Bump number of accumulated directions */
1685 314122 (*nvalid)++;
1686 }
1687 }
1688
1689 /* 2. Test N */
1690 /* If N point within IMAP boudaries ... */
1691
2/2
✓ Branch 5 → 6 taken 370031 times.
✓ Branch 5 → 8 taken 9779 times.
 379810 if(n >= 0){
1692 370031 iptr = imap + (n*mw) + mx;
1693 /* If valid direction ... */
1694
2/2
✓ Branch 6 → 7 taken 329085 times.
✓ Branch 6 → 8 taken 40946 times.
 370031 if(*iptr != INVALID_DIR){
1695 /* Accumulate cosine and sine components of the direction */
1696 329085 cospart += dir2rad->cos[*iptr];
1697 329085 sinpart += dir2rad->sin[*iptr];
1698 /* Bump number of accumulated directions */
1699 329085 (*nvalid)++;
1700 }
1701 }
1702
1703 /* 3. Test NE */
1704 /* If NE point within IMAP boudaries ... */
1705
2/2
✓ Branch 8 → 9 taken 359669 times.
✓ Branch 8 → 11 taken 20141 times.
 379810 if((e < mw) && (n >= 0)){
1706 359669 iptr = imap + (n*mw) + e;
1707 /* If valid direction ... */
1708
2/2
✓ Branch 9 → 10 taken 314166 times.
✓ Branch 9 → 11 taken 45503 times.
 359669 if(*iptr != INVALID_DIR){
1709 /* Accumulate cosine and sine components of the direction */
1710 314166 cospart += dir2rad->cos[*iptr];
1711 314166 sinpart += dir2rad->sin[*iptr];
1712 /* Bump number of accumulated directions */
1713 314166 (*nvalid)++;
1714 }
1715 }
1716
1717 /* 4. Test E */
1718 /* If E point within IMAP boudaries ... */
1719
2/2
✓ Branch 11 → 12 taken 369210 times.
✓ Branch 11 → 14 taken 10600 times.
 379810 if(e < mw){
1720 369210 iptr = imap + (my*mw) + e;
1721 /* If valid direction ... */
1722
2/2
✓ Branch 12 → 13 taken 326276 times.
✓ Branch 12 → 14 taken 42934 times.
 369210 if(*iptr != INVALID_DIR){
1723 /* Accumulate cosine and sine components of the direction */
1724 326276 cospart += dir2rad->cos[*iptr];
1725 326276 sinpart += dir2rad->sin[*iptr];
1726 /* Bump number of accumulated directions */
1727 326276 (*nvalid)++;
1728 }
1729 }
1730
1731 /* 5. Test SE */
1732 /* If SE point within IMAP boudaries ... */
1733
2/2
✓ Branch 14 → 15 taken 358123 times.
✓ Branch 14 → 17 taken 21687 times.
 379810 if((e < mw) && (s < mh)){
1734 358123 iptr = imap + (s*mw) + e;
1735 /* If valid direction ... */
1736
2/2
✓ Branch 15 → 16 taken 312447 times.
✓ Branch 15 → 17 taken 45676 times.
 358123 if(*iptr != INVALID_DIR){
1737 /* Accumulate cosine and sine components of the direction */
1738 312447 cospart += dir2rad->cos[*iptr];
1739 312447 sinpart += dir2rad->sin[*iptr];
1740 /* Bump number of accumulated directions */
1741 312447 (*nvalid)++;
1742 }
1743 }
1744
1745 /* 6. Test S */
1746 /* If S point within IMAP boudaries ... */
1747
2/2
✓ Branch 17 → 18 taken 368412 times.
✓ Branch 17 → 20 taken 11398 times.
 379810 if(s < mh){
1748 368412 iptr = imap + (s*mw) + mx;
1749 /* If valid direction ... */
1750
2/2
✓ Branch 18 → 19 taken 327368 times.
✓ Branch 18 → 20 taken 41044 times.
 368412 if(*iptr != INVALID_DIR){
1751 /* Accumulate cosine and sine components of the direction */
1752 327368 cospart += dir2rad->cos[*iptr];
1753 327368 sinpart += dir2rad->sin[*iptr];
1754 /* Bump number of accumulated directions */
1755 327368 (*nvalid)++;
1756 }
1757 }
1758
1759 /* 7. Test SW */
1760 /* If SW point within IMAP boudaries ... */
1761
2/2
✓ Branch 20 → 21 taken 357090 times.
✓ Branch 20 → 23 taken 22720 times.
 379810 if((w >= 0) && (s < mh)){
1762 357090 iptr = imap + (s*mw) + w;
1763 /* If valid direction ... */
1764
2/2
✓ Branch 21 → 22 taken 312679 times.
✓ Branch 21 → 23 taken 44411 times.
 357090 if(*iptr != INVALID_DIR){
1765 /* Accumulate cosine and sine components of the direction */
1766 312679 cospart += dir2rad->cos[*iptr];
1767 312679 sinpart += dir2rad->sin[*iptr];
1768 /* Bump number of accumulated directions */
1769 312679 (*nvalid)++;
1770 }
1771 }
1772
1773 /* 8. Test W */
1774 /* If W point within IMAP boudaries ... */
1775
2/2
✓ Branch 23 → 24 taken 368171 times.
✓ Branch 23 → 26 taken 11639 times.
 379810 if(w >= 0){
1776 368171 iptr = imap + (my*mw) + w;
1777 /* If valid direction ... */
1778
2/2
✓ Branch 24 → 25 taken 327698 times.
✓ Branch 24 → 26 taken 40473 times.
 368171 if(*iptr != INVALID_DIR){
1779 /* Accumulate cosine and sine components of the direction */
1780 327698 cospart += dir2rad->cos[*iptr];
1781 327698 sinpart += dir2rad->sin[*iptr];
1782 /* Bump number of accumulated directions */
1783 327698 (*nvalid)++;
1784 }
1785 }
1786
1787 /* If there were no neighbors found with valid direction ... */
1788
2/2
✓ Branch 26 → 27 taken 7854 times.
✓ Branch 26 → 28 taken 371956 times.
 379810 if(*nvalid == 0){
1789 /* Return INVALID direction. */
1790 7854 *dir_strength = 0;
1791 7854 *avrdir = INVALID_DIR;
1792 7854 return;
1793 }
1794
1795 /* Compute averages of accumulated cosine and sine direction components */
1796 371956 cospart /= (double)(*nvalid);
1797 371956 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 371956 *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 371956 *dir_strength = trunc_dbl_precision(*dir_strength, TRUNC_SCALE);
1806
1807 /* If the direction strength is not sufficiently high ... */
1808
2/2
✓ Branch 28 → 29 taken 13699 times.
✓ Branch 28 → 30 taken 358257 times.
 371956 if(*dir_strength < DIR_STRENGTH_MIN){
1809 /* Return INVALID direction. */
1810 13699 *dir_strength = 0;
1811 13699 *avrdir = INVALID_DIR;
1812 13699 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 358257 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 358257 pi2 = 2*M_PI;
1824 358257 theta += pi2;
1825 358257 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 358257 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 358257 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 30 → 31 not taken.
✓ Branch 30 → 32 taken 358257 times.
 358257 avr = trunc_dbl_precision(avr, TRUNC_SCALE);
1839
1/2
✗ Branch 33 → 34 not taken.
✓ Branch 33 → 35 taken 358257 times.
 358257 *avrdir = sround(avr);
1840
1841 /* Really do need to map values > NDIRS back onto [0..NDIRS) range. */
1842 358257 *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 57307 int num_valid_8nbrs(int *imap, const int mx, const int my,
1860 const int mw, const int mh)
1861 {
1862 57307 int e_ind, w_ind, n_ind, s_ind;
1863 57307 int nvalid;
1864
1865 /* Initialize VALID IMAP counter to zero. */
1866 57307 nvalid = 0;
1867
1868 /* Compute neighbor coordinates to current IMAP direction */
1869 57307 e_ind = mx+1; /* East index */
1870 57307 w_ind = mx-1; /* West index */
1871 57307 n_ind = my-1; /* North index */
1872 57307 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 2 → 3 taken 4042 times.
✓ Branch 2 → 4 taken 53265 times.
✓ Branch 4 → 3 taken 11457 times.
✓ Branch 4 → 5 taken 41808 times.
 57307 if((w_ind >= 0) && (n_ind >= 0) && (*(imap + (n_ind*mw) + w_ind) >= 0))
1877 /* Bump VALID counter. */
1878 57307 nvalid++;
1879
1880 /* 2. Test N IMAP value. */
1881
4/4
✓ Branch 5 → 6 taken 55255 times.
✓ Branch 5 → 8 taken 2052 times.
✓ Branch 6 → 7 taken 41658 times.
✓ Branch 6 → 8 taken 13597 times.
 57307 if((n_ind >= 0) && (*(imap + (n_ind*mw) + mx) >= 0))
1882 41658 nvalid++;
1883
1884 /* 3. Test NE IMAP value. */
1885
4/4
✓ Branch 8 → 9 taken 53432 times.
✓ Branch 8 → 11 taken 3875 times.
✓ Branch 9 → 10 taken 41942 times.
✓ Branch 9 → 11 taken 11490 times.
 57307 if((n_ind >= 0) && (e_ind < mw) && (*(imap + (n_ind*mw) + e_ind) >= 0))
1886 41942 nvalid++;
1887
1888 /* 4. Test E IMAP value. */
1889
4/4
✓ Branch 11 → 12 taken 55434 times.
✓ Branch 11 → 14 taken 1873 times.
✓ Branch 12 → 13 taken 41664 times.
✓ Branch 12 → 14 taken 13770 times.
 57307 if((e_ind < mw) && (*(imap + (my*mw) + e_ind) >= 0))
1890 41664 nvalid++;
1891
1892 /* 5. Test SE IMAP value. */
1893
4/4
✓ Branch 14 → 15 taken 53685 times.
✓ Branch 14 → 17 taken 3622 times.
✓ Branch 15 → 16 taken 42219 times.
✓ Branch 15 → 17 taken 11466 times.
 57307 if((e_ind < mw) && (s_ind < mh) && (*(imap + (s_ind*mw) + e_ind) >= 0))
1894 42219 nvalid++;
1895
1896 /* 6. Test S IMAP value. */
1897
4/4
✓ Branch 17 → 18 taken 55508 times.
✓ Branch 17 → 20 taken 1799 times.
✓ Branch 18 → 19 taken 41939 times.
✓ Branch 18 → 20 taken 13569 times.
 57307 if((s_ind < mh) && (*(imap + (s_ind*mw) + mx) >= 0))
1898 41939 nvalid++;
1899
1900 /* 7. Test SW IMAP value. */
1901
4/4
✓ Branch 20 → 21 taken 53518 times.
✓ Branch 20 → 23 taken 3789 times.
✓ Branch 21 → 22 taken 42100 times.
✓ Branch 21 → 23 taken 11418 times.
 57307 if((w_ind >= 0) && (s_ind < mh) && (*(imap + (s_ind*mw) + w_ind) >= 0))
1902 42100 nvalid++;
1903
1904 /* 8. Test W IMAP value. */
1905
4/4
✓ Branch 23 → 24 taken 55267 times.
✓ Branch 23 → 26 taken 2040 times.
✓ Branch 24 → 25 taken 41549 times.
✓ Branch 24 → 26 taken 13718 times.
 57307 if((w_ind >= 0) && (*(imap + (my*mw) + w_ind) >= 0))
1906 41549 nvalid++;
1907
1908 /* Return number of neighbors with VALID IMAP values. */
1909 57307 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 24 int vorticity(int *imap, const int mx, const int my,
1966 const int mw, const int mh, const int ndirs)
1967 {
1968 24 int e_ind, w_ind, n_ind, s_ind;
1969 24 int nw_val, n_val, ne_val, e_val, se_val, s_val, sw_val, w_val;
1970 24 int vmeasure;
1971
1972 /* Compute neighbor coordinates to current IMAP direction */
1973 24 e_ind = mx+1; /* East index */
1974 24 w_ind = mx-1; /* West index */
1975 24 n_ind = my-1; /* North index */
1976 24 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 2 → 3 taken 24 times.
✗ Branch 2 → 4 not taken.
 24 if((w_ind >= 0) && (n_ind >= 0))
1981 /* Set neighbor value to IMAP value. */
1982 24 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 4 → 5 taken 24 times.
✗ Branch 4 → 6 not taken.
 24 if(n_ind >= 0)
1989 24 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 6 → 7 taken 24 times.
✗ Branch 6 → 8 not taken.
 24 if((n_ind >= 0) && (e_ind < mw))
1995 24 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 8 → 9 taken 24 times.
✗ Branch 8 → 10 not taken.
 24 if(e_ind < mw)
2001 24 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 10 → 11 taken 24 times.
✗ Branch 10 → 12 not taken.
 24 if((e_ind < mw) && (s_ind < mh))
2007 24 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 12 → 13 taken 24 times.
✗ Branch 12 → 14 not taken.
 24 if(s_ind < mh)
2013 24 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 14 → 15 taken 24 times.
✗ Branch 14 → 16 not taken.
 24 if((w_ind >= 0) && (s_ind < mh))
2019 24 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 16 → 17 taken 24 times.
✗ Branch 16 → 18 not taken.
 24 if(w_ind >= 0)
2025 24 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 24 vmeasure = 0;
2034
2035 /* 1. NW & N */
2036 24 accum_nbr_vorticity(&vmeasure, nw_val, n_val, ndirs);
2037
2038 /* 2. N & NE */
2039 24 accum_nbr_vorticity(&vmeasure, n_val, ne_val, ndirs);
2040
2041 /* 3. NE & E */
2042 24 accum_nbr_vorticity(&vmeasure, ne_val, e_val, ndirs);
2043
2044 /* 4. E & SE */
2045 24 accum_nbr_vorticity(&vmeasure, e_val, se_val, ndirs);
2046
2047 /* 5. SE & S */
2048 24 accum_nbr_vorticity(&vmeasure, se_val, s_val, ndirs);
2049
2050 /* 6. S & SW */
2051 24 accum_nbr_vorticity(&vmeasure, s_val, sw_val, ndirs);
2052
2053 /* 7. SW & W */
2054 24 accum_nbr_vorticity(&vmeasure, sw_val, w_val, ndirs);
2055
2056 /* 8. W & NW */
2057 24 accum_nbr_vorticity(&vmeasure, w_val, nw_val, ndirs);
2058
2059 /* Return the accumulated vorticity measure. */
2060 24 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 192 void accum_nbr_vorticity(int *vmeasure, const int dir1, const int dir2,
2076 const int ndirs)
2077 {
2078 192 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 2 → 3 taken 110 times.
✓ Branch 2 → 9 taken 82 times.
✓ Branch 3 → 4 taken 94 times.
✓ Branch 3 → 9 taken 16 times.
 192 if((dir1 != dir2) && (dir1 >= 0)&&(dir2 >= 0)){
2084 /* Measure the clockwise distance from the first to the second */
2085 /* directions. */
2086 94 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 4 → 5 taken 42 times.
✓ Branch 4 → 6 taken 52 times.
 94 if(dist < 0)
2095 42 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 6 → 7 taken 39 times.
✓ Branch 6 → 8 taken 55 times.
 94 if(dist > (ndirs>>1))
2099 /* Decrement the vorticity measure. */
2100 39 (*vmeasure)--;
2101 else
2102 /* Otherwise, bump the vorticity measure. */
2103 55 (*vmeasure)++;
2104 }
2105 /* Otherwise both directions are either equal or */
2106 /* one or both directions are INVALID, so ignore. */
2107 192 }
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 41218 int curvature(int *imap, const int mx, const int my,
2126 const int mw, const int mh, const int ndirs)
2127 {
2128 41218 int *iptr;
2129 41218 int e_ind, w_ind, n_ind, s_ind;
2130 41218 int nw_val, n_val, ne_val, e_val, se_val, s_val, sw_val, w_val;
2131 41218 int cmeasure, dist;
2132
2133 /* Compute neighbor coordinates to current IMAP direction */
2134 41218 e_ind = mx+1; /* East index */
2135 41218 w_ind = mx-1; /* West index */
2136 41218 n_ind = my-1; /* North index */
2137 41218 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 2 → 3 taken 41218 times.
✗ Branch 2 → 4 not taken.
 41218 if((w_ind >= 0) && (n_ind >= 0))
2142 /* Set neighbor value to IMAP value. */
2143 41218 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 4 → 5 taken 41218 times.
✗ Branch 4 → 6 not taken.
 41218 if(n_ind >= 0)
2150 41218 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 6 → 7 taken 41218 times.
✗ Branch 6 → 8 not taken.
 41218 if((n_ind >= 0) && (e_ind < mw))
2156 41218 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 8 → 9 taken 41218 times.
✗ Branch 8 → 10 not taken.
 41218 if(e_ind < mw)
2162 41218 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 10 → 11 taken 41218 times.
✗ Branch 10 → 12 not taken.
 41218 if((e_ind < mw) && (s_ind < mh))
2168 41218 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 12 → 13 taken 41218 times.
✗ Branch 12 → 14 not taken.
 41218 if(s_ind < mh)
2174 41218 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 14 → 15 taken 41218 times.
✗ Branch 14 → 16 not taken.
 41218 if((w_ind >= 0) && (s_ind < mh))
2180 41218 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 16 → 17 taken 41218 times.
✗ Branch 16 → 18 not taken.
 41218 if(w_ind >= 0)
2186 41218 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 41218 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 41218 cmeasure = -1;
2199
2200 /* 1. With NW */
2201 /* Compute closest distance between neighboring directions. */
2202 41218 dist = closest_dir_dist(*iptr, nw_val, ndirs);
2203 /* Keep track of maximum. */
2204 41218 if(dist > cmeasure)
2205 cmeasure = dist;
2206
2207 /* 2. With N */
2208 41218 dist = closest_dir_dist(*iptr, n_val, ndirs);
2209 41218 if(dist > cmeasure)
2210 cmeasure = dist;
2211
2212 /* 3. With NE */
2213 41218 dist = closest_dir_dist(*iptr, ne_val, ndirs);
2214 41218 if(dist > cmeasure)
2215 cmeasure = dist;
2216
2217 /* 4. With E */
2218 41218 dist = closest_dir_dist(*iptr, e_val, ndirs);
2219 41218 if(dist > cmeasure)
2220 cmeasure = dist;
2221
2222 /* 5. With SE */
2223 41218 dist = closest_dir_dist(*iptr, se_val, ndirs);
2224 41218 if(dist > cmeasure)
2225 cmeasure = dist;
2226
2227 /* 6. With S */
2228 41218 dist = closest_dir_dist(*iptr, s_val, ndirs);
2229 41218 if(dist > cmeasure)
2230 cmeasure = dist;
2231
2232 /* 7. With SW */
2233 41218 dist = closest_dir_dist(*iptr, sw_val, ndirs);
2234 41218 if(dist > cmeasure)
2235 cmeasure = dist;
2236
2237 /* 8. With W */
2238 41218 dist = closest_dir_dist(*iptr, w_val, ndirs);
2239 41218 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 41218 return(cmeasure);
2245 }
2246