1 
2 /*
3 A* -------------------------------------------------------------------
4 B* This file contains source code for the PyMOL computer program
5 C* copyright 1998-2000 by Warren Lyford Delano of DeLano Scientific.
6 D* -------------------------------------------------------------------
7 E* It is unlawful to modify or remove this copyright notice.
8 F* -------------------------------------------------------------------
9 G* Please see the accompanying LICENSE file for further information.
10 H* -------------------------------------------------------------------
11 I* Additional authors of this source file include:
12 -*
13 -*
14 -*
15 Z* -------------------------------------------------------------------
16 */
17 #include"os_python.h"
18 
19 #include"os_predef.h"
20 #include"os_std.h"
21 #include"os_gl.h"
22 
23 #include"Base.h"
24 #include"MemoryDebug.h"
25 #include"OOMac.h"
26 #include"RepMesh.h"
27 #include"Map.h"
28 #include"Isosurf.h"
29 #include"Scene.h"
30 #include"Sphere.h"
31 #include"Setting.h"
32 #include"Color.h"
33 #include"main.h"
34 #include"PyMOLGlobals.h"
35 #include"Selector.h"
36 #include "ShaderMgr.h"
37 
38 typedef struct RepMesh {
39   Rep R;
40   pymol::vla<int> N;
41   int NTot;
42   pymol::vla<float> V;
43   float* VC;
44   int NDot;
45   float *Dot;
46   float Radius, Width;
47   int oneColorFlag;
48   int oneColor;
49   int *LastVisib;
50   int *LastColor;
51   float max_vdw;
52   int mesh_type;
53   CGO *shaderCGO;
54 } RepMesh;
55 
56 #include"ObjectMolecule.h"
57 
58 static
RepMeshFree(RepMesh * I)59 void RepMeshFree(RepMesh * I)
60 {
61   if (I->shaderCGO){
62     CGOFree(I->shaderCGO);
63     I->shaderCGO = 0;
64   }
65   FreeP(I->VC);
66   FreeP(I->LastColor);
67   FreeP(I->LastVisib);
68   OOFreeP(I);
69 }
70 
71 static
72 int RepMeshGetSolventDots(RepMesh * I, CoordSet * cs, float *min, float *max,
73                           float probe_radius);
74 
RepMeshCGOGenerate(RepMesh * I,RenderInfo * info)75 static int RepMeshCGOGenerate(RepMesh * I, RenderInfo * info)
76 {
77   PyMOLGlobals *G = I->R.G;
78   float *v = I->V.data();
79   float *vc = I->VC;
80   int *n = I->N.data();
81   int ok = true;
82   short use_shader;
83   short mesh_as_cylinders;
84   int c;
85   short dot_as_spheres = I->mesh_type==1 && SettingGet_i(G, I->R.cs->Setting, I->R.obj->Setting, cSetting_dot_as_spheres);
86   use_shader = SettingGetGlobal_b(G, cSetting_mesh_use_shader) &
87     SettingGetGlobal_b(G, cSetting_use_shaders);
88   mesh_as_cylinders = SettingGetGlobal_b(G, cSetting_render_as_cylinders) && SettingGetGlobal_b(G, cSetting_mesh_as_cylinders) && I->mesh_type!=1;
89 
90   ok &= CGOResetNormal(I->shaderCGO, true);
91 
92 #ifndef PURE_OPENGL_ES_2
93   int lighting =
94     SettingGet_i(G, I->R.cs->Setting, I->R.obj->Setting, cSetting_mesh_lighting);
95   if(!lighting) {
96     if(!use_shader && !info->line_lighting){
97       CGODisable(I->shaderCGO, GL_LIGHTING);
98     }
99   }
100 #endif
101   if (ok){
102     switch (I->mesh_type) {
103     case 0:
104       ok &= CGOSpecial(I->shaderCGO, LINEWIDTH_DYNAMIC_MESH);
105       break;
106     case 1:
107       ok &= CGOSpecial(I->shaderCGO, POINTSIZE_DYNAMIC_DOT_WIDTH);
108       break;
109     }
110   }
111 
112   ok &= CGOResetNormal(I->shaderCGO, false);
113 
114   switch (I->mesh_type) {
115   case 0:
116     if(n) {
117       if (ok){
118 	if(I->oneColorFlag) {
119 	  while(ok && *n) {
120 	    ok &= CGOColorv(I->shaderCGO, ColorGet(G, I->oneColor));
121 	    if (ok){
122 	      c = *(n++);
123 	      if (mesh_as_cylinders){
124 		float *origin, axis[3];
125 		for (; ok && c > 1; --c) {
126 		  origin = v;
127 		  v += 3;
128 		  axis[0] = v[0] - origin[0];
129 		  axis[1] = v[1] - origin[1];
130 		  axis[2] = v[2] - origin[2];
131 		  ok &= (bool)I->shaderCGO->add<cgo::draw::shadercylinder>(origin, axis, 1.f, 15);
132 		}
133 		if (c == 1) {
134 		  v += 3;
135 		}
136 	      } else {
137 		ok &= CGOBegin(I->shaderCGO, GL_LINE_STRIP);
138 		ok &= CGONormal(I->shaderCGO, 0.f,0.f,1.f);
139 		while(ok && c--) {
140 		  ok &= CGOVertexv(I->shaderCGO, v);
141 		  v += 3;
142 		}
143 		if (ok)
144 		  ok &= CGOEnd(I->shaderCGO);
145 	      }
146 	    }
147 	  }
148 	} else {
149 	  while(ok && *n) {
150 	    c = *(n++);
151 	    if (mesh_as_cylinders){
152 	      float *origin, axis[3], *color;
153 	      for (; ok && c > 1; --c) {
154 		ok &= CGOColorv(I->shaderCGO, vc);
155 		color = vc;
156 		origin = v;
157 		vc += 3;
158 		v += 3;
159 		axis[0] = v[0] - origin[0];
160 		axis[1] = v[1] - origin[1];
161 		axis[2] = v[2] - origin[2];
162 		if (*(color) != *(vc) || *(color+1) != *(vc+1) || *(color+2) != *(vc+2)){
163                   ok &= (bool)I->shaderCGO->add<cgo::draw::shadercylinder2ndcolor>(I->shaderCGO, origin, axis, 1.f, 15, vc);
164 		} else {
165                   ok &= (bool)I->shaderCGO->add<cgo::draw::shadercylinder>(origin, axis, 1.f, 15);
166 		}
167 	      }
168 	      if (c == 1) {
169 		v += 3;
170 		vc += 3;
171 	      }
172 	    } else {
173 	      ok &= CGOBegin(I->shaderCGO, GL_LINE_STRIP);
174 	      ok &= CGONormal(I->shaderCGO, 0.f,0.f,1.f);
175 	      while(ok && c--) {
176 		ok &= CGOColorv(I->shaderCGO, vc);
177 		vc += 3;
178 		if (ok){
179 		  ok &= CGOVertexv(I->shaderCGO, v);
180 		  v += 3;
181 		}
182 	      }
183 	      if (ok)
184 		ok &= CGOEnd(I->shaderCGO);
185 	    }
186 	  }
187 	}
188       }
189     }
190     break;
191   case 1:
192 #ifdef PURE_OPENGL_ES_2
193     /* TODO */
194 #else
195     glPointSize(SettingGet_f
196 		(G, I->R.cs->Setting, I->R.obj->Setting, cSetting_dot_width));
197 #endif
198     if(ok && n) {
199       if(I->oneColorFlag) {
200 	while(ok && *n) {
201 	  ok &= CGOColorv(I->shaderCGO, ColorGet(G, I->oneColor));
202 	  c = *(n++);
203 	  if (ok && !dot_as_spheres)
204 	    ok &= CGOBegin(I->shaderCGO, GL_POINTS);
205 	  while(ok && c--) {
206 	    if (dot_as_spheres)
207 	      ok &= CGOSphere(I->shaderCGO, v, 1.f);
208 	    else
209 	      ok &= CGOVertexv(I->shaderCGO, v);
210 	    v += 3;
211 	  }
212 	  if (ok && !dot_as_spheres)
213 	    ok &= CGOEnd(I->shaderCGO);
214 	}
215       } else {
216 	while(ok && *n) {
217 	  c = *(n++);
218 	  if (!dot_as_spheres)
219 	    ok &= CGOBegin(I->shaderCGO, GL_POINTS);
220 	  while(ok && c--) {
221 	    ok &= CGOColorv(I->shaderCGO, vc);
222 	    vc += 3;
223 	    if (ok){
224 	      if (dot_as_spheres)
225 		ok &= CGOSphere(I->shaderCGO, v, 1.f);
226 	      else
227 		ok &= CGOVertexv(I->shaderCGO, v);
228 	    }
229 	    v += 3;
230 	  }
231 	  if (!dot_as_spheres)
232 	    ok &= CGOEnd(I->shaderCGO);
233 	}
234       }
235     }
236     break;
237   }
238 
239   if (use_shader) {
240     if (ok){
241       CGO *convertcgo = NULL;
242       ok &= CGOStop(I->shaderCGO);
243       if (ok)
244 	convertcgo = CGOCombineBeginEnd(I->shaderCGO, 0);
245       CHECKOK(ok, convertcgo);
246       CGOFree(I->shaderCGO);
247       I->shaderCGO = convertcgo;
248       convertcgo = NULL;
249       if (ok){
250 	if (dot_as_spheres){
251 	  CGO *tmpCGO = CGONew(G);
252 	  if (ok) ok &= CGOEnable(tmpCGO, GL_SPHERE_SHADER);
253 	  if (ok) ok &= CGOSpecial(tmpCGO, DOTSIZE_WITH_SPHERESCALE);
254 	  convertcgo = CGOOptimizeSpheresToVBONonIndexedNoShader(I->shaderCGO,
255               CGO_BOUNDING_BOX_SZ + fsizeof<cgo::draw::sphere_buffers>() + 2);
256 	  if (ok) ok &= CGOAppendNoStop(tmpCGO, convertcgo);
257 	  if (ok) ok &= CGODisable(tmpCGO, GL_SPHERE_SHADER);
258 	  if (ok) ok &= CGOStop(tmpCGO);
259 	  CGOFreeWithoutVBOs(convertcgo);
260 	  convertcgo = tmpCGO;
261 	  convertcgo->use_shader = true;
262 	} else {
263 	  if (mesh_as_cylinders){
264 	    CGO *tmpCGO = CGONew(G);
265             ok &= CGOEnable(tmpCGO, GL_CYLINDER_SHADER);
266 	    if (ok) ok &= CGOSpecial(tmpCGO, MESH_WIDTH_FOR_SURFACES);
267             convertcgo = CGOConvertShaderCylindersToCylinderShader(I->shaderCGO,  tmpCGO);
268 	    if (ok) ok &= CGOAppendNoStop(tmpCGO, convertcgo);
269 	    if (ok) ok &= CGODisable(tmpCGO, GL_CYLINDER_SHADER);
270 	    if (ok) ok &= CGOStop(tmpCGO);
271 	    CGOFreeWithoutVBOs(convertcgo);
272 	    convertcgo = tmpCGO;
273 	    convertcgo->use_shader = convertcgo->has_draw_cylinder_buffers = true;
274 	  } else {
275 	    CGO *tmpCGO = CGONew(G);
276 	    if (ok) ok &= CGOEnable(tmpCGO, GL_DEFAULT_SHADER);
277 	    convertcgo = CGOOptimizeToVBONotIndexedNoShader(I->shaderCGO, 0);
278 	    if (ok) ok &= CGOAppendNoStop(tmpCGO, convertcgo);
279 	    if (ok) ok &= CGODisable(tmpCGO, GL_DEFAULT_SHADER);
280 	    if (ok) ok &= CGOStop(tmpCGO);
281 	    CGOFreeWithoutVBOs(convertcgo);
282 	    convertcgo = tmpCGO;
283 	  }
284 	}
285 	CHECKOK(ok, convertcgo);
286       }
287       if (convertcgo){
288 	convertcgo->use_shader = true;
289 	CGOFree(I->shaderCGO);
290 	I->shaderCGO = convertcgo;
291 	convertcgo = NULL;
292       }
293     }
294   }
295   return ok;
296 }
297 
298 
RepMeshRender(RepMesh * I,RenderInfo * info)299 static void RepMeshRender(RepMesh * I, RenderInfo * info)
300 {
301   CRay *ray = info->ray;
302   auto pick = info->pick;
303   PyMOLGlobals *G = I->R.G;
304   float *v = I->V.data();
305   float *vc = I->VC;
306   int *n = I->N.data();
307   int c;
308   const float *col = NULL;
309   float line_width = SceneGetDynamicLineWidth(info, I->Width);
310   short dot_as_spheres = I->mesh_type==1 && SettingGet_i(G, I->R.cs->Setting, I->R.obj->Setting, cSetting_dot_as_spheres);
311   int ok = true;
312 
313   if(ray) {
314     if(n) {
315       float radius;
316 
317       if(I->Radius <= 0.0F) {
318         radius = ray->PixelRadius * line_width / 2.0F;
319       } else {
320         radius = I->Radius;
321       }
322       /* looks like were missing some code here --
323 
324          what about mesh_type?
325 
326        */
327 
328       if(I->oneColorFlag)
329         col = ColorGet(G, I->oneColor);
330       ray->color3fv(ColorGet(G, I->R.obj->Color));
331       switch (I->mesh_type) {
332       case 0:
333         while(ok && *n) {
334           c = *(n++);
335           if(c--) {
336             vc += 3;
337             v += 3;
338             if(I->oneColorFlag) {
339               while(ok && c--) {
340                 ok &= ray->sausage3fv(v - 3, v, radius, col, col);
341                 v += 3;
342                 vc += 3;
343               }
344             } else {
345               while(ok && c--) {
346                 ok &= ray->sausage3fv(v - 3, v, radius, vc - 3, vc);
347                 v += 3;
348                 vc += 3;
349               }
350             }
351           }
352         }
353       case 1:
354         while(ok && *n) {
355           c = *(n++);
356           if(I->oneColorFlag) {
357             ray->color3fv(col);
358             while(ok && c--) {
359               ok &= ray->sphere3fv(v, radius);
360               v += 3;
361               vc += 3;
362             }
363           } else {
364             while(ok && c--) {
365               ray->color3fv(vc);
366               ok &= ray->sphere3fv(v, radius);
367               v += 3;
368               vc += 3;
369             }
370           }
371         }
372         break;
373       }
374     }
375   } else if(G->HaveGUI && G->ValidContext) {
376     if(pick) {
377       /* no picking meshes */
378     } else {
379       short use_shader, generate_shader_cgo = 0;
380       short mesh_as_cylinders ;
381       use_shader = SettingGetGlobal_b(G, cSetting_mesh_use_shader) &
382                    SettingGetGlobal_b(G, cSetting_use_shaders);
383       mesh_as_cylinders = SettingGetGlobal_b(G, cSetting_render_as_cylinders) && SettingGetGlobal_b(G, cSetting_mesh_as_cylinders) && I->mesh_type!=1;
384 
385       if (I->shaderCGO && !use_shader){
386 	CGOFree(I->shaderCGO);
387 	I->shaderCGO = 0;
388       }
389       if (I->shaderCGO && ((mesh_as_cylinders ^ I->shaderCGO->has_draw_cylinder_buffers) ||
390 			   (dot_as_spheres ^ I->shaderCGO->has_draw_sphere_buffers))){
391 	CGOFree(I->shaderCGO);
392 	I->shaderCGO = 0;
393       }
394 
395       if (use_shader){
396 	if (!I->shaderCGO){
397 	  I->shaderCGO = CGONew(G);
398 	  CHECKOK(ok, I->shaderCGO);
399 	  if (ok)
400 	    I->shaderCGO->use_shader = true;
401 	  generate_shader_cgo = 1;
402 	} else if (ok) {
403 	  CGORenderGL(I->shaderCGO, NULL, NULL, NULL, info, &I->R);
404 	  return;
405 	}
406       }
407 
408       if (ok){
409 	if (generate_shader_cgo){
410 	  ok &= RepMeshCGOGenerate(I, info);
411 	}
412       }
413 
414       int lighting =
415         SettingGet_i(G, I->R.cs->Setting, I->R.obj->Setting, cSetting_mesh_lighting);
416       if(!lighting) {
417         if(!info->line_lighting){
418 	  if (!use_shader && !generate_shader_cgo){
419 	    glDisable(GL_LIGHTING);
420 	  }
421 	}
422       }
423       if (!generate_shader_cgo){
424 	switch (I->mesh_type) {
425 	case 0:
426 	  if(info->width_scale_flag)
427 	    glLineWidth(line_width * info->width_scale);
428 	  else
429 	    glLineWidth(line_width);
430 	  break;
431 	case 1:
432 	  if(info->width_scale_flag)
433 	    glPointSize(SettingGet_f
434 			(G, I->R.cs->Setting, I->R.obj->Setting,
435 			 cSetting_dot_width) * info->width_scale);
436 	  else
437 	    glPointSize(SettingGet_f
438 			(G, I->R.cs->Setting, I->R.obj->Setting, cSetting_dot_width));
439 	  break;
440 	}
441       }
442 
443       if (ok){
444 	if (!generate_shader_cgo){
445 	  SceneResetNormal(G, false);
446 	}
447       }
448 
449       switch (I->mesh_type) {
450       case 0:
451 	if(n) {
452 	  if (!generate_shader_cgo){
453 	    if(I->oneColorFlag) {
454 	      while(*n) {
455 		glColor3fv(ColorGet(G, I->oneColor));
456 		c = *(n++);
457 		glBegin(GL_LINE_STRIP);
458 		while(c--) {
459 		  glVertex3fv(v);
460 		  v += 3;
461 		}
462 		glEnd();
463 	      }
464 	    } else {
465 	      while(*n) {
466 		c = *(n++);
467 		glBegin(GL_LINE_STRIP);
468 		while(c--) {
469 		  glColor3fv(vc);
470 		  vc += 3;
471 		  glVertex3fv(v);
472 		  v += 3;
473 		}
474 		glEnd();
475 	      }
476 	    }
477 	  }
478 	}
479 	break;
480       case 1:
481 	glPointSize(SettingGet_f
482 		    (G, I->R.cs->Setting, I->R.obj->Setting, cSetting_dot_width));
483 	if(ok && n) {
484 	  if (!generate_shader_cgo){
485 	    if(I->oneColorFlag) {
486 	      while(*n) {
487 		glColor3fv(ColorGet(G, I->oneColor));
488 		c = *(n++);
489 		glBegin(GL_POINTS);
490 		while(c--) {
491 		  glVertex3fv(v);
492 		  v += 3;
493 		}
494 		glEnd();
495 	      }
496 	    } else {
497 	      while(*n) {
498 		c = *(n++);
499 		glBegin(GL_POINTS);
500 		while(c--) {
501 		  glColor3fv(vc);
502 		  vc += 3;
503 		  glVertex3fv(v);
504 		  v += 3;
505 		}
506 		glEnd();
507 	      }
508 	    }
509 	  }
510 	}
511 	break;
512       }
513 
514       /* end of rendering, if using shaders, then render CGO */
515       if (use_shader) {
516 	if (ok){
517 	  {
518 	    const float *color;
519 	    color = ColorGet(G, I->R.obj->Color);
520 	    CGORenderGL(I->shaderCGO, color, NULL, NULL, info, &I->R);
521 	  }
522 	}
523       }
524 
525 #ifndef PURE_OPENGL_ES_2
526       if(!use_shader && !lighting)
527         glEnable(GL_LIGHTING);
528 #endif
529     }
530   }
531   if (!ok){
532     CGOFree(I->shaderCGO);
533     I->R.fInvalidate(&I->R, I->R.cs, cRepInvPurge);
534     I->R.cs->Active[cRepMesh] = false;
535   }
536 }
537 
538 static
RepMeshSameVis(RepMesh * I,CoordSet * cs)539 int RepMeshSameVis(RepMesh * I, CoordSet * cs)
540 {
541   int *lv, *lc;
542   int a;
543   AtomInfoType *ai;
544 
545   lv = I->LastVisib;
546   lc = I->LastColor;
547 
548   for(a = 0; a < cs->NIndex; a++) {
549     ai = cs->getAtomInfo(a);
550     if(*(lv++) != GET_BIT(ai->visRep, cRepMesh)) {
551       return false;
552     }
553     if(*(lc++) != ai->color) {
554       return false;
555     }
556   }
557   return true;
558 }
559 
560 static
RepMeshColor(RepMesh * I,CoordSet * cs)561 void RepMeshColor(RepMesh * I, CoordSet * cs)
562 {
563   PyMOLGlobals *G = cs->G;
564   MapType *map;
565   int a, i0, i, j, h, k, l, c1;
566   float *v0, *vc;
567   const float *c0;
568   int *lv, *lc;
569   int first_color;
570   ObjectMolecule *obj;
571   float probe_radius;
572   float dist, minDist;
573   int inclH;
574   int cullByFlag = false;
575   int mesh_mode;
576   int mesh_color;
577   AtomInfoType *ai2;
578   int state = I->R.context.state;
579 
580   obj = cs->Obj;
581 
582   probe_radius = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_solvent_radius);
583   mesh_color = SettingGet_color(G, cs->Setting, obj->Setting, cSetting_mesh_color);
584   mesh_mode = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_mode);
585   cullByFlag = (mesh_mode == cRepMesh_by_flags);
586   inclH = !(mesh_mode == cRepMesh_heavy_atoms);
587 
588   if(!I->LastVisib)
589     I->LastVisib = pymol::malloc<int>(cs->NIndex);
590   if(!I->LastColor)
591     I->LastColor = pymol::malloc<int>(cs->NIndex);
592   lv = I->LastVisib;
593   lc = I->LastColor;
594   for(a = 0; a < cs->NIndex; a++) {
595     ai2 = cs->getAtomInfo(a);
596     *(lv++) = GET_BIT(ai2->visRep, cRepMesh);
597     *(lc++) = ai2->color;
598   }
599 
600   if(I->mesh_type != 1) {
601     I->Width = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_mesh_width);
602     I->Radius = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_mesh_radius);
603   } else {
604     I->Width = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_dot_width);
605     I->Radius = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_dot_radius);
606   }
607 
608   if(I->NTot) {
609     obj = cs->Obj;
610     if(ColorCheckRamped(G, mesh_color)) {
611       I->oneColorFlag = false;
612     } else {
613       I->oneColorFlag = true;
614     }
615     first_color = -1;
616     if(!I->VC)
617       I->VC = pymol::malloc<float>(3 * I->NTot);
618     vc = I->VC;
619     /* now, assign colors to each point */
620     map = MapNew(G, I->max_vdw + probe_radius, cs->Coord, cs->NIndex, NULL);
621     if(map) {
622       MapSetupExpress(map);
623       for(a = 0; a < I->NTot; a++) {
624         AtomInfoType *ai0 = NULL;
625         c1 = 1;
626         minDist = FLT_MAX;
627         i0 = -1;
628         v0 = I->V.data() + 3 * a;
629         MapLocus(map, v0, &h, &k, &l);
630 
631         i = *(MapEStart(map, h, k, l));
632         if(i) {
633           j = map->EList[i++];
634           while(j >= 0) {
635             ai2 = obj->AtomInfo + cs->IdxToAtm[j];
636             if((inclH || (!ai2->isHydrogen())) &&
637                ((!cullByFlag) || (!(ai2->flags & cAtomFlag_ignore)))) {
638               dist = (float) diff3f(v0, cs->coordPtr(j)) - ai2->vdw;
639               if(dist < minDist) {
640                 i0 = j;
641                 ai0 = ai2;
642                 minDist = dist;
643               }
644             }
645             j = map->EList[i++];
646           }
647         }
648 
649         if(i0 >= 0) {
650           int at_mesh_color = AtomSettingGetWD(G, ai0, cSetting_mesh_color, mesh_color);
651 
652           if(at_mesh_color != -1) {
653             c1 = at_mesh_color;
654           } else {
655             c1 = ai0->color;
656           }
657           if(I->oneColorFlag) {
658             if(first_color >= 0) {
659               if(first_color != c1)
660                 I->oneColorFlag = false;
661             } else
662               first_color = c1;
663           }
664         }
665         /*
666            if(ColorCheckRamped(G,mesh_color)) {
667            c1 = mesh_color;
668            }
669          */
670 
671         if(ColorCheckRamped(G, c1)) {
672           I->oneColorFlag = false;
673           ColorGetRamped(G, c1, v0, vc, state);
674           vc += 3;
675         } else {
676           c0 = ColorGet(G, c1);
677           *(vc++) = *(c0++);
678           *(vc++) = *(c0++);
679           *(vc++) = *(c0++);
680         }
681       }
682       MapFree(map);
683     }
684     if(I->oneColorFlag) {
685       I->oneColor = first_color;
686     }
687     if (I->shaderCGO){
688       CGOFree(I->shaderCGO);
689       I->shaderCGO = 0;
690     }
691   }
692   /*
693      if(mesh_color>=0) {
694      I->oneColorFlag=1;
695      I->oneColor=mesh_color;
696      }
697    */
698 
699 }
700 
RepMeshNew(CoordSet * cs,int state)701 Rep *RepMeshNew(CoordSet * cs, int state)
702 {
703   PyMOLGlobals *G = cs->G;
704   ObjectMolecule *obj;
705   CoordSet *ccs;
706   int a, b, c, d, h, k, l, *n;
707   MapType *map = NULL, *smap = NULL;
708   /* grid */
709   Vector3f minE, maxE, sizeE;
710   float size;
711   int dims[3];
712   float gridSize;
713   Isofield *field;
714   Vector3f point;
715   float vLen, pVal, vdw;
716   int aNear;
717   float aLen;
718   int cur;
719   int meshFlag = false;
720   int inSolvFlag;
721   float probe_radius, probe_radius2;
722   float min_spacing;
723   int visFlag;
724   int mesh_mode;
725   int cullByFlag;
726   int inclH;
727   int solv_acc;
728   int mesh_type;
729   int mesh_skip;
730   int ok = true;
731 
732   AtomInfoType *ai1;
733 
734   OOAlloc(G, RepMesh);
735   CHECKOK (ok, I);
736   PRINTFD(G, FB_RepMesh)
737   " RepMeshNew-DEBUG: entered with coord-set %p\n", (void *) cs ENDFD;
738   if (ok){
739     obj = cs->Obj;
740     I->R.context.object = obj;
741     I->R.context.state = state;
742   }
743   if (ok){
744     probe_radius = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_solvent_radius);
745     probe_radius2 = probe_radius * probe_radius;
746     solv_acc = (SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_solvent));
747     mesh_type = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_type);
748     mesh_skip = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_skip);
749 
750     mesh_mode = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_mode);
751     cullByFlag = (mesh_mode == cRepMesh_by_flags);
752     inclH = !(mesh_mode == cRepMesh_heavy_atoms);
753   }
754   visFlag = false;
755   if(ok && (obj->RepVisCache & cRepMeshBit)){
756     for(a = 0; a < cs->NIndex; a++) {
757       ai1 = obj->AtomInfo + cs->IdxToAtm[a];
758       if((ai1->visRep & cRepMeshBit) &&
759          (inclH || (!ai1->isHydrogen())) &&
760          ((!cullByFlag) || (!(ai1->flags & (cAtomFlag_exfoliate | cAtomFlag_ignore))))) {
761         visFlag = true;
762         break;
763       }
764     }
765   }
766   if(!ok || !visFlag) {
767     OOFreeP(I);
768     return (NULL);              /* skip if no dots are visible */
769   }
770 
771   I->max_vdw = ObjectMoleculeGetMaxVDW(obj) + solv_acc * probe_radius;
772 
773   RepInit(G, &I->R);
774 
775   min_spacing = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_min_mesh_spacing);
776 
777   I->N = NULL;
778   I->NTot = 0;
779   I->V = NULL;
780   I->VC = NULL;
781   I->NDot = 0;
782   I->Dot = NULL;
783   I->R.fRender = (void (*)(struct Rep *, RenderInfo *)) RepMeshRender;
784   I->R.fFree = (void (*)(struct Rep *)) RepMeshFree;
785   I->R.obj = (CObject *) cs->Obj;
786   I->R.cs = cs;
787   I->R.fRecolor = (void (*)(struct Rep *, struct CoordSet *)) RepMeshColor;
788   I->R.fSameVis = (int (*)(struct Rep *, struct CoordSet *)) RepMeshSameVis;
789   I->LastVisib = NULL;
790   I->LastColor = NULL;
791   I->mesh_type = mesh_type;
792   I->Radius = SettingGet_f(G, cs->Setting, obj->Setting, cSetting_mesh_radius);
793   I->shaderCGO = 0;
794 
795   meshFlag = true;
796 
797   if(meshFlag) {
798     float trim_cutoff =
799       SettingGet_f(G, cs->Setting, obj->Setting, cSetting_mesh_cutoff);
800     int trim_flag = false;
801     float *trim_vla = NULL;
802     MapType *trim_map = NULL;
803 
804     int carve_state = 0;
805     int carve_flag = false;
806     float carve_cutoff =
807       SettingGet_f(G, cs->Setting, obj->Setting, cSetting_mesh_carve_cutoff);
808     const char *carve_selection = NULL;
809     float *carve_vla = NULL;
810     MapType *carve_map = NULL;
811 
812     int clear_state = 0;
813     int clear_flag = false;
814     float clear_cutoff =
815       SettingGet_f(G, cs->Setting, obj->Setting, cSetting_mesh_clear_cutoff);
816     const char *clear_selection = NULL;
817     float *clear_vla = NULL;
818     MapType *clear_map = NULL;
819 
820     int mesh_max = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_grid_max);
821     if(mesh_max < 1)
822       mesh_max = 1;
823 
824     if(trim_cutoff < 0.0F) {
825       trim_cutoff = I->max_vdw + probe_radius;
826     }
827 
828     if(trim_cutoff > 0.0F) {
829       int nc = 0;
830       trim_vla = VLAlloc(float, cs->NIndex * 3);
831       CHECKOK(ok, trim_vla);
832       for(c = 0; ok && c < cs->NIndex; c++) {
833         ai1 = obj->AtomInfo + cs->IdxToAtm[c];
834         if((ai1->visRep & cRepMeshBit) &&
835            (inclH || (!ai1->isHydrogen())) &&
836            ((!cullByFlag) ||
837             (!(ai1->flags & (cAtomFlag_ignore | cAtomFlag_exfoliate))))) {
838           VLACheck(trim_vla, float, nc * 3 + 2);
839 	  CHECKOK(ok, trim_vla);
840           if (ok) {
841             const float* src = cs->coordPtr(c);
842             float *dst = trim_vla + 3 * nc;
843             *(dst++) = *(src++);
844             *(dst++) = *(src++);
845             *(dst++) = *(src++);
846             nc++;
847           }
848         }
849       }
850       if(ok && nc) {
851         trim_map = MapNew(G, trim_cutoff, trim_vla, nc, NULL);
852 	CHECKOK(ok, trim_map);
853         if(ok) {
854           ok &= MapSetupExpress(trim_map);
855           trim_flag = true;
856         }
857       }
858     }
859 
860     if(ok && carve_cutoff > 0.0F) {
861       carve_state =
862         SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_carve_state) - 1;
863       carve_selection =
864         SettingGet_s(G, cs->Setting, obj->Setting, cSetting_mesh_carve_selection);
865       if(carve_selection)
866         carve_map = SelectorGetSpacialMapFromSeleCoord(G,
867                                                        SelectorIndexByName(G,
868                                                                            carve_selection),
869                                                        carve_state, carve_cutoff,
870                                                        &carve_vla);
871       if(carve_map) {
872         ok &= MapSetupExpress(carve_map);
873         carve_flag = true;
874       }
875     }
876     if(ok && clear_cutoff > 0.0F) {
877       clear_state =
878         SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_clear_state) - 1;
879       clear_selection =
880         SettingGet_s(G, cs->Setting, obj->Setting, cSetting_mesh_clear_selection);
881       if(clear_selection)
882         clear_map = SelectorGetSpacialMapFromSeleCoord(G,
883                                                        SelectorIndexByName(G,
884                                                                            clear_selection),
885                                                        clear_state, clear_cutoff,
886                                                        &clear_vla);
887       if(clear_map) {
888         ok &= MapSetupExpress(clear_map);
889         clear_flag = true;
890       }
891     }
892 
893     if (ok)
894       I->V = pymol::vla_take_ownership((float*)VLAMalloc(1000, sizeof(float), 9, false));
895     CHECKOK(ok, I->V);
896     if (ok)
897       I->N = pymol::vla_take_ownership((int*)VLAMalloc(100, sizeof(int), 9, false));
898     CHECKOK(ok, I->N);
899     if (ok)
900       I->N[0] = 0;
901 
902     if (ok){
903       for(c = 0; c < 3; c++) {
904 	minE[c] = FLT_MAX;
905 	maxE[c] = -FLT_MAX;
906       }
907       for(b = 0; b < obj->NCSet; b++) {
908 	ccs = obj->CSet[b];
909 	if(ccs) {
910 	  for(c = 0; c < ccs->NIndex; c++) {
911 	    ai1 = obj->AtomInfo + ccs->IdxToAtm[c];       /* WLD fixed 011218 */
912 	    if((ai1->visRep & cRepMeshBit) &&
913 	       (inclH || (!ai1->isHydrogen())) &&
914 	       ((!cullByFlag) ||
915 		(!(ai1->flags & (cAtomFlag_ignore | cAtomFlag_exfoliate))))) {
916 	      for(d = 0; d < 3; d++) {
917 		if(minE[d] > ccs->Coord[(3 * c) + d])
918 		  minE[d] = ccs->Coord[(3 * c) + d];
919 		if(maxE[d] < ccs->Coord[(3 * c) + d])
920 		  maxE[d] = ccs->Coord[(3 * c) + d];
921 	      }
922 	    }
923 	  }
924 	}
925       }
926       for(c = 0; c < 3; c++) {
927 	minE[c] -= (float) (I->max_vdw + 0.25F);
928 	maxE[c] += (float) (I->max_vdw + 0.25F);
929       }
930     }
931 
932     if (ok){
933       subtract3f(maxE, minE, sizeE);
934 
935       size = sizeE[0];
936       if(sizeE[1] > size)
937 	size = sizeE[1];
938       if(sizeE[2] > size)
939 	size = sizeE[2];
940 
941       gridSize = size / mesh_max; /* grid size is the largest axis divided by 25 */
942 
943       if(gridSize < min_spacing)
944 	gridSize = min_spacing;
945 
946       for(c = 0; c < 3; c++)
947 	dims[c] = (int) ((sizeE[c] / gridSize) + 1.5F);
948       field = new Isofield(G, dims);
949       CHECKOK(ok, field);
950     }
951 
952     if (ok){
953       for(a = 0; a < dims[0]; a++)
954 	for(b = 0; b < dims[1]; b++)
955 	  for(c = 0; c < dims[2]; c++)
956 	    F3(field->data, a, b, c) = 2.0;
957     }
958     OrthoBusyFast(G, 0, 1);
959     if(ok && !solv_acc)
960       ok &= RepMeshGetSolventDots(I, cs, minE, maxE, probe_radius);
961     if(ok) {
962       smap = MapNew(G, probe_radius, I->Dot, I->NDot, NULL);
963       map = MapNew(G, I->max_vdw + probe_radius, cs->Coord, cs->NIndex, NULL);
964     }
965     if(ok && map && smap) {
966       ok &= MapSetupExpress(smap);
967       if (ok)
968 	ok &= MapSetupExpress(map);
969       for(a = 0; ok && a < dims[0]; a++) {
970         OrthoBusyFast(G, dims[0] + a, dims[0] * 3);
971         point[0] = minE[0] + a * gridSize;
972         for(b = 0; ok && b < dims[1]; b++) {
973           point[1] = minE[1] + b * gridSize;
974           for(c = 0; ok && c < dims[2]; c++) {
975             float bestDist;
976             float dist2vdw;
977             float vdw_add = solv_acc * probe_radius;
978             point[2] = minE[2] + c * gridSize;
979             copy3f(point, F3Ptr(field->points, a, b, c));
980             aNear = -1;
981             bestDist = FLT_MAX;
982             aLen = FLT_MAX;
983             MapLocus(map, point, &h, &k, &l);
984             d = *(MapEStart(map, h, k, l));
985             if(d) {
986               cur = map->EList[d++];
987               while(ok && cur >= 0) {
988                 ai1 = obj->AtomInfo + cs->IdxToAtm[cur];
989                 if((inclH || (!ai1->isHydrogen())) &&
990                    ((!cullByFlag) || (!(ai1->flags & cAtomFlag_ignore)))) {
991                   vLen = (float) diff3f(point, cs->coordPtr(cur));
992                   dist2vdw = vLen - (ai1->vdw + vdw_add);
993                   if(dist2vdw < bestDist) {
994                     bestDist = dist2vdw;
995                     aLen = vLen;
996                     aNear = cur;
997                   }
998                 }
999                 cur = map->EList[d++];
1000 		ok &= !G->Interrupt;
1001               }
1002             }
1003             if(ok && aNear >= 0) {    /* near an atom... */
1004               pVal = aLen;      /* pVal is the distance from atom center */
1005               vdw = cs->Obj->AtomInfo[cs->IdxToAtm[aNear]].vdw + solv_acc * probe_radius;
1006               if((!solv_acc) && (pVal > vdw) && (pVal < (vdw + probe_radius))) {        /* point outside an atom */
1007                 inSolvFlag = false;
1008                 /* this point lies within a water radius of the atom, so
1009                    lets see if it is actually near a water */
1010                 aLen = FLT_MAX;
1011 
1012                 MapLocus(smap, point, &h, &k, &l);
1013                 d = *(MapEStart(smap, h, k, l));
1014                 if(d) {
1015                   cur = smap->EList[d++];
1016                   while(cur >= 0) {
1017                     vLen = diffsq3f(point, I->Dot + (cur * 3));
1018                     if(vLen < probe_radius2) {  /* within a water radius */
1019                       inSolvFlag = true;
1020                       break;
1021                     } else if(vLen < aLen) {    /* not within water radius */
1022                       aLen = vLen;      /* but nearest distance anyway */
1023                     }
1024                     cur = smap->EList[d++];
1025                   }
1026                 }
1027                 aLen = (float) sqrt1f(aLen);
1028                 if(inSolvFlag) {        /* point is inside a water, so a linear crossing point works fine */
1029                   pVal = pVal / vdw;
1030                 } else {        /* point is not inside a water, so we need to guestimate a value depending
1031                                  * on where the point is */
1032                   if(aLen < (2 * probe_radius)) {
1033                     pVal = 1.05F * (2 * probe_radius - aLen) / probe_radius;
1034                     /* within a radius, so assign based on water surface, but error on the side of exclusion */
1035 
1036                   } else {
1037                     pVal = 0.0; /* further than one water away... */
1038                   }
1039                 }
1040               } else {          /* either within atom or outside solvent shell */
1041                 pVal = pVal / vdw;
1042               }
1043               if(pVal < F3(field->data, a, b, c))
1044                 F3(field->data, a, b, c) = pVal;
1045             }
1046 	    ok &= !G->Interrupt;
1047           }
1048         }
1049 	ok &= !G->Interrupt;
1050 	if (!ok)
1051 	  break;
1052       }
1053     }
1054     MapFree(smap);
1055     MapFree(map);
1056     FreeP(I->Dot);
1057     OrthoBusyFast(G, 2, 3);
1058     if(ok) {
1059       ok &= IsosurfVolume(G, NULL, NULL, field, 1.0, I->N, I->V, NULL, mesh_type, mesh_skip,
1060                     1.0F);
1061     }
1062     DeleteP(field);
1063     if(ok && (I->N.data() && I->V.data() && (carve_flag || clear_flag || trim_flag))) {
1064       int cur_size = VLAGetSize(I->N);
1065       if((mesh_type == 0) && cur_size) {
1066         int *n = I->N.data();
1067         int *new_n = VLACalloc(int, cur_size);
1068         int new_size = 0;
1069         float *new_v = I->V.data();
1070         float *v = I->V.data();
1071 	CHECKOK(ok, new_n);
1072         while(ok && (c = *(n++))) {
1073           int new_c = 0;
1074           float *last_v = NULL;
1075           while(c--) {
1076             int a_keeper = false;
1077             if(trim_map) {
1078               int i = *(MapLocusEStart(trim_map, v));
1079               if(i) {
1080                 int j = trim_map->EList[i++];
1081                 while(j >= 0) {
1082                   float *v_targ = trim_vla + 3 * j;
1083                   if(within3f(v_targ, v, trim_cutoff)) {
1084                     a_keeper = true;
1085                     break;
1086                   }
1087                   j = trim_map->EList[i++];
1088                 }
1089               }
1090             } else {
1091               a_keeper = true;
1092             }
1093 
1094             if(a_keeper && carve_map) {
1095               int i = *(MapLocusEStart(carve_map, v));
1096               a_keeper = false;
1097               if(i) {
1098                 int j = carve_map->EList[i++];
1099                 while(j >= 0) {
1100                   float *v_targ = carve_vla + 3 * j;
1101                   if(within3f(v_targ, v, carve_cutoff)) {
1102                     a_keeper = true;
1103                     break;
1104                   }
1105                   j = carve_map->EList[i++];
1106                 }
1107               }
1108             }
1109 
1110             if(clear_map) {
1111               int i = *(MapLocusEStart(clear_map, v));
1112               if(i) {
1113                 int j = clear_map->EList[i++];
1114                 while(j >= 0) {
1115                   if(within3f(clear_vla + 3 * j, v, clear_cutoff)) {
1116                     a_keeper = false;
1117                     break;
1118                   }
1119                   j = clear_map->EList[i++];
1120                 }
1121               }
1122             }
1123             if(a_keeper) {
1124               if(new_c) {
1125                 new_c++;
1126                 *(new_v++) = *(v++);
1127                 *(new_v++) = *(v++);
1128                 *(new_v++) = *(v++);
1129               } else {
1130                 if(last_v) {
1131                   new_c += 2;
1132                   *(new_v++) = *(last_v++);
1133                   *(new_v++) = *(last_v++);
1134                   *(new_v++) = *(last_v++);
1135                   *(new_v++) = *(v++);
1136                   *(new_v++) = *(v++);
1137                   *(new_v++) = *(v++);
1138                   last_v = NULL;
1139                 } else {
1140                   last_v = v;
1141                   v += 3;
1142                 }
1143               }
1144             } else {
1145               last_v = NULL;
1146               v += 3;
1147               if(new_c) {
1148                 VLACheck(new_n, int, new_size + 1);     /* extends the zero count sentinel */
1149 		CHECKOK(ok, new_n);
1150 		if (ok){
1151 		  new_n[new_size] = new_c;
1152 		  new_size++;
1153 		  new_c = 0;
1154 		}
1155               }
1156             }
1157           }
1158           if(ok && new_c) {
1159             VLACheck(new_n, int, new_size + 1); /* extends the zero count sentinel */
1160 	    CHECKOK(ok, new_n);
1161             new_n[new_size] = new_c;
1162             new_size++;
1163             new_c = 0;
1164           }
1165         }
1166         I->N = pymol::vla_take_ownership(new_n);
1167       }
1168     }
1169     MapFree(trim_map);
1170     MapFree(carve_map);
1171     MapFree(clear_map);
1172     VLAFreeP(trim_vla);
1173     VLAFreeP(carve_vla);
1174     VLAFreeP(clear_vla);
1175     n = I->N.data();
1176     I->NTot = 0;
1177     if (ok){
1178       while(*n)
1179 	I->NTot += *(n++);
1180       RepMeshColor(I, cs);
1181     }
1182     OrthoBusyFast(G, 3, 4);
1183   }
1184   OrthoBusyFast(G, 4, 4);
1185   if(!ok) {
1186     RepMeshFree(I);
1187     I = NULL;
1188   }
1189   return (Rep *) I;
1190 }
1191 
RepMeshGetSolventDots(RepMesh * I,CoordSet * cs,float * min,float * max,float probe_radius)1192 int RepMeshGetSolventDots(RepMesh * I, CoordSet * cs, float *min, float *max,
1193                           float probe_radius)
1194 {
1195   PyMOLGlobals *G = cs->G;
1196   ObjectMolecule *obj = cs->Obj;
1197   int a, b, c, a1, a2, flag, i, h, k, l, j;
1198   int ok = true;
1199   float *v, *v0, vdw;
1200   MapType *map;
1201   int inFlag, *p, *dot_flag;
1202   SphereRec *sp = G->Sphere->Sphere[0];
1203   int cavity_cull;
1204   float probe_radius_plus;
1205   int dotCnt, maxCnt, maxDot = 0;
1206   int cnt;
1207   int inclH, mesh_mode, cullByFlag;
1208   AtomInfoType *ai1, *ai2;
1209   int ds = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_quality);
1210 
1211   if(ds < 0)
1212     ds = 0;
1213   if(ds > 4)
1214     ds = 4;
1215   sp = G->Sphere->Sphere[ds];
1216 
1217   cavity_cull = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_cavity_cull);
1218 
1219   mesh_mode = SettingGet_i(G, cs->Setting, obj->Setting, cSetting_mesh_mode);
1220   cullByFlag = (mesh_mode == cRepMesh_by_flags);
1221   inclH = !(mesh_mode == cRepMesh_heavy_atoms);
1222 
1223   I->Dot = pymol::malloc<float>(cs->NIndex * 3 * sp->nDot);
1224   ErrChkPtr(G, I->Dot);
1225 
1226   probe_radius_plus = probe_radius * 1.5F;
1227 
1228   I->NDot = 0;
1229   map = MapNew(G, I->max_vdw + probe_radius, cs->Coord, cs->NIndex, NULL);
1230   if(map) {
1231     MapSetupExpress(map);
1232     maxCnt = 0;
1233     v = I->Dot;
1234     for(a = 0; a < cs->NIndex; a++) {
1235 
1236       ai1 = obj->AtomInfo + cs->IdxToAtm[a];
1237       if((inclH || (!ai1->isHydrogen())) &&
1238          ((!cullByFlag) || (!(ai1->flags & (cAtomFlag_ignore))))) {
1239         OrthoBusyFast(G, a, cs->NIndex * 3);
1240         dotCnt = 0;
1241         a1 = cs->IdxToAtm[a];
1242         v0 = cs->coordPtr(a);
1243         vdw = cs->Obj->AtomInfo[a1].vdw + probe_radius;
1244         inFlag = true;
1245         for(c = 0; c < 3; c++) {
1246           if((min[c] - v0[c]) > vdw) {
1247             inFlag = false;
1248             break;
1249           };
1250           if((v0[c] - max[c]) > vdw) {
1251             inFlag = false;
1252             break;
1253           };
1254         }
1255         if(inFlag)
1256           for(b = 0; b < sp->nDot; b++) {
1257             v[0] = v0[0] + vdw * sp->dot[b][0];
1258             v[1] = v0[1] + vdw * sp->dot[b][1];
1259             v[2] = v0[2] + vdw * sp->dot[b][2];
1260             MapLocus(map, v, &h, &k, &l);
1261             flag = true;
1262             i = *(MapEStart(map, h, k, l));
1263             if(i) {
1264               j = map->EList[i++];
1265               while(j >= 0) {
1266 
1267                 ai2 = obj->AtomInfo + cs->IdxToAtm[j];
1268                 if((inclH || (!ai2->isHydrogen())) &&
1269                    ((!cullByFlag) || (!(ai2->flags & cAtomFlag_ignore))))
1270                   if(j != a) {
1271                     a2 = cs->IdxToAtm[j];
1272                     if(within3f
1273                        (cs->coordPtr(j), v, cs->Obj->AtomInfo[a2].vdw + probe_radius)) {
1274                       flag = false;
1275                       break;
1276                     }
1277                   }
1278                 j = map->EList[i++];
1279               }
1280             }
1281             if(flag) {
1282               dotCnt++;
1283               v += 3;
1284               I->NDot++;
1285             }
1286           }
1287         if(dotCnt > maxCnt) {
1288           maxCnt = dotCnt;
1289           maxDot = I->NDot - 1;
1290         }
1291       }
1292       ok &= !G->Interrupt;
1293       if(!ok) {
1294         break;
1295       }
1296     }
1297     MapFree(map);
1298   }
1299 
1300   if(ok && (cavity_cull > 0)) {
1301     dot_flag = pymol::malloc<int>(I->NDot);
1302     ErrChkPtr(G, dot_flag);
1303     for(a = 0; a < I->NDot; a++) {
1304       dot_flag[a] = 0;
1305     }
1306     dot_flag[maxDot] = 1;       /* this guarantees that we have a valid dot */
1307 
1308     map = MapNew(G, probe_radius_plus, I->Dot, I->NDot, NULL);
1309     if(map) {
1310       MapSetupExpress(map);
1311 
1312       flag = true;
1313       while(flag) {
1314         p = dot_flag;
1315         v = I->Dot;
1316 
1317         flag = false;
1318         for(a = 0; a < I->NDot; a++) {
1319           if(!dot_flag[a]) {
1320             cnt = 0;
1321             MapLocus(map, v, &h, &k, &l);
1322             i = *(MapEStart(map, h, k, l));
1323             if(i) {
1324               j = map->EList[i++];
1325               while(j >= 0) {
1326                 if(j != a) {
1327                   if(within3f(I->Dot + (3 * j), v, probe_radius_plus)) {
1328                     if(dot_flag[j]) {
1329                       *p = true;
1330                       flag = true;
1331                       break;
1332                     }
1333                     cnt++;
1334                     if(cnt > cavity_cull) {
1335                       *p = true;
1336                       flag = true;
1337                       break;
1338                     }
1339                   }
1340                 }
1341                 j = map->EList[i++];
1342               }
1343             }
1344           }
1345           v += 3;
1346           p++;
1347         }
1348 	ok &= !G->Interrupt;
1349 	if(!ok) {
1350 	  break;
1351 	}
1352       }
1353       MapFree(map);
1354     }
1355     v0 = I->Dot;
1356     v = I->Dot;
1357     p = dot_flag;
1358     c = I->NDot;
1359     I->NDot = 0;
1360     for(a = 0; a < c; a++) {
1361       if(*(p++)) {
1362         *(v0++) = *(v++);
1363         *(v0++) = *(v++);
1364         *(v0++) = *(v++);
1365         I->NDot++;
1366       } else {
1367         v += 3;
1368       }
1369     }
1370     FreeP(dot_flag);
1371   }
1372   if(!ok) {
1373     FreeP(I->Dot);
1374     I->NDot = 0;
1375   }
1376   return ok;
1377 }
1378