JunkBox_Lib  1.10.2
gmt.c
Go to the documentation of this file.
1 
8 #include "gmt.h"
9 #include "jbxl_state.h"
10 
11 
24 WSGraph 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 
90 WSGraph 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 
140 FSGraph 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 
188 WSGraph 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 
237 FSGraph 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 
285 WSGraph 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 
350 FSGraph 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 
414 WSGraph 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 
492 FSGraph 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 
569 WSGraph 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 
646 FSGraph 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 
722 WSGraph 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 
829 FSGraph 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 
935 VSGraph 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 
1012 VSGraph 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 
1088 WSGraph 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 
1162 FSGraph 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 
1235 VSGraph 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 
1336 VSGraph 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 
1411 WSGraph 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 
1462 WSGraph 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 
1521 WSGraph 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 
1558 FMask 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 
1627 WSGraph 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 
1701 WSGraph 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 
1781 WSGraph 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 
1888 WSGraph 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 
2021 int 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
Definition: gdata.h:70
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
Definition: gdata.h:113
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
Definition: gdata.h:56
ISGraph::state
int state
状態
Definition: gdata.h:60
ISGraph::gp
int * gp
グラフィックデータへのポインタ. xs*ys*zs*4Byte
Definition: gdata.h:61
VSGraph
Definition: gdata.h:84
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
Definition: gdata.h:42
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