source: CIVL/examples/modelbuilder/oneFilePhi3D.c

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 8.3 KB
Line 
1#include<stdio.h>
2#include <math.h>
3#include <stdlib.h>
4
5#define SPEED 1
6#define DEFAULT_BORDER_LOCATION -1
7#define DEFAULT_BORDER_DISTANCE 100000
8#define DEFAULT_INTERIOR_DISTANCE 90000
9#define max(a, b) ((a > b) ? a : b)
10#define min(a, b) ((a < b) ? a : b)
11
12struct phi_type {
13 int x, y, z;
14 double dx, dy, dz, F;
15 int *location;
16 double *distance;
17};
18
19typedef struct phi_type Phi;
20int max_x, max_xy;
21
22int main() {
23
24 void create_phi_function(Phi *);
25 void destroy_phi_function(Phi *);
26 void update_distance(Phi *, int);
27 void set_distance_negative_inside(Phi *, int);
28 void adjust_boundary(Phi *);
29 void run_fsm(Phi *, int);
30 void calc_dist_field(Phi *);
31 void fast_sweep(Phi *);
32 double solveEikonal(Phi *, int);
33
34
35 Phi *pf = (Phi *)malloc(sizeof(Phi));
36 if (!pf) {
37 printf("Error allocation memory for the phi function.\n");
38 exit(1);
39 }
40
41 create_phi_function(pf);
42
43 // print dimensions
44 printf("Dimensions:\n");
45 printf("x: %d\tdx: %f\n", pf->x, pf->dx);
46 printf("y: %d\tdy: %f\n", pf->y, pf->dy);
47 printf("z: %d\tdz: %f\n", pf->z, pf->dz);
48
49 // print size
50 printf("Size:\n");
51 printf("x: %d\tdx: %f\n", pf->x, pf->dx);
52 printf("y: %d\tdy: %f\n", pf->y, pf->dy);
53 printf("z: %d\tdz: %f\n", pf->z, pf->dz);
54
55 calc_dist_field(pf);
56
57 destroy_phi_function(pf);
58
59 free(pf);
60
61
62 return 0;
63
64}
65
66
67void create_phi_function(Phi *pf) {
68
69 // initialize fields for the phi function
70 pf->x = 255;
71 pf->dx = 4.0;
72 pf->y = 191;
73 pf->dy = 4.0;
74 pf->z = 127;
75 pf->dz = 1.0;
76 pf->F = 1.0;
77 printf("%lf",pf->F);
78 printf("***");
79 // allocate memory for location and distance arrays
80 int totalNodes = (pf->x + 2) * (pf->y + 2) * (pf->z + 2);
81 pf->location = (int *)malloc(sizeof(int) * totalNodes);
82 pf->distance = (double *)malloc(sizeof(double) * totalNodes);
83
84 }
85
86void destroy_phi_function(Phi *pf) {
87 free(pf->location);
88 free(pf->distance);
89}
90 void update_distance(Phi *pf, int totalNodes) {
91
92 int *l = &pf->location[0];
93 double *d = &pf->distance[0];
94
95 int i;
96 for (i = 0; i < totalNodes; i++) {
97 if (*l != DEFAULT_BORDER_LOCATION && *d != DEFAULT_BORDER_DISTANCE) {
98 //*d = (*l == 1 && *d == DEFAULT_BORDER_DISTANCE)
99 // ? -1
100 // : (*d > 0.0 || *d < 0.0) ? *d : DEFAULT_INTERIOR_DISTANCE;
101 *d = (*d > 0.0 || *d < 0.0) ? *d : DEFAULT_INTERIOR_DISTANCE;
102
103
104
105 }
106 l++;
107 d++;
108 }
109}
110
111 void set_distance_negative_inside(Phi *pf, int totalNodes) {
112
113 int *l = &pf->location[0];
114 double *d = &pf->distance[0];
115
116
117 int i;
118 for (i = 0; i < totalNodes; i++) {
119 if (*l != DEFAULT_BORDER_LOCATION && *d != DEFAULT_BORDER_DISTANCE) {
120 if (*l == 1) {
121 *d = -1;
122 }
123 }
124 l++;
125 d++;
126 }
127}
128
129 void adjust_boundary(Phi *pf) {
130
131 int x, y, z, i, j, k, xy;
132 x = pf->x + 2;
133 y = pf->y + 2;
134 z = pf->z + 2;
135 xy = x * y;
136
137 for (i = 0; i < z; i++) {
138 for (j = 0; j < y; j++) {
139 for (k = 0; k < x; k++) {
140 int I = i, J = j, K = k;
141 I = (i == z - 1) ? I - 1 : (!i) ? I + 1 : I;
142 J = (j == y - 1) ? J - 1 : (!j) ? J + 1 : J;
143 K = (k == x - 1) ? K - 1 : (!k) ? K + 1 : K;
144 if (i != I || j != J || k != K) {
145 pf->distance[i * xy + j * x + k] = pf->distance[I * xy + J * x + K];
146 }
147 }
148 }
149 }
150}
151
152double solveEikonal(Phi *pf, int index) {
153
154 double dist_new = 0;
155 double dist_old = pf->distance[index];
156
157 double dx = pf->dx, dy = pf->dy, dz = pf->dz;
158 double minX = min(pf->distance[index - 1], pf->distance[index + 1]);
159 double minY =
160 min(pf->distance[abs(index - max_x)], pf->distance[abs(index + max_x)]);
161 double minZ =
162 min(pf->distance[abs(index - max_xy)], pf->distance[abs(index + max_xy)]);
163 if(dx!= 4.000000 ||dy!=4.000000 || dz != 1.000000){
164 printf("%lf",dx);
165 printf(",");
166 printf("%lf",dy);
167 printf(",");
168 printf("%lf",dz);
169 printf(" ");
170 }
171
172 double m[] = { minX, minY, minZ };
173 double d[] = { dx, dy, dz };
174
175
176 // sort the mins
177 int i, j;
178 double tmp_m, tmp_d;
179 for (i = 1; i < 3; i++) {
180 for (j = 0; j < 3 - i; j++) {
181
182 if (m[j] > m[j + 1]) {
183 tmp_m = m[j];
184 tmp_d = d[j];
185 m[j] = m[j + 1];
186 d[j] = d[j + 1];
187 m[j + 1] = tmp_m;
188 d[j + 1] = tmp_d;
189 }
190 }
191
192 }
193
194 // simplifying the variables
195 double m_0 = m[0], m_1 = m[1], m_2 = m[2];
196 double d_0 = d[0], d_1 = d[1], d_2 = d[2];
197 double m2_0 = m_0 * m_0, m2_1 = m_1 * m_1, m2_2 = m_2 * m_2;
198 double d2_0 = d_0 * d_0, d2_1 = d_1 * d_1, d2_2 = d_2 * d_2;
199
200 if(d2_0==0 || d2_1 ==0)
201 {
202 printf("%lf",d2_0);
203 printf(",");
204 printf("%lf",d2_1);
205 }
206 dist_new = m_0 + d_0;
207 if (dist_new > m_1) {
208
209 double s = sqrt(-m2_0 + 2 * m_0 * m_1 - m2_1 + d2_0 + d2_1);
210 dist_new = (m_1 * d2_0 + m_0 * d2_1 + d_0 * d_1 * s) / (d2_0 + d2_1);
211
212 if (dist_new > m_2) {
213
214 double a =
215 sqrt(-m2_0 * d2_1 - m2_0 * d2_2 + 2 * m_0 * m_1 * d2_2 - m2_1 * d2_0 -
216 m2_1 * d2_2 + 2 * m_0 * m_2 * d2_1 - m2_2 * d2_0 - m2_2 * d2_1 +
217 2 * m_1 * m_2 * d2_0 + d2_0 * d2_1 + d2_0 * d2_2 + d2_1 * d2_2);
218
219 dist_new = (m_2 * d2_0 * d2_1 + m_1 * d2_0 * d2_2 + m_0 * d2_1 * d2_2 +
220 d_0 * d_1 * d_2 * a) /
221 (d2_0 * d2_1 + d2_0 * d2_2 + d2_1 * d2_2);
222 }
223 }
224
225
226
227
228 return min(dist_old, dist_new);
229}
230
231
232void fast_sweep(Phi *pf) {
233
234 int s, i, j, k, index;
235 // specifies the sweeping directions
236
237 int sweeps[8][3] = { { 1, 1, 1 },
238 { 0, 1, 0 },
239 { 0, 1, 1 },
240 { 1, 1, 0 },
241 { 0, 0, 0 },
242 { 1, 0, 1 },
243 { 1, 0, 0 },
244 { 0, 0, 1 } };
245
246 printf("Please wait sweeping.....\n");
247 for (s = 0; s < 8; ++s) {
248 // printf("Fast Sweeping start..... [%d/%d]\n", s, 7);
249
250 int iStart = (sweeps[s][0]) ? 1 : pf->z;
251 int iEnd = (sweeps[s][0]) ? pf->z + 1 : 0;
252
253 int jStart = (sweeps[s][1]) ? 1 : pf->y;
254 int jEnd = (sweeps[s][1]) ? pf->y + 1 : 0;
255
256 int kStart = (sweeps[s][2]) ? 1 : pf->x;
257 int kEnd = (sweeps[s][2]) ? pf->x + 1 : 0;
258
259 for (i = iStart; i != iEnd; i = (sweeps[s][0]) ? i + 1 : i - 1) {
260 for (j = jStart; j != jEnd; j = (sweeps[s][1]) ? j + 1 : j - 1) {
261 for (k = kStart; k != kEnd; k = (sweeps[s][2]) ? k + 1 : k - 1) {
262 index = i * max_xy + j * max_x + k;
263
264 //printf("%d",solveEikonal(pf,index));
265 pf->distance[index] = solveEikonal(pf, index);
266 //printf("%lf",pf->distance[index]);
267 //printf(",");
268 //printf(" ");
269
270
271
272 }
273 }
274
275 }
276 }
277
278
279 printf("Sweeping completed.......\n");
280}
281
282
283void run_fsm(Phi *pf, int iterations) {
284
285 max_x = pf->x + 2;
286 max_xy = max_x * (pf->y + 2);
287
288 double start, finish; // for timing
289 int itr = 0;
290 while (itr++ < iterations) {
291
292 // GET_TIME(start);
293
294 fast_sweep(pf);
295
296 // GET_TIME(finish);
297 //printf("Serial FSM time: %f s.\n", finish - start);
298 }
299}
300
301
302void calc_dist_field(Phi *pf) {
303
304 // get the total number of nodes
305 // in the grid
306 int totalNodes = (pf->x + 2) * (pf->y + 2) * (pf->z + 2);
307
308 // update the distance values
309 update_distance(pf, totalNodes);
310
311 // use the fast sweeping method
312 // to get the solution for the Eikonal Equation
313 int itr = 1;
314 run_fsm(pf, itr);
315
316 // set the distance values to negative
317 // for inside region
318 set_distance_negative_inside(pf, totalNodes);
319
320 adjust_boundary(pf);
321}
322
323
324
325
326
327
Note: See TracBrowser for help on using the repository browser.