1 /* functions that support interpolation in latitude-longitude projections */
2
3 #include <stdlib.h>
4 #include <grass/gis.h>
5 #include "main.h"
6
7 /************************************************************************/
8 /* This function initializes the search for nearest neighbors */
9 /* by locating the two data closest to the specified column in */
10 /* a linked list of row data */
11
first_west_LL(EW * ewptr,SHORT col)12 int first_west_LL(EW * ewptr, SHORT col)
13 {
14 if (ewptr->start == NULL) /* no data in this row */
15 ewptr->walive = ewptr->ealive = FALSE;
16 else if (ewptr->start == ewptr->start->prior) { /* single datum */
17 ewptr->west = ewptr->east = ewptr->start;
18 ewptr->walive = FALSE;
19 ewptr->ealive = TRUE;
20 }
21 else { /* two or more data in this row */
22 while (col > ewptr->start->x && ewptr->start->x < ewptr->start->next->x) /* start is west of col */
23 ewptr->start = ewptr->start->next;
24 ewptr->east = ewptr->start;
25 ewptr->west = ewptr->start->prior;
26 ewptr->walive = ewptr->ealive = TRUE;
27 }
28
29 return 0;
30 }
31
32
offset_distance_LL(SHORT offset)33 double offset_distance_LL(SHORT offset)
34 {
35 extern double *lat_diff;
36
37 return (*(lat_diff + abs((int)offset)));
38 }
39
40
41 /************************************************************************/
42 /* Return TRUE if search is exhausted both west and east in a row */
43
completed_row_LL(EW * ewptr)44 int completed_row_LL(EW * ewptr)
45 {
46 return (!ewptr->walive && !ewptr->ealive);
47 }
48
49
find_neighbors_LL(EW * ewptr,NEIGHBOR * nbr_head,SHORT row,SHORT col,int npoints,SHORT * neighbors)50 int find_neighbors_LL(EW * ewptr, NEIGHBOR * nbr_head, SHORT row, SHORT col,
51 int npoints, SHORT * neighbors)
52 {
53 MELEMENT **Mptr; /* double indirection !! */
54 int westward = 1; /* 1 if west of interpolation point */
55 double distance;
56 short *active; /* TRUE if active search in this direction */
57
58 active = &ewptr->walive; /* TRUE if searching west in this row */
59 Mptr = &ewptr->west; /* process search west first, then east */
60 do {
61 if (*active) {
62 distance = distance_LL(row, col, *Mptr);
63
64 if (*neighbors < npoints)
65 add_neighbor(Mptr, nbr_head, distance, ++(*neighbors));
66 else if (!replace_neighbor(Mptr, nbr_head, distance))
67 *active = FALSE; /* curtail search in this direction */
68
69 if (*active)
70 if (westward)
71 extend_west(ewptr);
72 else
73 extend_east(ewptr);
74 }
75
76 active = &ewptr->ealive;
77 Mptr = &ewptr->east;
78 } while (westward--); /* repeat loop for east and quit */
79
80 return 0;
81 }
82
83
84 /************************************************************************/
85 /* This function exhausts all possible nearest neighhbors */
86 /* within the row indexed by the ew search pointer */
87
exhaust_search_LL(EW * ewptr,NEIGHBOR * nbr_head,SHORT row,SHORT col)88 int exhaust_search_LL(EW * ewptr, NEIGHBOR * nbr_head, SHORT row, SHORT col)
89 {
90 double distance;
91
92 while (ewptr->walive) { /* search active */
93 distance = distance_LL(row, col, ewptr->west);
94
95 if (!replace_neighbor(&ewptr->west, nbr_head, distance))
96 break; /* curtail search in this direction */
97 else
98 extend_west(ewptr);
99 }
100
101 while (ewptr->ealive) { /* search active */
102 distance = distance_LL(row, col, ewptr->east);
103
104 if (!replace_neighbor(&ewptr->east, nbr_head, distance))
105 break; /* curtail search in this direction */
106 else
107 extend_east(ewptr);
108 }
109
110 return 0;
111 }
112
113
extend_west(EW * ewptr)114 int extend_west(EW * ewptr)
115 {
116 if (ewptr->west->prior != ewptr->east)
117 ewptr->west = ewptr->west->prior;
118 else
119 ewptr->walive = FALSE;
120
121 return 0;
122 }
123
124
extend_east(EW * ewptr)125 int extend_east(EW * ewptr)
126 {
127 if (ewptr->east->next != ewptr->west)
128 ewptr->east = ewptr->east->next;
129 else
130 ewptr->ealive = FALSE;
131
132 return 0;
133 }
134
135
distance_LL(SHORT row,SHORT col,MELEMENT * Mptr)136 double distance_LL(SHORT row, SHORT col, MELEMENT * Mptr)
137 {
138 extern double *rowlook, *collook;
139
140 /* use lookup tables to increase distance calculation efficiency */
141 LL_set_geodesic_distance(rowlook, row, Mptr->y);
142 return (LL_geodesic_distance(*(collook + abs((int)(col - Mptr->x)))));
143 }
144
145
146 /************************************************************************/
147 /* Lookup tables storing pre-processed latitude and longitude data */
148 /* are created for later use in selecting nearest neighbors */
149
LL_lookup_tables(SHORT nrows,SHORT ncols)150 int LL_lookup_tables(SHORT nrows, SHORT ncols)
151 {
152 extern double *rowlook, *collook, *lat_diff;
153 extern struct Cell_head window;
154 double *nextrow, *nextcol, *next_diff,
155 lon = 0., lat = window.north - (0.5 * window.ns_res);
156 SHORT i;
157
158 nextrow = rowlook = (double *)G_calloc(nrows, sizeof(double));
159 for (i = 0; i < nrows; i++, nextrow++, lat -= window.ns_res)
160 *nextrow = LL_set_geodesic_distance_lat(lat);
161
162 nextcol = collook = (double *)G_calloc(ncols, sizeof(double));
163 for (i = 0; i < ncols; i++, nextcol++, lon += window.ew_res)
164 *nextcol = set_sdlmr(lon);
165
166 /* compute distance between latitudes at same longitude */
167 next_diff = lat_diff = (double *)G_calloc(nrows, sizeof(double));
168 for (i = 0; i < nrows; i++, next_diff++) {
169 LL_set_geodesic_distance(rowlook, 0, i);
170 *next_diff = LL_geodesic_distance(0.);
171 }
172
173 return 0;
174 }
175