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