JunkBox_Lib 1.10.1
Loading...
Searching...
No Matches
gmt.c
Go to the documentation of this file.
1
8#include "gmt.h"
9#include "jbxl_state.h"
10
11
24WSGraph Laplacian(WSGraph vp, int mode)
25{
26 int i, j;
27 int da, db, dc, dd, de, df, dg, dh;
28 WSGraph lp;
29
30 lp = make_WSGraph(vp.xs, vp.ys, 1);
31 if (lp.gp==NULL) return lp;
32
33 if (vp.gp==NULL) {
34 lp.state = JBXL_GRAPH_NODATA_ERROR;
35 return lp;
36 }
37
38 if (mode==4) {
39 for (j=1; j<vp.ys-1; j++) {
40 for (i=1; i<vp.xs-1; i++) {
41 da = Px(vp, i+1, j) + Px(vp, i-1, j);
42 db = Px(vp, i, j);
43 dc = Px(vp, i, j+1) + Px(vp, i, j-1);
44 Px(lp, i, j) = da - 4*db + dc;
45 }
46 }
47 }
48
49 else if (mode==8) {
50 for (j=1; j<vp.ys-1; j++) {
51 for (i=1; i<vp.xs-1; i++) {
52 da = Px(vp, i+1, j) + Px(vp, i-1, j);
53 db = Px(vp, i, j+1) + Px(vp, i, j-1);
54 dc = Px(vp, i, j);
55 dd = Px(vp, i+1, j+1) + Px(vp, i-1, j+1);
56 de = Px(vp, i+1, j-1) + Px(vp, i-1, j-1);
57 Px(lp, i, j) = da + db - 8*dc + dd + de;
58 }
59 }
60 }
61
62 else {
63 for (j=2; j<vp.ys-2; j++) {
64 for (i=2; i<vp.xs-2; i++) {
65 da = Px(vp, i, j);
66 db = Px(vp, i+1, j) + Px(vp, i-1, j) + Px(vp, i, j+1) + Px(vp, i, j-1);
67 dc = Px(vp, i-1, j-2) + Px(vp, i, j-2) + Px(vp, i+1, j-2);
68 dd = Px(vp, i-1, j+2) + Px(vp, i, j+2) + Px(vp, i+1, j+2);
69 de = Px(vp, i-2, j-1) + Px(vp, i-2, j) + Px(vp, i-2, j+1);
70 df = Px(vp, i+2, j-1) + Px(vp, i+2, j) + Px(vp, i+2, j+1);
71 dg = Px(vp, i-2, j-2) + Px(vp, i+2, j-2);
72 dh = Px(vp, i-2, j+2) + Px(vp, i+2, j+2);
73 Px(lp, i, j) = (sWord)((-24*da-8*db+4*(dc+dd+de+df)+2*(dg+dh))/64);
74 }
75 }
76 }
77
78 return lp;
79}
80
81
90WSGraph xSobel(WSGraph vp)
91{
92 int i, j, k;
93 int da, db, dc, dd, de, nr;
94 WSGraph xp;
95
96 memset(&xp, 0, sizeof(WSGraph));
97 if (vp.gp==NULL) {
98 xp.state = JBXL_GRAPH_NODATA_ERROR;
99 return xp;
100 }
101 xp.state = JBXL_NORMAL;
102
103 if (vp.zs<=0) vp.zs = 1;
104 xp = make_WSGraph(vp.xs, vp.ys, vp.zs);
105 if (xp.gp==NULL) return xp;
106
107 for (k=0; k<vp.zs; k++) {
108 for (j=1; j<vp.ys-1; j++) {
109 for (i=1; i<vp.xs-1; i++) {
110 da = Vx(vp, i+1, j-1, k) - Vx(vp, i-1, j-1, k);
111 db = Vx(vp, i+1, j, k) - Vx(vp, i-1, j, k);
112 dc = Vx(vp, i+1, j+1, k) - Vx(vp, i-1, j+1, k);
113 if (k==0 || k==vp.zs-1) {
114 dd = de = 0;
115 nr = 8;
116 }
117 else {
118 dd = Vx(vp, i+1, j, k-1) - Vx(vp, i-1, j, k-1);
119 de = Vx(vp, i+1, j, k+1) - Vx(vp, i-1, j, k+1);
120 nr = 12;
121 }
122 Vx(xp, i, j, k) = (sWord)((da + 2*db + dc + dd + de)/nr);
123 }
124 }
125 }
126
127 return xp;
128}
129
130
140FSGraph fxSobel(FSGraph vp)
141{
142 int i, j, k;
143 double da, db, dc, dd, de, nr;
144 FSGraph xp;
145
146 memset(&xp, 0, sizeof(FSGraph));
147 if (vp.gp==NULL) {
148 xp.state = JBXL_GRAPH_NODATA_ERROR;
149 return xp;
150 }
151
152 if (vp.zs<=0) vp.zs = 1;
153 xp = make_FSGraph(vp.xs, vp.ys, vp.zs);
154 if (xp.gp==NULL) return xp;
155
156 for (k=0; k<vp.zs; k++) {
157 for (j=1; j<vp.ys-1; j++) {
158 for (i=1; i<vp.xs-1; i++) {
159 da = Vx(vp, i+1, j-1, k) - Vx(vp, i-1, j-1, k);
160 db = Vx(vp, i+1, j, k) - Vx(vp, i-1, j, k);
161 dc = Vx(vp, i+1, j+1, k) - Vx(vp, i-1, j+1, k);
162 if (k==0 || k==vp.zs-1) {
163 dd = de = 0.;
164 nr = 8.;
165 }
166 else {
167 dd = Vx(vp, i+1, j, k-1) - Vx(vp, i-1, j, k-1);
168 de = Vx(vp, i+1, j, k+1) - Vx(vp, i-1, j, k+1);
169 nr = 12.;
170 }
171 Vx(xp, i, j, k) = (da + 2.*db + dc + dd + de)/nr;
172 }
173 }
174 }
175
176 return xp;
177}
178
179
188WSGraph ySobel(WSGraph vp)
189{
190 int i, j, k;
191 int da, db, dc, dd, de, nr;
192 WSGraph xp;
193
194 memset(&xp, 0, sizeof(WSGraph));
195 if (vp.gp==NULL) {
196 xp.state = JBXL_GRAPH_NODATA_ERROR;
197 return xp;
198 }
199
200 if (vp.zs<=0) vp.zs = 1;
201 xp = make_WSGraph(vp.xs, vp.ys, vp.zs);
202 if (xp.gp==NULL) return xp;
203
204 for (k=0; k<vp.zs; k++) {
205 for (j=1; j<vp.ys-1; j++) {
206 for (i=1; i<vp.xs-1; i++) {
207 da = Vx(vp, i-1, j+1, k) - Vx(vp, i-1, j-1, k);
208 db = Vx(vp, i, j+1, k) - Vx(vp, i, j-1, k);
209 dc = Vx(vp, i+1, j+1, k) - Vx(vp, i+1, j-1, k);
210 if (k==0 || k==vp.zs-1) {
211 dd = de = 0;
212 nr = 8;
213 }
214 else {
215 dd = Vx(vp, i, j+1, k-1) - Vx(vp, i, j-1, k-1);
216 de = Vx(vp, i, j+1, k+1) - Vx(vp, i, j-1, k+1);
217 nr = 12;
218 }
219 Vx(xp, i, j, k) = (sWord)((da + 2*db + dc + dd + de)/nr);
220 }
221 }
222 }
223
224 return xp;
225}
226
227
237FSGraph fySobel(FSGraph vp)
238{
239 int i, j, k;
240 double da, db, dc, dd, de, nr;
241 FSGraph xp;
242
243 memset(&xp, 0, sizeof(FSGraph));
244 if (vp.gp==NULL) {
245 xp.state = JBXL_GRAPH_NODATA_ERROR;
246 return xp;
247 }
248
249 if (vp.zs<=0) vp.zs = 1;
250 xp = make_FSGraph(vp.xs, vp.ys, vp.zs);
251 if (xp.gp==NULL) return xp;
252
253 for (k=0; k<vp.zs; k++) {
254 for (j=1; j<vp.ys-1; j++) {
255 for (i=1; i<vp.xs-1; i++) {
256 da = Vx(vp, i-1, j+1, k) - Vx(vp, i-1, j-1, k);
257 db = Vx(vp, i, j+1, k) - Vx(vp, i, j-1, k);
258 dc = Vx(vp, i+1, j+1, k) - Vx(vp, i+1, j-1, k);
259 if (k==0 || k==vp.zs-1) {
260 dd = de = 0.;
261 nr = 8.;
262 }
263 else {
264 dd = Vx(vp, i, j+1, k-1) - Vx(vp, i, j-1, k-1);
265 de = Vx(vp, i, j+1, k+1) - Vx(vp, i, j-1, k+1);
266 nr = 12.;
267 }
268 Vx(xp, i, j, k) = (da + 2.*db + dc + dd + de)/nr;
269 }
270 }
271 }
272
273 return xp;
274}
275
276
285WSGraph zSobel(WSGraph vp)
286{
287 int i, j, k;
288 int da, db, dc, dd, de;
289 WSGraph xp;
290
291 memset(&xp, 0, sizeof(WSGraph));
292 if (vp.gp==NULL) {
293 xp.state = JBXL_GRAPH_NODATA_ERROR;
294 return xp;
295 }
296
297 if (vp.zs<=1) {
298 //fprintf(stderr,"ZSOBEL: no 3D data inputed.\n");
299 xp.state = JBXL_GRAPH_IVDARG_ERROR;
300 return xp;
301 }
302
303 xp = make_WSGraph(vp.xs, vp.ys, vp.zs);
304 if (xp.gp==NULL) return xp;
305
306 for (k=1; k<vp.zs-1; k++) {
307 for (j=1; j<vp.ys-1; j++) {
308 for (i=1; i<vp.xs-1; i++) {
309 da = Vx(vp, i-1, j, k+1) - Vx(vp, i-1, j, k-1);
310 db = Vx(vp, i+1, j, k+1) - Vx(vp, i+1, j, k-1);
311 dc = Vx(vp, i, j, k+1) - Vx(vp, i, j, k-1);
312 dd = Vx(vp, i, j-1, k+1) - Vx(vp, i, j-1, k-1);
313 de = Vx(vp, i, j+1, k+1) - Vx(vp, i, j+1, k-1);
314 Vx(xp, i, j, k) = (sWord)((da + db + 2*dc + dd + de)/12);
315 }
316 }
317 }
318
319 for (j=1; j<vp.ys-1; j++) {
320 for (i=1; i<vp.xs-1; i++) {
321 da = Vx(vp, i-1, j, 1) - Vx(vp, i-1, j, 0);
322 db = Vx(vp, i+1, j, 1) - Vx(vp, i+1, j, 0);
323 dc = Vx(vp, i, j, 1) - Vx(vp, i, j, 0);
324 dd = Vx(vp, i, j-1, 1) - Vx(vp, i, j-1, 0);
325 de = Vx(vp, i, j+1, 1) - Vx(vp, i, j+1, 0);
326 Vx(xp, i, j, 0) = (sWord)((da + db + 2*dc + dd + de)/12);
327
328 da = Vx(vp, i-1, j, vp.zs-1) - Vx(vp, i-1, j, vp.zs-2);
329 db = Vx(vp, i+1, j, vp.zs-1) - Vx(vp, i+1, j, vp.zs-2);
330 dc = Vx(vp, i, j, vp.zs-1) - Vx(vp, i, j, vp.zs-2);
331 dd = Vx(vp, i, j-1, vp.zs-1) - Vx(vp, i, j-1, vp.zs-2);
332 de = Vx(vp, i, j+1, vp.zs-1) - Vx(vp, i, j+1, vp.zs-2);
333 Vx(xp, i, j, vp.zs-1) = (sWord)((da + db + 2*dc + dd + de)/12);
334 }
335 }
336
337 return xp;
338}
339
340
350FSGraph fzSobel(FSGraph vp)
351{
352 int i, j, k;
353 double da, db, dc, dd, de;
354 FSGraph xp;
355
356 memset(&xp, 0, sizeof(FSGraph));
357 if (vp.gp==NULL) {
358 xp.state = JBXL_GRAPH_NODATA_ERROR;
359 return xp;
360 }
361
362 if (vp.zs<=1) {
363 //fprintf(stderr,"FZSOBEL: no 3D data inputed.\n");
364 xp.state = JBXL_GRAPH_IVDARG_ERROR;
365 return xp;
366 }
367
368 xp = make_FSGraph(vp.xs, vp.ys, vp.zs);
369 if (xp.gp==NULL) return xp;
370
371 for (k=1; k<vp.zs-1; k++) {
372 for (j=1; j<vp.ys-1; j++) {
373 for (i=1; i<vp.xs-1; i++) {
374 da = Vx(vp, i-1, j, k+1) - Vx(vp, i-1, j, k-1);
375 db = Vx(vp, i+1, j, k+1) - Vx(vp, i+1, j, k-1);
376 dc = Vx(vp, i, j, k+1) - Vx(vp, i, j, k-1);
377 dd = Vx(vp, i, j-1, k+1) - Vx(vp, i, j-1, k-1);
378 de = Vx(vp, i, j+1, k+1) - Vx(vp, i, j+1, k-1);
379 Vx(xp, i, j, k) = (da + db + 2.*dc + dd + de)/12.;
380 }
381 }
382 }
383
384 for (j=1; j<vp.ys-1; j++) {
385 for (i=1; i<vp.xs-1; i++) {
386 da = Vx(vp, i-1, j, 1) - Vx(vp, i-1, j, 0);
387 db = Vx(vp, i+1, j, 1) - Vx(vp, i+1, j, 0);
388 dc = Vx(vp, i, j, 1) - Vx(vp, i, j, 0);
389 dd = Vx(vp, i, j-1, 1) - Vx(vp, i, j-1, 0);
390 de = Vx(vp, i, j+1, 1) - Vx(vp, i, j+1, 0);
391 Vx(xp, i, j, 0) = (da + db + 2.*dc + dd + de)/12.;
392
393 da = Vx(vp, i-1, j, vp.zs-1) - Vx(vp, i-1, j, vp.zs-2);
394 db = Vx(vp, i+1, j, vp.zs-1) - Vx(vp, i+1, j, vp.zs-2);
395 dc = Vx(vp, i, j, vp.zs-1) - Vx(vp, i, j, vp.zs-2);
396 dd = Vx(vp, i, j-1, vp.zs-1) - Vx(vp, i, j-1, vp.zs-2);
397 de = Vx(vp, i, j+1, vp.zs-1) - Vx(vp, i, j+1, vp.zs-2);
398 Vx(xp, i, j, vp.zs-1) = (da + db + 2.*dc + dd + de)/12.;
399 }
400 }
401
402 return xp;
403}
404
405
414WSGraph xxSobel(WSGraph vp)
415{
416 int x, y, z, xs, ys, zs, cx, cy, cz, ps;
417 int da, db, dc, dd, de;
418 int df, dg, dh, di, dj, dk, dl, dm;
419 WSGraph px;
420
421 memset(&px, 0, sizeof(WSGraph));
422 if (vp.gp==NULL) {
423 px.state = JBXL_GRAPH_NODATA_ERROR;
424 return px;
425 }
426
427 xs = vp.xs;
428 ys = vp.ys;
429 zs = vp.zs;
430 ps = xs*ys;
431
432 if (zs<5) {
433 px = make_WSGraph(xs, ys, 1);
434 if (px.gp==NULL) return px;
435
436 for (y=2; y<ys-2; y++) {
437 cy = y*xs;
438 for (x=2; x<xs-2; x++) {
439 cx = cy + x;
440 da = vp.gp[cx-2*xs+2] - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2];
441 db = vp.gp[cx -xs+2] - 2*vp.gp[cx -xs] + vp.gp[cx -xs-2];
442 dc = vp.gp[cx +2] - 2*vp.gp[cx] + vp.gp[cx -2];
443 dd = vp.gp[cx +xs+2] - 2*vp.gp[cx +xs] + vp.gp[cx +xs-2];
444 de = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2];
445 px.gp[cx] = (sWord)((da + 4*db + 6*dc + 4*dd + de)/64);
446 }
447 }
448 }
449
450 else {
451 px = make_WSGraph(xs, ys, zs);
452 if (px.gp==NULL) return px;
453
454 for (z=2; z<zs-2; z++) {
455 cz = z*ps;
456 for (y=2; y<ys-2; y++) {
457 cy = cz + y*xs;
458 for (x=2; x<xs-2; x++) {
459 cx = cy + x;
460 da = vp.gp[cx +2] - 2*vp.gp[cx] + vp.gp[cx-2];
461 db = vp.gp[cx+xs+2] - 2*vp.gp[cx+xs] + vp.gp[cx+xs-2];
462 dc = vp.gp[cx-xs+2] - 2*vp.gp[cx-xs] + vp.gp[cx-xs-2];
463 dd = vp.gp[cx+ps+2] - 2*vp.gp[cx+ps] + vp.gp[cx+ps-2];
464 de = vp.gp[cx-ps+2] - 2*vp.gp[cx-ps] + vp.gp[cx-ps-2];
465 df = vp.gp[cx+xs+ps+2] - 2*vp.gp[cx+xs+ps] + vp.gp[cx+xs+ps-2];
466 dg = vp.gp[cx+xs-ps+2] - 2*vp.gp[cx+xs-ps] + vp.gp[cx+xs-ps-2];
467 dh = vp.gp[cx-xs+ps+2] - 2*vp.gp[cx-xs+ps] + vp.gp[cx-xs+ps-2];
468 di = vp.gp[cx-xs-ps+2] - 2*vp.gp[cx-xs-ps] + vp.gp[cx-xs-ps-2];
469 dj = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2];
470 dk = vp.gp[cx-2*xs+2] - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2];
471 dl = vp.gp[cx+2*ps+2] - 2*vp.gp[cx+2*ps] + vp.gp[cx+2*ps-2];
472 dm = vp.gp[cx-2*ps+2] - 2*vp.gp[cx-2*ps] + vp.gp[cx-2*ps-2];
473 px.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
474 }
475 }
476 }
477 }
478
479 return px;
480}
481
482
492FSGraph fxxSobel(FSGraph vp)
493{
494 int x, y, z, xs, ys, zs, cx, cy, cz, ps;
495 double da, db, dc, dd, de;
496 double df, dg, dh, di, dj, dk, dl, dm;
497 FSGraph px;
498
499 memset(&px, 0, sizeof(FSGraph));
500 if (vp.gp==NULL) {
501 px.state = JBXL_GRAPH_NODATA_ERROR;
502 return px;
503 }
504
505 xs = vp.xs;
506 ys = vp.ys;
507 zs = vp.zs;
508 ps = xs*ys;
509
510 if (zs<5) {
511 px = make_FSGraph(xs, ys, 1);
512 if (px.gp==NULL) return px;
513
514 for (y=2; y<ys-2; y++) {
515 cy = y*xs;
516 for (x=2; x<xs-2; x++) {
517 cx = cy + x;
518 da = vp.gp[cx-2*xs+2] - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2];
519 db = vp.gp[cx -xs+2] - 2*vp.gp[cx -xs] + vp.gp[cx -xs-2];
520 dc = vp.gp[cx +2] - 2*vp.gp[cx] + vp.gp[cx -2];
521 dd = vp.gp[cx +xs+2] - 2*vp.gp[cx +xs] + vp.gp[cx +xs-2];
522 de = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2];
523 px.gp[cx] = (da + 4*db + 6*dc + 4*dd + de)/64.;
524 }
525 }
526 }
527
528 else {
529 px = make_FSGraph(xs, ys, zs);
530 if (px.gp==NULL) return px;
531
532 for (z=2; z<zs-2; z++) {
533 cz = z*ps;
534 for (y=2; y<ys-2; y++) {
535 cy = cz + y*xs;
536 for (x=2; x<xs-2; x++) {
537 cx = cy + x;
538 da = vp.gp[cx +2] - 2*vp.gp[cx] + vp.gp[cx-2];
539 db = vp.gp[cx+xs+2] - 2*vp.gp[cx+xs] + vp.gp[cx+xs-2];
540 dc = vp.gp[cx-xs+2] - 2*vp.gp[cx-xs] + vp.gp[cx-xs-2];
541 dd = vp.gp[cx+ps+2] - 2*vp.gp[cx+ps] + vp.gp[cx+ps-2];
542 de = vp.gp[cx-ps+2] - 2*vp.gp[cx-ps] + vp.gp[cx-ps-2];
543 df = vp.gp[cx+xs+ps+2] - 2*vp.gp[cx+xs+ps] + vp.gp[cx+xs+ps-2];
544 dg = vp.gp[cx+xs-ps+2] - 2*vp.gp[cx+xs-ps] + vp.gp[cx+xs-ps-2];
545 dh = vp.gp[cx-xs+ps+2] - 2*vp.gp[cx-xs+ps] + vp.gp[cx-xs+ps-2];
546 di = vp.gp[cx-xs-ps+2] - 2*vp.gp[cx-xs-ps] + vp.gp[cx-xs-ps-2];
547 dj = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2];
548 dk = vp.gp[cx-2*xs+2] - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2];
549 dl = vp.gp[cx+2*ps+2] - 2*vp.gp[cx+2*ps] + vp.gp[cx+2*ps-2];
550 dm = vp.gp[cx-2*ps+2] - 2*vp.gp[cx-2*ps] + vp.gp[cx-2*ps-2];
551 px.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
552 }
553 }
554 }
555 }
556
557 return px;
558}
559
560
569WSGraph yySobel(WSGraph vp)
570{
571 int x, y, z, xs, ys, zs, cx, cy, cz, ps;
572 int da, db, dc, dd, de;
573 int df, dg, dh, di, dj, dk, dl, dm;
574 WSGraph py;
575
576 memset(&py, 0, sizeof(WSGraph));
577 if (vp.gp==NULL) {
578 py.state = JBXL_GRAPH_NODATA_ERROR;
579 return py;
580 }
581
582 xs = vp.xs;
583 ys = vp.ys;
584 zs = vp.zs;
585 ps = xs*ys;
586
587 if (zs<5) {
588 py = make_WSGraph(xs, ys, 1);
589 if (py.gp==NULL) return py;
590
591 for (y=2; y<ys-2; y++) {
592 cy = y*xs;
593 for (x=2; x<xs-2; x++) {
594 cx = cy + x;
595 da = vp.gp[cx+2*xs-2] - 2*vp.gp[cx-2] + vp.gp[cx-2*xs-2];
596 db = vp.gp[cx+2*xs-1] - 2*vp.gp[cx-1] + vp.gp[cx-2*xs-1];
597 dc = vp.gp[cx+2*xs] - 2*vp.gp[cx] + vp.gp[cx-2*xs];
598 dd = vp.gp[cx+2*xs+1] - 2*vp.gp[cx+1] + vp.gp[cx-2*xs+1];
599 de = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2] + vp.gp[cx-2*xs+2];
600 py.gp[cx] = (sWord)((da + 4*db + 6*dc + 4*dd + de)/64);
601 }
602 }
603 }
604 else {
605 py = make_WSGraph(xs, ys, zs);
606 if (py.gp==NULL) return py;
607
608 for (z=2; z<zs-2; z++) {
609 cz = z*ps;
610 for (y=2; y<ys-2; y++) {
611 cy = cz + y*xs;
612 for (x=2; x<xs-2; x++) {
613 cx = cy + x;
614 da = vp.gp[cx+2*xs] - 2*vp.gp[cx] + vp.gp[cx-2*xs];
615 db = vp.gp[cx+1+2*xs] - 2*vp.gp[cx+1] + vp.gp[cx+1-2*xs];
616 dc = vp.gp[cx-1+2*xs] - 2*vp.gp[cx-1] + vp.gp[cx-1-2*xs];
617 dd = vp.gp[cx+ps+2*xs] - 2*vp.gp[cx+ps] + vp.gp[cx+ps-2*xs];
618 de = vp.gp[cx-ps+2*xs] - 2*vp.gp[cx-ps] + vp.gp[cx-ps-2*xs];
619 df = vp.gp[cx+1+ps+2*xs] - 2*vp.gp[cx+1+ps] + vp.gp[cx+1+ps-2*xs];
620 dg = vp.gp[cx+1-ps+2*xs] - 2*vp.gp[cx+1-ps] + vp.gp[cx+1-ps-2*xs];
621 dh = vp.gp[cx-1+ps+2*xs] - 2*vp.gp[cx-1+ps] + vp.gp[cx-1+ps-2*xs];
622 di = vp.gp[cx-1-ps+2*xs] - 2*vp.gp[cx-1-ps] + vp.gp[cx-1-ps-2*xs];
623 dj = vp.gp[cx+2+2*xs] - 2*vp.gp[cx+2] + vp.gp[cx+2-2*xs];
624 dk = vp.gp[cx-2+2*xs] - 2*vp.gp[cx-2] + vp.gp[cx-2-2*xs];
625 dl = vp.gp[cx+2*ps+2*xs] - 2*vp.gp[cx+2*ps] + vp.gp[cx+2*ps-2*xs];
626 dm = vp.gp[cx-2*ps+2*xs] - 2*vp.gp[cx-2*ps] + vp.gp[cx-2*ps-2*xs];
627 py.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
628 }
629 }
630 }
631 }
632
633 return py;
634}
635
636
646FSGraph fyySobel(FSGraph vp)
647{
648 int x, y, z, xs, ys, zs, cx, cy, cz, ps;
649 double da, db, dc, dd, de;
650 double df, dg, dh, di, dj, dk, dl, dm;
651 FSGraph py;
652
653 memset(&py, 0, sizeof(FSGraph));
654 if (vp.gp==NULL) {
655 py.state = JBXL_GRAPH_NODATA_ERROR;
656 return py;
657 }
658
659 xs = vp.xs;
660 ys = vp.ys;
661 zs = vp.zs;
662 ps = xs*ys;
663
664 if (zs<5) {
665 py = make_FSGraph(xs, ys, 1);
666 if (py.gp==NULL) return py;
667
668 for (y=2; y<ys-2; y++) {
669 cy = y*xs;
670 for (x=2; x<xs-2; x++) {
671 cx = cy + x;
672 da = vp.gp[cx+2*xs-2] - 2*vp.gp[cx-2] + vp.gp[cx-2*xs-2];
673 db = vp.gp[cx+2*xs-1] - 2*vp.gp[cx-1] + vp.gp[cx-2*xs-1];
674 dc = vp.gp[cx+2*xs] - 2*vp.gp[cx] + vp.gp[cx-2*xs];
675 dd = vp.gp[cx+2*xs+1] - 2*vp.gp[cx+1] + vp.gp[cx-2*xs+1];
676 de = vp.gp[cx+2*xs+2] - 2*vp.gp[cx+2] + vp.gp[cx-2*xs+2];
677 py.gp[cx] = (da + 4*db + 6*dc + 4*dd + de)/64.;
678 }
679 }
680 }
681 else {
682 py = make_FSGraph(xs, ys, zs);
683 if (py.gp==NULL) return py;
684
685 for (z=2; z<zs-2; z++) {
686 cz = z*ps;
687 for (y=2; y<ys-2; y++) {
688 cy = cz + y*xs;
689 for (x=2; x<xs-2; x++) {
690 cx = cy + x;
691 da = vp.gp[cx+2*xs] - 2*vp.gp[cx] + vp.gp[cx-2*xs];
692 db = vp.gp[cx+1+2*xs] - 2*vp.gp[cx+1] + vp.gp[cx+1-2*xs];
693 dc = vp.gp[cx-1+2*xs] - 2*vp.gp[cx-1] + vp.gp[cx-1-2*xs];
694 dd = vp.gp[cx+ps+2*xs] - 2*vp.gp[cx+ps] + vp.gp[cx+ps-2*xs];
695 de = vp.gp[cx-ps+2*xs] - 2*vp.gp[cx-ps] + vp.gp[cx-ps-2*xs];
696 df = vp.gp[cx+1+ps+2*xs] - 2*vp.gp[cx+1+ps] + vp.gp[cx+1+ps-2*xs];
697 dg = vp.gp[cx+1-ps+2*xs] - 2*vp.gp[cx+1-ps] + vp.gp[cx+1-ps-2*xs];
698 dh = vp.gp[cx-1+ps+2*xs] - 2*vp.gp[cx-1+ps] + vp.gp[cx-1+ps-2*xs];
699 di = vp.gp[cx-1-ps+2*xs] - 2*vp.gp[cx-1-ps] + vp.gp[cx-1-ps-2*xs];
700 dj = vp.gp[cx+2+2*xs] - 2*vp.gp[cx+2] + vp.gp[cx+2-2*xs];
701 dk = vp.gp[cx-2+2*xs] - 2*vp.gp[cx-2] + vp.gp[cx-2-2*xs];
702 dl = vp.gp[cx+2*ps+2*xs] - 2*vp.gp[cx+2*ps] + vp.gp[cx+2*ps-2*xs];
703 dm = vp.gp[cx-2*ps+2*xs] - 2*vp.gp[cx-2*ps] + vp.gp[cx-2*ps-2*xs];
704 py.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
705 }
706 }
707 }
708 }
709
710 return py;
711}
712
713
722WSGraph zzSobel(WSGraph vp)
723{
724 int x, y, z, xs, ys, zs, cx, cy, cz, ps;
725 int da, db, dc, dd, de;
726 int df, dg, dh, di, dj, dk, dl, dm;
727 WSGraph pz;
728
729 memset(&pz, 0, sizeof(WSGraph));
730 if (vp.gp==NULL) {
731 pz.state = JBXL_GRAPH_NODATA_ERROR;
732 return pz;
733 }
734
735 if (vp.zs<5) {
736 //fprintf(stderr,"ZZSOBEL: no 3D data inputed.\n");
737 pz.state = JBXL_GRAPH_IVDARG_ERROR;
738 return pz;
739 }
740
741 xs = vp.xs;
742 ys = vp.ys;
743 zs = vp.zs;
744 ps = xs*ys;
745 pz = make_WSGraph(xs, ys, zs);
746 if (pz.gp==NULL) return pz;
747
748 for (z=2; z<zs-2; z++) {
749 cz = z*ps;
750 for (y=2; y<ys-2; y++) {
751 cy = cz + y*xs;
752 for (x=2; x<xs-2; x++) {
753 cx = cy + x;
754 da = vp.gp[cx +2*ps] - 2*vp.gp[cx] + vp.gp[cx -2*ps];
755 db = vp.gp[cx+1 +2*ps] - 2*vp.gp[cx+1] + vp.gp[cx+1 -2*ps];
756 dc = vp.gp[cx-1 +2*ps] - 2*vp.gp[cx-1] + vp.gp[cx-1 -2*ps];
757 dd = vp.gp[cx+xs+2*ps] - 2*vp.gp[cx+xs] + vp.gp[cx+xs-2*ps];
758 de = vp.gp[cx-xs+2*ps] - 2*vp.gp[cx-xs] + vp.gp[cx-xs-2*ps];
759 df = vp.gp[cx+1+xs+2*ps] - 2*vp.gp[cx+1+xs] + vp.gp[cx+1+xs-2*ps];
760 dg = vp.gp[cx+1-xs+2*ps] - 2*vp.gp[cx+1-xs] + vp.gp[cx+1-xs-2*ps];
761 dh = vp.gp[cx-1+xs+2*ps] - 2*vp.gp[cx-1+xs] + vp.gp[cx-1+xs-2*ps];
762 di = vp.gp[cx-1-xs+2*ps] - 2*vp.gp[cx-1-xs] + vp.gp[cx-1-xs-2*ps];
763 dj = vp.gp[cx+2 +2*ps] - 2*vp.gp[cx+2] + vp.gp[cx+2 -2*ps];
764 dk = vp.gp[cx-2 +2*ps] - 2*vp.gp[cx-2] + vp.gp[cx-2 -2*ps];
765 dl = vp.gp[cx+2*xs+2*ps] - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2*ps];
766 dm = vp.gp[cx-2*xs+2*ps] - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2*ps];
767 pz.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
768 }
769 }
770 }
771
772 cz = ps;
773 for (y=2; y<ys-2; y++) {
774 cy = cz + y*xs;
775 for (x=2; x<xs-2; x++) {
776 cx = cy + x;
777 da = vp.gp[cx +2*ps] - 2*vp.gp[cx];
778 db = vp.gp[cx+1 +2*ps] - 2*vp.gp[cx+1];
779 dc = vp.gp[cx-1 +2*ps] - 2*vp.gp[cx-1];
780 dd = vp.gp[cx+xs+2*ps] - 2*vp.gp[cx+xs];
781 de = vp.gp[cx-xs+2*ps] - 2*vp.gp[cx-xs];
782 df = vp.gp[cx+1+xs+2*ps] - 2*vp.gp[cx+1+xs];
783 dg = vp.gp[cx+1-xs+2*ps] - 2*vp.gp[cx+1-xs];
784 dh = vp.gp[cx-1+xs+2*ps] - 2*vp.gp[cx-1+xs];
785 di = vp.gp[cx-1-xs+2*ps] - 2*vp.gp[cx-1-xs];
786 dj = vp.gp[cx+2 +2*ps] - 2*vp.gp[cx+2];
787 dk = vp.gp[cx-2 +2*ps] - 2*vp.gp[cx-2];
788 dl = vp.gp[cx+2*xs+2*ps] - 2*vp.gp[cx+2*xs];
789 dm = vp.gp[cx-2*xs+2*ps] - 2*vp.gp[cx-2*xs];
790 pz.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
791 }
792 }
793
794 cz = (zs-2)*ps;
795 for (y=2; y<ys-2; y++) {
796 cy = cz + y*xs;
797 for (x=2; x<xs-2; x++) {
798 cx = cy + x;
799 da = - 2*vp.gp[cx] + vp.gp[cx -2*ps];
800 db = - 2*vp.gp[cx+1] + vp.gp[cx+1 -2*ps];
801 dc = - 2*vp.gp[cx-1] + vp.gp[cx-1 -2*ps];
802 dd = - 2*vp.gp[cx+xs] + vp.gp[cx+xs-2*ps];
803 de = - 2*vp.gp[cx-xs] + vp.gp[cx-xs-2*ps];
804 df = - 2*vp.gp[cx+1+xs] + vp.gp[cx+1+xs-2*ps];
805 dg = - 2*vp.gp[cx+1-xs] + vp.gp[cx+1-xs-2*ps];
806 dh = - 2*vp.gp[cx-1+xs] + vp.gp[cx-1+xs-2*ps];
807 di = - 2*vp.gp[cx-1-xs] + vp.gp[cx-1-xs-2*ps];
808 dj = - 2*vp.gp[cx+2] + vp.gp[cx+2 -2*ps];
809 dk = - 2*vp.gp[cx-2] + vp.gp[cx-2 -2*ps];
810 dl = - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2*ps];
811 dm = - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2*ps];
812 pz.gp[cx] = (sWord)((8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144);
813 }
814 }
815
816 return pz;
817}
818
819
829FSGraph fzzSobel(FSGraph vp)
830{
831 int x, y, z, xs, ys, zs, cx, cy, cz, ps;
832 double da, db, dc, dd, de;
833 double df, dg, dh, di, dj, dk, dl, dm;
834 FSGraph pz;
835
836 memset(&pz, 0, sizeof(FSGraph));
837 if (vp.gp==NULL) {
838 pz.state = JBXL_GRAPH_NODATA_ERROR;
839 return pz;
840 }
841
842 if (vp.zs<5) {
843 //fprintf(stderr,"FZZSOBEL: no 3D data inputed.\n");
844 pz.state = JBXL_GRAPH_IVDARG_ERROR;
845 return pz;
846 }
847
848 xs = vp.xs;
849 ys = vp.ys;
850 zs = vp.zs;
851 ps = xs*ys;
852 pz = make_FSGraph(xs, ys, zs);
853 if (pz.gp==NULL) return pz;
854
855 for (z=2; z<zs-2; z++) {
856 cz = z*ps;
857 for (y=2; y<ys-2; y++) {
858 cy = cz + y*xs;
859 for (x=2; x<xs-2; x++) {
860 cx = cy + x;
861 da = vp.gp[cx +2*ps] - 2*vp.gp[cx] + vp.gp[cx -2*ps];
862 db = vp.gp[cx+1 +2*ps] - 2*vp.gp[cx+1] + vp.gp[cx+1 -2*ps];
863 dc = vp.gp[cx-1 +2*ps] - 2*vp.gp[cx-1] + vp.gp[cx-1 -2*ps];
864 dd = vp.gp[cx +xs+2*ps] - 2*vp.gp[cx +xs] + vp.gp[cx +xs-2*ps];
865 de = vp.gp[cx -xs+2*ps] - 2*vp.gp[cx -xs] + vp.gp[cx -xs-2*ps];
866 df = vp.gp[cx+1+xs+2*ps] - 2*vp.gp[cx+1+xs] + vp.gp[cx+1+xs-2*ps];
867 dg = vp.gp[cx+1-xs+2*ps] - 2*vp.gp[cx+1-xs] + vp.gp[cx+1-xs-2*ps];
868 dh = vp.gp[cx-1+xs+2*ps] - 2*vp.gp[cx-1+xs] + vp.gp[cx-1+xs-2*ps];
869 di = vp.gp[cx-1-xs+2*ps] - 2*vp.gp[cx-1-xs] + vp.gp[cx-1-xs-2*ps];
870 dj = vp.gp[cx+2 +2*ps] - 2*vp.gp[cx+2] + vp.gp[cx+2 -2*ps];
871 dk = vp.gp[cx-2 +2*ps] - 2*vp.gp[cx-2] + vp.gp[cx-2 -2*ps];
872 dl = vp.gp[cx+2*xs+2*ps] - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2*ps];
873 dm = vp.gp[cx-2*xs+2*ps] - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2*ps];
874 pz.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
875 }
876 }
877 }
878
879 cz = ps;
880 for (y=2; y<ys-2; y++) {
881 cy = cz + y*xs;
882 for (x=2; x<xs-2; x++) {
883 cx = cy + x;
884 da = vp.gp[cx +2*ps] - 2*vp.gp[cx];
885 db = vp.gp[cx+1 +2*ps] - 2*vp.gp[cx+1];
886 dc = vp.gp[cx-1 +2*ps] - 2*vp.gp[cx-1];
887 dd = vp.gp[cx+xs+2*ps] - 2*vp.gp[cx+xs];
888 de = vp.gp[cx-xs+2*ps] - 2*vp.gp[cx-xs];
889 df = vp.gp[cx+1+xs+2*ps] - 2*vp.gp[cx+1+xs];
890 dg = vp.gp[cx+1-xs+2*ps] - 2*vp.gp[cx+1-xs];
891 dh = vp.gp[cx-1+xs+2*ps] - 2*vp.gp[cx-1+xs];
892 di = vp.gp[cx-1-xs+2*ps] - 2*vp.gp[cx-1-xs];
893 dj = vp.gp[cx+2 +2*ps] - 2*vp.gp[cx+2];
894 dk = vp.gp[cx-2 +2*ps] - 2*vp.gp[cx-2];
895 dl = vp.gp[cx+2*xs+2*ps] - 2*vp.gp[cx+2*xs];
896 dm = vp.gp[cx-2*xs+2*ps] - 2*vp.gp[cx-2*xs];
897 pz.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
898 }
899 }
900
901 cz = (zs-2)*ps;
902 for (y=2; y<ys-2; y++) {
903 cy = cz + y*xs;
904 for (x=2; x<xs-2; x++) {
905 cx = cy + x;
906 da = - 2*vp.gp[cx] + vp.gp[cx -2*ps];
907 db = - 2*vp.gp[cx+1] + vp.gp[cx+1 -2*ps];
908 dc = - 2*vp.gp[cx-1] + vp.gp[cx-1 -2*ps];
909 dd = - 2*vp.gp[cx+xs] + vp.gp[cx+xs-2*ps];
910 de = - 2*vp.gp[cx-xs] + vp.gp[cx-xs-2*ps];
911 df = - 2*vp.gp[cx+1+xs] + vp.gp[cx+1+xs-2*ps];
912 dg = - 2*vp.gp[cx+1-xs] + vp.gp[cx+1-xs-2*ps];
913 dh = - 2*vp.gp[cx-1+xs] + vp.gp[cx-1+xs-2*ps];
914 di = - 2*vp.gp[cx-1-xs] + vp.gp[cx-1-xs-2*ps];
915 dj = - 2*vp.gp[cx+2] + vp.gp[cx+2 -2*ps];
916 dk = - 2*vp.gp[cx-2] + vp.gp[cx-2 -2*ps];
917 dl = - 2*vp.gp[cx+2*xs] + vp.gp[cx+2*xs-2*ps];
918 dm = - 2*vp.gp[cx-2*xs] + vp.gp[cx-2*xs-2*ps];
919 pz.gp[cx] = (8*da+4*(db+dc+dd+de)+2*(df+dg+dh+di)+dj+dk+dl+dm)/144.;
920 }
921 }
922
923 return pz;
924}
925
926
935VSGraph vNabra(WSGraph vp) // Sobel
936{
937 int i, xs, ys, zs;
938 double xx, yy, zz;
939 WSGraph px, py, pz;
940 VSGraph pn;
941
942 memset(&pn, 0, sizeof(VSGraph));
943 if (vp.gp==NULL) {
944 pn.state = JBXL_GRAPH_NODATA_ERROR;
945 return pn;
946 }
947
948 xs = vp.xs;
949 ys = vp.ys;
950 zs = vp.zs;
951 pn = make_VSGraph(xs, ys, zs);
952 if (pn.gp==NULL) return pn;
953
954 px = xSobel(vp);
955 if (px.gp==NULL) {
956 free_VSGraph(&pn);
957 pn.state = px.state;
958 return pn;
959 }
960 py = ySobel(vp);
961 if (py.gp==NULL) {
962 free_VSGraph(&pn);
963 free(px.gp);
964 pn.state = py.state;
965 return pn;
966 }
967
968 if (vp.zs<3) {
969 for (i=0; i<xs*ys; i++) {
970 xx = px.gp[i];
971 yy = py.gp[i];
972 pn.gp[i] = set_vector(xx, yy, 0.0);
973 unit_vector(pn.gp[i]);
974 }
975 }
976 else {
977 pz = zSobel(vp);
978 if (pz.gp==NULL) {
979 free_VSGraph(&pn);
980 free(px.gp);
981 free(py.gp);
982 pn.state = pz.state;
983 return pn;
984 }
985
986 for (i=0; i<xs*ys*zs; i++) {
987 xx = px.gp[i];
988 yy = py.gp[i];
989 zz = pz.gp[i];
990 pn.gp[i] = set_vector(xx, yy, zz);
991 unit_vector(pn.gp[i]);
992 }
993 free(pz.gp);
994 }
995
996 free(px.gp);
997 free(py.gp);
998
999 return pn;
1000}
1001
1002
1012VSGraph vfNabra(FSGraph vp) // Sobel
1013{
1014 int i, xs, ys, zs;
1015 double xx, yy, zz;
1016 FSGraph px, py, pz;
1017 VSGraph pn;
1018
1019 memset(&pn, 0, sizeof(VSGraph));
1020 if (vp.gp==NULL) {
1021 pn.state = JBXL_GRAPH_NODATA_ERROR;
1022 return pn;
1023 }
1024
1025 xs = vp.xs;
1026 ys = vp.ys;
1027 zs = vp.zs;
1028 pn = make_VSGraph(xs, ys, zs);
1029 if (pn.gp==NULL) return pn;
1030
1031 px = fxSobel(vp);
1032 if (px.gp==NULL) {
1033 free_VSGraph(&pn);
1034 pn.state = px.state;
1035 return pn;
1036 }
1037 py = fySobel(vp);
1038 if (py.gp==NULL) {
1039 free_VSGraph(&pn);
1040 free(px.gp);
1041 pn.state = py.state;
1042 return pn;
1043 }
1044
1045 if (vp.zs<3) {
1046 for (i=0; i<xs*ys; i++) {
1047 xx = px.gp[i];
1048 yy = py.gp[i];
1049 pn.gp[i] = set_vector(xx, yy, 0.0);
1050 unit_vector(pn.gp[i]);
1051 }
1052 }
1053 else {
1054 pz = fzSobel(vp);
1055 if (pz.gp==NULL) {
1056 free_VSGraph(&pn);
1057 free(px.gp);
1058 free(py.gp);
1059 pn.state = pz.state;
1060 return pn;
1061 }
1062
1063 for (i=0; i<xs*ys*zs; i++) {
1064 xx = px.gp[i];
1065 yy = py.gp[i];
1066 zz = pz.gp[i];
1067 pn.gp[i] = set_vector(xx, yy, zz);
1068 unit_vector(pn.gp[i]);
1069 }
1070 free(pz.gp);
1071 }
1072
1073 free(px.gp);
1074 free(py.gp);
1075
1076 return pn;
1077}
1078
1079
1088WSGraph Nabra(WSGraph vp) // Sobel
1089{
1090 int i, xs, ys, zs;
1091 int xx, yy, zz;
1092 WSGraph px, py, pz, pn;
1093
1094 memset(&pn, 0, sizeof(WSGraph));
1095 if (vp.gp==NULL) {
1096 pn.state = JBXL_GRAPH_NODATA_ERROR;
1097 return pn;
1098 }
1099
1100 xs = vp.xs;
1101 ys = vp.ys;
1102 zs = vp.zs;
1103 pn = make_WSGraph(xs, ys, zs);
1104 if (pn.gp==NULL) return pn;
1105
1106 px = xSobel(vp);
1107 if (px.gp==NULL) {
1108 free_WSGraph(&pn);
1109 pn.state = px.state;
1110 return pn;
1111 }
1112 py = ySobel(vp);
1113 if (py.gp==NULL) {
1114 free_WSGraph(&pn);
1115 free(px.gp);
1116 pn.state = py.state;
1117 return pn;
1118 }
1119
1120 if (vp.zs<3) {
1121 for (i=0; i<xs*ys; i++) {
1122 xx = px.gp[i];
1123 yy = py.gp[i];
1124 pn.gp[i] = (sWord)sqrt(xx*xx + yy*yy);
1125 }
1126 }
1127 else {
1128 pz = zSobel(vp);
1129 if (pz.gp==NULL) {
1130 free_WSGraph(&pn);
1131 free(px.gp);
1132 free(py.gp);
1133 pn.state = pz.state;
1134 return pn;
1135 }
1136
1137 for (i=0; i<xs*ys*zs; i++) {
1138 xx = px.gp[i];
1139 yy = py.gp[i];
1140 zz = pz.gp[i];
1141 pn.gp[i] = (sWord)sqrt(xx*xx + yy*yy + zz*zz);
1142 }
1143 free(pz.gp);
1144 }
1145
1146 free(px.gp);
1147 free(py.gp);
1148
1149 return pn;
1150}
1151
1152
1162FSGraph fNabra(FSGraph vp) // Sobel
1163{
1164 int i, xs, ys, zs;
1165 double xx, yy, zz;
1166 FSGraph px, py, pz, pn;
1167
1168 memset(&pn, 0, sizeof(FSGraph));
1169 if (vp.gp==NULL) {
1170 pn.state = JBXL_GRAPH_NODATA_ERROR;
1171 return pn;
1172 }
1173
1174 xs = vp.xs;
1175 ys = vp.ys;
1176 zs = vp.zs;
1177 pn = make_FSGraph(xs, ys, zs);
1178 if (pn.gp==NULL) return pn;
1179
1180 px = fxSobel(vp);
1181 if (px.gp==NULL) {
1182 free_FSGraph(&pn);
1183 pn.state = px.state;
1184 return pn;
1185 }
1186 py = fySobel(vp);
1187 if (py.gp==NULL) {
1188 free_FSGraph(&pn);
1189 free(px.gp);
1190 pn.state = py.state;
1191 return pn;
1192 }
1193
1194 if (vp.zs<3) {
1195 for (i=0; i<xs*ys; i++) {
1196 xx = px.gp[i];
1197 yy = py.gp[i];
1198 pn.gp[i] = sqrt(xx*xx + yy*yy);
1199 }
1200 }
1201 else {
1202 pz = fzSobel(vp);
1203 if (pz.gp==NULL) {
1204 free_FSGraph(&pn);
1205 free(px.gp);
1206 free(py.gp);
1207 pn.state = pz.state;
1208 return pn;
1209 }
1210
1211 for (i=0; i<xs*ys*zs; i++) {
1212 xx = px.gp[i];
1213 yy = py.gp[i];
1214 zz = pz.gp[i];
1215 pn.gp[i] = sqrt(xx*xx + yy*yy + zz*zz);
1216 }
1217 free(pz.gp);
1218 }
1219
1220 free(px.gp);
1221 free(py.gp);
1222
1223 return pn;
1224}
1225
1226
1235VSGraph curvature3D(FSGraph vp)
1236{
1237 int i;
1238 double alph, beta, gamm, K, H;
1239 double fx, fy, fz, fxy, fyz, fzx, fxx, fyy, fzz, nb;
1240 FSGraph px, py, pz, pxy, pyz, pzx, pxx, pyy, pzz, nab;
1241 VSGraph pp;
1242
1243 memset(&pp, 0, sizeof(VSGraph));
1244 if (vp.gp==NULL) {
1245 pp.state = JBXL_GRAPH_NODATA_ERROR;
1246 return pp;
1247 }
1248
1249 if (vp.zs<5) {
1250 //fprintf(stderr,"CURVATURE3D: z dimension is < 5.\n");
1251 pp.state = JBXL_GRAPH_IVDARG_ERROR;
1252 return pp;
1253 }
1254
1255 pp = make_VSGraph(vp.xs, vp.ys, vp.zs);
1256 if (pp.gp==NULL) return pp;
1257
1258 nab = fNabra(vp);
1259 px = fxSobel(vp);
1260 py = fySobel(vp);
1261 pz = fzSobel(vp);
1262 pxy = fySobel(px);
1263 pyz = fzSobel(py);
1264 pzx = fxSobel(pz);
1265 pxx = fxxSobel(vp);
1266 pyy = fyySobel(vp);
1267 pzz = fzzSobel(vp);
1268
1269 if (nab.gp==NULL || px.gp==NULL || py.gp==NULL || pz.gp==NULL ||
1270 pxy.gp==NULL || pyz.gp==NULL || pzx.gp==NULL ||
1271 pxx.gp==NULL || pyy.gp==NULL || pzz.gp==NULL) {
1272 freeNull(px.gp);
1273 freeNull(py.gp);
1274 freeNull(pz.gp);
1275 freeNull(pxy.gp);
1276 freeNull(pyz.gp);
1277 freeNull(pzx.gp);
1278 freeNull(pxx.gp);
1279 freeNull(pyy.gp);
1280 freeNull(pzz.gp);
1281 freeNull(nab.gp);
1282 free_VSGraph(&pp);
1283 pp.state = JBXL_GRAPH_ERROR;
1284 return pp;
1285 }
1286
1287 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
1288 nb = nab.gp[i];
1289 fx = px.gp[i];
1290 fy = py.gp[i];
1291 fz = pz.gp[i];
1292 fxy = pxy.gp[i];
1293 fyz = pyz.gp[i];
1294 fzx = pzx.gp[i];
1295 fxx = pxx.gp[i];
1296 fyy = pyy.gp[i];
1297 fzz = pzz.gp[i];
1298
1299 if (nb*(fx*fx+fy*fy) !=0) {
1300 alph = (2*fx*fy*fxy - fx*fx*fyy - fy*fy*fxx)/(fx*fx+fy*fy);
1301 beta = (2*fz*(fx*fx+fy*fy)*(fx*fzx+fy*fyz) - 2*fx*fy*fz*fz*fxy
1302 - fx*fx*fz*fz*fxx - fy*fy*fz*fz*fyy
1303 - (fx*fx+fy*fy)*(fx*fx+fy*fy)*fzz)/(nb*nb*(fx*fx+fy*fy));
1304 gamm = ((fx*fx+fy*fy)*(fy*fzx-fx*fyz) + (fx*fx-fy*fy)*fz*fxy
1305 - fx*fy*fz*(fxx-fyy))/(nb*(fx*fx+fy*fy));
1306
1307 K = alph*beta - gamm*gamm;
1308 H = (alph + beta)/2;
1309 pp.gp[i] = set_vector(K, H, 0.0);
1310 }
1311 }
1312
1313 free(px.gp);
1314 free(py.gp);
1315 free(pz.gp);
1316 free(pxy.gp);
1317 free(pyz.gp);
1318 free(pzx.gp);
1319 free(pxx.gp);
1320 free(pyy.gp);
1321 free(pzz.gp);
1322 free(nab.gp);
1323
1324 return pp;
1325}
1326
1327
1336VSGraph curvature(FSGraph vp)
1337{
1338 int i;
1339 double K, H, d;
1340 double Ix, Ixx, Iy, Iyy, Ixy;
1341 FSGraph px, py, pxy, pxx, pyy;
1342 VSGraph pp;
1343
1344 memset(&pp, 0, sizeof(VSGraph));
1345 if (vp.gp==NULL) {
1346 pp.state = JBXL_GRAPH_NODATA_ERROR;
1347 return pp;
1348 }
1349
1350 if (vp.zs>1) {
1351 //fprintf(stderr,"CURVATURE: z dimension is > 1.\n");
1352 pp.state = JBXL_GRAPH_IVDARG_ERROR;
1353 return pp;
1354 }
1355
1356 pp = make_VSGraph(vp.xs, vp.ys, vp.zs);
1357 if (pp.gp==NULL) return pp;
1358
1359 px = fxSobel(vp);
1360 py = fySobel(vp);
1361 pxy = fySobel(px);
1362 pxx = fxxSobel(vp);
1363 pyy = fyySobel(vp);
1364
1365 if (px.gp==NULL||py.gp==NULL||pxy.gp==NULL||pxx.gp==NULL||pyy.gp==NULL) {
1366 freeNull(px.gp);
1367 freeNull(py.gp);
1368 freeNull(pxy.gp);
1369 freeNull(pxx.gp);
1370 freeNull(pyy.gp);
1371 free_VSGraph(&pp);
1372 pp.state = JBXL_GRAPH_ERROR;
1373 return pp;
1374 }
1375
1376 for (i=0; i<vp.xs*vp.ys; i++) {
1377 Ix = px.gp[i];
1378 Iy = py.gp[i];
1379 Ixy = pxy.gp[i];
1380 Ixx = pxx.gp[i];
1381 Iyy = pyy.gp[i];
1382 d = 1. + Ix*Ix + Iy*Iy;
1383
1384 K = (Ixx*Iyy-Ixy*Ixy)/(d*d);
1385 H = (Ixx+Ixx*Iy*Iy+Iyy+Iyy*Ix*Ix-2*Ix*Ixy*Iy)/(2.*d*sqrt(d));
1386 pp.gp[i] = set_vector(K, H, 0.0);
1387 }
1388
1389 free(px.gp);
1390 free(py.gp);
1391 free(pxy.gp);
1392 free(pxx.gp);
1393 free(pyy.gp);
1394
1395 return pp;
1396}
1397
1398
1411WSGraph curv2WSGraph(VSGraph xp)
1412{
1413 int i;
1414 WSGraph vp;
1415 double K, H;
1416
1417 memset(&vp, 0, sizeof(WSGraph));
1418 if (xp.gp==NULL) {
1419 vp.state = JBXL_GRAPH_NODATA_ERROR;
1420 return vp;
1421 }
1422
1423 vp = make_WSGraph(xp.xs, xp.ys, xp.zs);
1424 if (vp.gp==NULL) return vp;
1425
1426 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
1427 K = xp.gp[i].x;
1428 H = xp.gp[i].y;
1429 if (K>0 && H<0) vp.gp[i] = PEAK;
1430 else if (K>0 && H>0) vp.gp[i] = PIT;
1431 else if (K<0 && H<0) vp.gp[i] = SADDLE_RIDGE;
1432 else if (K<0 && H>0) vp.gp[i] = SADDLE_VALLEY;
1433 else if (K>0 && H==0) vp.gp[i] = NONE_SHAPE;
1434 else if (K<0 && H==0) vp.gp[i] = MINIMAL;
1435 else if (K==0 && H<0) vp.gp[i] = RIDGE;
1436 else if (K==0 && H>0) vp.gp[i] = VALLEY;
1437 else if (K==0 && H==0) vp.gp[i] = FLAT;
1438 }
1439
1440 return vp;
1441}
1442
1443
1462WSGraph WSCurve(WSGraph gx, int mode, int cc)
1463{
1464 int i;
1465 FSGraph wp;
1466 VSGraph xp;
1467 WSGraph vp;
1468
1469 memset(&vp, 0, sizeof(WSGraph));
1470 if (gx.gp==NULL) {
1471 vp.state = JBXL_GRAPH_NODATA_ERROR;
1472 return vp;
1473 }
1474
1475 wp = W2FSGraph(gx);
1476 if (gx.zs<5) xp = curvature(wp);
1477 else xp = curvature3D(wp);
1478
1479 vp = curv2WSGraph(xp);
1480 freeNull(wp.gp);
1481 freeNull(xp.gp);
1482 if (vp.gp==NULL) return vp;
1483
1484 if (mode==ALL) {
1485 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
1486 if (vp.gp[i]==FLAT) vp.gp[i] = 500;
1487 else if (vp.gp[i]==PIT) vp.gp[i] = 1000;
1488 else if (vp.gp[i]==SADDLE_RIDGE) vp.gp[i] = 1500;
1489 else if (vp.gp[i]==SADDLE_VALLEY) vp.gp[i] = 2000;
1490 else if (vp.gp[i]==NONE_SHAPE) vp.gp[i] = 2500;
1491 else if (vp.gp[i]==MINIMAL) vp.gp[i] = 3000;
1492 else if (vp.gp[i]==RIDGE) vp.gp[i] = 3500;
1493 else if (vp.gp[i]==VALLEY) vp.gp[i] = 4000;
1494 else if (vp.gp[i]==PEAK) vp.gp[i] = 4500;
1495 else vp.gp[i] = 0;
1496 }
1497 }
1498 else {
1499 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
1500 if ((vp.gp[i]&mode)!=0) vp.gp[i] = cc;
1501 else vp.gp[i] = 0;
1502 }
1503 }
1504 return vp;
1505}
1506
1507
1521WSGraph edge_enhance(WSGraph gd, int mode)
1522{
1523 int i;
1524 WSGraph la, vp;
1525
1526 memset(&vp, 0, sizeof(WSGraph));
1527 if (gd.gp==NULL) {
1528 vp.state = JBXL_GRAPH_NODATA_ERROR;
1529 return vp;
1530 }
1531
1532 la = Laplacian(gd, mode);
1533 if (la.gp==NULL) return la;
1534
1535 vp = make_WSGraph(gd.xs, gd.ys, 1);
1536 if (vp.gp==NULL) return vp;
1537
1538 for (i=0; i<vp.xs*vp.ys; i++) vp.gp[i] = gd.gp[i] - la.gp[i];
1539 return vp;
1540}
1541
1542
1543
1545// Filter
1546
1558FMask gauss_mask(double sig, int ms, int md)
1559{
1560 int xx, yy, zz, ns, cp, dx, ps, sw;
1561 double min, *fm;
1562 FMask mask;
1563
1564 mask.mode = 0;
1565 mask.msize = 0;
1566 mask.imask = NULL;
1567
1568 if (md<=2) { // 2D
1569 md = 2;
1570 ps = ms*ms;
1571 sw = 0;
1572 }
1573 else { // 3D
1574 md = 3;
1575 ps = ms*ms*ms;
1576 sw = 1;
1577 }
1578
1579 ns = ms/2;
1580 min = (double)SINTMAX;
1581 fm = (double*)malloc(ps*sizeof(double));
1582 mask.imask = (int*)malloc(ps*sizeof(int));
1583 if (fm==NULL || mask.imask==NULL) {
1584 if (fm!=NULL) free(fm);
1585 if (mask.imask!=NULL) free(mask.imask);
1586 memset(&mask, 0, sizeof(FMask));
1587 return mask;
1588 }
1589 memset(fm, 0, ps*sizeof(double));
1590 memset(mask.imask, 0, ps*sizeof(int));
1591
1592 for (zz=-ns*sw; zz<=ns*sw; zz++) {
1593 for (yy=-ns; yy<=ns; yy++) {
1594 for (xx=-ns; xx<=ns; xx++) {
1595 cp = (zz+ns)*ms*ms*sw + (yy+ns)*ms + (xx+ns);
1596 fm[cp] = exp(-(xx*xx+yy*yy+zz*zz)/(sig*sig));
1597 if (fm[cp]!=0.0) min = Min(min, fm[cp]);
1598 }
1599 }
1600 }
1601
1602 dx = 0;
1603 for (xx=0; xx<ps; xx++) {
1604 mask.imask[xx] = (int)(fm[xx]/min+0.5);
1605 dx += mask.imask[xx];
1606 }
1607
1608 mask.msize = ms;
1609 mask.nfact = dx;
1610 mask.mode = md;
1611
1612 free(fm);
1613 return mask;
1614}
1615
1616
1627WSGraph imask(WSGraph xp, FMask mask)
1628{
1629 int i, x, y, z, cx;
1630 int xx, yy, zz, cp, cw, sw;
1631 int kc, xc, xs, ps, pm, mz, zc;
1632 int ms, nf, min;
1633 double dd;
1634 WSGraph vp;
1635
1636 memset(&vp, 0, sizeof(WSGraph));
1637 if (xp.gp==NULL) {
1638 vp.state = JBXL_GRAPH_NODATA_ERROR;
1639 return vp;
1640 }
1641
1642 if (xp.zs<=1 && mask.mode>2) {
1643 //fprintf(stderr, "IMASK: mismach mask dimension %d %d\n", xp.zs, mask.mode);
1644 vp.state = JBXL_GRAPH_IVDARG_ERROR;
1645 return vp;
1646 }
1647
1648 nf = mask.nfact;
1649 ms = mask.msize;
1650 if (mask.mode==2) {
1651 sw = 0;
1652 pm = ms*ms;
1653 }
1654 else {
1655 sw = 1;
1656 pm = ms*ms*ms;
1657 }
1658
1659 mz = Min(ms, xp.zs);
1660 kc = pm/2;
1661 zc = mz/2;
1662 xc = ms/2;
1663 xs = xp.xs;
1664 ps = xp.xs*xp.ys;
1665
1666 min = SINTMAX;
1667 for (i=0; i<xp.xs*xp.ys; i++) min = Min(min, xp.gp[i]);
1668 vp = make_WSGraph(xp.xs, xp.ys, 1);
1669 if (vp.gp==NULL) return vp;
1670
1671 for (i=0; i<vp.xs*vp.ys; i++) vp.gp[i] = min;
1672
1673 z = xp.zs/2;
1674 for (y=xc; y<xp.ys-xc; y++)
1675 for (x=xc; x<xp.xs-xc; x++) {
1676 cx = z*ps + y*xs + x;
1677 dd = 0.0;
1678 for (zz=-zc*sw; zz<=zc*sw; zz++)
1679 for (yy=-xc; yy<=xc; yy++)
1680 for (xx=-xc; xx<=xc; xx++) {
1681 cp = kc + xx + yy*ms + zz*ms*ms;
1682 cw = cx + xx + yy*xs + zz*ps;
1683 dd = dd + (double)xp.gp[cw]*mask.imask[cp];
1684 }
1685 vp.gp[y*xs + x] = (sWord)(dd/nf);
1686 }
1687 return vp;
1688}
1689
1690
1701WSGraph median(WSGraph xp, int ms) /* 3D */
1702{
1703 int i, j, x, y, z;
1704 int xx, yy, zz, cw, ux, mz;
1705 int kc, xc, zc, xs, ps, cx;
1706 WSGraph vp;
1707 sWord *me;
1708
1709 memset(&vp, 0, sizeof(WSGraph));
1710 if (xp.gp==NULL) {
1711 vp.state = JBXL_GRAPH_NODATA_ERROR;
1712 return vp;
1713 }
1714
1715 mz = Min(ms, xp.zs);
1716 me = (sWord*)malloc(ms*ms*mz*sizeof(sWord));
1717 if (me==NULL) {
1718 vp.state = JBXL_GRAPH_MEMORY_ERROR;
1719 return vp;
1720 }
1721 memset(me, 0, ms*ms*mz*sizeof(sWord));
1722
1723 kc = ms*ms*mz/2;
1724 xc = ms/2;
1725 zc = mz/2;
1726 xs = xp.xs;
1727 ps = xp.xs*xp.ys;
1728 vp = make_WSGraph(vp.xs, vp.ys, vp.zs);
1729 if (vp.gp==NULL) {
1730 free(me);
1731 return vp;
1732 }
1733
1734 z = xp.zs/2;
1735 for(y=xc; y<xp.ys-xc; y++) {
1736 for(x=xc; x<xp.xs-xc; x++) {
1737 cx = z*ps + y*xs + x;
1738 i = 0;
1739 for (zz=-zc; zz<=zc; zz++) {
1740 for (yy=-xc; yy<=xc; yy++) {
1741 for (xx=-xc; xx<=xc; xx++) {
1742 cw = cx + xx + yy*xs + zz*ps;
1743 me[i++] = xp.gp[cw];
1744 }
1745 }
1746 }
1747
1748 for (i=0; i<ms*ms*mz-1; i++) {
1749 for (j=i+1; j<ms*ms*mz; j++) {
1750 if (me[i]<me[j]) {
1751 ux = me[i];
1752 me[i] = me[j];
1753 me[j] = ux;
1754 }
1755 }
1756 }
1757 vp.gp[cx-z*ps] = me[kc];
1758 }
1759 }
1760
1761 free(me);
1762 return vp;
1763}
1764
1765
1781WSGraph to2d(WSGraph gd, int mode)
1782{
1783 int i, j, k, psize;
1784 int cx, cy, cz, cw, xx, yy;
1785 WSGraph vp;
1786
1787 memset(&vp, 0, sizeof(WSGraph));
1788 if (gd.gp==NULL) {
1789 vp.state = JBXL_GRAPH_NODATA_ERROR;
1790 return vp;
1791 }
1792
1793 psize = gd.xs*gd.ys;
1794
1795 if (mode==TOP_VIEW) {
1796 vp = make_WSGraph(gd.xs, gd.ys, 1);
1797 if (vp.gp==NULL) return vp;
1798
1799 for (k=0; k<gd.zs; k++) {
1800 cz = k*psize;
1801 for (j=0; j<gd.ys; j++) {
1802 cy = j*gd.xs;
1803 for (i=0; i<gd.xs; i++) {
1804 cx = cz + cy + i;
1805 cw = cy + i;
1806 if (gd.gp[cx]!=0) vp.gp[cw] = Max(vp.gp[cw], gd.gp[cx]);
1807 }
1808 }
1809 }
1810 }
1811
1812 else if (mode==TOP_VIEW_DEPTH) {
1813 vp = make_WSGraph(gd.xs, gd.ys, 1);
1814 if (vp.gp==NULL) return vp;
1815
1816 for (k=0; k<gd.zs; k++) {
1817 cz = k*psize;
1818 for (j=0; j<gd.ys; j++) {
1819 cy = j*gd.xs;
1820 for (i=0; i<gd.xs; i++) {
1821 cx = cz + cy + i;
1822 cw = cy + i;
1823 if (gd.gp[cx]!=0) vp.gp[cw] = Max(vp.gp[cw], (gd.zs-k)+100);
1824 }
1825 }
1826 }
1827 }
1828
1829 else if (mode==SIDEX_VIEW) {
1830 vp = make_WSGraph(gd.ys, gd.zs, 1);
1831 if (vp.gp==NULL) return vp;
1832
1833 for (k=0; k<gd.zs; k++) {
1834 cz = k*psize;
1835 yy = k;
1836 for (j=0; j<gd.ys; j++) {
1837 cy = j*gd.xs;
1838 xx = gd.ys - 1 - j;
1839 for (i=0; i<gd.xs; i++) {
1840 cx = cz + cy + i;
1841 cw = yy*vp.xs + xx;
1842 if (gd.gp[cx]!=0) vp.gp[cw] = Max(vp.gp[cw], gd.gp[cx]);
1843 }
1844 }
1845 }
1846 }
1847
1848 else if (mode==SIDEY_VIEW) {
1849 vp = make_WSGraph(gd.xs, gd.zs, 1);
1850 if (vp.gp==NULL) return vp;
1851
1852 for (k=0; k<gd.zs; k++) {
1853 cz = k*psize;
1854 yy = k;
1855 for (j=0; j<gd.ys; j++) {
1856 cy = j*gd.xs;
1857 for (i=0; i<gd.xs; i++) {
1858 cx = cz + cy + i;
1859 xx = i;
1860 cw = yy*vp.xs + xx;
1861 if (gd.gp[cx]!=0) vp.gp[cw] = Max(vp.gp[cw], gd.gp[cx]);
1862 }
1863 }
1864 }
1865 }
1866
1867 else {
1868 memset(&vp, 0, sizeof(WSGraph));
1869 vp.state = JBXL_GRAPH_IVDARG_ERROR;
1870 //fprintf(stderr,"TO2D: unknown mode = %d.\n",mode);
1871 }
1872
1873 return vp;
1874}
1875
1876
1888WSGraph euclid_distance(WSGraph vp, int* rr, int bc)
1889{
1890 int i, j, k, l;
1891 int df, db, d, w;
1892 int rmax, rstart, rend;
1893 WSGraph wp;
1894 ISGraph pp, buff;
1895
1896 memset(&wp, 0, sizeof(WSGraph));
1897 if (vp.gp==NULL) {
1898 wp.state = JBXL_GRAPH_NODATA_ERROR;
1899 return wp;
1900 }
1901
1902 wp = make_WSGraph(vp.xs, vp.ys, vp.zs);
1903 if (wp.gp==NULL) return wp;
1904
1905 pp = make_ISGraph(vp.xs, vp.ys, vp.zs);
1906 if (pp.gp==NULL) {
1907 free_WSGraph(&wp);
1908 wp.state = pp.state;
1909 return wp;
1910 }
1911
1912 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
1913 if (vp.gp[i]>=bc) pp.gp[i] = 1;
1914 else pp.gp[i] = 0;
1915 }
1916
1917 for (k=1; k<=vp.zs; k++) {
1918 for (j=1; j<=vp.ys; j++) {
1919 df = vp.xs;
1920 for (i=1; i<=vp.xs; i++) {
1921 if (Vxt(pp, i, j, k)!=0) df = df + 1;
1922 else df = 0;
1923 Vxt(pp, i, j, k) = df*df;
1924 }
1925 }
1926 }
1927
1928 for (k=1; k<=vp.zs; k++) {
1929 for (j=1; j<=vp.ys; j++) {
1930 db = vp.xs;
1931 for (i=vp.xs; i>=1; i--) {
1932 if (Vxt(pp, i, j, k)!=0) db = db + 1;
1933 else db = 0;
1934 Vxt(pp, i, j, k) = Min(Vxt(pp, i, j, k), db*db);
1935 }
1936 }
1937 }
1938
1939 buff = make_ISGraph(vp.ys, 1, 1);
1940 for (k=1; k<=vp.zs; k++) {
1941 for (i=1; i<=vp.xs; i++) {
1942 for (j=1; j<=vp.ys; j++) {
1943 Lxt(buff, j) = Vxt(pp, i, j, k);
1944 }
1945 for (j=1; j<=vp.ys; j++) {
1946 d = Lxt(buff, j);
1947 if (d!=0) {
1948 rmax = (int)sqrt((double)d) + 1;
1949 rstart = Min(rmax, j-1);
1950 rend = Min(rmax, vp.ys-j);
1951 for (l=-rstart; l<=rend; l++) {
1952 w = Lxt(buff, j+l) + l*l;
1953 if (w<d) d = w;
1954 }
1955 }
1956 Vxt(pp, i, j, k) = d;
1957 }
1958 }
1959 }
1960 free(buff.gp);
1961
1962 *rr = 0;
1963 buff = make_ISGraph(vp.zs, 1, 1);
1964 for (j=1; j<=vp.ys; j++) {
1965 for (i=1; i<=vp.xs; i++) {
1966 for (k=1; k<=vp.zs; k++) {
1967 Lxt(buff, k) = Vxt(pp, i, j, k);
1968 }
1969 for (k=1; k<=vp.zs; k++) {
1970 d = Lxt(buff, k);
1971 if (d!=0) {
1972 rmax = (int)sqrt((double)d) + 1;
1973 rstart = Min(rmax, k-1);
1974 rend = Min(rmax, vp.zs-k);
1975 for (l=-rstart; l<=rend; l++) {
1976 w = Lxt(buff, k+l) + l*l;
1977 if (w<d) d = w;
1978 }
1979 *rr = Max(*rr, d);
1980 }
1981 Vxt(pp, i, j, k) = d;
1982 }
1983 }
1984 }
1985 free(buff.gp);
1986
1987 for (i=0; i<wp.xs*wp.ys*wp.zs; i++) {
1988 wp.gp[i] = (sWord)pp.gp[i];
1989 if (pp.gp[i]>32767) {
1990 fprintf(stderr,"EUCLID_DISTANCE: WARNING: distance is too long = %d!\n",pp.gp[i]);
1991 }
1992 }
1993 free(pp.gp);
1994
1995 return wp;
1996}
1997
1998
2021int out_round(WSGraph vp, int x, int y, IRBound* rb, int mode)
2022{
2023 int i, j, sp, cp, w, ll, ss;
2024 int xx, yy, vx, vy, ix, eflg=OFF;
2025 int r8[8]={-1, 1, -1, -1, 1, -1, 1, 1};
2026 int r4[8]={ 0, 1, -1, 0, 0, -1, 1, 0};
2027 int* cc;
2028
2029 if (vp.gp==NULL) return JBXL_GRAPH_IVDARG_ERROR;
2030
2031 i = y*vp.xs + x;
2032 rb->xmax = rb->xmin = x;
2033 rb->ymax = rb->ymin = y;
2034 sp = cp = i;
2035 ss = 0;
2036 ll = 0;
2037 vx = 1;
2038 vy = 0;
2039
2040 if (vp.gp[sp]==0 || sp==0) {
2041 //fprintf(stderr,"OUT_ROUND: irregular start point!! sp = %d\n",sp);
2042 return JBXL_GRAPH_IVDPARAM_ERROR;
2043 }
2044
2045 if (mode==8){
2046 ix = 8;
2047 cc = r8;
2048 }
2049 else if (mode==4) {
2050 ix = 4;
2051 cc = r4;
2052 }
2053 else {
2054 //fprintf(stderr,"OUT_ROUND: invalid mode = %d!!\n",mode);
2055 return JBXL_GRAPH_IVDMODE_ERROR;
2056 }
2057
2058 do {
2059 w = abs(vx)+abs(vy);
2060 xx = (vx*cc[0]+vy*cc[1])/w;
2061 yy = (vx*cc[2]+vy*cc[3])/w;
2062 for (j=1; j<=ix; j++) {
2063 if (vp.gp[cp+yy*vp.xs+xx]!=0) {
2064 vx = xx;
2065 vy = yy;
2066 cp = cp + yy*vp.xs + xx;
2067 xx = cp%vp.xs;
2068 yy = (cp-xx)/vp.xs;
2069 rb->xmax = Max(rb->xmax, xx);
2070 rb->ymax = Max(rb->ymax, yy);
2071 rb->xmin = Min(rb->xmin, xx);
2072 rb->ymin = Min(rb->ymin, yy);
2073 break;
2074 }
2075 else {
2076 if(sp==cp && xx==-1 && yy==0) {
2077 eflg = ON;
2078 break;
2079 }
2080 w = abs(xx)+abs(yy);
2081 vx = (xx*cc[4]+yy*cc[5])/w;
2082 vy = (xx*cc[6]+yy*cc[7])/w;
2083 xx = vx;
2084 yy = vy;
2085 }
2086 }
2087 ll++;
2088 if (abs(vx)+abs(vy)==2) ss++;
2089 //
2090 } while(eflg==OFF);
2091
2092 if (mode==4) ss = 0;
2093 (rb->xmax)++;
2094 (rb->ymax)++;
2095 (rb->xmin)--;
2096 (rb->ymin)--;
2097 rb->misc = ss;
2098
2099 return ll;
2100}
2101
Min
#define Min(x, y)
Definition common.h:250
OFF
#define OFF
Definition common.h:231
Max
#define Max(x, y)
Definition common.h:247
sWord
short sWord
2Byte
Definition common.h:335
SINTMAX
#define SINTMAX
Definition common.h:204
ON
#define ON
Definition common.h:230
make_WSGraph
WSGraph make_WSGraph(int xs, int ys, int zs)
Definition gdata.c:108
W2FSGraph
FSGraph W2FSGraph(WSGraph vp)
Definition gdata.c:618
make_VSGraph
VSGraph make_VSGraph(int xs, int ys, int zs)
Definition gdata.c:234
make_FSGraph
FSGraph make_FSGraph(int xs, int ys, int zs)
Definition gdata.c:150
make_ISGraph
ISGraph make_ISGraph(int xs, int ys, int zs)
Definition gdata.c:192
free_WSGraph
#define free_WSGraph(v)
Definition gdata.h:176
free_FSGraph
#define free_FSGraph(v)
Definition gdata.h:178
free_VSGraph
#define free_VSGraph(v)
Definition gdata.h:180
Lxt
#define Lxt(v, i)
座標を 1 から数える．
Definition gdata.h:199
Px
#define Px(v, i, j)
2次元画像データ vの (i, j) のデータを参照する．
Definition gdata.h:192
Vx
#define Vx(v, i, j, k)
3次元画像データ vの (i, j, k) のデータを参照する．
Definition gdata.h:193
Vxt
#define Vxt(v, i, j, k)
Definition gdata.h:201
yySobel
WSGraph yySobel(WSGraph vp)
Definition gmt.c:569
curv2WSGraph
WSGraph curv2WSGraph(VSGraph xp)
Definition gmt.c:1411
fzzSobel
FSGraph fzzSobel(FSGraph vp)
Definition gmt.c:829
ySobel
WSGraph ySobel(WSGraph vp)
Definition gmt.c:188
vNabra
VSGraph vNabra(WSGraph vp)
Definition gmt.c:935
curvature
VSGraph curvature(FSGraph vp)
Definition gmt.c:1336
xSobel
WSGraph xSobel(WSGraph vp)
Definition gmt.c:90
fySobel
FSGraph fySobel(FSGraph vp)
Definition gmt.c:237
fyySobel
FSGraph fyySobel(FSGraph vp)
Definition gmt.c:646
fxxSobel
FSGraph fxxSobel(FSGraph vp)
Definition gmt.c:492
gauss_mask
FMask gauss_mask(double sig, int ms, int md)
Definition gmt.c:1558
zSobel
WSGraph zSobel(WSGraph vp)
Definition gmt.c:285
fxSobel
FSGraph fxSobel(FSGraph vp)
Definition gmt.c:140
euclid_distance
WSGraph euclid_distance(WSGraph vp, int *rr, int bc)
Definition gmt.c:1888
median
WSGraph median(WSGraph xp, int ms)
Definition gmt.c:1701
fNabra
FSGraph fNabra(FSGraph vp)
Definition gmt.c:1162
to2d
WSGraph to2d(WSGraph gd, int mode)
Definition gmt.c:1781
out_round
int out_round(WSGraph vp, int x, int y, IRBound *rb, int mode)
Definition gmt.c:2021
Nabra
WSGraph Nabra(WSGraph vp)
Definition gmt.c:1088
fzSobel
FSGraph fzSobel(FSGraph vp)
Definition gmt.c:350
xxSobel
WSGraph xxSobel(WSGraph vp)
Definition gmt.c:414
curvature3D
VSGraph curvature3D(FSGraph vp)
Definition gmt.c:1235
vfNabra
VSGraph vfNabra(FSGraph vp)
Definition gmt.c:1012
zzSobel
WSGraph zzSobel(WSGraph vp)
Definition gmt.c:722
Laplacian
WSGraph Laplacian(WSGraph vp, int mode)
Definition gmt.c:24
imask
WSGraph imask(WSGraph xp, FMask mask)
Definition gmt.c:1627
WSCurve
WSGraph WSCurve(WSGraph gx, int mode, int cc)
Definition gmt.c:1462
edge_enhance
WSGraph edge_enhance(WSGraph gd, int mode)
Definition gmt.c:1521
gmt.h
グラフィック用数学ライブラリ＆フィルタ ヘッダ
SADDLE_RIDGE
#define SADDLE_RIDGE
Definition gmt.h:28
ALL
#define ALL
Definition gmt.h:24
SIDEY_VIEW
#define SIDEY_VIEW
Definition gmt.h:38
MINIMAL
#define MINIMAL
Definition gmt.h:30
RIDGE
#define RIDGE
Definition gmt.h:31
VALLEY
#define VALLEY
Definition gmt.h:32
TOP_VIEW
#define TOP_VIEW
Definition gmt.h:35
FLAT
#define FLAT
Definition gmt.h:33
SADDLE_VALLEY
#define SADDLE_VALLEY
Definition gmt.h:29
PEAK
#define PEAK
Definition gmt.h:26
NONE_SHAPE
#define NONE_SHAPE
Definition gmt.h:25
TOP_VIEW_DEPTH
#define TOP_VIEW_DEPTH
Definition gmt.h:39
PIT
#define PIT
Definition gmt.h:27
SIDEX_VIEW
#define SIDEX_VIEW
Definition gmt.h:37
jbxl_state.h
JunkBox_Lib 状態ヘッダ
JBXL_GRAPH_IVDARG_ERROR
#define JBXL_GRAPH_IVDARG_ERROR
無効な引数
Definition jbxl_state.h:178
JBXL_GRAPH_IVDPARAM_ERROR
#define JBXL_GRAPH_IVDPARAM_ERROR
無効なパラメータ
Definition jbxl_state.h:180
JBXL_GRAPH_NODATA_ERROR
#define JBXL_GRAPH_NODATA_ERROR
データが無い
Definition jbxl_state.h:171
JBXL_GRAPH_ERROR
#define JBXL_GRAPH_ERROR
GRAPH ライブラリーのエラー
Definition jbxl_state.h:167
JBXL_NORMAL
#define JBXL_NORMAL
正常
Definition jbxl_state.h:32
JBXL_GRAPH_IVDMODE_ERROR
#define JBXL_GRAPH_IVDMODE_ERROR
無効なモード
Definition jbxl_state.h:179
JBXL_GRAPH_MEMORY_ERROR
#define JBXL_GRAPH_MEMORY_ERROR
メモリエラー
Definition jbxl_state.h:170
unit_vector
vector unit_vector(vector a)
Definition matrix.c:18
set_vector
vector set_vector(double x, double y, double z)
Definition matrix.c:82
FMask
Definition gmt.h:16
FMask::msize
int msize
size of mask
Definition gmt.h:18
FMask::mode
int mode
2:2D, 3:3D
Definition gmt.h:17
FMask::imask
int * imask
pointer to mask
Definition gmt.h:20
FMask::nfact
int nfact
normalized facter
Definition gmt.h:19
FSGraph::zs
int zs
zサイズ. 4Byte. 2Dの場合は 1.
Definition gdata.h:73
FSGraph::xs
int xs
xサイズ. 4Byte.
Definition gdata.h:71
FSGraph::state
int state
状態
Definition gdata.h:74
FSGraph::gp
double * gp
グラフィックデータへのポインタ. xs*ys*zs*sizeof(double)
Definition gdata.h:75
FSGraph::ys
int ys
yサイズ. 4Byte.
Definition gdata.h:72
IRBound::xmax
int xmax
x軸境界の最大値．
Definition gdata.h:115
IRBound::misc
int misc
多目的用．
Definition gdata.h:120
IRBound::ymax
int ymax
y軸境界の最大値．
Definition gdata.h:117
IRBound::xmin
int xmin
x軸境界の最小値．
Definition gdata.h:114
IRBound::ymin
int ymin
y軸境界の最小値．
Definition gdata.h:116
ISGraph::state
int state
状態
Definition gdata.h:60
ISGraph::gp
int * gp
グラフィックデータへのポインタ. xs*ys*zs*4Byte
Definition gdata.h:61
VSGraph::zs
int zs
zサイズ. 4Byte. 2Dの場合は 1.
Definition gdata.h:87
VSGraph::xs
int xs
xサイズ. 4Byte.
Definition gdata.h:85
VSGraph::state
int state
状態
Definition gdata.h:88
VSGraph::gp
vector * gp
グラフィックデータへのポインタ. xs*ys*zs*sizeof(vector).
Definition gdata.h:89
VSGraph::ys
int ys
yサイズ. 4Byte.
Definition gdata.h:86
WSGraph::zs
int zs
zサイズ. 4Byte. 2Dの場合は 1.
Definition gdata.h:45
WSGraph::xs
int xs
xサイズ. 4Byte.
Definition gdata.h:43
WSGraph::state
int state
状態
Definition gdata.h:46
WSGraph::gp
sWord * gp
グラフィックデータへのポインタ. xs*ys*zs*2Byte.
Definition gdata.h:47
WSGraph::ys
int ys
yサイズ. 4Byte.
Definition gdata.h:44
vector::y
double y
y方向成分
Definition matrix.h:31
vector::x
double x
x方向成分
Definition matrix.h:30
freeNull
#define freeNull(p)
Definition tools.h:201