1
2
3 /*
4 A* -------------------------------------------------------------------
5 B* This file contains source code for the PyMOL computer program
6 C* Copyright (c) Schrodinger, LLC.
7 D* -------------------------------------------------------------------
8 E* It is unlawful to modify or remove this copyright notice.
9 F* -------------------------------------------------------------------
10 G* Please see the accompanying LICENSE file for further information.
11 H* -------------------------------------------------------------------
12 I* Additional authors of this source file include:
13 -*
14 -*
15 -*
16 Z* -------------------------------------------------------------------
17 */
18 #include"os_python.h"
19 #include"os_predef.h"
20 #include"os_std.h"
21
22 #include"OVRandom.h"
23 #include"OVContext.h"
24
25 #include"Isosurf.h"
26 #include"MemoryDebug.h"
27 #include"Err.h"
28 #include"Symmetry.h"
29 #include"Vector.h"
30 #include"Feedback.h"
31 #include"PConv.h"
32 #include"P.h"
33 #include"Util.h"
34
35 #define Trace_OFF
36
37 #define O3(field,P1,P2,P3,offs) Ffloat3(field,(P1)+offs[0],(P2)+offs[1],(P3)+offs[2])
38
39 #define O3Ptr(field,P1,P2,P3,offs) Ffloat3p(field,(P1)+offs[0],(P2)+offs[1],(P3)+offs[2])
40
41 #define O4(field,P1,P2,P3,P4,offs) Ffloat4(field,(P1)+offs[0],(P2)+offs[1],(P3)+offs[2],P4)
42
43 #define O4Ptr(field,P1,P2,P3,P4,offs) Ffloat4p(field,(P1)+offs[0],(P2)+offs[1],(P3)+offs[2],P4)
44
45 #define I3(field,P1,P2,P3) Fint3(field,P1,P2,P3)
46
47 #define I3Ptr(field,P1,P2,P3) Fint3p(field,P1,P2,P3)
48
49 #define I4(field,P1,P2,P3,P4) Fint4(field,P1,P2,P3,P4)
50
51 #define I4Ptr(field,P1,P2,P3,P4) Fint4p(field,P1,P2,P3,P4)
52
53 typedef struct PointType {
54 float Point[3];
55 int NLink;
56 struct PointType* Link[4];
57 } PointType;
58
59 #define EdgePtPtr(field,P2,P3,P4,P5) ((PointType*)Fvoid4p(field,P2,P3,P4,P5))
60
61 #define EdgePt(field,P2,P3,P4,P5) (*((PointType*)Fvoid4p(field,P2,P3,P4,P5)))
62
63 struct CIsosurf {
64 PyMOLGlobals *G;
65 CField *VertexCodes;
66 CField *ActiveEdges;
67 CField *Point;
68 int NLine;
69 int Skip;
70 int AbsDim[3], CurDim[3], CurOff[3];
71 int Max[3];
72 CField *Coord, *Data;
73 float Level;
74 int Code[256];
75
76 pymol::vla<int>* Num = nullptr;
77 int NSeg;
78 pymol::vla<float>* Line = nullptr;
79
80 };
81
82 static int IsosurfAlloc(PyMOLGlobals * G, CIsosurf * II);
83 static void IsosurfPurge(CIsosurf * II);
84 static int IsosurfCurrent(CIsosurf * II);
85 static int IsosurfCodeVertices(CIsosurf * II);
86 static int IsosurfFindActiveEdges(CIsosurf * II);
87 static int IsosurfFindLines(CIsosurf * II);
88 static int IsosurfDrawLines(CIsosurf * II);
89 static void IsosurfCode(CIsosurf * II, const char *bits1, const char *bits2);
90 static int IsosurfDrawPoints(CIsosurf * II);
91 static int IsosurfPoints(CIsosurf * II);
92 static int IsosurfGradients(PyMOLGlobals * G, CSetting * set1, CSetting * set2,
93 CIsosurf * II, Isofield * field,
94 int *range, float min_level, float max_level);
95
96 #define IsosurfSubSize 64
97
_IsosurfFree(CIsosurf * I)98 static void _IsosurfFree(CIsosurf * I)
99 {
100 FreeP(I);
101 }
102
IsosurfFree(PyMOLGlobals * G)103 void IsosurfFree(PyMOLGlobals * G)
104 {
105 _IsosurfFree(G->Isosurf);
106 G->Isosurf = NULL;
107 }
108
109
110 /*===========================================================================*/
IsosurfAsPyList(PyMOLGlobals * G,Isofield * field)111 PyObject *IsosurfAsPyList(PyMOLGlobals * G, Isofield * field)
112 {
113 PyObject *result = NULL;
114
115 result = PyList_New(4);
116
117 PyList_SetItem(result, 0, PConvIntArrayToPyList(field->dimensions, 3));
118 PyList_SetItem(result, 1, PyInt_FromLong(field->save_points));
119 PyList_SetItem(result, 2, FieldAsPyList(G, field->data.get()));
120 if(field->save_points)
121 PyList_SetItem(result, 3, FieldAsPyList(G, field->points.get()));
122 else
123 PyList_SetItem(result, 3, PConvAutoNone(NULL));
124 return (PConvAutoNone(result));
125 }
126
127
128 /*===========================================================================*/
129 inline
IsosurfInterpolate(CIsosurf * I,float * v1,float * l1,float * v2,float * l2,float * pt)130 static void IsosurfInterpolate(CIsosurf * I, float *v1, float *l1, float *v2, float *l2,
131 float *pt)
132 {
133 float ratio;
134 ratio = (I->Level - *l1) / (*l2 - *l1);
135 pt[0] = v1[0] + (v2[0] - v1[0]) * ratio;
136 pt[1] = v1[1] + (v2[1] - v1[1]) * ratio;
137 pt[2] = v1[2] + (v2[2] - v1[2]) * ratio;
138 }
139
140
141 /*===========================================================================*/
IsosurfNewFromPyList(PyMOLGlobals * G,PyObject * list)142 Isofield *IsosurfNewFromPyList(PyMOLGlobals * G, PyObject * list)
143 {
144 int ok = true;
145 int dim4[4];
146 int a;
147
148 Isofield *result = NULL;
149 if(ok)
150 ok = (list != NULL);
151 if(ok)
152 ok = PyList_Check(list);
153 /* TO ENABLE BACKWARDS COMPATIBILITY...
154 Always check ll when adding new PyList_GetItem's */
155 if(ok)
156 ok = (result = new Isofield()) != nullptr;
157 if(ok)
158 ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list, 0), result->dimensions, 3);
159 if(ok)
160 ok = PConvPyIntToInt(PyList_GetItem(list, 1), &result->save_points);
161 if(ok){
162 result->data.reset(FieldNewFromPyList_From_List(G, list, 2));
163 ok = result->data != nullptr;
164 }
165 if(ok) {
166 if(result->save_points) {
167 result->points.reset(FieldNewFromPyList_From_List(G, list, 3));
168 ok = result->points != nullptr;
169 }
170 else {
171 for(a = 0; a < 3; a++)
172 dim4[a] = result->dimensions[a];
173 dim4[3] = 3;
174 result->points.reset(CField::make<float>(G, dim4, 4));
175 ok = result->points != nullptr;
176 }
177 }
178 if(!ok) {
179 DeleteP(result);
180 }
181 return (result);
182 }
183
184 /*===========================================================================*/
IsofieldComputeGradients(PyMOLGlobals * G,Isofield * field)185 void IsofieldComputeGradients(PyMOLGlobals * G, Isofield * field)
186 {
187 int dim[4];
188 int a, b, c;
189 CField *data = field->data.get();
190
191 if(!field->gradients) {
192
193 /* compute gradients relative to grid axis spacing */
194
195 for(a = 0; a < 3; a++)
196 dim[a] = field->dimensions[a];
197 dim[3] = 3;
198 field->gradients.reset(CField::make<float>(G, dim, 4));
199 auto gradients = field->gradients.get();
200 dim[3] = 3;
201
202 /* bulk internal */
203
204 for(a = 1; a < (dim[0] - 1); a++) {
205 for(b = 1; b < (dim[1] - 1); b++) {
206 for(c = 1; c < (dim[2] - 1); c++) {
207 F4(gradients, a, b, c, 0) =
208 (F3(data, a + 1, b, c) - F3(data, a - 1, b, c)) / 2.0F;
209 F4(gradients, a, b, c, 1) =
210 (F3(data, a, b + 1, c) - F3(data, a, b - 1, c)) / 2.0F;
211 F4(gradients, a, b, c, 2) =
212 (F3(data, a, b, c + 1) - F3(data, a, b, c - 1)) / 2.0F;
213 }
214 }
215 }
216
217 for(a = 0; a < dim[0]; a += (dim[0] - 1)) {
218
219 /* 'a' faces */
220 for(b = 1; b < (dim[1] - 1); b++) {
221 for(c = 1; c < (dim[2] - 1); c++) {
222 if(!a) {
223 F4(gradients, a, b, c, 0) = (F3(data, a + 1, b, c) - F3(data, a, b, c));
224 } else {
225 F4(gradients, a, b, c, 0) = (F3(data, a, b, c) - F3(data, a - 1, b, c));
226 }
227 F4(gradients, a, b, c, 1) =
228 (F3(data, a, b + 1, c) - F3(data, a, b - 1, c)) / 2.0F;
229 F4(gradients, a, b, c, 2) =
230 (F3(data, a, b, c + 1) - F3(data, a, b, c - 1)) / 2.0F;
231 }
232 }
233
234 /* 'c' edges and all eight corners */
235 for(b = 0; b < dim[1]; b += (dim[1] - 1)) {
236 for(c = 0; c < dim[2]; c++) {
237
238 if(!a) {
239 F4(gradients, a, b, c, 0) = (F3(data, a + 1, b, c) - F3(data, a, b, c));
240 } else {
241 F4(gradients, a, b, c, 0) = (F3(data, a, b, c) - F3(data, a - 1, b, c));
242 }
243
244 if(!b) {
245 F4(gradients, a, b, c, 1) = (F3(data, a, b + 1, c) - F3(data, a, b, c));
246 } else {
247 F4(gradients, a, b, c, 1) = (F3(data, a, b, c) - F3(data, a, b - 1, c));
248 }
249
250 if(!c) {
251 F4(gradients, a, b, c, 2) = (F3(data, a, b, c + 1) - F3(data, a, b, c));
252 } else if(c < (dim[2] - 1)) {
253 F4(gradients, a, b, c, 2) =
254 (F3(data, a, b, c + 1) - F3(data, a, b, c - 1)) / 2.0F;
255 } else {
256 F4(gradients, a, b, c, 2) = (F3(data, a, b, c) - F3(data, a, b, c - 1));
257 }
258 }
259 }
260
261 /* 'b' edges */
262 for(c = 0; c < dim[2]; c += (dim[2] - 1)) {
263 for(b = 1; b < (dim[1] - 1); b++) {
264 if(!a) {
265 F4(gradients, a, b, c, 0) = (F3(data, a + 1, b, c) - F3(data, a, b, c));
266 } else {
267 F4(gradients, a, b, c, 0) = (F3(data, a, b, c) - F3(data, a - 1, b, c));
268 }
269
270 F4(gradients, a, b, c, 1) =
271 (F3(data, a, b + 1, c) - F3(data, a, b - 1, c)) / 2.0F;
272
273 if(!c) {
274 F4(gradients, a, b, c, 2) = (F3(data, a, b, c + 1) - F3(data, a, b, c));
275 } else if(c < (dim[2] - 1)) {
276 F4(gradients, a, b, c, 2) =
277 (F3(data, a, b, c + 1) - F3(data, a, b, c - 1)) / 2.0F;
278 } else {
279 F4(gradients, a, b, c, 2) = (F3(data, a, b, c) - F3(data, a, b, c - 1));
280 }
281 }
282 }
283 }
284
285 for(b = 0; b < dim[1]; b += (dim[1] - 1)) {
286
287 for(a = 1; a < (dim[0] - 1); a++) {
288
289 /* 'b' faces */
290
291 for(c = 1; c < (dim[2] - 1); c++) {
292 F4(gradients, a, b, c, 0) =
293 (F3(data, a + 1, b, c) - F3(data, a - 1, b, c)) / 2.0F;
294 if(!b) {
295 F4(gradients, a, b, c, 1) = (F3(data, a, b + 1, c) - F3(data, a, b, c));
296 } else {
297 F4(gradients, a, b, c, 1) = (F3(data, a, b, c) - F3(data, a, b - 1, c));
298 }
299 F4(gradients, a, b, c, 2) =
300 (F3(data, a, b, c + 1) - F3(data, a, b, c - 1)) / 2.0F;
301 }
302
303 /* 'a' edges */
304
305 for(c = 0; c < dim[2]; c += (dim[2] - 1)) {
306 F4(gradients, a, b, c, 0) =
307 (F3(data, a + 1, b, c) - F3(data, a - 1, b, c)) / 2.0F;
308 if(!b) {
309 F4(gradients, a, b, c, 1) = (F3(data, a, b + 1, c) - F3(data, a, b, c));
310 } else {
311 F4(gradients, a, b, c, 1) = (F3(data, a, b, c) - F3(data, a, b - 1, c));
312 }
313 if(!c) {
314 F4(gradients, a, b, c, 2) = (F3(data, a, b, c + 1) - F3(data, a, b, c));
315 } else {
316 F4(gradients, a, b, c, 2) = (F3(data, a, b, c) - F3(data, a, b, c - 1));
317 }
318 }
319 }
320 }
321
322 /* 'c' faces */
323
324 for(c = 0; c < dim[2]; c += (dim[2] - 1)) {
325 for(a = 1; a < (dim[0] - 1); a++) {
326 for(b = 1; b < (dim[1] - 1); b++) {
327 F4(gradients, a, b, c, 0) =
328 (F3(data, a + 1, b, c) - F3(data, a - 1, b, c)) / 2.0F;
329 F4(gradients, a, b, c, 1) =
330 (F3(data, a, b + 1, c) - F3(data, a, b - 1, c)) / 2.0F;
331 if(!c) {
332 F4(gradients, a, b, c, 2) = (F3(data, a, b, c + 1) - F3(data, a, b, c));
333 } else {
334 F4(gradients, a, b, c, 2) = (F3(data, a, b, c) - F3(data, a, b, c - 1));
335 }
336 }
337 }
338 }
339 }
340 }
341
342
343 /*===========================================================================*/
Isofield(PyMOLGlobals * G,const int * const dims)344 Isofield::Isofield(PyMOLGlobals * G, const int * const dims)
345 {
346 int dim4[4];
347 std::copy_n(dims, 3, dim4);
348 dim4[3] = 3;
349
350 /* Warning: ...FromPyList also allocs and inits from the heap */
351
352 data.reset(CField::make<float>(G, dims, 3));
353 points.reset(CField::make<float>(G, dim4, 4));
354 std::copy_n(dims, 3, dimensions);
355 }
356
357 /*===========================================================================*/
IsosurfCode(CIsosurf * II,const char * bits1,const char * bits2)358 static void IsosurfCode(CIsosurf * II, const char *bits1, const char *bits2)
359 {
360 CIsosurf *I = II;
361 int c;
362 int b;
363 int sum1, sum2;
364
365 c = 0;
366 while(bits1[c])
367 c++;
368 c--;
369 sum1 = 0;
370 b = 1;
371 while(c >= 0) {
372 if(bits1[c] == '1')
373 sum1 = sum1 + b;
374 b = b + b;
375 c--;
376 }
377
378 c = 0;
379 while(bits2[c])
380 c++;
381 c--;
382 sum2 = 0;
383 b = 1;
384 while(c >= 0) {
385 if(bits2[c] == '1')
386 sum2 = sum2 + b;
387 b = b + b;
388 c--;
389 }
390
391 I->Code[sum1] = sum2;
392 #ifdef Trace
393 printf("IsosurfCode: %s (%i) -> %s (%i)\n", bits1, sum1, bits2, sum2);
394 #endif
395 }
396
397
398 /*===========================================================================*/
IsosurfNew(PyMOLGlobals * G)399 static CIsosurf *IsosurfNew(PyMOLGlobals * G)
400 {
401 int c;
402 CIsosurf *I = pymol::calloc<CIsosurf>(1);
403 I->G = G;
404 I->VertexCodes = NULL;
405 I->ActiveEdges = NULL;
406 I->Point = NULL;
407 I->Line = NULL;
408 I->Skip = 0;
409 for(c = 0; c < 256; c++)
410 I->Code[c] = -1;
411
412
413 /*___
414 | / |
415 |/ |
416 |___|
417 32
418 */
419 IsosurfCode(I, "10000010", "100000");
420 IsosurfCode(I, "01000001", "100000");
421
422
423 /*___
424 | \ |
425 | \|
426 |___|
427 16
428 */
429 IsosurfCode(I, "10010000", "010000");
430 IsosurfCode(I, "01100000", "010000");
431
432
433 /*___
434 | |
435 | /|
436 |_/_|
437 8
438 */
439 IsosurfCode(I, "00101000", "001000");
440 IsosurfCode(I, "00010100", "001000");
441
442
443 /*___
444 | |
445 |\ |
446 |_\_|
447 4
448 */
449 IsosurfCode(I, "00001001", "000100");
450 IsosurfCode(I, "00000110", "000100");
451
452
453 /*___
454 | \ |
455 |\ \|
456 |_\_|
457 16+4=20
458 */
459
460 IsosurfCode(I, "01101001", "010100");
461
462
463 /*___
464 | / |
465 |/ /|
466 |_/_|
467 32+8=40
468 */
469 IsosurfCode(I, "10010110", "101000");
470
471
472 /*___
473 | | |
474 | | |
475 |_|_|
476 2
477 */
478 IsosurfCode(I, "10001000", "000010");
479 IsosurfCode(I, "01000100", "000010");
480
481
482 /*___
483 | |
484 |---|
485 |___|
486 1
487 */
488 IsosurfCode(I, "00100010", "000001");
489 IsosurfCode(I, "00010001", "000001");
490
491 return (I);
492 }
493
IsosurfInit(PyMOLGlobals * G)494 int IsosurfInit(PyMOLGlobals * G)
495 {
496 G->Isosurf = IsosurfNew(G);
497 return 1;
498 }
499
500
501 /*===========================================================================*/
IsosurfExpand(Isofield * field1,Isofield * field2,CCrystal * cryst,CSymmetry * sym,int * range)502 int IsosurfExpand(Isofield * field1, Isofield * field2, CCrystal * cryst,
503 CSymmetry * sym, int *range)
504 {
505 float rmn[3], rmx[3];
506 float imn[3], imx[3];
507 float fstep[3], frange[3];
508 int field1max[3];
509 int expanded = false;
510 int missing = false;
511
512 field1max[0] = field1->dimensions[0] - 1;
513 field1max[1] = field1->dimensions[1] - 1;
514 field1max[2] = field1->dimensions[2] - 1;
515 {
516 int a;
517 for(a = 0; a < 3; a++) {
518 rmn[a] = F4(field1->points, 0, 0, 0, a);
519 rmx[a] = F4(field1->points, field1max[0], field1max[1], field1max[2], a);
520 }
521 }
522
523 /* get min/max extents of map1 in fractional space */
524
525 transform33f3f(cryst->RealToFrac, rmn, imn);
526 transform33f3f(cryst->RealToFrac, rmx, imx);
527
528 /* compute step size */
529
530 subtract3f(imx, imn, frange);
531
532 fstep[0] = frange[0] / field1max[0];
533 fstep[1] = frange[1] / field1max[1];
534 fstep[2] = frange[2] / field1max[2];
535
536 /* compute coordinate points for second field */
537
538 if (SymmetryAttemptGeneration(sym)) {
539 int i, j, k;
540 int i_stop, j_stop, k_stop;
541 float frac[3];
542 int nMat = sym->getNSymMat();
543
544 i_stop = field2->dimensions[0];
545 j_stop = field2->dimensions[1];
546 k_stop = field2->dimensions[2];
547 for(i = 0; i < i_stop; i++) {
548 frac[0] = imn[0] + fstep[0] * (i + range[0]);
549 for(j = 0; j < j_stop; j++) {
550 frac[1] = imn[1] + fstep[1] * (j + range[1]);
551 for(k = 0; k < k_stop; k++) {
552 float average = 0.0F;
553 float extrapolate_average = 0.0F;
554 int cnt = 0;
555 int extrapolate_cnt = 0;
556
557 /* first compute the coordinate */
558
559 float *ptr = F4Ptr(field2->points, i, j, k, 0);
560 frac[2] = imn[2] + fstep[2] * (k + range[2]);
561 transform33f3f(cryst->FracToReal, frac, ptr);
562
563 /* then compute the value at the coordinate */
564
565 for(int n = nMat - 1; n >= 0; n--) {
566 const float *matrix = sym->getSymMat(n);
567 float test_frac[3];
568
569 transform44f3f(matrix, frac, test_frac);
570
571 /* we're assuming that the identity matrix appears in the list */
572
573 test_frac[0] -= imn[0];
574 test_frac[1] -= imn[1];
575 test_frac[2] -= imn[2];
576
577 test_frac[0] -= (int) floor(test_frac[0] + R_SMALL4);
578 test_frac[1] -= (int) floor(test_frac[1] + R_SMALL4);
579 test_frac[2] -= (int) floor(test_frac[2] + R_SMALL4);
580
581 {
582 int a, b, c;
583 float x, y, z;
584 a = (int) (test_frac[0] / fstep[0]);
585 b = (int) (test_frac[1] / fstep[1]);
586 c = (int) (test_frac[2] / fstep[2]);
587 x = (test_frac[0] / fstep[0]) - a;
588 y = (test_frac[1] / fstep[1]) - b;
589 z = (test_frac[2] / fstep[2]) - c;
590
591 if((a >= 0) && (b >= 0) && (c >= 0) &&
592 (a <= (field1max[0] + 1)) && (b <= (field1max[1] + 1))
593 && (c <= (field1max[2] + 1))) {
594 while(a >= field1max[0]) {
595 a--;
596 x += 1.0F;
597 }
598 while(b >= field1max[1]) {
599 b--;
600 y += 1.0F;
601 }
602 while(c >= field1max[2]) {
603 c--;
604 z += 1.0F;
605 }
606
607 {
608 const float sloppy_1 = 1.0F + R_SMALL4;
609 if((x <= sloppy_1) && (y <= sloppy_1) && (z <= sloppy_1)) {
610 if(!expanded) {
611 if((matrix[0] != 1.0F) || /* not identity matrix */
612 (matrix[5] != 1.0F) ||
613 (matrix[10] != 1.0F) || (matrix[15] != 1.0F) ||
614 /* and not inside source map */
615 ((imn[0] - frac[0]) > R_SMALL4)
616 || ((frac[0] - imx[0]) > R_SMALL4)
617 || ((imn[1] - frac[1]) > R_SMALL4)
618 || ((frac[1] - imx[1]) > R_SMALL4)
619 || ((imn[2] - frac[2]) > R_SMALL4)
620 || ((frac[2] - imx[2]) > R_SMALL4)) {
621 expanded = true;
622 }
623 }
624 /* found at least one point obtained through symmetry or priodicity */
625 if(x > 1.0F)
626 x = 1.0F;
627 if(y > 1.0F)
628 y = 1.0F;
629 if(z > 1.0F)
630 z = 1.0F;
631 average += FieldInterpolatef(field1->data.get(), a, b, c, x, y, z);
632 cnt++;
633 } else {
634 /* allow 1 cell of extrapolation -- this saves us
635 when someone issues map_double and is then technically
636 missing a plane of data at the edge of the cell */
637 if(((x-1.0F) < sloppy_1) && ((y-1.0F) < sloppy_1) && ((z-1.0F) < sloppy_1)) {
638 if(x > 1.0F)
639 x = 1.0F;
640 if(y > 1.0F)
641 y = 1.0F;
642 if(z > 1.0F)
643 z = 1.0F;
644 extrapolate_average += FieldInterpolatef(field1->data.get(), a, b, c, x, y, z);
645 extrapolate_cnt++;
646 }
647 }
648 }
649 }
650 }
651 }
652 if(cnt) {
653 F3(field2->data, i, j, k) = average / cnt;
654 } else if(extrapolate_cnt) {
655 F3(field2->data, i, j, k) = extrapolate_average / extrapolate_cnt;
656 } else {
657 missing = true;
658 F3(field2->data, i, j, k) = 0.0F; /* complain? */
659 }
660 }
661 }
662 }
663 }
664 if(expanded) {
665 if(missing) {
666 return -1;
667 } else {
668 return 1;
669 }
670 } else {
671 return 0;
672 }
673 }
674
IsosurfGetRange(PyMOLGlobals * G,Isofield * field,CCrystal * cryst,float * mn,float * mx,int * range,int clamp)675 int IsosurfGetRange(PyMOLGlobals * G, Isofield * field,
676 CCrystal * cryst, float *mn, float *mx, int *range, int clamp)
677 {
678 float rmn[3], rmx[3];
679 float imn[3], imx[3];
680 float mix[24], imix[24];
681 int a, b;
682 int clamped = false; /* clamped? */
683 PRINTFD(G, FB_Isosurface)
684 " IsosurfGetRange: entered mn: %4.2f %4.2f %4.2f mx: %4.2f %4.2f %4.2f\n",
685 mn[0], mn[1], mn[2], mx[0], mx[1], mx[2]
686 ENDFD;
687
688 for(a = 0; a < 3; a++) {
689 rmn[a] = F4(field->points, 0, 0, 0, a);
690 rmx[a] = F4(field->points,
691 field->dimensions[0] - 1,
692 field->dimensions[1] - 1, field->dimensions[2] - 1, a);
693 }
694
695 /* get min/max extents of map in fractional space */
696
697 transform33f3f(cryst->RealToFrac, rmn, imn);
698 transform33f3f(cryst->RealToFrac, rmx, imx);
699
700 mix[0] = mn[0];
701 mix[1] = mn[1];
702 mix[2] = mn[2];
703
704 mix[3] = mx[0];
705 mix[4] = mn[1];
706 mix[5] = mn[2];
707
708 mix[6] = mn[0];
709 mix[7] = mx[1];
710 mix[8] = mn[2];
711
712 mix[9] = mn[0];
713 mix[10] = mn[1];
714 mix[11] = mx[2];
715
716 mix[12] = mx[0];
717 mix[13] = mx[1];
718 mix[14] = mn[2];
719
720 mix[15] = mx[0];
721 mix[16] = mn[1];
722 mix[17] = mx[2];
723
724 mix[18] = mn[0];
725 mix[19] = mx[1];
726 mix[20] = mx[2];
727
728 mix[21] = mx[0];
729 mix[22] = mx[1];
730 mix[23] = mx[2];
731
732 /* compute min/max of query in fractional space */
733
734 for(b = 0; b < 8; b++) {
735 transform33f3f(cryst->RealToFrac, mix + 3 * b, imix + 3 * b);
736 }
737
738 for(a = 0; a < 3; a++) {
739 if(imx[a] != imn[a]) { /* protect against div by zero */
740 int b;
741 int mini = 0, maxi = 0, tst_min, tst_max;
742 float cur;
743 for(b = 0; b < 8; b++) {
744 cur =
745 ((field->dimensions[a] - 1) * (imix[a + 3 * b] - imn[a]) / (imx[a] - imn[a]));
746 tst_min = (int) floor(cur);
747 tst_max = ((int) ceil(cur)) + 1;
748
749 if(!b) {
750 mini = tst_min;
751 maxi = tst_max;
752 } else {
753 if(mini > tst_min)
754 mini = tst_min;
755 if(maxi <= tst_max)
756 maxi = tst_max;
757 }
758 }
759
760 range[a] = mini;
761
762 range[a + 3] = maxi;
763 } else {
764 range[a] = 0;
765 range[a + 3] = 1;
766 }
767 if(range[a] < 0) {
768 if(clamp)
769 range[a] = 0;
770 clamped = true;
771 }
772 if(range[a] > field->dimensions[a]) {
773 if(clamp)
774 range[a] = field->dimensions[a];
775 clamped = true;
776 }
777 if(range[a + 3] < 0) {
778 if(clamp)
779 range[a + 3] = 0;
780 clamped = true;
781 }
782 if(range[a + 3] > field->dimensions[a]) {
783 if(clamp)
784 range[a + 3] = field->dimensions[a];
785 clamped = true;
786 }
787 }
788 PRINTFD(G, FB_Isosurface)
789 " IsosurfGetRange: returning range: %d %d %d %d %d %d\n",
790 range[0], range[1], range[2], range[3], range[4], range[5]
791 ENDFD;
792 return clamped;
793 }
794
795
796 /*===========================================================================*/
IsosurfVolume(PyMOLGlobals * G,CSetting * set1,CSetting * set2,Isofield * field,float level,pymol::vla<int> & num,pymol::vla<float> & vert,int * range,int mode,int skip,float alt_level)797 int IsosurfVolume(PyMOLGlobals* G, CSetting* set1, CSetting* set2,
798 Isofield* field, float level, pymol::vla<int>& num, pymol::vla<float>& vert,
799 int* range, int mode, int skip, float alt_level)
800 {
801 int ok = true;
802 CIsosurf *I;
803 if(PIsGlutThread()) {
804 I = G->Isosurf;
805 } else {
806 I = IsosurfNew(G);
807 }
808 CHECKOK(ok, I);
809 {
810 int Steps[3];
811 int c, i, j, k;
812 int x, y, z;
813 int range_store[6];
814 I->Num = std::addressof(num);
815 I->Line = std::addressof(vert);
816 I->Skip = skip;
817 if(range) {
818 for(c = 0; c < 3; c++) {
819 I->AbsDim[c] = field->dimensions[c];
820 I->CurDim[c] = IsosurfSubSize + 1;
821 Steps[c] = ((range[3 + c] - range[c]) - 2) / IsosurfSubSize + 1;
822 }
823 } else {
824 range = range_store;
825 for(c = 0; c < 3; c++) {
826 range[c] = 0;
827 range[3 + c] = field->dimensions[c];
828 I->AbsDim[c] = field->dimensions[c];
829 I->CurDim[c] = IsosurfSubSize + 1;
830 Steps[c] = (I->AbsDim[c] - 2) / IsosurfSubSize + 1;
831 }
832 }
833
834 I->Coord = field->points.get();
835 I->Data = field->data.get();
836 I->Level = level;
837 if(ok)
838 ok = IsosurfAlloc(G, I);
839
840 I->NLine = 0;
841 I->NSeg = 0;
842 I->Num->check(I->NSeg);
843 (*I->Num)[I->NSeg] = I->NLine;
844
845 if(ok) {
846 switch (mode) {
847 case 3:
848 ok = IsosurfGradients(G, set1, set2, I, field, range, level, alt_level);
849 IsosurfPurge(I);
850 break;
851 default:
852 for(i = 0; i < Steps[0]; i++) {
853 for(j = 0; j < Steps[1]; j++) {
854 for(k = 0; k < Steps[2]; k++) {
855 if(ok) {
856 I->CurOff[0] = IsosurfSubSize * i;
857 I->CurOff[1] = IsosurfSubSize * j;
858 I->CurOff[2] = IsosurfSubSize * k;
859 for(c = 0; c < 3; c++)
860 I->CurOff[c] += range[c];
861 for(c = 0; c < 3; c++) {
862 I->Max[c] = range[3 + c] - I->CurOff[c];
863 if(I->Max[c] > (IsosurfSubSize + 1))
864 I->Max[c] = (IsosurfSubSize + 1);
865 }
866 if(!(i || j || k)) {
867 for(x = 0; x < I->Max[0]; x++)
868 for(y = 0; y < I->Max[1]; y++)
869 for(z = 0; z < I->Max[2]; z++)
870 for(c = 0; c < 3; c++)
871 EdgePt(I->Point, x, y, z, c).NLink = 0;
872 }
873 #ifdef Trace
874 for(c = 0; c < 3; c++)
875 printf(" IsosurfVolume: c: %i CurOff[c]: %i Max[c] %i\n", c,
876 I->CurOff[c], I->Max[c]);
877 #endif
878
879 if(ok)
880 switch (mode) {
881 case 0: /* standard mode - want lines */
882 ok = IsosurfCurrent(I);
883 break;
884 case 1: /* point mode - just want points on the isosurface */
885 ok = IsosurfPoints(I);
886 break;
887 case 2:
888 /* reserved */
889 break;
890 }
891 if(G->Interrupt) {
892 ok = false;
893 }
894 }
895 }
896 }
897 }
898 IsosurfPurge(I);
899 break;
900 }
901 }
902
903 if(mode) {
904 PRINTFB(G, FB_Isomesh, FB_Blather)
905 " IsosurfVolume: Surface generated using %d dots.\n", I->NLine ENDFB(G);
906 } else {
907 PRINTFB(G, FB_Isomesh, FB_Blather)
908 " IsosurfVolume: Surface generated using %d lines.\n", I->NLine ENDFB(G);
909 }
910
911 if(!ok) {
912 I->NLine = 0;
913 I->NSeg = 0;
914 }
915 /* shrinks sizes for more efficient RAM usage */
916
917 I->Line->resize(I->NLine * 3);
918 I->Num->resize(I->NSeg + 1);
919 (*I->Num)[I->NSeg] = 0; /* important - must terminate the segment list */
920
921 if(!PIsGlutThread()) {
922 _IsosurfFree(I);
923 }
924 }
925 return (ok);
926 }
927
928
929 /*===========================================================================*/
IsosurfAlloc(PyMOLGlobals * G,CIsosurf * II)930 static int IsosurfAlloc(PyMOLGlobals * G, CIsosurf * II)
931 {
932 CIsosurf *I = II;
933
934 int ok = true;
935 int dim4[4];
936 int a;
937 for(a = 0; a < 3; a++)
938 dim4[a] = I->CurDim[a];
939 dim4[3] = 3;
940
941 I->VertexCodes = CField::make<int>(G, I->CurDim, 3);
942 I->ActiveEdges = CField::make<int>(G, dim4, 4);
943 I->Point = CField::make<PointType>(G, dim4, 4);
944 if(!(I->VertexCodes && I->ActiveEdges && I->Point)) {
945 IsosurfPurge(I);
946 ok = false;
947 }
948 #ifdef Trace
949 printf(" IsosurfAlloc: ok: %i\n", ok);
950 #endif
951 return (ok);
952 }
953
954
955 /*===========================================================================*/
IsosurfPurge(CIsosurf * II)956 static void IsosurfPurge(CIsosurf * II)
957 {
958 CIsosurf *I = II;
959 DeleteP(I->VertexCodes);
960 DeleteP(I->ActiveEdges);
961 DeleteP(I->Point);
962 }
963
964
965 /*===========================================================================*/
IsosurfCurrent(CIsosurf * II)966 static int IsosurfCurrent(CIsosurf * II)
967 {
968 CIsosurf *I = II;
969 int ok = true;
970 if(IsosurfCodeVertices(I)) {
971 if(ok)
972 ok = IsosurfFindActiveEdges(I);
973 if(ok)
974 ok = IsosurfFindLines(I);
975 if(ok)
976 ok = IsosurfDrawLines(I);
977 }
978 return (ok);
979 }
980
981
982 /*===========================================================================*/
IsosurfPoints(CIsosurf * II)983 static int IsosurfPoints(CIsosurf * II)
984 {
985 CIsosurf *I = II;
986 int ok = true;
987 if(IsosurfCodeVertices(I)) {
988 if(ok)
989 ok = IsosurfFindActiveEdges(I);
990 if(ok)
991 ok = IsosurfDrawPoints(I);
992 }
993 return (ok);
994 }
995
996
997 /*===========================================================================*/
998
IsosurfGradients(PyMOLGlobals * G,CSetting * set1,CSetting * set2,CIsosurf * II,Isofield * field,int * range,float min_level,float max_level)999 static int IsosurfGradients(PyMOLGlobals * G, CSetting * set1, CSetting * set2,
1000 CIsosurf * II, Isofield * field,
1001 int *range, float min_level, float max_level)
1002 {
1003 CIsosurf *I = II;
1004 int ok = true;
1005
1006 /* use local copies for performance reasons */
1007
1008 int n_seg = I->NSeg;
1009 int n_line = I->NLine;
1010 auto& i_line = *I->Line;
1011 auto i_data = I->Data;
1012 auto& i_num = *I->Num;
1013
1014 /* get cascaded state, object, or global settings */
1015
1016 int spacing = SettingGet_i(G, set1, set2, cSetting_gradient_spacing);
1017 float step_size = SettingGet_f(G, set1, set2, cSetting_gradient_step_size);
1018 float max_walk = SettingGet_f(G, set1, set2, cSetting_gradient_max_length);
1019 float min_walk = SettingGet_f(G, set1, set2, cSetting_gradient_min_length);
1020 float min_slope = SettingGet_f(G, set1, set2, cSetting_gradient_min_slope);
1021 float min_dot = SettingGet_f(G, set1, set2, cSetting_gradient_normal_min_dot);
1022 float symmetry = SettingGet_f(G, set1, set2, cSetting_gradient_symmetry);
1023
1024 int symmetry_flag = false; /* are we searching for symmetric segments? */
1025 if(symmetry != 0.0F)
1026 symmetry_flag = true; /* (very slow process) */
1027 if(symmetry > 1.0F)
1028 symmetry = 1.0F / symmetry;
1029
1030 /* clamp dangerous parameters */
1031
1032 if(step_size < 0.01F)
1033 step_size = 0.01F;
1034
1035 if(min_slope < 0.00001F)
1036 min_slope = 0.00001F;
1037
1038 /* make sure we have gradients available for map */
1039
1040 if(!field->gradients)
1041 IsofieldComputeGradients(G, field);
1042
1043 /* and that map has a minimum size */
1044
1045 if(field->gradients) {
1046
1047 /* locals for performance */
1048
1049 CField *gradients = field->gradients.get();
1050 CField *points = field->points.get();
1051
1052 /* flags marking excluded regions to avoid (currently wasteful) */
1053 int *flag = NULL;
1054
1055 /* variable length array for recording segment paths */
1056 int *active_cell = VLAlloc(int, 1000);
1057
1058 /* ordered list of coordinates for processing */
1059 int *order = NULL;
1060
1061 int range_size; /* total points in region being drawn */
1062 int range_dim[3]; /* dimension of drawn region */
1063 int flag_stride[3]; /* stride values for flag array */
1064
1065 range_dim[0] = (range[3] - range[0]);
1066 range_dim[1] = (range[4] - range[1]);
1067 range_dim[2] = (range[5] - range[2]);
1068
1069 range_size = range_dim[0] * range_dim[1] * range_dim[2];
1070 if (ok)
1071 flag = pymol::calloc<int>(range_size);
1072 CHECKOK(ok, flag);
1073 flag_stride[0] = 1;
1074 flag_stride[1] = range_dim[0];
1075 flag_stride[2] = range_dim[0] * range_dim[1];
1076
1077 if (ok)
1078 order = pymol::calloc<int>(3 * range_size);
1079 CHECKOK(ok, order);
1080 if(order && flag && (range_dim[0] > 1) && (range_dim[1] > 1) && (range_dim[2] > 1)) {
1081
1082 {
1083 /* compute approximate cell spacing */
1084
1085 float average_cell_axis_dist;
1086 float *pos[4];
1087 pos[0] = Ffloat4p(I->Coord, 0, 0, 0, 0);
1088 pos[1] = Ffloat4p(I->Coord, 1, 0, 0, 0);
1089 pos[2] = Ffloat4p(I->Coord, 0, 1, 0, 0);
1090 pos[3] = Ffloat4p(I->Coord, 0, 0, 1, 0);
1091
1092 average_cell_axis_dist = (float) ((diff3f(pos[0], pos[1]) +
1093 diff3f(pos[0], pos[2]) +
1094 diff3f(pos[0], pos[3])) / 3.0);
1095
1096 /* scale parameters into cell units */
1097
1098 max_walk /= average_cell_axis_dist;
1099 min_walk /= average_cell_axis_dist;
1100 step_size /= average_cell_axis_dist;
1101 min_slope *= average_cell_axis_dist;
1102 }
1103
1104 {
1105 /* generate randomized list of cell coordinates */
1106
1107 /* always use same seed for same volume */
1108 OVRandom *my_rand = OVRandom_NewBySeed(G->Context->heap, range_size);
1109 {
1110 /* fill */
1111 int i, j, k, *p = order;
1112 for(k = range[2]; k < range[5]; k++) {
1113 for(j = range[1]; j < range[4]; j++) {
1114 for(i = range[0]; i < range[3]; i++) {
1115 p[0] = i;
1116 p[1] = j;
1117 p[2] = k;
1118 p += 3;
1119 }
1120 }
1121 }
1122 }
1123 {
1124 /* shuffle */
1125 int a;
1126 for(a = 0; a < range_size; a++) {
1127 int *p = order + 3 * (int) (range_size * OVRandom_Get_float64_exc1(my_rand));
1128 int *q = order + 3 * (int) (range_size * OVRandom_Get_float64_exc1(my_rand));
1129 int t0 = p[0], t1 = p[1], t2 = p[2];
1130 p[0] = q[0];
1131 p[1] = q[1];
1132 p[2] = q[2];
1133 q[0] = t0;
1134 q[1] = t1;
1135 q[2] = t2;
1136 }
1137 }
1138 OVRandom_Del(my_rand);
1139 }
1140
1141 {
1142 /* now draw our lines */
1143
1144 int a;
1145 int *start_locus = order;
1146 float prev_grad_normal[3] = { 0.0F, 0.0F, 0.0F };
1147 for(a = 0; a < range_size; a++) {
1148 int n_active_cell = 0; /* how many cells have we traversed */
1149 float walk = max_walk; /* distance remaining to travel */
1150
1151 int abort_n_line = n_line; /* for backtracking */
1152 int abort_n_seg = n_seg;
1153
1154 int pass; /* the pass are we on */
1155
1156 float symmetry_max = FLT_MIN; /* if we're trying to get symmetric segments */
1157 float symmetry_min = FLT_MAX;
1158
1159 for(pass = 0; pass < 2; pass++) { /* one pass down the gradient, one up */
1160
1161 int have_prev = false; /* flag & storage for previous gradient & locus */
1162 int *prev_locus = NULL;
1163
1164 int locus[3]; /* what cell are we in? */
1165 float fract[3] = { 0.0F, 0.0F, 0.0F }; /* where in the cell are we? */
1166 int n_vert = 0;
1167
1168 locus[0] = start_locus[0];
1169 locus[1] = start_locus[1];
1170 locus[2] = start_locus[2];
1171
1172 for(;;) {
1173
1174 /* master segment extension loop, using "break" to exit */
1175
1176 {
1177 /* normalize locus and fract before each new step */
1178 int done = false;
1179 int b;
1180 for(b = 0; b < 3; b++) {
1181
1182 while(fract[b] < 0.0F) { /* force fract >= 0.0 */
1183 fract[b] += 1.0F;
1184 locus[b]--;
1185 }
1186 while(fract[b] >= 1.0F) { /* force fract into [0.0-1.0) */
1187 fract[b] -= 1.0F;
1188 locus[b]++;
1189 }
1190 while(locus[b] > (range[b + 3] - 2)) { /* above range? done */
1191 if(fract[b] <= 0.0F) {
1192 locus[b]--;
1193 fract[b] += 1.0F;
1194 if(locus[b] < range[b]) { /* below range? done */
1195 done = true;
1196 break;
1197 }
1198 } else {
1199 done = true;
1200 break;
1201 }
1202 }
1203 while(locus[b] < range[b]) { /* below range? done */
1204 if(fract[b] > 1.0F) {
1205 locus[b]++;
1206 fract[b] -= 1.0F;
1207 if(locus[b] > (range[b + 3] - 2)) { /* above range? done */
1208 done = true;
1209 break;
1210 }
1211 } else {
1212 done = true;
1213 break;
1214 }
1215 }
1216
1217 }
1218 if(done)
1219 break;
1220 }
1221
1222 /* have we moved cells? */
1223
1224 if((!have_prev) || (have_prev && ((locus[0] != prev_locus[0]) ||
1225 (locus[1] != prev_locus[1]) ||
1226 (locus[2] != prev_locus[2])))) {
1227 /* above: prev_locus may be NULL, so relying upon shortcut logic eval */
1228
1229 /* stop if we hit a flagged cell (flag always in lower corner) */
1230
1231 if(*(flag + (((locus[0] - range[0]) * flag_stride[0]) +
1232 ((locus[1] - range[1]) * flag_stride[1]) +
1233 ((locus[2] - range[2]) * flag_stride[2])))) {
1234 break;
1235 }
1236
1237 }
1238
1239 {
1240 /* stop if level exceeds desired ranges */
1241 float level = FieldInterpolatef(i_data, locus[0], locus[1], locus[2],
1242 fract[0], fract[1], fract[2]);
1243 if((level < min_level) || (level > max_level))
1244 break;
1245 if(symmetry_flag) {
1246 if(symmetry_min > level)
1247 symmetry_min = level;
1248 if(symmetry_max < level)
1249 symmetry_max = level;
1250 }
1251 }
1252
1253 {
1254 /* interpolate gradient relative to grid */
1255
1256 float interp_gradient[3];
1257
1258 FieldInterpolate3f(gradients, locus, fract, interp_gradient);
1259
1260 if(length3f(interp_gradient) < min_slope) {
1261 /* if region is too flat, then bail */
1262 break;
1263 }
1264
1265 {
1266 /* add a line vertex at this point */
1267
1268 {
1269 float *f;
1270 VLACheck(i_line, float, n_line * 3 + 2);
1271 f = i_line + (n_line * 3);
1272 FieldInterpolate3f(points, locus, fract, f);
1273 n_line++;
1274 n_vert++;
1275 }
1276
1277 /* record locus for subsequent oblation */
1278
1279 if((!have_prev) || (have_prev && ((locus[0] != prev_locus[0]) ||
1280 (locus[1] != prev_locus[1]) ||
1281 (locus[2] != prev_locus[2])))) {
1282
1283 VLACheck(active_cell, int, n_active_cell * 3 + 2);
1284 {
1285 int *xrd = active_cell + (n_active_cell * 3);
1286 xrd[0] = locus[0];
1287 xrd[1] = locus[1];
1288 xrd[2] = locus[2];
1289 n_active_cell++;
1290 prev_locus = xrd; /* warning: volatile pointer */
1291 }
1292 }
1293
1294 /* adjust length of gradient vector */
1295
1296 normalize3f(interp_gradient);
1297
1298 /* make sure gradient isn't too divergent to take another step */
1299
1300 if(have_prev) {
1301 if(dot_product3f(interp_gradient, prev_grad_normal) < min_dot)
1302 break;
1303 }
1304
1305 /* take another step */
1306
1307 copy3f(interp_gradient, prev_grad_normal);
1308
1309 /* scale and flip sign */
1310
1311 if(pass) {
1312 scale3f(interp_gradient, -step_size, interp_gradient);
1313 } else {
1314 scale3f(interp_gradient, step_size, interp_gradient);
1315 }
1316
1317 /* record progress */
1318
1319 walk -= step_size;
1320
1321 /* leave if max_walk is reached */
1322
1323 if(walk < 0.0F) {
1324 break;
1325 } else {
1326 /* otherwise move */
1327 add3f(interp_gradient, fract, fract);
1328 }
1329 have_prev = true;
1330 }
1331 }
1332 } /* for */
1333
1334 if(n_vert < 2) { /* quash isolated vertices */
1335 if(n_vert) {
1336 n_line = i_num[n_seg];
1337 }
1338 } else if(n_line != i_num[n_seg]) { /* or retain count of new line segment */
1339 VLACheck(i_num, int, n_seg + 1);
1340 i_num[n_seg] = n_line - i_num[n_seg];
1341 n_seg++;
1342 i_num[n_seg] = n_line;
1343 }
1344
1345 }
1346 {
1347 int abort_segment = false;
1348 if(symmetry_flag) {
1349 if((symmetry_max * symmetry_min) >= 0.0F) /* abort if not both +/- pot. sampled */
1350 abort_segment = true;
1351 else {
1352 float symmetry_ratio = (float) (fabs(symmetry_max) / fabs(symmetry_min));
1353 if(symmetry_ratio > 1.0F)
1354 symmetry_ratio = 1.0F / symmetry_ratio;
1355 if(symmetry_ratio < symmetry) /* abort if +/- weren't close enough in magnitude */
1356 abort_segment = true;
1357 }
1358 }
1359 if((max_walk - walk) < min_walk) { /* ignore too-short segments */
1360 abort_segment = true;
1361 }
1362
1363 if(abort_segment) {
1364 n_seg = abort_n_seg;
1365 n_line = abort_n_line;
1366 i_num[n_seg] = n_line;
1367 } else {
1368 /* otherwise, keep line and oblate neighborhood */
1369
1370 int *ac = active_cell;
1371 int b;
1372 int cutoff_sq = spacing * spacing;
1373 for(b = 0; b < n_active_cell; b++) {
1374 int ii = ac[0], jj = ac[1], kk = ac[2];
1375 int i0 = ii - spacing;
1376 int j0 = jj - spacing;
1377 int k0 = kk - spacing;
1378
1379 int i1 = ii + spacing + 1;
1380 int j1 = jj + spacing + 1;
1381 int k1 = kk + spacing + 1;
1382
1383 if(i0 < range[0])
1384 i0 = range[0];
1385 if(i1 >= range[3])
1386 i1 = range[3] - 1;
1387 if(j0 < range[1])
1388 j0 = range[1];
1389 if(j1 >= range[4])
1390 j1 = range[4] - 1;
1391 if(k0 < range[2])
1392 k0 = range[2];
1393 if(k1 >= range[5])
1394 k1 = range[5] - 1;
1395
1396 {
1397 int i, j, k;
1398 int *flag1 = flag + (((i0 - range[0]) * flag_stride[0]) +
1399 ((j0 - range[1]) * flag_stride[1]) +
1400 ((k0 - range[2]) * flag_stride[2]));
1401
1402 /* highly optimized spherical flag-fill routine */
1403
1404 for(k = k0; k < k1; k++) {
1405 int *flag2 = flag1;
1406 int kk_sq = (kk - k);
1407 kk_sq = kk_sq * kk_sq;
1408
1409 for(j = j0; j < j1; j++) {
1410 int *flag3 = flag2;
1411 int jj_sq = (jj - j);
1412 jj_sq = (jj_sq * jj_sq) + kk_sq;
1413
1414 if(!(jj_sq > cutoff_sq)) {
1415 for(i = i0; i < i1; i++) {
1416 if(!*flag3) {
1417 int tot_sq = (ii - i);
1418 tot_sq = (tot_sq * tot_sq) + jj_sq;
1419 if(!(tot_sq > cutoff_sq)) {
1420 *flag3 = true;
1421 }
1422 }
1423 flag3++;
1424 } /* for i */
1425 }
1426 flag2 += flag_stride[1];
1427 } /* for j */
1428 flag1 += flag_stride[2];
1429 } /* for k */
1430 }
1431
1432 ac += 3; /* advance to next active cell */
1433 } /* for b in active_cell */
1434 }
1435 }
1436 start_locus += 3;
1437 } /* for a in range_size */
1438 }
1439
1440 }
1441 /* purge memory */
1442 VLAFreeP(active_cell);
1443 FreeP(order);
1444 FreeP(flag);
1445 }
1446
1447 /* restore modified local copies */
1448 I->NLine = n_line;
1449 I->NSeg = n_seg;
1450 return (ok);
1451 }
1452
1453
1454 /*===========================================================================*/
IsosurfDrawPoints(CIsosurf * II)1455 static int IsosurfDrawPoints(CIsosurf * II)
1456 {
1457 CIsosurf *I = II;
1458 float *a, *b;
1459 int i, j, k;
1460 int ok = true;
1461
1462 if(ok) {
1463 for(i = 0; i < (I->Max[0] - 1); i++) {
1464 for(j = 0; j < I->Max[1]; j++) {
1465 for(k = 0; k < I->Max[2]; k++) {
1466 if((I3(I->VertexCodes, i, j, k)) && (!I3(I->VertexCodes, i + 1, j, k))) {
1467 IsosurfInterpolate(I,
1468 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
1469 O3Ptr(I->Data, i, j, k, I->CurOff),
1470 O4Ptr(I->Coord, i + 1, j, k, 0, I->CurOff),
1471 O3Ptr(I->Data, i + 1, j, k, I->CurOff),
1472 &(EdgePt(I->Point, i, j, k, 0).Point[0]));
1473
1474 I->Line->check(I->NLine * 3 + 2);
1475 a = I->Line->data() + (I->NLine * 3);
1476 b = &(EdgePt(I->Point, i, j, k, 0).Point[0]);
1477 *(a++) = *(b++);
1478 *(a++) = *(b++);
1479 *a = *b;
1480 I->NLine++;
1481 } else if(!(I3(I->VertexCodes, i, j, k)) && (I3(I->VertexCodes, i + 1, j, k))) {
1482 IsosurfInterpolate(I,
1483 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
1484 O3Ptr(I->Data, i, j, k, I->CurOff),
1485 O4Ptr(I->Coord, i + 1, j, k, 0, I->CurOff),
1486 O3Ptr(I->Data, i + 1, j, k, I->CurOff),
1487 &(EdgePt(I->Point, i, j, k, 0).Point[0]));
1488
1489 I->Line->check(I->NLine * 3 + 2);
1490 a = I->Line->data() + (I->NLine * 3);
1491 b = &(EdgePt(I->Point, i, j, k, 0).Point[0]);
1492 *(a++) = *(b++);
1493 *(a++) = *(b++);
1494 *a = *b;
1495 I->NLine++;
1496 } else
1497 I4(I->ActiveEdges, i, j, k, 0) = 0;
1498 }
1499 }
1500 if(I->G->Interrupt) {
1501 ok = false;
1502 break;
1503 }
1504 }
1505 }
1506
1507 if(ok) {
1508 for(i = 0; i < I->Max[0]; i++) {
1509 for(j = 0; j < (I->Max[1] - 1); j++) {
1510 for(k = 0; k < I->Max[2]; k++) {
1511 if((I3(I->VertexCodes, i, j, k)) && (!I3(I->VertexCodes, i, j + 1, k))) {
1512 I4(I->ActiveEdges, i, j, k, 1) = 2;
1513 IsosurfInterpolate(I,
1514 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
1515 O3Ptr(I->Data, i, j, k, I->CurOff),
1516 O4Ptr(I->Coord, i, j + 1, k, 0, I->CurOff),
1517 O3Ptr(I->Data, i, j + 1, k, I->CurOff),
1518 &(EdgePt(I->Point, i, j, k, 1).Point[0]));
1519
1520 I->Line->check(I->NLine * 3 + 2);
1521 a = I->Line->data() + (I->NLine * 3);
1522 b = &(EdgePt(I->Point, i, j, k, 1).Point[0]);
1523 *(a++) = *(b++);
1524 *(a++) = *(b++);
1525 *a = *b;
1526 I->NLine++;
1527
1528 } else if(!(I3(I->VertexCodes, i, j, k)) && (I3(I->VertexCodes, i, j + 1, k))) {
1529 IsosurfInterpolate(I,
1530 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
1531 O3Ptr(I->Data, i, j, k, I->CurOff),
1532 O4Ptr(I->Coord, i, j + 1, k, 0, I->CurOff),
1533 O3Ptr(I->Data, i, j + 1, k, I->CurOff),
1534 &(EdgePt(I->Point, i, j, k, 1).Point[0]));
1535
1536 I->Line->check(I->NLine * 3 + 2);
1537 a = I->Line->data() + (I->NLine * 3);
1538 b = &(EdgePt(I->Point, i, j, k, 1).Point[0]);
1539 *(a++) = *(b++);
1540 *(a++) = *(b++);
1541 *a = *b;
1542 I->NLine++;
1543
1544 }
1545 }
1546 }
1547 if(I->G->Interrupt) {
1548 ok = false;
1549 break;
1550 }
1551 }
1552 }
1553
1554 for(i = 0; i < I->Max[0]; i++) {
1555 for(j = 0; j < I->Max[1]; j++) {
1556 for(k = 0; k < (I->Max[2] - 1); k++) {
1557 if((I3(I->VertexCodes, i, j, k)) && (!I3(I->VertexCodes, i, j, k + 1))) {
1558 IsosurfInterpolate(I,
1559 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
1560 O3Ptr(I->Data, i, j, k, I->CurOff),
1561 O4Ptr(I->Coord, i, j, k + 1, 0, I->CurOff),
1562 O3Ptr(I->Data, i, j, k + 1, I->CurOff),
1563 &(EdgePt(I->Point, i, j, k, 2).Point[0]));
1564
1565 I->Line->check(I->NLine * 3 + 2);
1566 a = I->Line->data() + (I->NLine * 3);
1567 b = &(EdgePt(I->Point, i, j, k, 2).Point[0]);
1568 *(a++) = *(b++);
1569 *(a++) = *(b++);
1570 *a = *b;
1571 I->NLine++;
1572
1573 } else if(!(I3(I->VertexCodes, i, j, k)) && (I3(I->VertexCodes, i, j, k + 1))) {
1574 IsosurfInterpolate(I,
1575 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
1576 O3Ptr(I->Data, i, j, k, I->CurOff),
1577 O4Ptr(I->Coord, i, j, k + 1, 0, I->CurOff),
1578 O3Ptr(I->Data, i, j, k + 1, I->CurOff),
1579 &(EdgePt(I->Point, i, j, k, 2).Point[0]));
1580
1581 I->Line->check(I->NLine * 3 + 2);
1582 a = I->Line->data() + (I->NLine * 3);
1583 b = &(EdgePt(I->Point, i, j, k, 2).Point[0]);
1584 *(a++) = *(b++);
1585 *(a++) = *(b++);
1586 *a = *b;
1587 I->NLine++;
1588
1589 }
1590 }
1591 if(I->G->Interrupt) {
1592 ok = false;
1593 break;
1594 }
1595 }
1596 }
1597 if(ok) {
1598 if(I->NLine != (*I->Num)[I->NSeg]) { /* any new points? */
1599 I->Num->check(I->NSeg + 1);
1600 (*I->Num)[I->NSeg] = I->NLine - (*I->Num)[I->NSeg];
1601 I->NSeg++;
1602 (*I->Num)[I->NSeg] = I->NLine;
1603 }
1604 }
1605 return (ok);
1606 }
1607
1608
1609 /*===========================================================================*/
IsosurfDrawLines(CIsosurf * II)1610 static int IsosurfDrawLines(CIsosurf * II)
1611 {
1612 CIsosurf *I = II;
1613 int c, i, j, k;
1614 float *a, *b;
1615 int ok = true;
1616 PointType *Cur, *Start, *Next;
1617 int MaxLinks, MaxL, Cnt;
1618 int NLink;
1619 #ifdef Trace
1620 int LCount = 0;
1621 #endif
1622
1623 for(i = 0; i < I->Max[0]; i++) {
1624 for(j = 0; j < I->Max[1]; j++) {
1625 for(k = 0; k < I->Max[2]; k++) {
1626 for(c = 0; c < 3; c++) {
1627 Start = EdgePtPtr(I->Point, i, j, k, c);
1628 while(Start->NLink) {
1629 Cur = Start;
1630 I->Line->check(I->NLine * 3 + 2);
1631 a = I->Line->data() + (I->NLine * 3);
1632 b = Cur->Point;
1633 *(a++) = *(b++);
1634 *(a++) = *(b++);
1635 *a = *b;
1636 I->NLine++;
1637
1638 while(Cur) {
1639 if(Cur->NLink) {
1640 Cur->NLink--;
1641 NLink = Cur->NLink;
1642 /* Choose point which has most links */
1643 MaxL = NLink;
1644 MaxLinks = Cur->Link[MaxL]->NLink;
1645 Cnt = MaxL - 1;
1646 while(Cnt >= 0) {
1647 if((Cur->Link[Cnt]->NLink) > MaxLinks) {
1648 MaxL = Cnt;
1649 MaxLinks = Cur->Link[Cnt]->NLink;
1650 }
1651 Cnt--;
1652 }
1653 Next = Cur->Link[MaxL];
1654 if(MaxL != NLink)
1655 Cur->Link[MaxL] = Cur->Link[NLink];
1656 /* Remove double link */
1657 Next->NLink--;
1658 NLink = Next->NLink;
1659 Cnt = NLink;
1660 while(Cnt >= 0) {
1661 if(Next->Link[Cnt] == Cur)
1662 break;
1663 else
1664 Cnt--;
1665 }
1666 if(Cnt >= 0) {
1667 if(Cnt != NLink)
1668 Next->Link[Cnt] = Next->Link[NLink];
1669 }
1670 #ifdef Trace
1671 else
1672 printf(" error: IsosurfDrawLines: can't find double link\n");
1673 #endif
1674
1675 Cur = Next;
1676 I->Line->check(I->NLine * 3 + 2);
1677 a = I->Line->data() + (I->NLine * 3);
1678 b = Cur->Point;
1679 *(a++) = *(b++);
1680 *(a++) = *(b++);
1681 *a = *b;
1682 I->NLine++;
1683 } else {
1684 #ifdef Trace
1685 LCount++;
1686 #endif
1687 Cur = NULL;
1688 if(I->NLine != (*I->Num)[I->NSeg]) { /* any new lines? */
1689 I->Num->check(I->NSeg + 1);
1690 (*I->Num)[I->NSeg] = I->NLine - (*I->Num)[I->NSeg];
1691 I->NSeg++;
1692 I->Num->check(I->NSeg);
1693 (*I->Num)[I->NSeg] = I->NLine;
1694 }
1695 }
1696 }
1697 }
1698 }
1699 }
1700 }
1701 if(I->G->Interrupt) {
1702 ok = false;
1703 break;
1704 }
1705 }
1706 #ifdef Trace
1707 if(ok)
1708 printf(" DrawLineCount: %i\n", LCount);
1709 #endif
1710 return (ok);
1711 }
1712
1713
1714 /*===========================================================================*/
IsosurfFindLines(CIsosurf * II)1715 static int IsosurfFindLines(CIsosurf * II)
1716 {
1717 CIsosurf *I = II;
1718 int i, j, k, ip1, jp1, kp1;
1719 int ok = true;
1720 int index, cod;
1721 int Max0m1, Max1m1, Max2m1;
1722 int skip = I->Skip;
1723 #ifdef Trace
1724 int LCount = 0;
1725 #endif
1726 PointType *p1, *p2;
1727
1728 Max0m1 = I->Max[0] - 1;
1729 Max1m1 = I->Max[1] - 1;
1730 Max2m1 = I->Max[2] - 1;
1731 for(i = 0; i < I->Max[0]; i++) {
1732 for(j = 0; j < I->Max[1]; j++) {
1733 for(k = 0; k < I->Max[2]; k++) {
1734 ip1 = i + 1;
1735 jp1 = j + 1;
1736 kp1 = k + 1;
1737 if((j < Max1m1) && (k < Max2m1) && ((!skip) || !(i % skip))) { /* i-plane */
1738 index = I4(I->ActiveEdges, i, j, k, 1) << 2;
1739 index = (index + I4(I->ActiveEdges, i, jp1, k, 2)) << 2;
1740 index = (index + I4(I->ActiveEdges, i, j, kp1, 1)) << 2;
1741 index = index + I4(I->ActiveEdges, i, j, k, 2);
1742 if(index) {
1743 cod = I->Code[index];
1744 #ifdef Trace
1745 if(index && (cod < 0))
1746 printf("IsosurfFindLines: bad index: %i \n", index);
1747 #endif
1748 while(cod > 0) {
1749 p1 = NULL;
1750 p2 = NULL;
1751 switch (cod) {
1752 case 40:
1753 case 32:
1754 cod = cod - 32;
1755 p1 = EdgePtPtr(I->Point, i, j, k, 1);
1756 p2 = EdgePtPtr(I->Point, i, j, k, 2);
1757 break;
1758 case 20:
1759 case 16:
1760 cod = cod - 16;
1761 p1 = EdgePtPtr(I->Point, i, j, k, 1);
1762 p2 = EdgePtPtr(I->Point, i, jp1, k, 2);
1763 break;
1764 case 8:
1765 cod = cod - 8;
1766 p1 = EdgePtPtr(I->Point, i, j, kp1, 1);
1767 p2 = EdgePtPtr(I->Point, i, jp1, k, 2);
1768 break;
1769 case 4:
1770 cod = cod - 4;
1771 p1 = EdgePtPtr(I->Point, i, j, kp1, 1);
1772 p2 = EdgePtPtr(I->Point, i, j, k, 2);
1773 break;
1774 case 2:
1775 cod = cod - 2;
1776 p1 = EdgePtPtr(I->Point, i, j, k, 1);
1777 p2 = EdgePtPtr(I->Point, i, j, kp1, 1);
1778 break;
1779 case 1:
1780 cod = cod - 1;
1781 p1 = EdgePtPtr(I->Point, i, j, k, 2);
1782 p2 = EdgePtPtr(I->Point, i, jp1, k, 2);
1783 break;
1784 default:
1785 cod = 0;
1786 p1 = NULL;
1787 p2 = NULL;
1788 break;
1789 }
1790 if(p1 && p2) {
1791 p1->Link[p1->NLink] = p2;
1792 p1->NLink++;
1793 p2->Link[p2->NLink] = p1;
1794 p2->NLink++;
1795 #ifdef Trace
1796 LCount++;
1797 #endif
1798 }
1799 }
1800 }
1801 }
1802 if((i < Max0m1) && (j < Max1m1) && ((!skip) || !(k % skip))) { /* k-plane */
1803 index = I4(I->ActiveEdges, i, j, k, 0) << 2;
1804 index = (index + I4(I->ActiveEdges, ip1, j, k, 1)) << 2;
1805 index = (index + I4(I->ActiveEdges, i, jp1, k, 0)) << 2;
1806 index = index + I4(I->ActiveEdges, i, j, k, 1);
1807 if(index) {
1808 cod = I->Code[index];
1809 #ifdef Trace
1810 if(index && (cod < 0))
1811 printf("IsosurfFindLines: bad index: %i \n", index);
1812 #endif
1813 while(cod > 0) {
1814 switch (cod) {
1815 case 40:
1816 case 32:
1817 cod = cod - 32;
1818 p1 = EdgePtPtr(I->Point, i, j, k, 0);
1819 p2 = EdgePtPtr(I->Point, i, j, k, 1);
1820 break;
1821 case 20:
1822 case 16:
1823 cod = cod - 16;
1824 p1 = EdgePtPtr(I->Point, i, j, k, 0);
1825 p2 = EdgePtPtr(I->Point, ip1, j, k, 1);
1826 break;
1827 case 8:
1828 cod = cod - 8;
1829 p1 = EdgePtPtr(I->Point, i, jp1, k, 0);
1830 p2 = EdgePtPtr(I->Point, ip1, j, k, 1);
1831 break;
1832 case 4:
1833 cod = cod - 4;
1834 p1 = EdgePtPtr(I->Point, i, jp1, k, 0);
1835 p2 = EdgePtPtr(I->Point, i, j, k, 1);
1836 break;
1837 case 2:
1838 cod = cod - 2;
1839 p1 = EdgePtPtr(I->Point, i, j, k, 0);
1840 p2 = EdgePtPtr(I->Point, i, jp1, k, 0);
1841 break;
1842 case 1:
1843 cod = cod - 1;
1844 p1 = EdgePtPtr(I->Point, i, j, k, 1);
1845 p2 = EdgePtPtr(I->Point, ip1, j, k, 1);
1846 break;
1847 default:
1848 cod = 0;
1849 p1 = NULL;
1850 p2 = NULL;
1851 break;
1852 }
1853 if(p1 && p2) {
1854 p1->Link[p1->NLink] = p2;
1855 p1->NLink++;
1856 p2->Link[p2->NLink] = p1;
1857 p2->NLink++;
1858 #ifdef Trace
1859 LCount++;
1860 #endif
1861 }
1862 }
1863 }
1864 }
1865 if((i < Max0m1) && (k < Max2m1) && ((!skip) || !(j % skip))) { /* j-plane */
1866 index = I4(I->ActiveEdges, i, j, k, 0) << 2;
1867 index = (index + I4(I->ActiveEdges, ip1, j, k, 2)) << 2;
1868 index = (index + I4(I->ActiveEdges, i, j, kp1, 0)) << 2;
1869 index = index + I4(I->ActiveEdges, i, j, k, 2);
1870 if(index) {
1871 cod = I->Code[index];
1872 #ifdef Trace
1873 if(index && (cod < 0))
1874 printf("IsosurfFindLines: bad index: %i \n", index);
1875 #endif
1876 while(cod > 0) {
1877 switch (cod) {
1878 case 40:
1879 case 32:
1880 cod = cod - 32;
1881 p1 = EdgePtPtr(I->Point, i, j, k, 0);
1882 p2 = EdgePtPtr(I->Point, i, j, k, 2);
1883 break;
1884 case 20:
1885 case 16:
1886 cod = cod - 16;
1887 p1 = EdgePtPtr(I->Point, i, j, k, 0);
1888 p2 = EdgePtPtr(I->Point, ip1, j, k, 2);
1889 break;
1890 case 8:
1891 cod = cod - 8;
1892 p1 = EdgePtPtr(I->Point, i, j, k + 1, 0);
1893 p2 = EdgePtPtr(I->Point, ip1, j, k, 2);
1894 break;
1895 case 4:
1896 cod = cod - 4;
1897 p1 = EdgePtPtr(I->Point, i, j, kp1, 0);
1898 p2 = EdgePtPtr(I->Point, i, j, k, 2);
1899 break;
1900 case 2:
1901 cod = cod - 2;
1902 p1 = EdgePtPtr(I->Point, i, j, k, 0);
1903 p2 = EdgePtPtr(I->Point, i, j, kp1, 0);
1904 break;
1905 case 1:
1906 cod = cod - 1;
1907 p1 = EdgePtPtr(I->Point, i, j, k, 2);
1908 p2 = EdgePtPtr(I->Point, ip1, j, k, 2);
1909 break;
1910 default:
1911 cod = 0;
1912 p1 = NULL;
1913 p2 = NULL;
1914 break;
1915 }
1916 if(p1 && p2) {
1917 p1->Link[p1->NLink] = p2;
1918 p1->NLink++;
1919 p2->Link[p2->NLink] = p1;
1920 p2->NLink++;
1921 #ifdef Trace
1922 LCount++;
1923 #endif
1924 }
1925 }
1926 }
1927 }
1928 }
1929 }
1930 if(I->G->Interrupt) {
1931 ok = false;
1932 break;
1933 }
1934
1935 }
1936 #ifdef Trace
1937 printf(" IsosurfFindLines: %i lines found\n", LCount);
1938 #endif
1939 return (ok);
1940 }
1941
1942
1943 /*===========================================================================*/
IsosurfFindActiveEdges(CIsosurf * II)1944 static int IsosurfFindActiveEdges(CIsosurf * II)
1945 {
1946 CIsosurf *I = II;
1947 int i, j, k;
1948 int ok = true;
1949 #ifdef Trace
1950 int ECount = 0;
1951 #endif
1952
1953 if(ok) {
1954 for(i = 0; i < (I->Max[0] - 1); i++) {
1955 for(j = 0; j < I->Max[1]; j++) {
1956 for(k = 0; k < I->Max[2]; k++) {
1957 if((I3(I->VertexCodes, i, j, k)) && (!I3(I->VertexCodes, i + 1, j, k))) {
1958 #ifdef Trace
1959 ECount++;
1960 #endif
1961 I4(I->ActiveEdges, i, j, k, 0) = 2;
1962 IsosurfInterpolate(I,
1963 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
1964 O3Ptr(I->Data, i, j, k, I->CurOff),
1965 O4Ptr(I->Coord, i + 1, j, k, 0, I->CurOff),
1966 O3Ptr(I->Data, i + 1, j, k, I->CurOff),
1967 &(EdgePt(I->Point, i, j, k, 0).Point[0]));
1968 } else if(!(I3(I->VertexCodes, i, j, k)) && (I3(I->VertexCodes, i + 1, j, k))) {
1969 #ifdef Trace
1970 ECount++;
1971 #endif
1972 I4(I->ActiveEdges, i, j, k, 0) = 1;
1973 IsosurfInterpolate(I,
1974 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
1975 O3Ptr(I->Data, i, j, k, I->CurOff),
1976 O4Ptr(I->Coord, i + 1, j, k, 0, I->CurOff),
1977 O3Ptr(I->Data, i + 1, j, k, I->CurOff),
1978 &(EdgePt(I->Point, i, j, k, 0).Point[0]));
1979 } else
1980 I4(I->ActiveEdges, i, j, k, 0) = 0;
1981 }
1982 }
1983 if(I->G->Interrupt) {
1984 ok = false;
1985 break;
1986 }
1987 }
1988 }
1989 if(ok) {
1990 for(i = 0; i < I->Max[0]; i++) {
1991 for(j = 0; j < (I->Max[1] - 1); j++) {
1992 for(k = 0; k < I->Max[2]; k++) {
1993 if((I3(I->VertexCodes, i, j, k)) && (!I3(I->VertexCodes, i, j + 1, k))) {
1994 #ifdef Trace
1995 ECount++;
1996 #endif
1997 I4(I->ActiveEdges, i, j, k, 1) = 2;
1998 IsosurfInterpolate(I,
1999 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
2000 O3Ptr(I->Data, i, j, k, I->CurOff),
2001 O4Ptr(I->Coord, i, j + 1, k, 0, I->CurOff),
2002 O3Ptr(I->Data, i, j + 1, k, I->CurOff),
2003 &(EdgePt(I->Point, i, j, k, 1).Point[0]));
2004 } else if(!(I3(I->VertexCodes, i, j, k)) && (I3(I->VertexCodes, i, j + 1, k))) {
2005 #ifdef Trace
2006 ECount++;
2007 #endif
2008 I4(I->ActiveEdges, i, j, k, 1) = 1;
2009 IsosurfInterpolate(I,
2010 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
2011 O3Ptr(I->Data, i, j, k, I->CurOff),
2012 O4Ptr(I->Coord, i, j + 1, k, 0, I->CurOff),
2013 O3Ptr(I->Data, i, j + 1, k, I->CurOff),
2014 &(EdgePt(I->Point, i, j, k, 1).Point[0]));
2015 } else {
2016 I4(I->ActiveEdges, i, j, k, 1) = 0;
2017 }
2018 }
2019 }
2020 if(I->G->Interrupt) {
2021 ok = false;
2022 break;
2023 }
2024 }
2025 }
2026 if(ok) {
2027 for(i = 0; i < I->Max[0]; i++) {
2028 for(j = 0; j < I->Max[1]; j++) {
2029 for(k = 0; k < (I->Max[2] - 1); k++) {
2030 if((I3(I->VertexCodes, i, j, k)) && (!I3(I->VertexCodes, i, j, k + 1))) {
2031 #ifdef Trace
2032 ECount++;
2033 #endif
2034 I4(I->ActiveEdges, i, j, k, 2) = 2;
2035 IsosurfInterpolate(I,
2036 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
2037 O3Ptr(I->Data, i, j, k, I->CurOff),
2038 O4Ptr(I->Coord, i, j, k + 1, 0, I->CurOff),
2039 O3Ptr(I->Data, i, j, k + 1, I->CurOff),
2040 &(EdgePt(I->Point, i, j, k, 2).Point[0]));
2041 } else if(!(I3(I->VertexCodes, i, j, k)) && (I3(I->VertexCodes, i, j, k + 1))) {
2042 #ifdef Trace
2043 ECount++;
2044 #endif
2045 I4(I->ActiveEdges, i, j, k, 2) = 1;
2046 IsosurfInterpolate(I,
2047 O4Ptr(I->Coord, i, j, k, 0, I->CurOff),
2048 O3Ptr(I->Data, i, j, k, I->CurOff),
2049 O4Ptr(I->Coord, i, j, k + 1, 0, I->CurOff),
2050 O3Ptr(I->Data, i, j, k + 1, I->CurOff),
2051 &(EdgePt(I->Point, i, j, k, 2).Point[0]));
2052 } else {
2053 I4(I->ActiveEdges, i, j, k, 2) = 0;
2054 }
2055 }
2056 }
2057 if(I->G->Interrupt) {
2058 ok = false;
2059 break;
2060 }
2061 }
2062 }
2063 #ifdef Trace
2064 printf(" IsosurfFindActiveEdges: %i active edges found\n", ECount);
2065 #endif
2066 return (ok);
2067 }
2068
2069
2070 /*===========================================================================*/
IsosurfCodeVertices(CIsosurf * II)2071 static int IsosurfCodeVertices(CIsosurf * II)
2072 {
2073 CIsosurf *I = II;
2074 int i, j, k;
2075 int VCount = 0;
2076 int ok = true;
2077 for(i = 0; i < I->Max[0]; i++) {
2078 for(j = 0; j < I->Max[1]; j++) {
2079 for(k = 0; k < I->Max[2]; k++) {
2080 if((O3(I->Data, i, j, k, I->CurOff) > I->Level)) {
2081 I3(I->VertexCodes, i, j, k) = 1;
2082 VCount++;
2083 } else {
2084 I3(I->VertexCodes, i, j, k) = 0;
2085 }
2086 }
2087 }
2088 if(I->G->Interrupt) {
2089 ok = false;
2090 break;
2091 }
2092 }
2093 #ifdef Trace
2094 printf(" IsosurfCodeVertices: %i of %i vertices above level\n", VCount,
2095 I->CurDim[0] * I->CurDim[1] * I->CurDim[2]);
2096 #endif
2097 if(!ok)
2098 VCount = 0;
2099 return (VCount);
2100 }
2101
2102 /*
2103 * corner: output buffer of size 8 * 3
2104 */
IsofieldGetCorners(PyMOLGlobals * G,Isofield * field,float * corner)2105 void IsofieldGetCorners(PyMOLGlobals * G, Isofield * field, float * corner) {
2106 CField* points = field->points.get();
2107 for(int a = 0; a < 8; a++) {
2108 int i = (a & 1) ? (points->dim[0] - 1) : 0;
2109 int j = (a & 2) ? (points->dim[1] - 1) : 0;
2110 int k = (a & 4) ? (points->dim[2] - 1) : 0;
2111 memcpy(corner + a * 3, F3Ptr(points, i, j, k), 3 * sizeof(float));
2112 }
2113 }
2114