-
Notifications
You must be signed in to change notification settings - Fork 0
/
rolling_coin_smoothing.c
executable file
·152 lines (128 loc) · 5.79 KB
/
rolling_coin_smoothing.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
#include "bathymetrictools.h"
/*
* This file contains:
* - Navigationally safe Rolling Coin surface manipulation (2.5D) functions
* - Coin is defined by structured datatype "Coin"
*/
/*
* Surface manipulation (~smoothing) function
* - Iterates over cells
* - Modifies the surface
* - Memory management and no data handling
*/
void coinRollSurface(struct FloatSurface *src, struct Coin *penny) {
printf("Rolling Coin..");
fflush(stdout);
const int radius = penny->radius; // Valid indexes of coin are normally [-radius, radius]
int *limits = calloc(4, sizeof(int)); // An array to hold valid coin index ranges for special cases
const float nodata = src->nodata;
const float placeholder = -999999.0;
float shoalest;
// Create new temporary float** (2D) array:
float **temp = createFloatArray(src->cols, src->rows);
// Initialize all cells to an elevation of 10 000 (meters):
for (int row = 0; row < src->rows; row ) {
for (int col = 0; col < src->cols; col ) {
temp[row][col] = 10000.0;
}
}
// Iterate over depth model cells and smooth surface:
for (int row = 0; row < src->rows; row ) {
for (int col = 0; col < src->cols; col ) {
// Get coin index limits:
getCoinIndexRange(src->rows, src->cols, row, col, penny->radius, limits);
shoalest = placeholder; // Reset shoalest value
shoalest = getShoalestDepthOnCoin(src, penny, row, col); // Get actual shoalest value
if (fabs(shoalest - nodata) > EPSILON) { // Depth values found (!= nodata)
for (int row_coin = limits[0]; row_coin <= limits[2]; row_coin ) { // Use valid Coin Row indexes (coin may be partial)
for (int col_coin = limits[1]; col_coin <= limits[3]; col_coin ) { // Use valid Coin Col indexes (coin may be partial)
if (penny->array[row_coin radius][col_coin radius] == TRUE) { // On the coin
if (temp[row row_coin][col col_coin] > shoalest) { // If current depth of cell is shoaler than shoalest
temp[row row_coin][col col_coin] = shoalest; // "Press" depth to coin area
}
}
}
}
}
}
}
// Restore original nodata (safety first):
for (int row = 0; row < src->rows; row ) {
for (int col = 0; col < src->cols; col ) {
if ((fabs(src->array[row][col] - nodata) < EPSILON)) { // (value == nodata)
temp[row][col] = nodata;
}
}
}
// Free memory allocated for temp array:
float **destruct = src->array; // Store original data array pointer
src->array = temp; // Replace original data array with smoothed data array
freeFloatArray(destruct, src->rows); // Free original data array
free(limits); // Free index range list memory
printf("Done\n");
fflush(stdout);
}
/*
* Checks coin area, finds shoalest depth (maximum elevation)
* - Edge effects: uses partial coin when necessary
* - Returns shoalest depth on neighborhood (coin area) or 'No data'
*/
float getShoalestDepthOnCoin(struct FloatSurface *src, struct Coin *penny, const int row_index, const int col_index) {
const float placeholder = -999999.0;
float ret = placeholder;
int *limits;
// Get coin index range:
limits = calloc(4, sizeof(int));
getCoinIndexRange(src->rows, src->cols, row_index, col_index, penny->radius, limits);
// Check coin area using coin index ranges based on current cell:
for (int row_coin = limits[0]; row_coin <= limits[2]; row_coin ) {
for (int col_coin = limits[1]; col_coin <= limits[3]; col_coin ) {
if (penny->array[penny->radius row_coin][penny->radius col_coin] == TRUE) { // If cell is on coin --> check
if (src->array[row_index row_coin][col_index col_coin] > ret) {
if (fabs(src->array[row_index row_coin][col_index col_coin] - src->nodata) > EPSILON) { // != NO DATA
ret = src->array[row_index row_coin][col_index col_coin];
}
}
}
}
}
free(limits); // Free index range list memory
if (fabs(ret - src->nodata) < EPSILON) {
ret = src->nodata;
}
return ret;
}
/*
* Calculates the allowed range of coin cell indexes
* - Edge effects are handled by using partial coin
* - Partial coin is used near data edges and in data corners
*
* Allowed index ranges are saved to a list passed as a parameter.
* List contents are as follows [list index]:
* [0]: Min coin row index
* [1]: Min coin column index
* [2]: Max coin row index
* [3]: Max coin column index
*/
void getCoinIndexRange(const int rows, const int cols, const int current_row, const int current_col, const int coin_radius, int *index_ranges) {
if (current_row < coin_radius) {
index_ranges[0] = -current_row;
} else {
index_ranges[0] = -coin_radius;
}
if (current_col < coin_radius) {
index_ranges[1] = -current_col;
} else {
index_ranges[1] = -coin_radius;
}
if ((rows - 1 - current_row) < coin_radius) {
index_ranges[2] = rows - 1 - current_row;
} else {
index_ranges[2] = coin_radius;
}
if ((cols - 1 - current_col) < coin_radius) {
index_ranges[3] = cols - 1 - current_col;
} else {
index_ranges[3] = coin_radius;
}
}