JunkBox_Lib++ (for Windows) 1.10.1
Loading...
Searching...
No Matches
Gmt.h
Go to the documentation of this file.
1#ifndef __JBXL_CPP_GRAPHIC_MATH_H_
2#define __JBXL_CPP_GRAPHIC_MATH_H_
3
10#include "Gdata.h"
11
12
13//
14namespace jbxl {
15
16
17/*
18point() は遅いので使用しない．
19
20template <typename R, typename T> MSGraph<R> Laplacian(MSGraph<T> vp, int mode=0);
21
22template <typename R, typename T> MSGraph<R> xSobel(MSGraph<T> vp);
23template <typename R, typename T> MSGraph<R> ySobel(MSGraph<T> vp);
24template <typename R, typename T> MSGraph<R> zSobel(MSGraph<T> vp);
25template <typename R, typename T> MSGraph<R> xxSobel(MSGraph<T> vp);
26template <typename R, typename T> MSGraph<R> yySobel(MSGraph<T> vp);
27template <typename R, typename T> MSGraph<R> zzSobel(MSGraph<T> vp);
28
29template <typename R, typename T> MSGraph<Vector<R> > vNabla(MSGraph<T> vp);
30template <typename R, typename T> MSGraph<R> Nabla(MSGraph<T> vp);
31
32template <typename R, typename T> MSGraph<R> edgeEnhance(MSGraph<T> gd, int mode=0);
33template <typename T> MSGraph<T> medianFilter(MSGraph<T> xp, int ms=3);
34template <typename T> MSGraph<int> euclidDistance(MSGraph<T> vp, int bc, int& rr);
35
36template <typename R, typename T> MSGraph<R> MSMaskFilter(MSGraph<R> vp, MSGraph<T> filter, int abs=FALSE)
37*/
38
39
40
42// テンプレート定義
43
58template <typename R, typename T> MSGraph<R> Laplacian(MSGraph<T> vp, int mode=0)
59{
60 int i, j;
61 int nx, ny, xs, xs2;
62 R da, db, dc, dd, de, df, dg, dh;
63 MSGraph<R> lp;
64
65 lp.mimicry(vp);
66 lp.base = lp.zero = 0;
67
68 if (lp.isNull()) {
69 DEBUG_MODE PRINT_MESG("LAPLACIAN: No More Memory!!!\n");
70 lp.state = JBXL_GRAPH_MEMORY_ERROR;
71 return lp;
72 }
73
74 xs = vp.xs;
75 xs2 = 2*xs;
76
77 if (mode==4) {
78 for (j=1; j<vp.ys-1; j++) {
79 ny = j*vp.xs;
80 for (i=1; i<vp.xs-1; i++) {
81 nx = ny + i;
82 da = vp.gp[nx+1] + vp.gp[nx-1];
83 db = vp.gp[nx];
84 dc = vp.gp[nx+xs] + vp.gp[nx-xs];
85 lp.gp[nx] = (R)(da - 4.*db + dc);
86 //da = vp.point(i+1, j) + vp.point(i-1, j);
87 //db = vp.point(i, j);
88 //dc = vp.point(i, j+1) + vp.point(i, j-1);
89 //lp.point(i, j) = (R)(da - 4.*db + dc);
90 }
91 }
92 }
93
94 else if (mode==8) {
95 for (j=1; j<vp.ys-1; j++) {
96 ny = j*vp.xs;
97 for (i=1; i<vp.xs-1; i++) {
98 nx = ny + i;
99 da = vp.gp[nx+1] + vp.gp[nx-1];
100 db = vp.gp[nx+xs] + vp.gp[nx-xs];
101 dc = vp.gp[nx];
102 dd = vp.gp[nx+1+xs] + vp.gp[nx-1+xs];
103 de = vp.gp[nx+1-xs] + vp.gp[nx-1-xs];
104 lp.gp[nx] = (R)(da + db - 8.*dc + dd + de);
105 //da = vp.point(i+1, j) + vp.point(i-1, j);
106 //db = vp.point(i, j+1) + vp.point(i, j-1);
107 //dc = vp.point(i, j);
108 //dd = vp.point(i+1, j+1) + vp.point(i-1, j+1);
109 //de = vp.point(i+1, j-1) + vp.point(i-1, j-1);
110 //lp.point(i, j) = (R)(da + db - 8.*dc + dd + de);
111 }
112 }
113 }
114
115 else {
116 for (j=2; j<vp.ys-2; j++) {
117 ny = j*vp.xs;
118 for (i=2; i<vp.xs-2; i++) {
119 nx = ny + i;
120 da = vp.gp[nx];
121 db = vp.gp[nx+1] + vp.gp[nx-1] + vp.gp[nx+xs] + vp.gp[nx-xs];
122 dc = vp.gp[nx-1-xs2] + vp.gp[nx-xs2] + vp.gp[nx+1-xs2];
123 dd = vp.gp[nx-1+xs2] + vp.gp[nx+xs2] + vp.gp[nx+1+xs2];
124 de = vp.gp[nx-2-xs ] + vp.gp[nx-2] + vp.gp[nx-2+xs];
125 df = vp.gp[nx+2-xs ] + vp.gp[nx+2] + vp.gp[nx+2+xs];
126 dg = vp.gp[nx-2-xs2] + vp.gp[nx+2-xs2];
127 dh = vp.gp[nx-2+xs2] + vp.gp[nx+2+xs2];
128 lp.gp[nx] = (R)((-12.*da - 4.*db + 2.*(dc+dd+de+df) + dg + dh)/32.);
129
130 //da = vp.point(i, j);
131 //db = vp.point(i+1, j) + vp.point(i-1, j) + vp.point(i, j+1) + vp.point(i, j-1);
132 //dc = vp.point(i-1, j-2) + vp.point(i, j-2) + vp.point(i+1, j-2);
133 //dd = vp.point(i-1, j+2) + vp.point(i, j+2) + vp.point(i+1, j+2);
134 //de = vp.point(i-2, j-1) + vp.point(i-2, j) + vp.point(i-2, j+1);
135 //df = vp.point(i+2, j-1) + vp.point(i+2, j) + vp.point(i+2, j+1);
136 //dg = vp.point(i-2, j-2) + vp.point(i+2, j-2);
137 //dh = vp.point(i-2, j+2) + vp.point(i+2, j+2);
138 //lp.point(i, j) = (R)((-12.*da - 4.*db + 2.*(dc+dd+de+df) + dg + dh)/32.);
139 }
140 }
141 }
142 return lp;
143}
144
145
154template <typename R, typename T> MSGraph<R> xSobel(MSGraph<T> vp)
155{
156 int i, j, k;
157 int pl, nx, ny, nz;
158 R da, db, dc, dd, de, nr;
159 MSGraph<R> xp;
160
161 xp.mimicry(vp);
162 xp.base = xp.zero = 0;
163
164 if (xp.isNull()) {
165 DEBUG_MODE PRINT_MESG("XSOBEL: No More Memory!!!\n");
166 xp.state = JBXL_GRAPH_MEMORY_ERROR;
167 return xp;
168 }
169
170 // カウンタ
171 CVCounter* vcounter = NULL;
172 if (vp.zs>2) vcounter = GetUsableGlobalCounter();
173 if (vcounter!=NULL) vcounter->SetMax(vp.zs);
174
175 pl = vp.xs*vp.ys;
176 for (k=0; k<vp.zs; k++) {
177 nz = k*pl;
178 for (j=1; j<vp.ys-1; j++) {
179 ny = nz + j*vp.xs;
180 for (i=1; i<vp.xs-1; i++) {
181 nx = ny + i;
182 da = vp.gp[nx+1-vp.xs] - vp.gp[nx-1-vp.xs]; // 1/4 * (da+dc)
183 db = vp.gp[nx+1] - vp.gp[nx-1]; // 1/2 * db
184 dc = vp.gp[nx+1+vp.xs] - vp.gp[nx-1+vp.xs];
185 //da = vp.point(i+1, j-1, k) - vp.point(i-1, j-1, k); // 1/4 * (da+dc)
186 //db = vp.point(i+1, j, k) - vp.point(i-1, j, k); // 1/2 * db
187 //dc = vp.point(i+1, j+1, k) - vp.point(i-1, j+1, k);
188
189 if (k==0 || k==vp.zs-1) {
190 dd = de = 0;
191 nr = 8;
192 }
193 else {
194 dd = vp.gp[nx+1-pl] - vp.gp[nx-1-pl]; // 1/4 * (dd+de)
195 de = vp.gp[nx+1+pl] - vp.gp[nx-1+pl];
196 //dd = vp.point(i+1, j, k-1) - vp.point(i-1, j, k-1); // 1/4 * (dd+de)
197 //de = vp.point(i+1, j, k+1) - vp.point(i-1, j, k+1);
198 nr = 12;
199 }
200
201 xp.gp[nx] = (R)((da + 2.*db + dc + dd + de)/nr);
202 }
203 }
204 if (vcounter!=NULL) vcounter->StepIt();
205 }
206
207 return xp;
208}
209
210
219template<typename R, typename T> MSGraph<R> ySobel(MSGraph<T> vp)
220{
221 int i, j, k;
222 int pl, nx, ny, nz;
223 R da, db, dc, dd, de, nr;
224 MSGraph<R> xp;
225
226 xp.mimicry(vp);
227 xp.base = xp.zero = 0;
228
229 if (xp.isNull()) {
230 DEBUG_MODE PRINT_MESG("YSOBEL: No More Memory!!\n");
231 xp.state = JBXL_GRAPH_MEMORY_ERROR;
232 return xp;
233 }
234
235 // カウンタ
236 CVCounter* vcounter = NULL;
237 if (vp.zs>2) vcounter = GetUsableGlobalCounter();
238 if (vcounter!=NULL) vcounter->SetMax(vp.zs);
239
240 pl = vp.xs*vp.ys;
241 for (k=0; k<vp.zs; k++) {
242 nz = k*pl;
243 for (j=1; j<vp.ys-1; j++) {
244 ny = nz + j*vp.xs;
245 for (i=1; i<vp.xs-1; i++) {
246 nx = ny + i;
247 da = vp.gp[nx-1+vp.xs] - vp.gp[nx-1-vp.xs];
248 db = vp.gp[nx +vp.xs] - vp.gp[nx -vp.xs];
249 dc = vp.gp[nx+1+vp.xs] - vp.gp[nx+1-vp.xs];
250 //da = vp.point(i-1, j+1, k) - vp.point(i-1, j-1, k);
251 //db = vp.point(i , j+1, k) - vp.point(i, j-1, k);
252 //dc = vp.point(i+1, j+1, k) - vp.point(i+1, j-1, k);
253
254 if (k==0 || k==vp.zs-1) {
255 dd = de = 0;
256 nr = 8;
257 }
258 else {
259 dd = vp.gp[nx+vp.xs-pl] - vp.gp[nx-vp.xs-pl];
260 de = vp.gp[nx+vp.xs+pl] - vp.gp[nx-vp.xs+pl];
261 //dd = vp.point(i, j+1, k-1) - vp.point(i, j-1, k-1);
262 //de = vp.point(i, j+1, k+1) - vp.point(i, j-1, k+1);
263 nr = 12;
264 }
265
266 xp.gp[nx] = (R)((da + 2.*db + dc + dd + de)/nr);
267 }
268 }
269 if (vcounter!=NULL) vcounter->StepIt();
270 }
271
272 return xp;
273}
274
275
284template<typename R, typename T> MSGraph<R> zSobel(MSGraph<T> vp)
285{
286 int i, j, k;
287 int pl, nx, ny, nz;
288 R da, db, dc, dd, de;
289 MSGraph<R> xp;
290
291 xp.mimicry(vp);
292 xp.base = xp.zero = 0;
293
294 if (xp.isNull()) {
295 DEBUG_MODE PRINT_MESG("ZSOBEL: No More Memory!!\n");
296 xp.state = JBXL_GRAPH_MEMORY_ERROR;
297 return xp;
298 }
299 if (vp.zs<2) {
300 xp.state = JBXL_GRAPH_NODATA_ERROR;
301 return xp; // 0 のグラフィックデータを返す
302 }
303
304 // カウンタ
305 CVCounter* vcounter = NULL;
306 if (vp.zs>2) vcounter = GetUsableGlobalCounter();
307 if (vcounter!=NULL) vcounter->SetMax(vp.zs-1);
308
309 pl = vp.xs*vp.ys;
310 for (k=1; k<vp.zs-1; k++) {
311 nz = k*pl;
312 for (j=1; j<vp.ys-1; j++) {
313 ny = nz + j*vp.xs;
314 for (i=1; i<vp.xs-1; i++) {
315 nx = ny +i;
316 da = vp.gp[nx-1+pl] - vp.gp[nx-1-pl];
317 db = vp.gp[nx+1+pl] - vp.gp[nx+1-pl];
318 dc = vp.gp[nx +pl] - vp.gp[nx -pl];
319 dd = vp.gp[nx-vp.xs+pl] - vp.gp[nx-vp.xs-pl];
320 de = vp.gp[nx+vp.xs+pl] - vp.gp[nx+vp.xs-pl];
321 //da = vp.point(i-1, j, k+1) - vp.point(i-1, j, k-1);
322 //db = vp.point(i+1, j, k+1) - vp.point(i+1, j, k-1);
323 //dc = vp.point(i, j, k+1) - vp.point(i, j, k-1);
324 //dd = vp.point(i, j-1, k+1) - vp.point(i, j-1, k-1);
325 //de = vp.point(i, j+1, k+1) - vp.point(i, j+1, k-1);
326 xp.gp[nx] = (R)((da + db + 2.*dc + dd + de)/12.);
327 }
328 }
329 if (vcounter!=NULL) vcounter->StepIt();
330 }
331
332 // k==0 and k==vp.zs-1
333 nz = (vp.zs-1)*pl;
334 for (j=1; j<vp.ys-1; j++) {
335 ny = j*vp.xs;
336 for (i=1; i<vp.xs-1; i++) {
337 nx = ny + i;
338 da = vp.gp[nx];
339 db = vp.gp[nx+pl];
340 dc = vp.gp[nx+1 +pl] + vp.gp[nx-1 +pl];
341 dd = vp.gp[nx+vp.xs+pl] + vp.gp[nx-vp.xs+pl];
342 //da = vp.point(i, j, 0);
343 //db = vp.point(i, j, 1);
344 //dc = vp.point(i+1, j, 1) + vp.point(i-1, j, 1);
345 //dd = vp.point(i, j+1, 1) + vp.point(i, j-1, 1);
346 xp.gp[nx] = (R)((2.*db + dc + dd)/6. - da);
347 }
348
349 ny = ny + nz;
350 for (i=1; i<vp.xs-1; i++) {
351 nx = ny + i;
352 da = vp.gp[nx];
353 db = vp.gp[nx-pl];
354 dc = vp.gp[nx+1 -pl] + vp.gp[nx-1 -pl];
355 dd = vp.gp[nx+vp.xs-pl] + vp.gp[nx-vp.xs-pl];
356 //da = vp.point(i, j, vp.zs-1);
357 //db = vp.point(i, j, vp.zs-2);
358 //dc = vp.point(i+1, j, vp.zs-2) + vp.point(i-1, j, vp.zs-2);
359 //dd = vp.point(i, j+1, vp.zs-2) + vp.point(i, j-1, vp.zs-2);
360 xp.gp[nx] = (R)(da - (2.*db + dc + dd)/6.);
361 }
362 }
363 if (vcounter!=NULL) vcounter->PutFill();
364
365 return xp;
366}
367
368
377template<typename R, typename T> MSGraph<R> xxSobel(MSGraph<T> vp)
378{
379 int i, j, k;
380 int pl, nx, ny, nz, pl2, xs, xs2;
381 R da, db, dc, dd, de;
382 R df, dg, dh, di, dj, dk, dl, dm, nr;
383 MSGraph<R> xp;
384
385 xp.mimicry(vp);
386 xp.base = xp.zero = 0;
387
388 if (xp.isNull()) {
389 DEBUG_MODE PRINT_MESG("XXSOBEL: No More Memory!!\n");
390 xp.state = JBXL_GRAPH_MEMORY_ERROR;
391 return xp;
392 }
393
394 // カウンタ
395 CVCounter* vcounter = NULL;
396 if (vp.zs>2) vcounter = GetUsableGlobalCounter();
397 if (vcounter!=NULL) vcounter->SetMax(vp.zs);
398
399 pl = vp.xs*vp.ys;
400 pl2 = 2*pl;
401 xs = vp.xs;
402 xs2 = 2*vp.xs;
403
404 for (k=0; k<vp.zs; k++) {
405 nz = k*pl;
406 for (j=2; j<vp.ys-2; j++) {
407 ny = nz + j*vp.xs;
408 for (i=2; i<vp.xs-2; i++) {
409 nx = ny + i;
410 da = vp.gp[nx+2-xs2] - 2*vp.gp[nx-xs2] + vp.gp[nx-2-xs2];
411 db = vp.gp[nx+2-xs ] - 2*vp.gp[nx-xs] + vp.gp[nx-2-xs];
412 dc = vp.gp[nx+2] - 2*vp.gp[nx] + vp.gp[nx-2];
413 dd = vp.gp[nx+2+xs] - 2*vp.gp[nx+xs] + vp.gp[nx-2+xs];
414 de = vp.gp[nx+2+xs2] - 2*vp.gp[nx+xs2] + vp.gp[nx-2+xs2];
415 //da = vp.point(i+2, j-2, k) - 2*vp.point(i, j-2, k) + vp.point(i-2, j-2, k);
416 //db = vp.point(i+2, j-1, k) - 2*vp.point(i, j-1, k) + vp.point(i-2, j-1, k);
417 //dc = vp.point(i+2, j, k) - 2*vp.point(i, j, k) + vp.point(i-2, j, k);
418 //dd = vp.point(i+2, j+1, k) - 2*vp.point(i, j+1, k) + vp.point(i-2, j+1, k);
419 //de = vp.point(i+2, j+2, k) - 2*vp.point(i, j+2, k) + vp.point(i-2, j+2, k);
420
421 if (k==0 || k==vp.zs-1) {
422 dc = (R)(6.*dc);
423 df = dg = dh = di = dj = dk = dl = dm = 0;
424 nr = 64;
425 }
426 else {
427 dc = (R)(8.*dc);
428 df = vp.gp[nx+2-xs-pl] - 2*vp.gp[nx-xs-pl] + vp.gp[nx-2-xs-pl];
429 dg = vp.gp[nx+2 -pl] - 2*vp.gp[nx -pl] + vp.gp[nx-2 -pl];
430 dh = vp.gp[nx+2+xs-pl] - 2*vp.gp[nx+xs-pl] + vp.gp[nx-2+xs-pl];
431 di = vp.gp[nx+2-xs+pl] - 2*vp.gp[nx-xs+pl] + vp.gp[nx-2-xs+pl];
432 dj = vp.gp[nx+2 +pl] - 2*vp.gp[nx +pl] + vp.gp[nx-2 +pl];
433 dk = vp.gp[nx+2+xs+pl] - 2*vp.gp[nx+xs+pl] + vp.gp[nx-2+xs+pl];
434 //df = vp.point(i+2, j-1, k-1) - 2*vp.point(i, j-1, k-1) + vp.point(i-2, j-1, k-1);
435 //dg = vp.point(i+2, j, k-1) - 2*vp.point(i, j , k-1) + vp.point(i-2, j, k-1);
436 //dh = vp.point(i+2, j+1, k-1) - 2*vp.point(i, j+1, k-1) + vp.point(i-2, j+1, k-1);
437 //di = vp.point(i+2, j-1, k+1) - 2*vp.point(i, j-1, k+1) + vp.point(i-2, j-1, k+1);
438 //dj = vp.point(i+2, j, k+1) - 2*vp.point(i, j, k+1) + vp.point(i-2, j, k+1);
439 //dk = vp.point(i+2, j+1, k+1) - 2*vp.point(i, j+1, k+1) + vp.point(i-2, j+1, k+1);
440
441 if (k==1 || k==vp.zs-2) {
442 dl = dm = 0;
443 nr = 136;
444 }
445 else {
446 dl = vp.gp[nx+2-pl2] - 2*vp.gp[nx-pl2] + vp.gp[nx-2-pl2];
447 dm = vp.gp[nx+2+pl2] - 2*vp.gp[nx+pl2] + vp.gp[nx-2+pl2];\
448 //dl = vp.point(i+2, j, k-2) - 2*vp.point(i, j, k-2) + vp.point(i-2, j, k-2);
449 //dm = vp.point(i+2, j, k+2) - 2*vp.point(i, j, k+2) + vp.point(i-2, j, k+2);
450 nr = 144;
451 }
452 }
453 xp.gp[nx] = (R)((dc + 4.*(db+dd+dg+dj) + 2.*(df+dh+di+dk) + da+de+dl+dm)/nr);
454 }
455 }
456 if (vcounter!=NULL) vcounter->StepIt();
457 }
458 return xp;
459}
460
461
470template<typename R, typename T> MSGraph<R> yySobel(MSGraph<T> vp)
471{
472 int i, j, k;
473 int pl, nx, ny, nz, pl2, xs, xs2;
474 R da, db, dc, dd, de;
475 R df, dg, dh, di, dj, dk, dl, dm, nr;
476 MSGraph<R> xp;
477
478 xp.mimicry(vp);
479 xp.base = xp.zero = 0;
480
481 if (xp.isNull()) {
482 DEBUG_MODE PRINT_MESG("YYSOBEL: No More Memory!!\n");
483 xp.state = JBXL_GRAPH_MEMORY_ERROR;
484 return xp;
485 }
486
487 // カウンタ
488 CVCounter* vcounter = NULL;
489 if (vp.zs>2) vcounter = GetUsableGlobalCounter();
490 if (vcounter!=NULL) vcounter->SetMax(vp.zs);
491
492 pl = vp.xs*vp.ys;
493 pl2 = 2*pl;
494 xs = vp.xs;
495 xs2 = 2*vp.xs;
496
497 for (k=0; k<vp.zs; k++) {
498 nz = k*pl;
499 for (j=2; j<vp.ys-2; j++) {
500 ny = nz + j*vp.xs;
501 for (i=2; i<vp.xs-2; i++) {
502 nx = ny + i;
503 da = vp.gp[nx-2+xs2] - 2*vp.gp[nx-2] + vp.gp[nx-2-xs2];
504 db = vp.gp[nx-1+xs2] - 2*vp.gp[nx-1] + vp.gp[nx-1-xs2];
505 dc = vp.gp[nx +xs2] - 2*vp.gp[nx] + vp.gp[nx -xs2];
506 dd = vp.gp[nx+1+xs2] - 2*vp.gp[nx+1] + vp.gp[nx+1-xs2];
507 de = vp.gp[nx+2+xs2] - 2*vp.gp[nx+2] + vp.gp[nx+2-xs2];
508 //da = vp.point(i-2, j+2, k) - 2*vp.point(i-2, j, k) + vp.point(i-2, j-2, k);
509 //db = vp.point(i-1, j+2, k) - 2*vp.point(i-1, j, k) + vp.point(i-1, j-2, k);
510 //dc = vp.point(i, j+2, k) - 2*vp.point(i, j, k) + vp.point(i, j-2, k);
511 //dd = vp.point(i+1, j+2, k) - 2*vp.point(i+1, j, k) + vp.point(i+1, j-2, k);
512 //de = vp.point(i+2, j+2, k) - 2*vp.point(i+2, j, k) + vp.point(i+2, j-2, k);
513
514 if (k==0 || k==vp.zs-1) {
515 dc = (R)(6.*dc);
516 df = dg = dh = di = dj = dk = dl = dm = 0;
517 nr = 64;
518 }
519 else {
520 dc = (R)(8.*dc);
521 df = vp.gp[nx-1+xs2-pl] - 2*vp.gp[nx-1-pl] + vp.gp[nx-1-xs2-pl];
522 dg = vp.gp[nx +xs2-pl] - 2*vp.gp[nx -pl] + vp.gp[nx -xs2-pl];
523 dh = vp.gp[nx+1+xs2-pl] - 2*vp.gp[nx+1-pl] + vp.gp[nx+1-xs2-pl];
524 di = vp.gp[nx-1+xs2+pl] - 2*vp.gp[nx-1+pl] + vp.gp[nx-1-xs2+pl];
525 dj = vp.gp[nx +xs2+pl] - 2*vp.gp[nx +pl] + vp.gp[nx -xs2+pl];
526 dk = vp.gp[nx+1+xs2+pl] - 2*vp.gp[nx+1+pl] + vp.gp[nx+1-xs2+pl];
527 //df = vp.point(i-1, j+2, k-1) - 2*vp.point(i-1, j, k-1) + vp.point(i-1, j-2, k-1);
528 //dg = vp.point(i, j+2, k-1) - 2*vp.point(i, j, k-1) + vp.point(i, j-2, k-1);
529 //dh = vp.point(i+1, j+2, k-1) - 2*vp.point(i+1, j, k-1) + vp.point(i+1, j-2, k-1);
530 //di = vp.point(i-1, j+2, k+1) - 2*vp.point(i-1, j, k+1) + vp.point(i-1, j-2, k+1);
531 //dj = vp.point(i, j+2, k+1) - 2*vp.point(i, j, k+1) + vp.point(i, j-2, k+1);
532 //dk = vp.point(i+1, j+2, k+1) - 2*vp.point(i+1, j, k+1) + vp.point(i+1, j-2, k+1);
533
534 if (k==1 || k==vp.zs-2) {
535 dl = dm = 0;
536 nr = 136;
537 }
538 else {
539 dl = vp.gp[nx+xs2-pl2] - 2*vp.gp[nx-pl2] + vp.gp[nx-xs2-pl2];
540 dm = vp.gp[nx+xs2+pl2] - 2*vp.gp[nx+pl2] + vp.gp[nx-xs2+pl2];
541 //dl = vp.point(i, j+2, k-2) - 2*vp.point(i, j, k-2) + vp.point(i, j-2, k-2);
542 //dm = vp.point(i, j+2, k+2) - 2*vp.point(i, j, k+2) + vp.point(i, j-2, k+2);
543 nr = 144;
544 }
545 }
546 xp.gp[nx] = (R)((dc + 4.*(db+dd+dg+dj) + 2.*(df+dh+di+dk) + da+de+dl+dm)/nr);
547 }
548 }
549 if (vcounter!=NULL) vcounter->StepIt();
550 }
551
552 return xp;
553}
554
555
564template<typename R, typename T> MSGraph<R> zzSobel(MSGraph<T> vp)
565{
566 int i, j, k;
567 R da, db, dc, dd, de;
568 R df, dg, dh, di, dj, dk, dl, dm;
569 MSGraph<R> pp, xp;
570
571 if (vp.zs<2) { // 0 のグラフィックデータを返す
572 pp.mimicry(vp);
573 pp.state = JBXL_GRAPH_NODATA_ERROR;
574 return pp;
575 }
576
577 // カウンタ
578 CVCounter* vcounter = NULL;
579 CVCounter* ccounter = NULL;
580 if (vp.zs>2) vcounter = GetUsableGlobalCounter();
581 if (vcounter!=NULL) {
582 vcounter->SetMax(200);
583 ccounter = vcounter->MakeChildCounter(100);
584 SetGlobalCounter(ccounter);
585 }
586
587 pp = zSobel<R>(vp);
588
589 if (vcounter!=NULL) {
590 vcounter->DeleteChildCounter();
591 ccounter = vcounter->MakeChildCounter(100);
592 SetGlobalCounter(ccounter);
593 }
594
595 if (!pp.isNull()) {
596 xp = zSobel<R>(pp);
597 pp.free();
598 }
599 else xp = pp;
600
601 if (vcounter!=NULL) {
602 vcounter->DeleteChildCounter();
603 SetGlobalCounter(vcounter);
604 vcounter->PutFill();
605 }
606
607/* if (vp.zs<5) return xp;
608 for (k=2; k<vp.zs-2; k++) {
609 for (j=2; j<vp.ys-2; j++) {
610 for (i=2; i<vp.xs-2; i++) {
611 da = vp.point(i, j, k+2) - 2*vp.point(i, j, k) + vp.point(i, j, k-2);
612 db = vp.point(i+1, j, k+2) - 2*vp.point(i+1, j, k) + vp.point(i+1, j, k-2);
613 dc = vp.point(i-1, j, k+2) - 2*vp.point(i-1, j, k) + vp.point(i-1, j, k-2);
614 dd = vp.point(i, j+1, k+2) - 2*vp.point(i, j+1, k) + vp.point(i, j+1, k-2);
615 de = vp.point(i, j-1, k+2) - 2*vp.point(i, j-1, k) + vp.point(i, j-1, k-2);
616 df = vp.point(i+1, j+1, k+2) - 2*vp.point(i+1, j+1, k) + vp.point(i+1, j+1, k-2);
617 dg = vp.point(i+1, j-1, k+2) - 2*vp.point(i+1, j-1, k) + vp.point(i+1, j-1, k-2);
618 dh = vp.point(i-1, j+1, k+2) - 2*vp.point(i-1, j+1, k) + vp.point(i-1, j+1, k-2);
619 di = vp.point(i-1, j-1, k+2) - 2*vp.point(i-1, j-1, k) + vp.point(i-1, j-1, k-2);
620 dj = vp.point(i+2, j, k+2) - 2*vp.point(i+2, j, k) + vp.point(i+2, j, k-2);
621 dk = vp.point(i-2, j, k+2) - 2*vp.point(i-2, j, k) + vp.point(i-2, j, k-2);
622 dl = vp.point(i, j+2, k+2) - 2*vp.point(i, j+2, k) + vp.point(i, j+2, k-2);
623 dm = vp.point(i, j-2, k+2) - 2*vp.point(i, j-2, k) + vp.point(i, j-2, k-2);
624 xp.point(i, j, k) = (R)((8.*da + 4.*(db+dc+dd+de) + 2.*(df+dg+dh+di) +dj+dk+dl+dm)/144.);
625 }
626 }
627 }
628*/
629 return xp;
630}
631
632
641template <typename R, typename T> MSGraph<Vector<R> > vNabla(MSGraph<T> vp)
642{
643 int i;
644 MSGraph<R> px, py, pz;
645 MSGraph<Vector<R> > nv;
646
647 //MSGraph<Vector<R> > nv(vp.xs, vp.ys, vp.zs);
648 nv.xs = vp.xs;
649 nv.ys = vp.ys;
650 nv.zs = vp.zs;
651 nv.zero.set(vp.zero, vp.zero, vp.zero);
652 nv.base.set(vp.base, vp.base, vp.base);
653 nv.RZxy = vp.RZxy;
654 nv.rbound = vp.rbound;
655
656 nv.gp = (Vector<R>*)malloc(sizeof(Vector<R>)*nv.xs*nv.ys*nv.zs);
657 if (nv.isNull()) {
658 DEBUG_MODE PRINT_MESG("vNabla: No More Memory!!\n");
659 nv.state = JBXL_GRAPH_MEMORY_ERROR;
660 return nv;
661 }
662 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
663 nv.gp[i] = nv.base;
664 }
665
666 px = xSobel<R>(vp);
667 if (px.gp==NULL) {
668 nv.state = px.state;
669 return nv;
670 }
671
672 py = ySobel<R>(vp);
673 if (py.gp==NULL) {
674 px.free();
675 nv.state = py.state;
676 return nv;
677 }
678
679 pz = zSobel<R>(vp); // 2Dなら 0が入る
680 if (pz.gp==NULL) {
681 px.free();
682 py.free();
683 nv.state = pz.state;
684 return nv;
685 }
686
687 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
688 (nv.gp[i])->set_Vector(px.gp[i], py.gp[i], pz.gp[i]);
689 }
690
691 px.free();
692 py.free();
693 pz.free();
694
695 return nv;
696}
697
698
707template <typename R, typename T> MSGraph<R> Nabla(MSGraph<T> vp)
708{
709 int i;
710 R xx, yy, zz;
711 MSGraph<R> px, py, pz, nv;
712
713 nv.mimicry(vp);
714 if (nv.isNull()) {
715 DEBUG_MODE PRINT_MESG("Nabla: No More Memory!!\n");
716 nv.state = JBXL_GRAPH_MEMORY_ERROR;
717 return nv;
718 }
719
720 px = xSobel<R>(vp);
721 if (px.gp==NULL) {
722 nv.state = px.state;
723 return nv;
724 }
725
726 py = ySobel<R>(vp);
727 if (py.gp==NULL) {
728 px.free();
729 nv.state = py.state;
730 return nv;
731 }
732
733 pz = zSobel<R>(vp);
734 if (pz.gp==NULL) {
735 px.free();
736 py.free();
737 nv.state = pz.state;
738 return nv;
739 }
740
741 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
742 xx = px.gp[i];
743 yy = py.gp[i];
744 zz = pz.gp[i];
745 nv.gp[i] = (R)sqrt((double)xx*xx + yy*yy + zz*zz);
746 }
747
748 px.free();
749 py.free();
750 pz.free();
751
752 return nv;
753}
754
755
768template <typename R, typename T> MSGraph<R> edgeEnhance(MSGraph<T> gd, int mode=0)
769{
770 int i;
771 MSGraph<R> la, vp;
772
773 vp.mimicry(gd);
774 if (vp.isNull()) {
775 DEBUG_MODE PRINT_MESG("edgeEnhance: No More Memory!!\n");
776 vp.state = JBXL_GRAPH_MEMORY_ERROR;
777 return vp;
778 }
779
780 la = Laplacian<R>(gd, mode);
781 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
782 vp.gp[i] = gd.gp[i] - la.gp[i];
783 }
784 la.free();
785
786 return vp;
787}
788
789
800template <typename T> MSGraph<T> medianFilter(MSGraph<T> xp, int ms=3)
801{
802 int i, j, x, y, z, cx;
803 int xx, yy, zz, cw, ux, mz;
804 int kc, xc, zc, xs, ps;
805 T* me;
806 MSGraph<T> vp;
807
808 mz = Min(ms, xp.zs);
809 me = (T*)malloc(ms*ms*mz*sizeof(T));
810 if (me!=NULL) memset(me, 0, ms*ms*mz*sizeof(T));
811
812 vp.mimicry(xp);
813 if (vp.isNull()) {
814 free(me);
815 DEBUG_MODE PRINT_MESG("medianFilter: No More Memory!!\n");
816 vp.state = JBXL_GRAPH_MEMORY_ERROR;
817 return vp;
818 }
819
820 kc = ms*ms*mz/2;
821 xc = ms/2;
822 zc = mz/2;
823 xs = xp.xs;
824 ps = xp.xs*xp.ys;
825 z = xp.zs/2;
826 for(y=xc; y<xp.ys-xc; y++)
827 for(x=xc; x<xp.xs-xc; x++) {
828 cx = z*ps + y*xs + x;
829 i = 0;
830 for (zz=-zc; zz<=zc; zz++)
831 for (yy=-xc; yy<=xc; yy++)
832 for (xx=-xc; xx<=xc; xx++) {
833 cw = cx + xx + yy*xs + zz*ps;
834 me[i++] = xp.gp[cw];
835 }
836 for (i=0; i<ms*ms*mz-1; i++)
837 for (j=i+1; j<ms*ms*mz; j++) {
838 if (me[i]<me[j]) {
839 ux = me[i];
840 me[i] = me[j];
841 me[j] = ux;
842 }
843 }
844 vp.gp[cx-z*ps] = me[kc];
845 }
846
847 free(me);
848 return vp;
849}
850
851
864template <typename T> MSGraph<int> euclidDistance(MSGraph<T> vp, int bc, int& rr)
865{
866 int i, j, k, l, df, d, w;
867 int rmax, rstart, rend;
868 int nx, ny, nz, pl;
869
870 rr = -1;
871 MSGraph<int> pp(vp.xs, vp.ys, vp.zs, (int)vp.zero, (int)vp.base, vp.RZxy);
872 if (pp.isNull()) {
873 DEBUG_MODE PRINT_MESG("euclidDistance: No More Memory!! E1\n");
874 pp.state = JBXL_GRAPH_MEMORY_ERROR;
875 return pp;
876 }
877
878 for (i=0; i<vp.xs*vp.ys*vp.zs; i++) {
879 if (vp.gp[i]>=bc) pp.gp[i] = 1;
880 else pp.gp[i] = 0;
881 }
882
883 pl = vp.xs*vp.ys;
884
885 for (k=0; k<vp.zs; k++) {
886 nz = k*pl;
887 for (j=0; j<vp.ys; j++) {
888 df = vp.xs;
889 ny = nz + j*vp.xs;
890 for (i=0; i<vp.xs; i++) {
891 nx = ny + i;
892 if (pp.gp[nx]!=0) df = df + 1;
893 else df = 0;
894 pp.gp[nx] = df*df;
895 }
896 }
897 }
898
899 for (k=0; k<vp.zs; k++) {
900 nz = k*pl;
901 for (j=0; j<vp.ys; j++) {
902 df = vp.xs;
903 ny = nz + j*vp.xs;
904 for (i=vp.xs-1; i>=0; i--) {
905 nx = ny + i;
906 if (pp.gp[nx]!=0) df = df + 1;
907 else df = 0;
908 pp.gp[nx] = Min(pp.gp[nx], df*df);
909 }
910 }
911 }
912
913 rmax = Max(vp.ys, vp.zs);
914 MSGraph<int> buf(rmax);
915 if (buf.isNull()) {
916 pp.free();
917 DEBUG_MODE PRINT_MESG("euclidDistance: No More Memory!! E2\n");
918 pp.state = JBXL_GRAPH_MEMORY_ERROR;
919 return pp;
920 }
921
922 for (k=0; k<vp.zs; k++) {
923 nz = k*pl;
924 for (i=0; i<vp.xs; i++) {
925 nx = nz + i;
926 for (j=0; j<vp.ys; j++) buf.gp[j] = pp.gp[nx+j*vp.xs];
927 for (j=0; j<vp.ys; j++) {
928 ny = nx + j*vp.xs;
929 d = buf.gp[j];
930 if (d!=0) {
931 rmax = (int)sqrt((double)d) + 1;
932 rstart = Min(rmax, j);
933 rend = Min(rmax, vp.ys-j-1);
934 for (l=-rstart; l<=rend; l++) {
935 w = buf.gp[j+l] + l*l;
936 if (w<d) d = w;
937 }
938 }
939 pp.gp[ny] = d;
940 }
941 }
942 }
943 buf.clear();
944
945 rr = 0;
946 for (j=0; j<vp.ys; j++) {
947 ny = j*vp.xs;
948 for (i=0; i<vp.xs; i++) {
949 nx = ny + i;
950 for (k=0; k<vp.zs; k++) buf.gp[k] = pp.gp[nx+k*pl];
951 for (k=0; k<vp.zs; k++) {
952 nz = nx + k*pl;
953 d = buf.gp[k];
954 if (d!=0) {
955 rmax = (int)sqrt((double)d) + 1;
956 rstart = Min(rmax, k);
957 rend = Min(rmax, vp.zs-k-1);
958 for (l=-rstart; l<=rend; l++) {
959 w = buf.gp[k+l] + l*l;
960 if (w<d) d = w;
961 }
962 rr = Max(rr, d);
963 }
964 pp.gp[nz] = d;
965 }
966 }
967 }
968 buf.free();
969
970 return pp;
971}
972
973
974
976// フィルター
977//
978
979#define FILTER_NON 0
980#define FILTER_ABS 1
981#define FILTER_MINMAX 2
982#define FILTER_NORM 3
983
984
985//
986// mode:
987//
988template <typename R, typename T> MSGraph<R> MSMaskFilter(MSGraph<R> vp, MSGraph<T> filter, int mode=FILTER_NON)
989{
990 MSGraph<R> xp;
991
992 if (vp.xs<filter.xs || vp.ys<filter.ys || vp.zs<filter.zs) {
993 DEBUG_MODE PRINT_MESG("MSMaskFilter: Error: mismach filter dimension!!\n");
994 xp.state = JBXL_GRAPH_NODATA_ERROR;
995 return xp;
996 }
997 if (filter.norm==0.0) {
998 DEBUG_MODE PRINT_MESG("MSMaskFilter: Error: norm of filter is zero!!\n");
999 xp.state = JBXL_GRAPH_NODATA_ERROR;
1000 return xp;
1001 }
1002
1003 xp.mimicry(vp);
1004
1005 int xs = filter.xs/2;
1006 int ys = filter.ys/2;
1007 int zs = filter.zs/2;
1008
1009 for (int k=zs; k<xp.zs-zs; k++) {
1010 for (int j=ys; j<xp.ys-ys; j++) {
1011 for (int i=xs; i<xp.xs-xs; i++) {
1012
1013 T conv = (T)0;
1014 for (int n=-zs; n<=zs; n++) {
1015 for (int m=-ys; m<=ys; m++) {
1016 for (int l=-xs; l<=xs; l++) {
1017 conv += filter.point(xs+l, ys+m, zs+n) * vp.point(i+l, j+m, k+n);
1018 }
1019 }
1020 }
1021
1022 R pt = (R)(conv/filter.norm);
1023
1024 if (mode==FILTER_ABS && pt<(R)0) pt = -pt;
1025 else if (mode==FILTER_MINMAX) {
1026 if (pt<(R)vp.min) pt = (R)vp.min;
1027 else if (pt>(R)vp.max) pt = (R)vp.max;
1028 }
1029 xp.point(i, j, k) = pt;
1030 }
1031 }
1032 }
1033
1034 xp.get_minmax();
1035 if (mode==FILTER_NORM && xp.max!=xp.min) {
1036 for (int i=0; i<xp.xs*xp.ys*xp.zs; i++) {
1037 xp.gp[i] = (R)(((T)(xp.gp[i]-xp.min)*(vp.max-vp.min))/(xp.max-xp.min) + vp.min);
1038 }
1039 xp.get_minmax();
1040 }
1041
1042 return xp;
1043}
1044
1045
1046} // namespace
1047
1048
1049#endif
1050
Gdata.h
グラフィックデータ定義用ヘッダ
FILTER_NON
#define FILTER_NON
何もしない
Definition Gmt.h:979
FILTER_MINMAX
#define FILTER_MINMAX
元のデータの範囲に限定
Definition Gmt.h:981
FILTER_NORM
#define FILTER_NORM
元のデータの範囲に伸張
Definition Gmt.h:982
FILTER_ABS
#define FILTER_ABS
絶対値
Definition Gmt.h:980
jbxl::CVCounter::StepIt
virtual void StepIt(int n=1)
カウンタのメモリを増やす
Definition ClassBox.h:171
jbxl::CVCounter::SetMax
virtual void SetMax(int m)
カウンタの最大値（最終目標）を設定
Definition ClassBox.h:161
jbxl::CVCounter::DeleteChildCounter
virtual void DeleteChildCounter()
子カウンタの削除（有効領域の無効化）
Definition ClassBox.h:179
jbxl::CVCounter::MakeChildCounter
virtual CVCounter * MakeChildCounter(int n)
子カウンタの作成（有効領域を再定義）
Definition ClassBox.h:178
jbxl::CVCounter::PutFill
virtual void PutFill()
取り敢えずの目標（最短目標）までカウンタを進める．
Definition ClassBox.h:164
jbxl::MSGraph::clear
void clear(T v)
全空間を画素値 v にする
Definition Gdata.h:352
jbxl::MSGraph::gp
T * gp
グラフィックデータへのポインタ.
Definition Gdata.h:81
jbxl::MSGraph::base
T base
画措置の底上げの値．
Definition Gdata.h:83
jbxl::MSGraph::point
T & point(int x, int y=0, int z=0)
座標(x,y,z)の画素値の参照
Definition Gdata.h:114
jbxl::MSGraph::zs
int zs
zサイズ. 4Byte. 2Dの場合は 1.
Definition Gdata.h:80
jbxl::MSGraph::get_minmax
void get_minmax(void)
min, max を獲得
Definition Gdata.h:298
jbxl::MSGraph::max
T max
画素値の最大値
Definition Gdata.h:87
jbxl::MSGraph::min
T min
画素値の最小値
Definition Gdata.h:88
jbxl::MSGraph::norm
double norm
規格化定数．フィルタのときに使用．
Definition Gdata.h:93
jbxl::MSGraph::xs
int xs
xサイズ. 4Byte.
Definition Gdata.h:78
jbxl::MSGraph::state
int state
エラー制御
Definition Gdata.h:90
jbxl::MSGraph::RZxy
double RZxy
Z軸の歪．Z軸の間隔を 1とした XY軸の間隔．(X or Y)/Z.
Definition Gdata.h:92
jbxl::MSGraph::zero
T zero
画素値のゼロ位．
Definition Gdata.h:82
jbxl::MSGraph::free
void free(void)
グラフィックデータを開放する
Definition Gdata.h:378
jbxl::MSGraph::isNull
bool isNull(void)
グラフィックデータを持っていないか？
Definition Gdata.h:194
jbxl::MSGraph::mimicry
void mimicry(MSGraph< R > s)
Definition Gdata.h:134
jbxl::MSGraph::rbound
RBound< int > rbound
画像の境界情報
Definition Gdata.h:91
jbxl::MSGraph::ys
int ys
yサイズ. 4Byte.
Definition Gdata.h:79
Min
#define Min(x, y)
Definition common.h:250
Max
#define Max(x, y)
Definition common.h:247
JBXL_GRAPH_NODATA_ERROR
#define JBXL_GRAPH_NODATA_ERROR
データが無い
Definition jbxl_state.h:171
JBXL_GRAPH_MEMORY_ERROR
#define JBXL_GRAPH_MEMORY_ERROR
メモリエラー
Definition jbxl_state.h:170
jbxl
Definition Brep.h:29
jbxl::Laplacian
MSGraph< R > Laplacian(MSGraph< T > vp, int mode=0)
Definition Gmt.h:58
jbxl::GetUsableGlobalCounter
CVCounter * GetUsableGlobalCounter()
現在有効なグローバルカウンタを得る．（子カウンタを得るかもしれない）
Definition ClassBox.h:185
jbxl::euclidDistance
MSGraph< int > euclidDistance(MSGraph< T > vp, int bc, int &rr)
Definition Gmt.h:864
jbxl::xSobel
MSGraph< R > xSobel(MSGraph< T > vp)
Definition Gmt.h:154
jbxl::xxSobel
MSGraph< R > xxSobel(MSGraph< T > vp)
Definition Gmt.h:377
jbxl::edgeEnhance
MSGraph< R > edgeEnhance(MSGraph< T > gd, int mode=0)
Definition Gmt.h:768
jbxl::ySobel
MSGraph< R > ySobel(MSGraph< T > vp)
Definition Gmt.h:219
jbxl::zzSobel
MSGraph< R > zzSobel(MSGraph< T > vp)
Definition Gmt.h:564
jbxl::SetGlobalCounter
void SetGlobalCounter(CVCounter *counter)
グローバルカウンタのセット
Definition ClassBox.h:192
jbxl::MSMaskFilter
MSGraph< R > MSMaskFilter(MSGraph< R > vp, MSGraph< T > filter, int mode=FILTER_NON)
Definition Gmt.h:988
jbxl::zSobel
MSGraph< R > zSobel(MSGraph< T > vp)
Definition Gmt.h:284
jbxl::vNabla
MSGraph< Vector< R > > vNabla(MSGraph< T > vp)
Definition Gmt.h:641
jbxl::medianFilter
MSGraph< T > medianFilter(MSGraph< T > xp, int ms=3)
Definition Gmt.h:800
jbxl::Nabla
MSGraph< R > Nabla(MSGraph< T > vp)
Definition Gmt.h:707
jbxl::yySobel
MSGraph< R > yySobel(MSGraph< T > vp)
Definition Gmt.h:470
PRINT_MESG
#define PRINT_MESG(...)
環境依存用の出力関数．MS Windows用は未実装
Definition tools.h:469
DEBUG_MODE
#define DEBUG_MODE
Definition tools.h:502