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 -* Larry Coopet (various optimizations)
14 -*
15 -*
16 Z* ------------------------------------------------------------------- */
17 #include"os_python.h"
18 #include"os_predef.h"
19 #include"os_std.h"
20
21 #include"MemoryDebug.h"
22 #include"Err.h"
23 #include"OOMac.h"
24 #include"Map.h"
25 #include"Setting.h"
26 #include"Feedback.h"
27 #include"MemoryCache.h"
28 #include"Base.h"
29
30 static MapType *_MapNew(PyMOLGlobals * G, float range, const float *vert, int nVert,
31 const float *extent, const int *flag, int group_id, int block_id);
32
MapGetDiv(MapType * I)33 float MapGetDiv(MapType * I)
34 {
35 return I->Div;
36 }
37
~MapType()38 MapType::~MapType()
39 {
40 auto I = this;
41 {
42 CacheFreeP(I->G, I->Head, I->group_id, I->block_base + cCache_map_head_offset, false);
43 CacheFreeP(I->G, I->Link, I->group_id, I->block_base + cCache_map_link_offset, false);
44 CacheFreeP(I->G, I->EHead, I->group_id, I->block_base + cCache_map_ehead_offset,
45 false);
46 CacheFreeP(I->G, I->EMask, I->group_id, I->block_base + cCache_map_emask_offset,
47 false);
48 VLACacheFreeP(I->G, I->EList, I->group_id, I->block_base + cCache_map_elist_offset,
49 false);
50 }
51 }
52
MapCacheInit(MapCache * M,MapType * I,int group_id,int block_base)53 int MapCacheInit(MapCache * M, MapType * I, int group_id, int block_base)
54 {
55 PyMOLGlobals *G = I->G;
56 int ok = true;
57
58 M->G = G;
59 M->block_base = I->block_base;
60 M->Cache =
61 CacheCalloc(G, int, I->NVert, group_id, block_base + cCache_map_cache_offset);
62 CHECKOK(ok, M->Cache);
63 if (ok)
64 M->CacheLink =
65 CacheAlloc(G, int, I->NVert, group_id, block_base + cCache_map_cache_link_offset);
66 CHECKOK(ok, M->CacheLink);
67 M->CacheStart = -1;
68 return ok;
69 /* p=M->Cache;
70 for(a=0;a<I->NVert;a++)
71 *(p++) = 0; */
72 }
73
MapCacheReset(MapCache * M)74 void MapCacheReset(MapCache * M)
75 {
76 int i = M->CacheStart;
77 int *cachep = M->Cache;
78 int *clinkp = M->CacheLink;
79 int i1 = 0, i2 = 0, i3 = 0, i4 = 0, ii;
80 while(i >= 0) { /* believe it or not, unrolling gives us almost 10%!!! */
81 ii = clinkp[i];
82 i1 = i;
83 i = ii;
84 if(ii >= 0) {
85 ii = clinkp[ii];
86 i2 = i;
87 i = ii;
88 }
89 cachep[i1] = 0; /* this doesn't look safe, but it is. i1-i4 are always valid indices */
90 if(ii >= 0) {
91 ii = clinkp[ii];
92 i3 = i;
93 i = ii;
94 }
95 cachep[i2] = 0; /* this doesn't look safe, but it is. i1-i4 are always valid indices */
96 if(ii >= 0) {
97 ii = clinkp[ii];
98 i4 = i;
99 i = ii;
100 }
101 cachep[i3] = 0; /* this doesn't look safe, but it is. i1-i4 are always valid indices */
102 cachep[i4] = 0; /* this doesn't look safe, but it is. i1-i4 are always valid indices */
103 }
104 M->CacheStart = -1;
105 }
106
MapCacheFree(MapCache * M,int group_id,int block_base)107 void MapCacheFree(MapCache * M, int group_id, int block_base)
108 {
109 CacheFreeP(M->G, M->Cache, group_id, block_base + cCache_map_cache_offset, false);
110 CacheFreeP(M->G, M->CacheLink, group_id, block_base + cCache_map_cache_link_offset,
111 false);
112 }
113
114 #define MapSafety 0.01F
115
MapInside(MapType * I,const float * v,int * a,int * b,int * c)116 int MapInside(MapType * I, const float *v, int *a, int *b, int *c)
117 { /* special version for ray-tracing */
118 int atmp, btmp, ctmp;
119 float iDiv = I->recipDiv;
120
121 atmp = (int) ((v[0] - I->Min[0]) * iDiv) + MapBorder;
122 btmp = (int) ((v[1] - I->Min[1]) * iDiv) + MapBorder;
123 ctmp = (int) ((v[2] - I->Min[2]) * iDiv) + MapBorder;
124
125 if(atmp < I->iMin[0]) {
126 if((I->iMin[0] - atmp) > 3)
127 return (-1);
128 else
129 atmp = I->iMin[0];
130 } else if(atmp > I->iMax[0]) {
131 if((atmp - I->iMax[0]) > 3)
132 return (-1);
133 else
134 atmp = I->iMax[0];
135 }
136
137 if(btmp < I->iMin[1]) {
138 if((I->iMin[1] - btmp) > 3)
139 return (-1);
140 else
141 btmp = I->iMin[1];
142 } else if(btmp > I->iMax[1]) {
143 if((btmp - I->iMax[1]) > 3)
144 return (-1);
145 else
146 btmp = I->iMax[1];
147 }
148
149 /* printf("%d %d %d\n",ctmp,I->iMin[2],I->iMax[2]); */
150 if(ctmp < I->iMin[2]) {
151 if((I->iMin[2] - ctmp) > 3)
152 return (-1);
153 else
154 ctmp = I->iMin[2];
155 } else if(ctmp > I->iMax[2]) {
156 if((ctmp - I->iMax[2]) > 3)
157 return (0); /* just keep searching... */
158 else
159 ctmp = I->iMax[2];
160 }
161
162 /* printf("%d %d %d %d %d %d %d %d %d\n",
163 atmp,I->iMin[0],I->iMax[0],
164 btmp,I->iMin[1],I->iMax[1],
165 ctmp,I->iMin[2],I->iMax[2]);
166 */
167
168 if(!*(MapEStart(I, atmp, btmp, ctmp)))
169 return (0);
170
171 *a = atmp;
172 *b = btmp;
173 *c = ctmp;
174
175 return (1);
176 }
177
MapInsideXY(MapType * I,const float * v,int * a,int * b,int * c)178 int MapInsideXY(MapType * I, const float *v, int *a, int *b, int *c)
179 { /* special version for ray-tracing */
180 int atmp, btmp, ctmp;
181 const float iDiv = I->recipDiv;
182
183 atmp = (int) ((v[0] - I->Min[0]) * iDiv) + MapBorder;
184 btmp = (int) ((v[1] - I->Min[1]) * iDiv) + MapBorder;
185 ctmp = (int) ((v[2] - I->Min[2]) * iDiv) + MapBorder + 1;
186
187 if(atmp < I->iMin[0]) {
188 if((I->iMin[0] - atmp) > 1)
189 return (false);
190 else
191 atmp = I->iMin[0];
192 } else if(atmp > I->iMax[0]) {
193 if((atmp - I->iMax[0]) > 1)
194 return (false);
195 else
196 atmp = I->iMax[0];
197 }
198
199 if(btmp < I->iMin[1]) {
200 if((I->iMin[1] - btmp) > 1)
201 return (false);
202 else
203 btmp = I->iMin[1];
204 } else if(btmp > I->iMax[1]) {
205 if((btmp - I->iMax[1]) > 1)
206 return (false);
207 else
208 btmp = I->iMax[1];
209 }
210
211 if(!*(I->EMask + I->Dim[1] * atmp + btmp))
212 return (false);
213
214 if(ctmp < I->iMin[2])
215 ctmp = I->iMin[2];
216 else if(ctmp > I->iMax[2])
217 ctmp = I->iMax[2];
218
219 *a = atmp;
220 *b = btmp;
221 *c = ctmp;
222
223 return (true);
224 }
225
226 #define ELIST_GROW_FACTOR 3
227
MapSetupExpressXY(MapType * I,int n_vert,int negative_start)228 int MapSetupExpressXY(MapType * I, int n_vert, int negative_start)
229 { /* setup a list of XY neighbors for each square */
230 PyMOLGlobals *G = I->G;
231 int n, a, b, c, flag;
232 int d, e, i;
233 unsigned int mapSize;
234 int st, dim2;
235 int n_alloc = n_vert * 15; /* emprical est. */
236 int ok = true;
237
238 PRINTFD(G, FB_Map)
239 " MapSetupExpressXY-Debug: entered.\n" ENDFD;
240
241 mapSize = I->Dim[0] * I->Dim[1] * I->Dim[2];
242 I->EHead =
243 CacheCalloc(G, int, mapSize, I->group_id, I->block_base + cCache_map_ehead_offset);
244 CHECKOK(ok, I->EHead);
245 if (ok)
246 I->EList = (int*) VLACacheMalloc(G, n_alloc, sizeof(int), ELIST_GROW_FACTOR, 0,
247 I->group_id, I->block_base + cCache_map_elist_offset);
248 CHECKOK(ok, I->EList);
249 if (ok)
250 I->EMask = CacheCalloc(G, int, I->Dim[0] * I->Dim[1],
251 I->group_id, I->block_base + cCache_map_emask_offset);
252 CHECKOK(ok, I->EMask);
253
254 n = 1;
255 dim2 = I->Dim[2];
256
257 for(a = I->iMin[0]; ok && a <= I->iMax[0]; a++) {
258 for(b = I->iMin[1]; ok && b <= I->iMax[1]; b++) {
259 for(c = I->iMin[2]; ok && c <= I->iMax[2]; c++) { /* a better alternative exists... */
260 int *iPtr1 = (I->Head + ((a - 1) * I->D1D2) + ((b - 1) * dim2) + c);
261
262 st = n;
263 flag = false;
264
265 for(d = a - 1; d <= a + 1; d++) {
266 /*int *iPtr2 = (I->Head + (d * I->D1D2) + ((b-1)*dim2) + c); */
267 int *iPtr2 = iPtr1;
268
269 for(e = b - 1; e <= b + 1; e++) {
270 /*i = *MapFirst(I,d,e,c); */
271 i = *iPtr2;
272 if(i >= 0) {
273 flag = true;
274 while(i >= 0) {
275
276 VLACacheCheck(G, I->EList, int, n, I->group_id,
277 I->block_base + cCache_map_elist_offset);
278 CHECKOK(ok, I->EList);
279 I->EList[n] = i;
280 n++;
281 i = MapNext(I, i);
282 }
283 }
284
285 iPtr2 += dim2;
286 }
287
288 iPtr1 += I->D1D2;
289 }
290
291 if(ok && flag) {
292 *(I->EMask + I->Dim[1] * a + b) = true;
293 *(MapEStart(I, a, b, c)) = negative_start ? -st : st;
294 VLACacheCheck(G, I->EList, int, n, I->group_id,
295 I->block_base + cCache_map_elist_offset);
296 CHECKOK(ok, I->EList);
297 I->EList[n] = -1;
298 n++;
299 }
300 }
301 }
302 }
303
304 PRINTFB(G, FB_Map, FB_Blather)
305 " MapSetupExpressXY: %d rows in express table\n", n ENDFB(G);
306
307 if (ok){
308 I->NEElem = n;
309 VLACacheSize(G, I->EList, int, I->NEElem, I->group_id,
310 I->block_base + cCache_map_elist_offset);
311 CHECKOK(ok, I->EList);
312 }
313 PRINTFD(G, FB_Map)
314 " MapSetupExpressXY-Debug: leaving...\n" ENDFD;
315 return ok;
316 }
317
MapSetupExpressXYVert(MapType * I,float * vert,int n_vert,int negative_start)318 int MapSetupExpressXYVert(MapType * I, float *vert, int n_vert, int negative_start)
319 { /* setup a list of XY neighbors for each square */
320 PyMOLGlobals *G = I->G;
321 int h, n, a, b, c;
322 int j, k, dim2;
323 float *v;
324 int *eBase, *hBase;
325 int n_alloc = n_vert * 15; /* emprical est. */
326 int ok = true;
327
328 PRINTFD(G, FB_Map)
329 " MapSetupExpressXYVert-Debug: entered n_vert = %d negative_start = %d\n", n_vert,
330 negative_start ENDFD;
331
332 /*mapSize = I->Dim[0]*I->Dim[1]*I->Dim[2]; */
333 I->EHead = CacheCalloc(G, int, I->Dim[0] * I->Dim[1] * I->Dim[2],
334 I->group_id, I->block_base + cCache_map_ehead_offset);
335 CHECKOK(ok, I->EHead);
336 if (ok)
337 I->EMask = CacheCalloc(G, int, I->Dim[0] * I->Dim[1],
338 I->group_id, I->block_base + cCache_map_emask_offset);
339 CHECKOK(ok, I->EMask);
340 if (ok)
341 I->EList = (int*) VLACacheMalloc(G, n_alloc, sizeof(int), ELIST_GROW_FACTOR, 0, I->group_id, I->block_base + cCache_map_elist_offset); /* autozero */
342 CHECKOK(ok, I->EList);
343
344 n = 1;
345 v = vert;
346 dim2 = I->Dim[2];
347
348 for(h = 0; h < n_vert; h++) {
349 MapLocus(I, v, &j, &k, &c);
350
351 eBase = I->EHead + ((j - 1) * I->D1D2) + ((k - 1) * dim2) + c;
352 hBase = I->Head + (((j - 1) - 1) * I->D1D2);
353
354 for(a = j - 1; ok && a <= j + 1; a++) {
355 int *ePtr1 = eBase;
356
357 for(b = k - 1; ok && b <= k + 1; b++) {
358 if(*ePtr1 == 0) { /* if we haven't yet expanded this voxel... */
359 int st = n;
360 int flag = false;
361 int *hPtr1 = hBase + dim2 * (b - 1) + (c - 1);
362
363 int d, e, f;
364
365 for(d = a - 1; ok && d <= a + 1; d++) {
366 int *hPtr2 = hPtr1;
367 for(e = b - 1; ok && e <= b + 1; e++) {
368 int *hPtr3 = hPtr2;
369 for(f = c - 1; ok && f <= c + 1; f++) {
370 /* register int i = *MapFirst(I,d,e,f); */
371 int i = *hPtr3;
372
373 if(i > -1) {
374 flag = true;
375 while(ok && i > -1) {
376 VLACacheCheck(G, I->EList, int, n, I->group_id,
377 I->block_base + cCache_map_elist_offset);
378 CHECKOK(ok, I->EList);
379 I->EList[n] = i;
380 n++;
381 i = MapNext(I, i);
382 }
383 }
384 hPtr3++;
385 }
386 hPtr2 += dim2;
387 }
388 hPtr1 += I->D1D2;
389 }
390
391 if(flag) {
392 *(I->EMask + I->Dim[1] * a + b) = true;
393 *(MapEStart(I, a, b, c)) = negative_start ? -st : st;
394 VLACacheCheck(G, I->EList, int, n, I->group_id,
395 I->block_base + cCache_map_elist_offset);
396 CHECKOK(ok, I->EList);
397 I->EList[n] = -1;
398 n++;
399 }
400 }
401
402 ePtr1 += dim2;
403 }
404
405 eBase += I->D1D2;
406 hBase += I->D1D2;
407 }
408
409 v += 3;
410 }
411
412 PRINTFB(G, FB_Map, FB_Blather)
413 " MapSetupExpressXYVert: %d rows in express table\n", n ENDFB(G);
414
415 if (ok){
416 I->NEElem = n;
417 VLACacheSize(G, I->EList, int, I->NEElem, I->group_id,
418 I->block_base + cCache_map_elist_offset);
419 CHECKOK(ok, I->EList);
420 }
421 PRINTFD(G, FB_Map)
422 " MapSetupExpressXYVert-Debug: leaving...\n" ENDFD;
423 return ok;
424 }
425
MapSetupExpressPerp(MapType * I,const float * vert,float front,int nVertHint,int negative_start,const int * spanner)426 int MapSetupExpressPerp(MapType * I, const float *vert, float front, int nVertHint,
427 int negative_start, const int *spanner)
428 {
429 PyMOLGlobals *G = I->G;
430 int n = 0;
431 int a, b, c, i;
432
433 unsigned int mapSize;
434 int st;
435 int n_alloc = nVertHint * 15; /* emprical est. */
436 int ok = true;
437
438 int iMin0 = I->iMin[0];
439 int iMin1 = I->iMin[1];
440 int iMax0 = I->iMax[0];
441 int iMax1 = I->iMax[1];
442 float iDiv = I->recipDiv;
443 float min0 = I->Min[0] * iDiv;
444 float min1 = I->Min[1] * iDiv;
445 float base0, base1;
446 float perp_factor, premult;
447 int *emask, dim1, *link, *ptr1, *ptr2;
448
449 PRINTFD(G, FB_Map)
450 " MapSetupExpress-Debug: entered.\n" ENDFD;
451
452 mapSize = I->Dim[0] * I->Dim[1] * I->Dim[2];
453 I->EHead = CacheCalloc(G, int, mapSize,
454 I->group_id, I->block_base + cCache_map_ehead_offset);
455 CHECKOK(ok, I->EHead);
456 if (ok)
457 I->EList = (int*) VLACacheMalloc(G, n_alloc, sizeof(int), ELIST_GROW_FACTOR, 0,
458 I->group_id, I->block_base + cCache_map_elist_offset);
459 CHECKOK(ok, I->EList);
460 if (ok)
461 I->EMask = CacheCalloc(G, int, I->Dim[0] * I->Dim[1],
462 I->group_id, I->block_base + cCache_map_emask_offset);
463 CHECKOK(ok, I->EMask);
464
465 emask = I->EMask;
466 dim1 = I->Dim[1];
467 link = I->Link;
468 premult = -front * iDiv;
469
470 n = 1;
471 for(a = (iMin0 - 1); ok && a <= (iMax0 + 1); a++)
472 for(b = (iMin1 - 1); ok && b <= (iMax1 + 1); b++)
473 for(c = (I->iMin[2] - 1); ok && c <= (I->iMax[2] + 1); c++) {
474
475 int d, e, f;
476
477 /* compute a "shadow" mask for all vertices */
478
479 i = *MapFirst(I, a, b, c);
480 while(i >= 0) {
481 const float* v0 = vert + 3 * i;
482 perp_factor = premult / v0[2];
483 base0 = v0[0] * perp_factor;
484 base1 = v0[1] * perp_factor;
485
486 d = (int) (base0 - min0);
487 e = (int) (base1 - min1);
488
489 d += MapBorder;
490 e += MapBorder;
491
492 if(d < iMin0) {
493 d = iMin0;
494 } else if(d > iMax0) {
495 d = iMax0;
496 }
497
498 if(e < iMin1) {
499 e = iMin1;
500 } else if(e > iMax1) {
501 e = iMax1;
502 }
503 i = link[i];
504 ptr2 = (ptr1 = emask + dim1 * (d - 1) + (e - 1));
505 *(ptr2++) = true;
506 *(ptr2++) = true;
507 *(ptr2++) = true;
508 ptr2 = (ptr1 += dim1);
509 *(ptr2++) = true;
510 *(ptr2++) = true;
511 *(ptr2++) = true;
512 ptr2 = (ptr1 += dim1);
513 *(ptr2++) = true;
514 *(ptr2++) = true;
515 *(ptr2++) = true;
516 }
517
518 {
519 const int am1 = a - 1, ap1 = a + 1, bm1 = b - 1, bp1 = b + 1, cm1 = c - 1, cp1 =
520 c + 1;
521 const int dim2 = I->Dim[2];
522 int flag = false;
523 int *hPtr1 = I->Head + ((am1) * I->D1D2) + ((bm1) * dim2) + cm1;
524 st = n;
525 for(d = am1; ok && d <= ap1; d++) {
526 int *hPtr2 = hPtr1;
527 for(e = bm1; ok && e <= bp1; e++) {
528 int *hPtr3 = hPtr2;
529 for(f = cm1; ok && f <= cp1; f++) {
530 i = *(hPtr3++);
531 /* i=*MapFirst(I,d,e,f); */
532 if(i >= 0) {
533 flag = true;
534 while(ok && i >= 0) {
535 if((!spanner) || (f == c) || spanner[i]) {
536 /* for non-voxel-spanners, only spread in the XY plane (memory use ~ 9X instead of 27X -- a big difference!) */
537 VLACacheCheck(G, I->EList, int, n, I->group_id,
538 I->block_base + cCache_map_elist_offset);
539 CHECKOK(ok, I->EList);
540 I->EList[n] = i;
541 n++;
542 }
543 i = link[i];
544 }
545 }
546 }
547 hPtr2 += dim2;
548 }
549 hPtr1 += I->D1D2;
550 }
551 if(ok && flag) {
552 *(MapEStart(I, a, b, c)) = negative_start ? -st : st;
553 VLACacheCheck(G, I->EList, int, n, I->group_id,
554 I->block_base + cCache_map_elist_offset);
555 CHECKOK(ok, I->EList);
556 I->EList[n] = -1;
557 n++;
558 }
559 }
560 }
561 PRINTFB(G, FB_Map, FB_Blather)
562 " MapSetupExpressPerp: %d rows in express table \n", n ENDFB(G);
563 if (ok){
564 I->NEElem = n;
565 VLACacheSize(G, I->EList, int, I->NEElem, I->group_id,
566 I->block_base + cCache_map_elist_offset);
567 CHECKOK(ok, I->EList);
568 }
569 PRINTFD(G, FB_Map)
570 " MapSetupExpress-Debug: leaving...n=%d\n", n ENDFD;
571 return ok;
572 }
573
MapSetupExpress(MapType * I)574 int MapSetupExpress(MapType * I)
575 { /* setup a list of neighbors for each square */
576 PyMOLGlobals *G = I->G;
577 int n = 0;
578 int c, d, e, f, i, cm1, cp2, D1D2 = I->D1D2, D2 = I->Dim[2];
579 int mx2 = I->iMax[2];
580 int *link = I->Link;
581 int st, flag;
582 int *i_ptr3, *i_ptr4, *i_ptr5;
583 int *e_list = NULL;
584 int mx0 = I->iMax[0], mx1 = I->iMax[1], a, am1, ap2, *i_ptr1, b, bm1, bp2, *i_ptr2;
585 unsigned int mapSize;
586 int ok = true;
587
588 PRINTFD(G, FB_Map)
589 " MapSetupExpress-Debug: entered.\n" ENDFD;
590
591 mapSize = I->Dim[0] * I->Dim[1] * I->Dim[2];
592 I->EHead =
593 CacheCalloc(G, int, mapSize, group_id, I->block_base + cCache_map_ehead_offset);
594 CHECKOK(ok, I->EHead);
595 if (ok)
596 e_list = (int*) VLACacheMalloc(G, 1000, sizeof(int), 5, 0, group_id, block_offset);
597 CHECKOK(ok, e_list);
598
599 n = 1;
600 for(a = (I->iMin[0] - 1); ok && a <= mx0; a++) {
601 am1 = a - 1;
602 ap2 = a + 2;
603 i_ptr1 = I->Head + am1 * D1D2;
604 for(b = (I->iMin[1] - 1); ok && b <= mx1; b++) {
605 bm1 = b - 1;
606 bp2 = b + 2;
607 i_ptr2 = i_ptr1 + bm1 * D2;
608 for(c = (I->iMin[2] - 1); ok && c <= mx2; c++) {
609 st = n;
610 cm1 = c - 1;
611 cp2 = c + 2;
612 flag = false;
613 i_ptr5 = (i_ptr4 = (i_ptr3 = i_ptr2 + cm1));
614
615 for(d = am1; ok && d < ap2; d++) {
616 for(e = bm1; ok && e < bp2; e++) {
617 for(f = cm1; ok && f < cp2; f++) {
618 /*i=*MapFirst(I,d,e,f); */
619 if((i = *(i_ptr5++)) >= 0) {
620 flag = true;
621 do {
622 VLACacheCheck(G, e_list, int, n, group_id, block_offset);
623 CHECKOK(ok, e_list);
624 if (ok)
625 e_list[n++] = i;
626 /*i=MapNext(I,i); */
627 } while(ok && (i = link[i]) >= 0);
628 }
629 ok &= !G->Interrupt;
630 }
631 if (ok)
632 i_ptr5 = (i_ptr4 += D2);
633 }
634 if (ok)
635 i_ptr5 = (i_ptr4 = (i_ptr3 += D1D2));
636 }
637 if (ok){
638 if(flag) {
639 *(MapEStart(I, a, b, c)) = st;
640 VLACacheCheck(G, e_list, int, n, group_id, block_offset);
641 CHECKOK(ok, e_list);
642 e_list[n] = -1;
643 n++;
644 } else {
645 *(MapEStart(I, a, b, c)) = 0;
646 }
647 }
648 }
649 }
650 }
651 if (ok){
652 I->EList = e_list;
653 I->NEElem = n;
654 VLACacheSize(G, I->EList, int, I->NEElem, group_id, block_offset);
655 CHECKOK(ok, I->EList);
656 }
657 PRINTFD(G, FB_Map)
658 " MapSetupExpress-Debug: leaving...n=%d\n", n ENDFD;
659 return ok;
660 }
661
MapLocus(MapType * I,const float * v,int * a,int * b,int * c)662 void MapLocus(MapType * I, const float *v, int *a, int *b, int *c)
663 {
664 int at, bt, ct;
665 float invDiv = I->recipDiv;
666
667 at = (int) ((v[0] - I->Min[0]) * invDiv) + MapBorder;
668 bt = (int) ((v[1] - I->Min[1]) * invDiv) + MapBorder;
669 ct = (int) ((v[2] - I->Min[2]) * invDiv) + MapBorder;
670
671 /* range checking... */
672 if(at < I->iMin[0])
673 at = I->iMin[0];
674 else if(at > I->iMax[0])
675 at = I->iMax[0];
676
677 if(bt < I->iMin[1])
678 bt = I->iMin[1];
679 else if(bt > I->iMax[1])
680 bt = I->iMax[1];
681
682 if(ct < I->iMin[2])
683 ct = I->iMin[2];
684 else if(ct > I->iMax[2])
685 ct = I->iMax[2];
686
687 *a = at;
688 *b = bt;
689 *c = ct;
690 }
691
MapLocusEStart(MapType * I,const float * v)692 int *MapLocusEStart(MapType * I, const float *v)
693 {
694 int a, b, c;
695 float invDiv = I->recipDiv;
696 a = (int) (((v[0] - I->Min[0]) * invDiv) + MapBorder);
697 b = (int) (((v[1] - I->Min[1]) * invDiv) + MapBorder);
698 c = (int) (((v[2] - I->Min[2]) * invDiv) + MapBorder);
699 if(a < I->iMin[0])
700 a = I->iMin[0];
701 else if(a > I->iMax[0])
702 a = I->iMax[0];
703 if(b < I->iMin[1])
704 b = I->iMin[1];
705 else if(b > I->iMax[1])
706 b = I->iMax[1];
707 if(c < I->iMin[2])
708 c = I->iMin[2];
709 else if(c > I->iMax[2])
710 c = I->iMax[2];
711 return (I->EHead + ((a) * I->D1D2) + ((b) * I->Dim[2]) + (c));
712 }
713
MapExclLocus(MapType * I,const float * v,int * a,int * b,int * c)714 int MapExclLocus(MapType * I, const float *v, int *a, int *b, int *c)
715 {
716 float invDiv = I->recipDiv;
717
718 *a = (int) (((v[0] - I->Min[0]) * invDiv) + MapBorder);
719 if(*a < I->iMin[0])
720 return (0);
721 else if(*a > I->iMax[0])
722 return (0);
723 *b = (int) (((v[1] - I->Min[1]) * invDiv) + MapBorder);
724 if(*b < I->iMin[1])
725 return (0);
726 else if(*b > I->iMax[1])
727 return (0);
728 *c = (int) (((v[2] - I->Min[2]) * invDiv) + MapBorder);
729 if(*c < I->iMin[2])
730 return (0);
731 else if(*c > I->iMax[2])
732 return (0);
733 return (1);
734 }
735
736 /**
737 * Return EList start index for points in proximity to `v`.
738 * Return 0 if `v` is outside the grid or there are no points.
739 */
MapExclLocusEStart(MapType * map,const float * v)740 int MapExclLocusEStart(MapType* map, const float* v)
741 {
742 int h, k, l;
743 if (!MapExclLocus(map, v, &h, &k, &l))
744 return 0;
745 return *(MapEStart(map, h, k, l));
746 }
747
MapGetSeparation(PyMOLGlobals * G,float range,const float * mx,const float * mn,float * diagonal)748 float MapGetSeparation(PyMOLGlobals * G, float range, const float *mx, const float *mn,
749 float *diagonal)
750 {
751 float maxSize;
752 float size, maxSubDiv, divSize, subDiv[3];
753 float maxCubed, subDivCubed;
754 int a;
755 maxSize = SettingGetGlobal_i(G, cSetting_hash_max);
756
757 maxCubed = maxSize * maxSize * maxSize;
758
759 /* find longest axis: diagonal = max-min, for
760 * each axis and find the largest */
761 subtract3f(mx, mn, diagonal);
762 diagonal[0] = (float) fabs(diagonal[0]);
763 diagonal[1] = (float) fabs(diagonal[1]);
764 diagonal[2] = (float) fabs(diagonal[2]);
765 /* find largest */
766 size = diagonal[0];
767 if(diagonal[1] > size)
768 size = diagonal[1];
769 if(diagonal[2] > size)
770 size = diagonal[2];
771 /* err check size and diagonal */
772 if(size == 0.0) {
773 diagonal[0] = 1.0;
774 diagonal[1] = 1.0;
775 diagonal[2] = 1.0;
776 size = 1.0;
777 }
778
779 /* compute maximum number of subdivisions */
780 maxSubDiv = (float) (size / (range + MapSafety));
781 if(maxSubDiv < 1.0F)
782 maxSubDiv = 1.0F;
783
784 /* find resulting divSize */
785 divSize = size / maxSubDiv;
786 if(divSize < MapSafety)
787 divSize = MapSafety;
788 for(a = 0; a < 3; a++) {
789 subDiv[a] = (float) ((int) ((diagonal[a] / divSize) + 0.5F));
790 subDiv[a] = (subDiv[a] < 1.0F) ? 1.0F : subDiv[a];
791 }
792 subDivCubed = subDiv[0] * subDiv[1] * subDiv[2];
793
794 if(subDivCubed > maxCubed) {
795 divSize = (float) (divSize / pow(maxCubed / subDivCubed, 0.33333F));
796 } else if(subDivCubed < maxCubed) {
797 divSize = (float) (divSize * pow(subDivCubed / maxCubed, 0.33333F));
798 }
799
800 if(divSize < (range + MapSafety))
801 divSize = range + MapSafety;
802
803 if(Feedback(G, FB_Map, FB_Debugging)) {
804 PRINTF
805 " MapGetSeparation: range %8.3f divSize %8.3f size %8.3f\n", range, divSize, size
806 ENDF(G);
807 /* dump3f(mx,"mx");
808 dump3f(mn,"mn");
809 dump3f(diagonal,"diagonal");
810 printf("%8.3f\n",divSize);
811 printf("divSize %8.3f\n",divSize);
812 */
813 }
814 return (divSize);
815 }
816
MapNew(PyMOLGlobals * G,float range,const float * vert,int nVert,const float * extent)817 MapType *MapNew(PyMOLGlobals * G, float range, const float *vert, int nVert, const float *extent)
818 {
819 return (_MapNew(G, range, vert, nVert, extent, NULL, -1, 0));
820 }
821
MapNewCached(PyMOLGlobals * G,float range,const float * vert,int nVert,const float * extent,int group_id,int block_id)822 MapType *MapNewCached(PyMOLGlobals * G, float range, const float *vert, int nVert,
823 const float *extent, int group_id, int block_id)
824 {
825 return (_MapNew(G, range, vert, nVert, extent, NULL, group_id, block_id));
826 }
827
MapNewFlagged(PyMOLGlobals * G,float range,const float * vert,int nVert,const float * extent,const int * flag)828 MapType *MapNewFlagged(PyMOLGlobals * G, float range, const float *vert, int nVert,
829 const float *extent, const int *flag)
830 {
831 return (_MapNew(G, range, vert, nVert, extent, flag, -1, 0));
832 }
833
_MapNew(PyMOLGlobals * G,float range,const float * vert,int nVert,const float * extent,const int * flag,int group_id,int block_base)834 static MapType *_MapNew(PyMOLGlobals * G, float range, const float *vert, int nVert,
835 const float *extent, const int *flag, int group_id, int block_base)
836 {
837 int a, c;
838 int mapSize;
839 int h, k, l;
840 int *list;
841 const float *v;
842 int firstFlag;
843 Vector3f diagonal;
844 int ok = true;
845
846 auto I = new MapType();
847 PRINTFD(G, FB_Map)
848 " MapNew-Debug: entered.\n" ENDFD;
849 CHECKOK(ok, I);
850 if (!ok){
851 return NULL;
852 }
853 /* Initialize */
854 I->G = G;
855 I->group_id = group_id;
856 I->block_base = block_base;
857
858 /* initialize an empty cache for the map */
859 I->Link = CacheAlloc(G, int, nVert, group_id, block_base + cCache_map_link_offset);
860 CHECKOK(ok, I->Link);
861 if (!ok){
862 MapFree(I);
863 return NULL;
864 }
865 for(a = 0; a < nVert; a++)
866 I->Link[a] = -1;
867
868 /* map extents; set if valid, otherwise determine based on the flagged vertices */
869 if(extent) {
870 /* valid, so copy */
871 I->Min[0] = extent[0];
872 I->Max[0] = extent[1];
873 I->Min[1] = extent[2];
874 I->Max[1] = extent[3];
875 I->Min[2] = extent[4];
876 I->Max[2] = extent[5];
877 } else {
878 /* blank, so determine */
879 I->Min[0] = 0.0F;
880 I->Max[0] = 0.0F;
881 I->Min[1] = 0.0F;
882 I->Max[1] = 0.0F;
883 I->Min[2] = 0.0F;
884 I->Max[2] = 0.0F;
885
886 /* flag is an array of ints, one per vertex to signify inclusion for consideration in this map */
887 if(flag) {
888 firstFlag = true;
889 v = vert;
890 for(a = 0; a < nVert; a++) {
891 /* if we consider this vertex */
892 if(flag[a]) {
893 /* first-time setup*/
894 if(firstFlag) {
895 for(c = 0; c < 3; c++) {
896 I->Min[c] = v[c];
897 I->Max[c] = v[c];
898 }
899 firstFlag = false;
900 } else {
901 /* min/max extents, over all vertices */
902 for(c = 0; c < 3; c++) {
903 if(I->Min[c] > v[c])
904 I->Min[c] = v[c];
905 if(I->Max[c] < v[c])
906 I->Max[c] = v[c];
907 }
908 }
909 }
910 v += 3;
911 }
912 } else {
913 /* no flag: do all vertices in the list */
914 if(nVert) {
915 v = vert;
916 for(c = 0; c < 3; c++) {
917 I->Min[c] = v[c];
918 I->Max[c] = v[c];
919 }
920 v += 3;
921 for(a = 1; a < nVert; a++) {
922 for(c = 0; c < 3; c++) {
923 if(I->Min[c] > v[c])
924 I->Min[c] = v[c];
925 if(I->Max[c] < v[c])
926 I->Max[c] = v[c];
927 }
928 v += 3;
929 }
930 }
931 }
932 }
933
934 /* sanity check */
935 for(a = 0; a < 3; a++) {
936 if(I->Min[a] > I->Max[a]) {
937 std::swap(I->Max[a], I->Min[a]);
938 }
939
940 // empirical limit to avoid crash in PYMOL-3002
941 const float SANITY_LIMIT = 1e10;
942 if(I->Min[a] < -SANITY_LIMIT) {
943 PRINTFB(G, FB_Map, FB_Warnings)
944 " %s-Warning: clamping Min %e -> %e\n", __FUNCTION__,
945 I->Min[a], -SANITY_LIMIT ENDFB(G);
946 I->Min[a] = -SANITY_LIMIT;
947 }
948 if(I->Max[a] > SANITY_LIMIT) {
949 PRINTFB(G, FB_Map, FB_Warnings)
950 " %s-Warning: clamping Max %e -> %e\n", __FUNCTION__,
951 I->Max[a], SANITY_LIMIT ENDFB(G);
952 I->Max[a] = SANITY_LIMIT;
953 }
954 }
955
956 if(Feedback(G, FB_Map, FB_Debugging)) {
957 printf(" MapSetup: %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
958 I->Min[0], I->Min[1], I->Min[2], I->Max[0], I->Max[1], I->Max[2]);
959 }
960 /* interesting */
961 for(c = 0; c < 3; c++) {
962 I->Min[c] -= MapSafety;
963 I->Max[c] += MapSafety;
964 }
965 /* pad the boundaries by "range" */
966 if(range < 0.0) { /* negative range is a flag to expand edges using "range". */
967 range = -range;
968 for(c = 0; c < 3; c++) {
969 I->Min[c] -= range;
970 I->Max[c] += range;
971 }
972 }
973
974 /* compute final box size ..................... */
975 I->Div = MapGetSeparation(G, range, I->Max, I->Min, diagonal);
976 I->recipDiv = 1.0F / (I->Div); /* cache this */
977
978 /* add borders to avoid special edge cases */
979 I->Dim[0] = (int) ((diagonal[0] / I->Div) + 1 + (2 * MapBorder));
980 I->Dim[1] = (int) ((diagonal[1] / I->Div) + 1 + (2 * MapBorder));
981 I->Dim[2] = (int) ((diagonal[2] / I->Div) + 1 + (2 * MapBorder));
982
983 if(Feedback(G, FB_Map, FB_Debugging)) {
984 printf(" MapSetup: nVert: %d\n", nVert);
985 printf(" MapSetup: I->Div: %8.3f\n", I->Div);
986 printf(" MapSetup: %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f\n",
987 I->Min[0], I->Min[1], I->Min[2], I->Max[0], I->Max[1], I->Max[2]);
988 printf(" MapSetup: %8d %8d %8d\n", I->Dim[0], I->Dim[1], I->Dim[2]);
989 }
990
991 I->D1D2 = I->Dim[1] * I->Dim[2];
992
993 I->iMin[0] = MapBorder;
994 I->iMin[1] = MapBorder;
995 I->iMin[2] = MapBorder;
996
997 I->iMax[0] = I->Dim[0] - (1 + MapBorder);
998 I->iMax[1] = I->Dim[1] - (1 + MapBorder);
999 I->iMax[2] = I->Dim[2] - (1 + MapBorder);
1000
1001 /* compute size and allocate */
1002 mapSize = I->Dim[0] * I->Dim[1] * I->Dim[2];
1003 I->Head = CacheAlloc(G, int, mapSize, group_id, block_base + cCache_map_head_offset);
1004 CHECKOK(ok, I->Head);
1005 if (!ok){
1006 MapFree(I);
1007 return NULL;
1008 }
1009 /* initialize */
1010 /* for(a=0;a<I->Dim[0];a++)
1011 for(b=0;b<I->Dim[1];b++)
1012 for(c=0;c<I->Dim[2];c++)
1013 *(MapFirst(I,a,b,c))=-1; */
1014
1015 /* Trick for fast clearing to -1! */
1016 memset(I->Head, 0xFF, mapSize * sizeof(int));
1017
1018 I->NVert = nVert;
1019
1020 PRINTFD(G, FB_Map)
1021 " MapNew-Debug: creating 3D hash...\n" ENDFD;
1022
1023 /* create 3-D hash of the vertices */
1024 if(flag) {
1025 v = vert;
1026 for(a = 0; a < nVert; a++) {
1027 if(flag[a])
1028 if(MapExclLocus(I, v, &h, &k, &l)) {
1029 list = MapFirst(I, h, k, l);
1030 I->Link[a] = *list;
1031 *list = a; /*add to top of list */
1032 }
1033 v += 3;
1034 }
1035 } else {
1036 v = vert;
1037 for(a = 0; a < nVert; a++) {
1038 if(MapExclLocus(I, v, &h, &k, &l)) {
1039 list = MapFirst(I, h, k, l);
1040 /* printf("LINK %d %d %d %d %5.2f %5.2f %5.2f\n",a,h,k,l,v[0],v[1],v[2]); */
1041 I->Link[a] = *list;
1042 *list = a; /*add to top of list */
1043 }
1044 v += 3;
1045 }
1046 }
1047
1048 PRINTFD(G, FB_Map)
1049 " MapNew-Debug: leaving...\n" ENDFD;
1050
1051 return (I);
1052 }
1053