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
18 #include <stdint.h>
19 #include <algorithm>
20
21 #include"os_python.h"
22 #include"os_numpy.h"
23 #include"os_std.h"
24 #include"os_gl.h"
25
26 #include"OOMac.h"
27 #include"ObjectMap.h"
28 #include"Base.h"
29 #include"MemoryDebug.h"
30 #include"Map.h"
31 #include"Parse.h"
32 #include"Isosurf.h"
33 #include"Vector.h"
34 #include"Color.h"
35 #include"main.h"
36 #include"Scene.h"
37 #include"PConv.h"
38 #include"Word.h"
39 #include"Vector.h"
40 #include"PyMOLGlobals.h"
41 #include"Matrix.h"
42 #include"Util.h"
43 #include"ShaderMgr.h"
44 #include"CGO.h"
45 #include"File.h"
46 #include"Executive.h"
47 #include"Field.h"
48
49 #define n_space_group_numbers 231
50 static const char * space_group_numbers[] = {
51 "",
52 "P 1", "P -1", "P 2", "P 21", "C 2",
53 "P M", "P C", "C M", "C C", "P 2/M",
54 "P 21/M", "C 2/M", "P 2/C", "P 21/C", "C 2/C",
55 "P 2 2 2", "P 2 2 21", "P 21 21 2", "P 21 21 21", "C 2 2 21",
56 "C 2 2 2", "F 2 2 2", "I 2 2 2", "I 21 21 21", "P M M 2",
57 "P M C 21", "P C C 2", "P M A 2", "P C A 21", "P N C 2",
58 "P M N 21", "P B A 2", "P N A 21", "P N N 2", "C M M 2",
59 "C M C 21", "C C C 2", "A M M 2", "A B M 2", "A M A 2",
60 "A B A 2", "F M M 2", "F D D 2", "I M M 2", "I B A 2",
61 "I M A 2", "P M M M", "P N N N", "P C C M", "P B A N",
62 "P M M A", "P N N A", "P M N A", "P C C A", "P B A M",
63 "P C C N", "P B C M", "P N N M", "P M M N", "P B C N",
64 "P B C A", "P N M A", "C M C M", "C M C A", "C M M M",
65 "C C C M", "C M M A", "C C C A", "F M M M", "F D D D",
66 "I M M M", "I B A M", "I B C A", "I M M A", "P 4",
67 "P 41", "P 42", "P 43", "I 4", "I 41",
68 "P -4", "I -4", "P 4/M", "P 42/M", "P 4/N",
69 "P 42/N", "I 4/M", "I 41/A", "P 4 2 2", "P 4 21 2",
70 "P 41 2 2", "P 41 21 2", "P 42 2 2", "P 42 21 2", "P 43 2 2",
71 "P 43 21 2", "I 4 2 2", "I 41 2 2", "P 4 M M", "P 4 B M",
72 "P 42 C M", "P 42 N M", "P 4 C C", "P 4 N C", "P 42 M C",
73 "P 42 B C", "I 4 M M", "I 4 C M", "I 41 M D", "I 41 C D",
74 "P -4 2 M", "P -4 2 C", "P -4 21 M", "P -4 21 C", "P -4 M 2",
75 "P -4 C 2", "P -4 B 2", "P -4 N 2", "I -4 M 2", "I -4 C 2",
76 "I -4 2 M", "I -4 2 D", "P 4/M M M", "P 4/M C C", "P 4/N B M",
77 "P 4/N N C", "P 4/M B M", "P 4/M N C", "P 4/N M M", "P 4/N C C",
78 "P 42/M M C", "P 42/M C M", "P 42/N B C", "P 42/N N M", "P 42/M B C",
79 "P 42/M N M", "P 42/N M C", "P 42/N C M", "I 4/M M M", "I 4/M C M",
80 "I 41/A M D", "I 41/A C D", "P 3", "P 31", "P 32",
81 "R 3", "P -3", "R -3", "P 3 1 2", "P 3 2 1",
82 "P 31 1 2", "P 31 2 1", "P 32 1 2", "P 32 2 1", "R 3 2",
83 "P 3 M 1", "P 3 1 M", "P 3 C 1", "P 3 1 C", "R 3 M",
84 "R 3 C", "P -3 1 M", "P -3 1 C", "P -3 M 1", "P -3 C 1",
85 "R -3 M", "R -3 C", "P 6", "P 61", "P 65",
86 "P 62", "P 64", "P 63", "P -6", "P 6/M",
87 "P 63/M", "P 6 2 2", "P 61 2 2", "P 65 2 2", "P 62 2 2",
88 "P 64 2 2", "P 63 2 2", "P 6 M M", "P 6 C C", "P 63 C M",
89 "P 63 M C", "P -6 M 2", "P -6 C 2", "P -6 2 M", "P -6 2 C",
90 "P 6/M M M", "P 6/M C C", "P 63/M C M", "P 63/M M C", "P 2 3",
91 "F 2 3", "I 2 3", "P 21 3", "I 21 3", "P M -3",
92 "P N -3", "F M -3", "F D -3", "I M -3", "P A -3",
93 "I A -3", "P 4 3 2", "P 42 3 2", "F 4 3 2", "F 41 3 2",
94 "I 4 3 2", "P 43 3 2", "P 41 3 2", "I 41 3 2", "P -4 3 M",
95 "F -4 3 M", "I -4 3 M", "P -4 3 N", "F -4 3 C", "I -4 3 D",
96 "P M -3 M", "P N -3 N", "P M -3 N", "P N -3 M", "F M -3 M",
97 "F M -3 C", "F D -3 M", "F D -3 C", "I M -3 M", "I A -3 D",
98 };
99
100 /* MapState::ValidXtal -- determines whether the MapState's Xtal type passed in is valid
101 * PARAMS
102 * ms, MapState
103 * RETURNS
104 * true/false
105 */
ObjectMapStateValidXtal(ObjectMapState * ms)106 int ObjectMapStateValidXtal(ObjectMapState * ms)
107 {
108 if(ms && ms->Active) {
109 switch (ms->MapSource) {
110 /* these are the only formats which are known to contain xtal
111 information */
112 case cMapSourceCrystallographic:
113 case cMapSourceCCP4:
114 case cMapSourceBRIX:
115 case cMapSourceGRD:
116 return true;
117 break;
118
119 }
120 }
121 return false;
122 }
123
124 /* Map::ValidXtal -- detemines whether the Map is valid
125 * PARAMS
126 * I, Map
127 * state, map's state
128 * RETURNS
129 * true/false
130 */
ObjectMapValidXtal(ObjectMap * I,int state)131 int ObjectMapValidXtal(ObjectMap * I, int state)
132 {
133 if((state >= 0) && (state < I->State.size())) {
134 ObjectMapState *ms = &I->State[state];
135 return ObjectMapStateValidXtal(ms);
136 }
137 return false;
138 }
139
140 /* MapState::GetExcludedStats --
141 * PARARMS
142 * G, usual PyMOL globals
143 * ms, MapState
144 * vert_vla, variable length array of vertices
145 * beyond, radius of exclusion
146 * within, radius of inlcusion
147 * level, map level
148 *
149 * RETURNS
150 * number of exluded pts?
151 */
ObjectMapStateGetExcludedStats(PyMOLGlobals * G,ObjectMapState * ms,float * vert_vla,float beyond,float within,float * level)152 int ObjectMapStateGetExcludedStats(PyMOLGlobals * G, ObjectMapState * ms, float *vert_vla,
153 float beyond, float within, float *level)
154 {
155 double sum = 0.0, sumsq = 0.0;
156 float mean, stdev;
157 int cnt = 0;
158 int list_size;
159 float cutoff = beyond;
160 MapType *voxelmap = NULL;
161
162 /* size of the VLA */
163 if(vert_vla) {
164 list_size = VLAGetSize(vert_vla) / 3;
165 } else {
166 list_size = 0;
167 }
168 if(cutoff < within)
169 cutoff = within;
170
171 /* make a new map from the VLA .............. */
172 if(list_size)
173 voxelmap = MapNew(G, -cutoff, vert_vla, list_size, NULL);
174
175 if(voxelmap || (!list_size)) {
176 int a, b, c;
177 int h, k, l, i, j;
178 const int *fdim = ms->FDim;
179 int within_flag, within_default = false;
180 int beyond_flag;
181
182 const Isofield *field = ms->Field.get();
183 if(list_size)
184 MapSetupExpress(voxelmap);
185
186 within_flag = true;
187 beyond_flag = true;
188
189 if(within < R_SMALL4)
190 within_default = true;
191 for(c = 0; c < fdim[2]; c++) {
192 for(b = 0; b < fdim[1]; b++) {
193 for(a = 0; a < fdim[0]; a++) {
194 if(list_size) {
195 within_flag = within_default;
196 beyond_flag = true;
197
198 const float* v = F4Ptr(field->points, a, b, c, 0);
199
200 MapLocus(voxelmap, v, &h, &k, &l);
201 i = *(MapEStart(voxelmap, h, k, l));
202 if(i) {
203 j = voxelmap->EList[i++];
204 while(j >= 0) {
205 if(!within_flag) {
206 if(within3f(vert_vla + 3 * j, v, within)) {
207 within_flag = true;
208 }
209 }
210 if(within3f(vert_vla + 3 * j, v, beyond)) {
211 beyond_flag = false;
212 break;
213 }
214 j = voxelmap->EList[i++];
215 }
216 }
217 }
218
219 if(within_flag && beyond_flag) { /* point isn't too close to any vertex */
220 const float f_val = F3(field->data, a, b, c);
221 sum += f_val;
222 sumsq += (f_val * f_val);
223 cnt++;
224 }
225 }
226 }
227 }
228 if(voxelmap)
229 MapFree(voxelmap);
230 }
231 if(cnt) {
232 mean = (float) (sum / cnt);
233 stdev = (float) sqrt1d((sumsq - (sum * sum / cnt)) / (cnt));
234 level[1] = mean;
235 level[0] = mean - stdev;
236 level[2] = mean + stdev;
237 }
238 return cnt;
239 }
240
241 /* MapState::ObjectMapStateGetDataRange -- get min/max from the scalar field
242 * PARAMS
243 * G
244 * usual PyMOLGlobals
245 * ms
246 * I
247 * RETURNS
248 * npts, but stores range in min/max
249 */
ObjectMapStateGetDataRange(PyMOLGlobals * G,ObjectMapState * ms,float * min,float * max)250 int ObjectMapStateGetDataRange(PyMOLGlobals * G, ObjectMapState * ms, float *min,
251 float *max)
252 {
253 float max_val = 0.0F, min_val = 0.0F;
254 CField *data = ms->Field->data.get();
255 int cnt = data->dim[0] * data->dim[1] * data->dim[2];
256 float *raw_data = (float *) data->data.data();
257 if(cnt) {
258 int a;
259 min_val = (max_val = *(raw_data++));
260 for(a = 1; a < cnt; a++) {
261 double f_val = *(raw_data++);
262 if(min_val > f_val)
263 min_val = f_val;
264 if(max_val < f_val)
265 max_val = f_val;
266 }
267 }
268 *min = min_val;
269 *max = max_val;
270 return cnt;
271 }
272
273 /* MapState::ObjectMapStateGetHistogram -- compute a map histogram
274 *
275 * INPUT PARAMS
276 *
277 * h_points - number of histogram points
278 * limit - limit the data to (mean - limit * stdev, mean + limit * stdev)
279 * if limit <= 0, don't trim the histogram
280 * min_arg, max_arg - limit the data, has priority over "limit" param;
281 * Unused if min_arg==max_arg
282 *
283 * OUTPUT PARAMS
284 *
285 * histogram - output histogram buffer. first four values are
286 * minimum, maximum, mean and stdev, followed by n_points
287 * of non-normalized histogram counts.
288 */
ObjectMapStateGetHistogram(PyMOLGlobals * G,ObjectMapState * ms,int n_points,float limit,float * histogram,float min_arg,float max_arg)289 int ObjectMapStateGetHistogram(PyMOLGlobals * G, ObjectMapState * ms,
290 int n_points, float limit, float *histogram,
291 float min_arg, float max_arg)
292 {
293 float max_val = 0.0f, min_val = 0.0f;
294 float sum = 0.0f, sumsq = 0.0f;
295 float min_his, max_his, irange, mean, stdev;
296 int pos;
297 CField *data = ms->Field->data.get();
298 int cnt = data->dim[0] * data->dim[1] * data->dim[2];
299 float *raw_data = (float *) data->data.data();
300 if(cnt) {
301 int a;
302
303 // compute min/max/mean/stdev
304 sum = min_val = (max_val = *(raw_data++));
305 sumsq = sum*sum;
306 for(a = 1; a < cnt; a++) {
307 double f_val = *(raw_data++);
308 if(min_val > f_val)
309 min_val = f_val;
310 if(max_val < f_val)
311 max_val = f_val;
312 sum += f_val;
313 sumsq += f_val*f_val;
314 }
315 mean = (float) (sum / cnt);
316 stdev = (float) sqrt1d((sumsq - (sum * sum / cnt)) / (cnt));
317
318 // adjust min/max to limit
319 if (min_arg != max_arg) {
320 min_his = min_arg;
321 max_his = max_arg;
322 } else if (limit > 0.0F) {
323 min_his = mean - limit*stdev;
324 if (min_his < min_val)
325 min_his = min_val;
326 max_his = mean + limit*stdev;
327 if (max_his > max_val)
328 max_his = max_val;
329 } else {
330 min_his = min_val;
331 max_his = max_val;
332 }
333
334 // Compute the histogram
335 if(n_points > 0) {
336 irange = (float)(n_points-1) / (max_his - min_his);
337 for (a = 0; a < n_points; a++)
338 histogram[a+4] = 0.0f;
339 raw_data = (float *) data->data.data();
340 for (a = 0; a < cnt; a++) {
341 double f_val = *(raw_data++);
342 pos = (int)(irange * (f_val-min_his));
343 if (pos >= 0 && pos < n_points) {
344 histogram[pos+4] += 1.0;
345 }
346 }
347 }
348 histogram[0] = min_his;
349 histogram[1] = max_his;
350 histogram[2] = mean;
351 histogram[3] = stdev;
352 } else {
353 histogram[0] = 0.0;
354 histogram[1] = 1.0;
355 histogram[2] = 1.0;
356 histogram[3] = 1.0;
357 }
358 return cnt;
359 }
360
361 /**
362 * @see ObjectMapStateInterpolate for parameter description
363 * @param state Object state (can be -2 for current state)
364 * @return True if the map state exists and `result` has been written to
365 */
ObjectMapInterpolate(ObjectMap * I,int state,const float * array,float * result,int * flag,int n)366 int ObjectMapInterpolate(ObjectMap * I, int state, const float *array, float *result, int *flag,
367 int n)
368 {
369 int ok = false;
370 float txf_buffer[3];
371 float *txf = txf_buffer;
372
373 ObjectMapState *ms = ObjectMapGetState(I, state);
374
375 if(ms) {
376 double *matrix = ObjectStateGetInvMatrix(ms);
377
378 if(matrix) {
379 /* we have to back-transform points */
380 if(n > 1) {
381 txf = pymol::malloc<float>(3 * n);
382 }
383
384 const float *src = array;
385 float *dst = txf;
386 array = txf;
387
388 for (int nn = n; nn--; src += 3, dst += 3) {
389 transform44d3f(matrix, src, dst);
390 }
391 }
392
393 ok = ObjectMapStateInterpolate(ms, array, result, flag, n);
394 ok = true;
395 }
396
397 if(txf != txf_buffer)
398 FreeP(txf);
399
400 return (ok);
401 }
402
ObjectMapStateTrim(PyMOLGlobals * G,ObjectMapState * ms,float * mn,float * mx,int quiet)403 static void ObjectMapStateTrim(PyMOLGlobals * G, ObjectMapState * ms,
404 float *mn, float *mx, int quiet)
405 {
406 int div[3];
407 int min[3];
408 int fdim[4];
409 int new_min[3], new_max[3], new_fdim[3];
410 int a, b, c, d, e, f;
411 float *vt, *vr;
412 float v[3];
413 float grid[3];
414 Isofield *field;
415 float orig_size = 1.0F;
416 float new_size = 1.0F;
417
418 if(ObjectMapStateValidXtal(ms)) {
419 float tst[3], frac_tst[3];
420 float frac_mn[3];
421 float frac_mx[3];
422 int hit_flag = false;
423
424 /* compute the limiting box extents in fractional space */
425
426 for(a = 0; a < 8; a++) {
427 tst[0] = (a & 0x1) ? mn[0] : mx[0];
428 tst[1] = (a & 0x2) ? mn[1] : mx[1];
429 tst[2] = (a & 0x4) ? mn[2] : mx[2];
430 transform33f3f(ms->Symmetry->Crystal.RealToFrac, tst, frac_tst);
431 if(!a) {
432 copy3f(frac_tst, frac_mn);
433 copy3f(frac_tst, frac_mx);
434 } else {
435 for(b = 0; b < 3; b++) {
436 frac_mn[b] = (frac_mn[b] > frac_tst[b]) ? frac_tst[b] : frac_mn[b];
437 frac_mx[b] = (frac_mx[b] < frac_tst[b]) ? frac_tst[b] : frac_mx[b];
438 }
439 }
440 }
441 for(a = 0; a < 3; a++) {
442 div[a] = ms->Div[a];
443 min[a] = ms->Min[a];
444 fdim[a] = ms->FDim[a];
445 }
446 fdim[3] = 3;
447
448 {
449 int first_flag[3] = { false, false, false };
450 for(d = 0; d < 3; d++) {
451 int tst_min, tst_max;
452 float v_min, v_max;
453 for(c = 0; c < (fdim[d] - 1); c++) {
454 tst_min = c + min[d];
455 tst_max = tst_min + 1;
456 v_min = tst_min / ((float) div[d]);
457 v_max = tst_max / ((float) div[d]);
458 if(((v_min >= frac_mn[d]) && (v_min <= frac_mx[d])) ||
459 ((v_max >= frac_mn[d]) && (v_max <= frac_mx[d]))) {
460 if(!first_flag[d]) {
461 first_flag[d] = true;
462 new_min[d] = tst_min;
463 new_max[d] = tst_max;
464 } else {
465 new_min[d] = (new_min[d] > tst_min) ? tst_min : new_min[d];
466 new_max[d] = (new_max[d] < tst_max) ? tst_max : new_max[d];
467 }
468 }
469 }
470 new_fdim[d] = (new_max[d] - new_min[d]) + 1;
471 }
472 hit_flag = first_flag[0] && first_flag[1] & first_flag[2];
473 }
474
475 if(hit_flag)
476 hit_flag = (new_fdim[0] != fdim[0]) || (new_fdim[1] != fdim[1]) ||
477 (new_fdim[2] != fdim[2]);
478
479 if(hit_flag) {
480 orig_size = fdim[0] * fdim[1] * fdim[2];
481 new_size = new_fdim[0] * new_fdim[1] * new_fdim[2];
482
483 field = new Isofield(G, new_fdim);
484 field->save_points = ms->Field->save_points;
485
486 for(c = 0; c < new_fdim[2]; c++) {
487 f = c + (new_min[2] - min[2]);
488 for(b = 0; b < new_fdim[1]; b++) {
489 e = b + (new_min[1] - min[1]);
490 for(a = 0; a < new_fdim[0]; a++) {
491 d = a + (new_min[0] - min[0]);
492 vt = F4Ptr(field->points, a, b, c, 0);
493 vr = F4Ptr(ms->Field->points, d, e, f, 0);
494 copy3f(vr, vt);
495 F3(field->data, a, b, c) = F3(ms->Field->data, d, e, f);
496 }
497 }
498 }
499 for(a = 0; a < 3; a++) {
500 ms->Min[a] = new_min[a];
501 ms->Max[a] = new_max[a];
502 ms->FDim[a] = new_fdim[a];
503 }
504 ms->Field.reset(field);
505
506 /* compute new extents */
507 v[2] = (ms->Min[2]) / ((float) ms->Div[2]);
508 v[1] = (ms->Min[1]) / ((float) ms->Div[1]);
509 v[0] = (ms->Min[0]) / ((float) ms->Div[0]);
510
511 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMin);
512
513 v[2] = ((ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
514 v[1] = ((ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
515 v[0] = ((ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
516
517 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMax);
518
519 /* new corner */
520 {
521 float vv[3];
522 d = 0;
523 for(c = 0; c < ms->FDim[2]; c += (ms->FDim[2] - 1)) {
524 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
525 for(b = 0; b < ms->FDim[1]; b += (ms->FDim[1] - 1)) {
526 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
527 for(a = 0; a < ms->FDim[0]; a += (ms->FDim[0] - 1)) {
528 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
529 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vv);
530 copy3f(vv, ms->Corner + 3 * d);
531 d++;
532 }
533 }
534 }
535 }
536 }
537 } else { /* not a crystal map */
538 int hit_flag = false;
539 float *origin = ms->Origin.data();
540
541 for(a = 0; a < 3; a++) {
542 min[a] = ms->Min[a];
543 grid[a] = ms->Grid[a];
544 fdim[a] = ms->FDim[a];
545 }
546 fdim[3] = 3;
547
548 {
549 int first_flag[3] = { false, false, false };
550 for(d = 0; d < 3; d++) {
551 int tst_min, tst_max;
552 float v_min, v_max;
553 for(c = 0; c < (fdim[d] - 1); c++) {
554 tst_min = c + min[d];
555 tst_max = tst_min + 1;
556 v_min = origin[d] + grid[d] * (tst_min);
557 v_max = origin[d] + grid[d] * (tst_max);
558 if(((v_min >= mn[d]) && (v_min <= mx[d])) ||
559 ((v_max >= mn[d]) && (v_max <= mx[d]))) {
560 if(!first_flag[d]) {
561 first_flag[d] = true;
562 hit_flag = true;
563 new_min[d] = tst_min;
564 new_max[d] = tst_max;
565 } else {
566 new_min[d] = (new_min[d] > tst_min) ? tst_min : new_min[d];
567 new_max[d] = (new_max[d] < tst_max) ? tst_max : new_max[d];
568 }
569 }
570 }
571 new_fdim[d] = (new_max[d] - new_min[d]) + 1;
572 }
573 hit_flag = first_flag[0] && first_flag[1] & first_flag[2];
574 }
575
576 if(hit_flag)
577 hit_flag = (new_fdim[0] != fdim[0]) || (new_fdim[1] != fdim[1]) ||
578 (new_fdim[2] != fdim[2]);
579
580 if(hit_flag) {
581
582 orig_size = fdim[0] * fdim[1] * fdim[2];
583 new_size = new_fdim[0] * new_fdim[1] * new_fdim[2];
584
585 field = new Isofield(G, new_fdim);
586 field->save_points = ms->Field->save_points;
587
588 for(c = 0; c < new_fdim[2]; c++) {
589 f = c + (new_min[2] - min[2]);
590 for(b = 0; b < new_fdim[1]; b++) {
591 e = b + (new_min[1] - min[1]);
592 for(a = 0; a < new_fdim[0]; a++) {
593 d = a + (new_min[0] - min[0]);
594 vt = F4Ptr(field->points, a, b, c, 0);
595 vr = F4Ptr(ms->Field->points, d, e, f, 0);
596 copy3f(vr, vt);
597 F3(field->data, a, b, c) = F3(ms->Field->data, d, e, f);
598 }
599 }
600 }
601 for(a = 0; a < 3; a++) {
602 ms->Min[a] = new_min[a];
603 ms->Max[a] = new_max[a];
604 ms->FDim[a] = new_fdim[a];
605 if(!ms->Dim.empty())
606 ms->Dim[a] = new_fdim[a];
607 }
608
609 ms->Field.reset(field);
610
611 for(e = 0; e < 3; e++) {
612 ms->ExtentMin[e] = ms->Origin[e] + ms->Grid[e] * ms->Min[e];
613 ms->ExtentMax[e] = ms->Origin[e] + ms->Grid[e] * ms->Max[e];
614 }
615
616 d = 0;
617 for(c = 0; c < ms->FDim[2]; c += ms->FDim[2] - 1) {
618 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
619
620 for(b = 0; b < ms->FDim[1]; b += ms->FDim[1] - 1) {
621 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
622
623 for(a = 0; a < ms->FDim[0]; a += ms->FDim[0] - 1) {
624 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
625 copy3f(v, ms->Corner + 3 * d);
626 d++;
627 }
628 }
629 }
630 }
631 }
632 if(!quiet) {
633 PRINTFB(G, FB_ObjectMap, FB_Actions)
634 " ObjectMap: Map volume reduced by %2.0f%%.\n",
635 (100 * (orig_size - new_size)) / orig_size ENDFB(G);
636 }
637 }
638
ObjectMapStateDouble(PyMOLGlobals * G,ObjectMapState * ms)639 static void ObjectMapStateDouble(PyMOLGlobals * G, ObjectMapState * ms)
640 {
641 int div[3];
642 int min[3];
643 int max[3];
644 int fdim[4];
645 int a, b, c;
646 float v[3], vr[3];
647 float *vt;
648 float x, y, z;
649 float grid[3];
650
651 Isofield *field;
652
653 if(ObjectMapStateValidXtal(ms)) {
654 for(a = 0; a < 3; a++) {
655 div[a] = ms->Div[a] * 2;
656 min[a] = ms->Min[a] * 2;
657 max[a] = ms->Max[a] * 2;
658 fdim[a] = ms->FDim[a] * 2 - 1;
659 }
660 fdim[3] = 3;
661 field = new Isofield(G, fdim);
662 field->save_points = ms->Field->save_points;
663 for(c = 0; c < fdim[2]; c++) {
664 v[2] = (c + min[2]) / ((float) div[2]);
665 z = (c & 0x1) ? 0.5F : 0.0F;
666 for(b = 0; b < fdim[1]; b++) {
667 v[1] = (b + min[1]) / ((float) div[1]);
668 y = (b & 0x1) ? 0.5F : 0.0F;
669 for(a = 0; a < fdim[0]; a++) {
670 v[0] = (a + min[0]) / ((float) div[0]);
671 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
672 x = (a & 0x1) ? 0.5F : 0.0F;
673 vt = F4Ptr(field->points, a, b, c, 0);
674 copy3f(vr, vt);
675 if((a & 0x1) || (b & 0x1) || (c & 0x1)) {
676 F3(field->data, a, b, c) = FieldInterpolatef(ms->Field->data.get(),
677 a / 2, b / 2, c / 2, x, y, z);
678 } else {
679 F3(field->data, a, b, c) = F3(ms->Field->data, a / 2, b / 2, c / 2);
680 }
681 }
682 }
683 }
684 for(a = 0; a < 3; a++) {
685 ms->Min[a] = min[a];
686 ms->Max[a] = max[a];
687 ms->FDim[a] = fdim[a];
688 ms->Div[a] = div[a];
689 }
690
691 ms->Field.reset(field);
692 } else {
693 for(a = 0; a < 3; a++) {
694 grid[a] = ms->Grid[a] / 2.0F;
695 min[a] = ms->Min[a] * 2;
696 max[a] = ms->Max[a] * 2;
697 fdim[a] = ms->FDim[a] * 2 - 1;
698 }
699 fdim[3] = 3;
700
701 field = new Isofield(G, fdim);
702 field->save_points = ms->Field->save_points;
703
704 for(c = 0; c < fdim[2]; c++) {
705 v[2] = ms->Origin[2] + grid[2] * (c + min[2]);
706 z = (c & 0x1) ? 0.5F : 0.0F;
707 for(b = 0; b < fdim[1]; b++) {
708 v[1] = ms->Origin[1] + grid[1] * (b + min[1]);
709 y = (b & 0x1) ? 0.5F : 0.0F;
710 for(a = 0; a < fdim[0]; a++) {
711 v[0] = ms->Origin[0] + grid[0] * (a + min[0]);
712 x = (a & 0x1) ? 0.5F : 0.0F;
713 vt = F4Ptr(field->points, a, b, c, 0);
714 copy3f(v, vt);
715 if((a & 0x1) || (b & 0x1) || (c & 0x1)) {
716 F3(field->data, a, b, c) = FieldInterpolatef(ms->Field->data.get(),
717 a / 2, b / 2, c / 2, x, y, z);
718 } else {
719 F3(field->data, a, b, c) = F3(ms->Field->data, a / 2, b / 2, c / 2);
720 }
721 }
722 }
723 }
724 for(a = 0; a < 3; a++) {
725 ms->Min[a] = min[a];
726 ms->Max[a] = max[a];
727 ms->FDim[a] = fdim[a];
728 if(!ms->Dim.empty())
729 ms->Dim[a] = fdim[a];
730 if(!ms->Grid.empty())
731 ms->Grid[a] = grid[a];
732 }
733 ms->Field.reset(field);
734 }
735 }
736
ObjectMapStateHalve(PyMOLGlobals * G,ObjectMapState * ms,int smooth)737 static void ObjectMapStateHalve(PyMOLGlobals * G, ObjectMapState * ms, int smooth)
738 {
739 int div[3];
740 int min[3];
741 int max[3];
742 int fdim[4];
743 int a, b, c;
744 float v[3], vr[3];
745 float *vt;
746 float x, y, z;
747 float grid[3];
748
749 Isofield *field;
750
751 if(ObjectMapStateValidXtal(ms)) {
752 int *old_div, *old_min, *old_max;
753 int a_2, b_2, c_2;
754
755 for(a = 0; a < 3; a++) {
756 div[a] = ms->Div[a] / 2;
757 /* note that when Div is odd, a one-cell extrapolation will
758 occur in order to preserve map size and get us onto a even
759 number of sampling intervals for the cell */
760
761 min[a] = ms->Min[a] / 2;
762 max[a] = ms->Max[a] / 2;
763 while((min[a] * 2) < ms->Min[a])
764 min[a]++;
765 while((max[a] * 2) > ms->Max[a])
766 max[a]--;
767
768 fdim[a] = (max[a] - min[a]) + 1;
769 }
770 fdim[3] = 3;
771 old_div = ms->Div;
772 old_min = ms->Min;
773 old_max = ms->Max;
774
775 if(smooth)
776 FieldSmooth3f(ms->Field->data.get());
777
778 field = new Isofield(G, fdim);
779 field->save_points = ms->Field->save_points;
780
781 /*
782 printf("old_div %d %d %d\n",old_div[0],old_div[1],old_div[2]);
783 printf("old_min %d %d %d\n",old_min[0],old_min[1],old_min[2]);
784 printf("min %d %d %d\n",min[0],min[1],min[2]);
785 printf("%d %d %d\n",ms->FDim[0],ms->FDim[1],ms->FDim[2]);
786 */
787
788 for(c = 0; c < fdim[2]; c++) {
789 v[2] = (c + min[2]) / ((float) div[2]);
790 c_2 = (c + min[2]) * 2 - old_min[2];
791 z = (v[2] - ((c_2 + old_min[2]) / (float) old_div[2])) * old_div[2];
792 if(c_2 >= old_max[2]) {
793 c_2 = old_max[2] - 1;
794 z = (v[2] - ((c_2 + old_min[2]) / (float) old_div[2])) * old_div[2];
795 }
796 for(b = 0; b < fdim[1]; b++) {
797 v[1] = (b + min[1]) / ((float) div[1]);
798 b_2 = (b + min[1]) * 2 - old_min[1];
799 y = (v[1] - ((b_2 + old_min[1]) / (float) old_div[1])) * old_div[1];
800 if(b_2 >= old_max[1]) {
801 b_2 = old_max[1] - 1;
802 y = (v[1] - ((b_2 + old_min[1]) / (float) old_div[1])) * old_div[1];
803 }
804 for(a = 0; a < fdim[0]; a++) {
805 v[0] = (a + min[0]) / ((float) div[0]);
806 a_2 = (a + min[0]) * 2 - old_min[0];
807 x = (v[0] - ((a_2 + old_min[0]) / (float) old_div[0])) * old_div[0];
808 if(a_2 >= old_max[0]) {
809 a_2 = old_max[0] - 1;
810 x = (v[0] - ((a_2 + old_min[0]) / (float) old_div[0])) * old_div[0];
811 }
812 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
813 vt = F4Ptr(field->points, a, b, c, 0);
814 copy3f(vr, vt);
815 F3(field->data, a, b, c) = FieldInterpolatef(ms->Field->data.get(),
816 a_2, b_2, c_2, x, y, z);
817 }
818 }
819 }
820 for(a = 0; a < 3; a++) {
821 ms->Min[a] = min[a];
822 ms->Max[a] = max[a];
823 ms->FDim[a] = fdim[a];
824 ms->Div[a] = div[a];
825 }
826
827 ms->Field.reset(field);
828
829 /* compute new extents */
830 v[2] = (ms->Min[2]) / ((float) ms->Div[2]);
831 v[1] = (ms->Min[1]) / ((float) ms->Div[1]);
832 v[0] = (ms->Min[0]) / ((float) ms->Div[0]);
833
834 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMin);
835
836 v[2] = ((ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
837 v[1] = ((ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
838 v[0] = ((ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
839
840 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMax);
841
842 /* new corner */
843 {
844 float vv[3];
845 int d = 0;
846
847 for(c = 0; c < ms->FDim[2]; c += (ms->FDim[2] - 1)) {
848 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
849 for(b = 0; b < ms->FDim[1]; b += (ms->FDim[1] - 1)) {
850 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
851 for(a = 0; a < ms->FDim[0]; a += (ms->FDim[0] - 1)) {
852 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
853 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vv);
854 copy3f(vv, ms->Corner + 3 * d);
855 d++;
856 }
857 }
858 }
859 }
860 } else {
861 for(a = 0; a < 3; a++) {
862 grid[a] = ms->Grid[a] * 2.0F;
863 min[a] = ms->Min[a] / 2;
864 max[a] = ms->Max[a] / 2;
865 fdim[a] = (ms->FDim[a] + 1) / 2;
866 }
867 fdim[3] = 3;
868
869 field = new Isofield(G, fdim);
870 field->save_points = ms->Field->save_points;
871
872 for(c = 0; c < fdim[2]; c++) {
873 v[2] = ms->Origin[2] + grid[2] * (c + min[2]);
874 for(b = 0; b < fdim[1]; b++) {
875 v[1] = ms->Origin[1] + grid[1] * (b + min[1]);
876 for(a = 0; a < fdim[0]; a++) {
877 v[0] = ms->Origin[0] + grid[0] * (a + min[0]);
878 vt = F4Ptr(field->points, a, b, c, 0);
879 copy3f(v, vt);
880 F3(field->data, a, b, c) = F3(ms->Field->data, a * 2, b * 2, c * 2);
881 }
882 }
883 }
884 for(a = 0; a < 3; a++) {
885 ms->Min[a] = min[a];
886 ms->Max[a] = max[a];
887 ms->FDim[a] = fdim[a];
888 if(!ms->Dim.empty())
889 ms->Dim[a] = fdim[a];
890 if(!ms->Grid.empty())
891 ms->Grid[a] = grid[a];
892 }
893 ms->Field.reset(field);
894
895 }
896 }
897
ObjectMapTrim(ObjectMap * I,int state,float * mn,float * mx,int quiet)898 pymol::Result<> ObjectMapTrim(
899 ObjectMap* I, int state, float* mn, float* mx, int quiet)
900 {
901 int update = false;
902
903 /* TO DO: convert mn and mx into map local coordinates if map itself is transformed... */
904
905 if(state < 0) {
906 for(auto& ms : I->State) {
907 if(ms.Active) {
908 ObjectMapStateTrim(I->G, &ms, mn, mx, quiet);
909 update = true;
910 }
911 }
912 } else if((state >= 0) && (state < I->State.size()) && (I->State[state].Active)) {
913 ObjectMapStateTrim(I->G, &I->State[state], mn, mx, quiet);
914 } else {
915 return pymol::make_error("Invalid state.");
916 }
917 if(update)
918 ObjectMapUpdateExtents(I);
919 return {};
920 }
921
ObjectMapDouble(ObjectMap * I,int state)922 pymol::Result<> ObjectMapDouble(ObjectMap* I, int state)
923 {
924 if(state < 0) {
925 for(auto& state : I->State) {
926 if(state.Active)
927 ObjectMapStateDouble(I->G, &state);
928 }
929 } else if((state >= 0) && (state < I->State.size()) && (I->State[state].Active)) {
930 ObjectMapStateDouble(I->G, &I->State[state]);
931 } else {
932 return pymol::make_error("Invalidate state.");
933 }
934 return {};
935 }
936
ObjectMapHalve(ObjectMap * I,int state,int smooth)937 pymol::Result<> ObjectMapHalve(ObjectMap * I, int state, int smooth)
938 {
939 if(state < 0) {
940 for(auto& state : I->State) {
941 if(state.Active)
942 ObjectMapStateHalve(I->G, &state, smooth);
943 }
944 } else if((state >= 0) && (state < I->State.size()) && (I->State[state].Active)) {
945 ObjectMapStateHalve(I->G, &I->State[state], smooth);
946 } else {
947 return pymol::make_error("Invalidate state.");
948 }
949 ObjectMapUpdateExtents(I);
950
951 return {};
952 }
953
ObjectMapStateContainsPoint(ObjectMapState * ms,float * point)954 int ObjectMapStateContainsPoint(ObjectMapState * ms, float *point)
955 {
956 int result = false;
957 float x, y, z;
958 int x_floor, y_floor, z_floor;
959 int x_ceil, y_ceil, z_ceil;
960
961 if(ObjectMapStateValidXtal(ms)) {
962 float frac[3];
963
964 transform33f3f(ms->Symmetry->Crystal.RealToFrac, point, frac);
965
966 x = (ms->Div[0] * frac[0]);
967 y = (ms->Div[1] * frac[1]);
968 z = (ms->Div[2] * frac[2]);
969 x_floor = floor(x);
970 x_ceil = ceil(x);
971 y_floor = floor(y);
972 y_ceil = ceil(y);
973 z_floor = floor(z);
974 z_ceil = ceil(z);
975
976 if((x_floor >= ms->Min[0]) && (x_ceil <= ms->Max[0]) &&
977 (y_floor >= ms->Min[1]) && (y_ceil <= ms->Max[1]) &&
978 (z_floor >= ms->Min[2]) && (z_ceil <= ms->Max[2]))
979 result = true;
980 } else {
981 x = (point[0] - ms->Origin[0]) / ms->Grid[0];
982 y = (point[1] - ms->Origin[1]) / ms->Grid[1];
983 z = (point[2] - ms->Origin[2]) / ms->Grid[2];
984
985 x_floor = floor(x);
986 x_ceil = ceil(x);
987 y_floor = floor(y);
988 y_ceil = ceil(y);
989 z_floor = floor(z);
990 z_ceil = ceil(z);
991
992 if((x_floor >= ms->Min[0]) && (x_ceil <= ms->Max[0]) &&
993 (y_floor >= ms->Min[1]) && (y_ceil <= ms->Max[1]) &&
994 (z_floor >= ms->Min[2]) && (z_ceil <= ms->Max[2]))
995 result = true;
996
997 if((x >= ms->Min[0]) && (x <= ms->Max[0]) &&
998 (y >= ms->Min[1]) && (y <= ms->Max[1]) && (z >= ms->Min[2]) && (z <= ms->Max[2]))
999 result = true;
1000 }
1001 return (result);
1002 }
1003
1004 /**
1005 * @param array Coordinate array of length `3*n`
1006 * @param[out] result Array of length `n` which will be populated with map
1007 * values at coorinate positions (linear interpolated between grid points)
1008 * @param[out] flag Array of length `n` which will be populated with booleans
1009 * indicating if points were within map bounds (optional, can be NULL)
1010 * @return False if any coordinate was out of bounds
1011 */
ObjectMapStateInterpolate(ObjectMapState * ms,const float * array,float * result,int * flag,int n)1012 int ObjectMapStateInterpolate(ObjectMapState * ms, const float *array, float *result, int *flag,
1013 int n)
1014 {
1015 int ok = true;
1016 const float *inp;
1017 int a, b, c;
1018 float x, y, z;
1019 inp = array;
1020
1021 if(ObjectMapStateValidXtal(ms)) {
1022 float frac[3];
1023
1024 while(n--) {
1025 /* get the fractional coordinate */
1026 transform33f3f(ms->Symmetry->Crystal.RealToFrac, inp, frac);
1027
1028 inp += 3;
1029
1030 /* compute the effective lattice offset as a function of cell spacing */
1031
1032 x = (ms->Div[0] * frac[0]);
1033 y = (ms->Div[1] * frac[1]);
1034 z = (ms->Div[2] * frac[2]);
1035
1036 /* now separate the integral and fractional parts for interpolation */
1037
1038 a = (int) floor(x + R_SMALL8);
1039 b = (int) floor(y + R_SMALL8);
1040 c = (int) floor(z + R_SMALL8);
1041 x -= a;
1042 y -= b;
1043 z -= c;
1044
1045 if(flag)
1046 *flag = 1;
1047 /* wrap into the map */
1048
1049 if(a < ms->Min[0]) {
1050 if(x < 0.99F) {
1051 ok = false;
1052 if(flag)
1053 *flag = 0;
1054 }
1055 x = 0.0F;
1056 a = ms->Min[0];
1057 } else if(a >= ms->FDim[0] + ms->Min[0] - 1) {
1058 if(x > 0.01F) {
1059 ok = false;
1060 if(flag)
1061 *flag = 0;
1062 }
1063 x = 0.0F;
1064 a = ms->FDim[0] + ms->Min[0] - 1;
1065 }
1066
1067 if(b < ms->Min[1]) {
1068 if(y < 0.99F) {
1069 ok = false;
1070 if(flag)
1071 *flag = 0;
1072 }
1073 y = 0.0F;
1074 b = ms->Min[1];
1075 } else if(b >= ms->FDim[1] + ms->Min[1] - 1) {
1076 if(y > 0.01F) {
1077 ok = false;
1078 if(flag)
1079 *flag = 0;
1080 }
1081 y = 0.0F;
1082 b = ms->FDim[1] + ms->Min[1] - 1;
1083 }
1084
1085 if(c < ms->Min[2]) {
1086 if(z < 0.99F) {
1087 ok = false;
1088 if(flag)
1089 *flag = 0;
1090 }
1091 z = 0.0F;
1092 c = ms->Min[2];
1093 } else if(c >= ms->FDim[2] + ms->Min[2] - 1) {
1094 if(z > 0.01) {
1095 ok = false;
1096 if(flag)
1097 *flag = 0;
1098 }
1099 z = 0.0F;
1100 c = ms->FDim[2] + ms->Min[2] - 1;
1101 }
1102 /* printf("%d %d %d %8.3f %8.3f %8.3f\n",a,b,c,x,y,z); */
1103 *(result++) = FieldInterpolatef(ms->Field->data.get(),
1104 a - ms->Min[0],
1105 b - ms->Min[1], c - ms->Min[2], x, y, z);
1106 if(flag)
1107 flag++;
1108
1109 }
1110 } else {
1111
1112 while(n--) {
1113
1114 x = (inp[0] - ms->Origin[0]) / ms->Grid[0];
1115 y = (inp[1] - ms->Origin[1]) / ms->Grid[1];
1116 z = (inp[2] - ms->Origin[2]) / ms->Grid[2];
1117 inp += 3;
1118
1119 a = (int) floor(x + R_SMALL8);
1120 b = (int) floor(y + R_SMALL8);
1121 c = (int) floor(z + R_SMALL8);
1122 x -= a;
1123 y -= b;
1124 z -= c;
1125
1126 if(flag)
1127 *flag = 1;
1128 if(a < ms->Min[0]) {
1129 x = 0.0F;
1130 a = ms->Min[0];
1131 ok = false;
1132 if(flag)
1133 *flag = 0;
1134 } else if(a >= ms->Max[0]) {
1135 x = 1.0F;
1136 a = ms->Max[0] - 1;
1137 ok = false;
1138 if(flag)
1139 *flag = 0;
1140 }
1141
1142 if(b < ms->Min[1]) {
1143 y = 0.0F;
1144 b = ms->Min[1];
1145 ok = false;
1146 if(flag)
1147 *flag = 0;
1148 } else if(b >= ms->Max[1]) {
1149 y = 1.0F;
1150 b = ms->Max[1] - 1;
1151 ok = false;
1152 if(flag)
1153 *flag = 0;
1154 }
1155
1156 if(c < ms->Min[2]) {
1157 z = 0.0F;
1158 c = ms->Min[2];
1159 ok = false;
1160 if(flag)
1161 *flag = 0;
1162 } else if(c >= ms->Max[2]) {
1163 z = 1.0F;
1164 c = ms->Max[2] - 1;
1165 ok = false;
1166 if(flag)
1167 *flag = 0;
1168 }
1169 /* printf("%d %d %d %8.3f %8.3f %8.3f\n",a,b,c,x,y,z); */
1170 *(result++) = FieldInterpolatef(ms->Field->data.get(),
1171 a - ms->Min[0],
1172 b - ms->Min[1], c - ms->Min[2], x, y, z);
1173 if(flag)
1174 flag++;
1175 }
1176 }
1177 return (ok);
1178 }
1179
1180 #ifndef _PYMOL_NOPY
1181 static int ObjectMapNumPyArrayToMapState(PyMOLGlobals * G, ObjectMapState * I,
1182 PyObject * ary, int quiet);
1183 #endif
1184
ObjectMapRegeneratePoints(ObjectMap * om)1185 void ObjectMapRegeneratePoints(ObjectMap * om){
1186 int i;
1187 for (i=0; i<om->State.size();i++){
1188 ObjectMapStateRegeneratePoints(&om->State[i]);
1189 }
1190 }
1191
ObjectMapStateRegeneratePoints(ObjectMapState * ms)1192 void ObjectMapStateRegeneratePoints(ObjectMapState * ms)
1193 {
1194 int a, b, c, e;
1195 float v[3], vr[3];
1196
1197 if(ObjectMapStateValidXtal(ms)) {
1198 for(c = 0; c < ms->FDim[2]; c++) {
1199 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
1200 for(b = 0; b < ms->FDim[1]; b++) {
1201 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
1202 for(a = 0; a < ms->FDim[0]; a++) {
1203 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
1204 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
1205 for(e = 0; e < 3; e++)
1206 F4(ms->Field->points, a, b, c, e) = vr[e];
1207 }
1208 }
1209 }
1210 } else {
1211 for(c = 0; c < ms->FDim[2]; c++) {
1212 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
1213 for(b = 0; b < ms->FDim[1]; b++) {
1214 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
1215 for(a = 0; a < ms->FDim[0]; a++) {
1216 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
1217 for(e = 0; e < 3; e++) {
1218 F4(ms->Field->points, a, b, c, e) = v[e];
1219 }
1220 }
1221 }
1222 }
1223 }
1224 }
1225
ObjectMapStateAsPyList(ObjectMapState * I)1226 static PyObject *ObjectMapStateAsPyList(ObjectMapState * I)
1227 {
1228 PyObject *result = NULL;
1229
1230 result = PyList_New(16);
1231 PyList_SetItem(result, 0, PyInt_FromLong(I->Active));
1232 if(I->Symmetry) {
1233 PyList_SetItem(result, 1, SymmetryAsPyList(I->Symmetry.get()));
1234 } else {
1235 PyList_SetItem(result, 1, PConvAutoNone(Py_None));
1236 }
1237 if(!I->Origin.empty()) {
1238 PyList_SetItem(result, 2, PConvFloatArrayToPyList(I->Origin.data(), 3));
1239 } else {
1240 PyList_SetItem(result, 2, PConvAutoNone(Py_None));
1241 }
1242 if(!I->Range.empty()) {
1243 PyList_SetItem(result, 3, PConvFloatArrayToPyList(I->Range.data(), 3));
1244 } else {
1245 PyList_SetItem(result, 3, PConvAutoNone(Py_None));
1246 }
1247 if(!I->Dim.empty()) {
1248 PyList_SetItem(result, 4, PConvIntArrayToPyList(I->Dim.data(), 3));
1249 } else {
1250 PyList_SetItem(result, 4, PConvAutoNone(Py_None));
1251 }
1252 if(!I->Grid.empty()) {
1253 PyList_SetItem(result, 5, PConvFloatArrayToPyList(I->Grid.data(), 3));
1254 } else {
1255 PyList_SetItem(result, 5, PConvAutoNone(Py_None));
1256 }
1257 PyList_SetItem(result, 6, PConvFloatArrayToPyList(I->Corner, 24));
1258 PyList_SetItem(result, 7, PConvFloatArrayToPyList(I->ExtentMin, 3));
1259 PyList_SetItem(result, 8, PConvFloatArrayToPyList(I->ExtentMax, 3));
1260 PyList_SetItem(result, 9, PyInt_FromLong(I->MapSource));
1261
1262 PyList_SetItem(result, 10, PConvIntArrayToPyList(I->Div, 3));
1263 PyList_SetItem(result, 11, PConvIntArrayToPyList(I->Min, 3));
1264 PyList_SetItem(result, 12, PConvIntArrayToPyList(I->Max, 3));
1265 PyList_SetItem(result, 13, PConvIntArrayToPyList(I->FDim, 4));
1266
1267 PyList_SetItem(result, 14, IsosurfAsPyList(I->G, I->Field.get()));
1268 PyList_SetItem(result, 15, ObjectStateAsPyList(I));
1269 return (PConvAutoNone(result));
1270 }
1271
ObjectMapAllStatesAsPyList(ObjectMap * I)1272 static PyObject *ObjectMapAllStatesAsPyList(ObjectMap * I)
1273 {
1274 PyObject *result = NULL;
1275 int a;
1276 result = PyList_New(I->State.size());
1277 for(a = 0; a < I->State.size(); a++) {
1278 if(I->State[a].Active) {
1279 PyList_SetItem(result, a, ObjectMapStateAsPyList(&I->State[a]));
1280 } else {
1281 PyList_SetItem(result, a, PConvAutoNone(NULL));
1282 }
1283 }
1284 return (PConvAutoNone(result));
1285
1286 }
1287
ObjectMapStateCopy(const ObjectMapState * src,ObjectMapState * I)1288 static int ObjectMapStateCopy(const ObjectMapState * src, ObjectMapState * I)
1289 {
1290 int ok = true;
1291 if(ok) {
1292 I->Active = src->Active;
1293 if(I->Active) {
1294 I->Symmetry = src->Symmetry;
1295 I->Origin = src->Origin;
1296 I->Range = src->Range;
1297 I->Grid = src->Grid;
1298 I->Dim = src->Dim;
1299 std::copy_n(src->Corner, 24, I->Corner);
1300
1301 copy3f(src->ExtentMin, I->ExtentMin);
1302 copy3f(src->ExtentMax, I->ExtentMax);
1303
1304 I->MapSource = src->MapSource;
1305
1306 copy3f(src->Div, I->Div);
1307 copy3f(src->Min, I->Min);
1308 copy3f(src->Max, I->Max);
1309 copy3f(src->FDim, I->FDim);
1310
1311 I->Field = src->Field;
1312 if(ok)
1313 ObjectMapStateRegeneratePoints(I);
1314 }
1315 }
1316 return (ok);
1317 }
1318
ObjectMapState(const ObjectMapState & src)1319 ObjectMapState::ObjectMapState(const ObjectMapState& src) : CObjectState(src)
1320 {
1321 ObjectMapStateCopy(&src, this);
1322 }
1323
operator =(const ObjectMapState & src)1324 ObjectMapState& ObjectMapState::operator=(const ObjectMapState& src)
1325 {
1326 CObjectState::operator=(src);
1327 ObjectMapStateCopy(&src, this);
1328 return *this;
1329 }
1330
ObjectMapStateFromPyList(PyMOLGlobals * G,ObjectMapState * I,PyObject * list)1331 static int ObjectMapStateFromPyList(PyMOLGlobals * G, ObjectMapState * I, PyObject * list)
1332 {
1333 int ok = true;
1334 int ll = 0;
1335 PyObject *tmp;
1336 if(ok)
1337 ok = (list != NULL);
1338 if(ok) {
1339 if(!PyList_Check(list))
1340 I->Active = false;
1341 else {
1342 if(ok)
1343 ok = PyList_Check(list);
1344 if(ok)
1345 ll = PyList_Size(list);
1346 /* TO SUPPORT BACKWARDS COMPATIBILITY...
1347 Always check ll when adding new PyList_GetItem's */
1348
1349 if(ok)
1350 ok = PConvPyIntToInt(PyList_GetItem(list, 0), &I->Active);
1351 if(ok) {
1352 tmp = PyList_GetItem(list, 1);
1353 if(tmp == Py_None)
1354 I->Symmetry = NULL;
1355 else {
1356 I->Symmetry.reset(SymmetryNewFromPyList(G, tmp));
1357 ok = I->Symmetry != nullptr;
1358 }
1359 CPythonVal_Free(tmp);
1360 }
1361 if(ok) {
1362 tmp = PyList_GetItem(list, 2);
1363 if(CPythonVal_IsNone(tmp))
1364 I->Origin.clear();
1365 else
1366 ok = PConvFromPyObject(G, tmp, I->Origin);
1367 CPythonVal_Free(tmp);
1368 }
1369 if(ok) {
1370 tmp = PyList_GetItem(list, 3);
1371 if(CPythonVal_IsNone(tmp))
1372 I->Range.clear();
1373 else
1374 ok = PConvFromPyObject(G, tmp, I->Range);
1375 CPythonVal_Free(tmp);
1376 }
1377 if(ok) {
1378 tmp = PyList_GetItem(list, 4);
1379 if(CPythonVal_IsNone(tmp))
1380 I->Dim.clear();
1381 else
1382 ok = PConvFromPyObject(G, tmp, I->Dim);
1383 CPythonVal_Free(tmp);
1384 }
1385 if(ok) {
1386 tmp = PyList_GetItem(list, 5);
1387 if(CPythonVal_IsNone(tmp))
1388 I->Grid.clear();
1389 else
1390 ok = PConvFromPyObject(G, tmp, I->Grid);
1391 CPythonVal_Free(tmp);
1392 }
1393 if(ok)
1394 ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list, 6), I->Corner, 24);
1395 if(ok)
1396 ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list, 7), I->ExtentMin, 3);
1397 if(ok)
1398 ok = PConvPyListToFloatArrayInPlace(PyList_GetItem(list, 8), I->ExtentMax, 3);
1399 if(ok)
1400 ok = PConvPyIntToInt(PyList_GetItem(list, 9), &I->MapSource);
1401 if(ok)
1402 ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list, 10), I->Div, 3);
1403 if(ok)
1404 ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list, 11), I->Min, 3);
1405 if(ok)
1406 ok = PConvPyListToIntArrayInPlace(PyList_GetItem(list, 12), I->Max, 3);
1407 if(ok)
1408 ok = CPythonVal_PConvPyListToIntArrayInPlace_From_List(G, list, 13, I->FDim, 4);
1409 if(ok){
1410 tmp = CPythonVal_PyList_GetItem(G, list, 14);
1411 I->Field.reset(IsosurfNewFromPyList(G, tmp));
1412 ok = I->Field != nullptr;
1413 CPythonVal_Free(tmp);
1414 }
1415 if(ok && (ll > 15)){
1416 tmp = CPythonVal_PyList_GetItem(G, list, 15);
1417 ok = ObjectStateFromPyList(G, tmp, I);
1418 CPythonVal_Free(tmp);
1419 }
1420 if(ok)
1421 ObjectMapStateRegeneratePoints(I);
1422 }
1423 }
1424 return (ok);
1425 }
1426
ObjectMapAllStatesFromPyList(ObjectMap * I,PyObject * list)1427 static int ObjectMapAllStatesFromPyList(ObjectMap * I, PyObject * list)
1428 {
1429 int ok = true;
1430 int a;
1431 if(ok)
1432 ok = PyList_Check(list);
1433 if(ok) {
1434 I->State.resize(PyList_Size(list), ObjectMapState(I->G));
1435 for(a = 0; a < I->State.size(); a++) {
1436 CPythonVal *val = CPythonVal_PyList_GetItem(I->G, list, a);
1437 ok = ObjectMapStateFromPyList(I->G, &I->State[a], val);
1438 CPythonVal_Free(val);
1439 if(!ok)
1440 break;
1441 }
1442 }
1443 return (ok);
1444 }
1445
ObjectMapAsPyList(ObjectMap * I)1446 PyObject *ObjectMapAsPyList(ObjectMap * I)
1447 {
1448 PyObject *result = NULL;
1449
1450 result = PyList_New(3);
1451 PyList_SetItem(result, 0, ObjectAsPyList(I));
1452 PyList_SetItem(result, 1, PyInt_FromLong(I->State.size()));
1453 PyList_SetItem(result, 2, ObjectMapAllStatesAsPyList(I));
1454
1455 return (PConvAutoNone(result));
1456 }
1457
ObjectMapNewFromPyList(PyMOLGlobals * G,PyObject * list,ObjectMap ** result)1458 int ObjectMapNewFromPyList(PyMOLGlobals * G, PyObject * list, ObjectMap ** result)
1459 {
1460 int ok = true;
1461 ObjectMap *I = NULL;
1462 (*result) = NULL;
1463
1464 if(ok)
1465 ok = (list != NULL);
1466 if(ok)
1467 ok = PyList_Check(list);
1468 /* TO SUPPORT BACKWARDS COMPATIBILITY...
1469 Always check ll when adding new PyList_GetItem's */
1470 I = new ObjectMap(G);
1471 if(ok)
1472 ok = (I != NULL);
1473
1474 if(ok){
1475 auto *val = PyList_GetItem(list, 0);
1476 ok = ObjectFromPyList(G, val, I);
1477 }
1478 if(ok){
1479 CPythonVal *val = CPythonVal_PyList_GetItem(G, list, 2);
1480 ok = ObjectMapAllStatesFromPyList(I, val);
1481 CPythonVal_Free(val);
1482 }
1483 if(ok) {
1484 (*result) = I;
1485 ObjectMapUpdateExtents(I);
1486 } else {
1487 /* cleanup? */
1488 }
1489
1490 return (ok);
1491 }
1492
ObjectMapNewCopy(PyMOLGlobals * G,const ObjectMap * src,ObjectMap ** result,int source_state,int target_state)1493 int ObjectMapNewCopy(PyMOLGlobals * G, const ObjectMap * src, ObjectMap ** result,
1494 int source_state, int target_state)
1495 {
1496 int ok = true;
1497 ObjectMap *I = NULL;
1498 I = new ObjectMap(G);
1499 if(ok)
1500 ok = (I != NULL);
1501 if(ok)
1502 ok = ObjectCopyHeader(I, src);
1503 if(ok) {
1504 if(source_state == -1) { /* all states */
1505 int state;
1506 VecCheckEmplace(I->State, I->State.size(), I->G);
1507 for(state = 0; state < src->State.size(); state++) {
1508 I->State[state] = src->State[state];
1509 }
1510 } else {
1511 if(target_state < 0)
1512 target_state = 0;
1513 if(source_state < 0)
1514 source_state = 0;
1515 VecCheckEmplace(I->State, target_state, G);
1516 if(source_state < src->State.size()) {
1517 I->State[target_state] = src->State[source_state];
1518 } else {
1519 ok = false;
1520 /* to do */
1521 }
1522 }
1523 }
1524 if(ok)
1525 *result = I;
1526 return ok;
1527 }
1528
ObjectMapStatePrime(ObjectMap * I,int state)1529 ObjectMapState *ObjectMapStatePrime(ObjectMap * I, int state)
1530 {
1531 if(state < 0)
1532 state = I->State.size();
1533 if(I->State.size() <= state) {
1534 VecCheckEmplace(I->State, state, I->G);
1535 }
1536 return &I->State[state];
1537 }
1538
ObjectMapUpdateExtents(ObjectMap * I)1539 void ObjectMapUpdateExtents(ObjectMap * I)
1540 {
1541 int a;
1542 float *min_ext, *max_ext;
1543 float tr_min[3], tr_max[3];
1544 I->ExtentFlag = false;
1545
1546 for(a = 0; a < I->State.size(); a++) {
1547 ObjectMapState *ms = &I->State[a];
1548 if(ms->Active) {
1549 if(!ms->Matrix.empty()) {
1550 transform44d3f(ms->Matrix.data(), ms->ExtentMin, tr_min);
1551 transform44d3f(ms->Matrix.data(), ms->ExtentMax, tr_max);
1552 {
1553 float tmp;
1554 int a;
1555 for(a = 0; a < 3; a++)
1556 if(tr_min[a] > tr_max[a]) {
1557 tmp = tr_min[a];
1558 tr_min[a] = tr_max[a];
1559 tr_max[a] = tmp;
1560 }
1561 }
1562 min_ext = tr_min;
1563 max_ext = tr_max;
1564 } else {
1565 min_ext = ms->ExtentMin;
1566 max_ext = ms->ExtentMax;
1567 }
1568
1569 if(!I->ExtentFlag) {
1570 copy3f(min_ext, I->ExtentMin);
1571 copy3f(max_ext, I->ExtentMax);
1572 I->ExtentFlag = true;
1573 } else {
1574 min3f(min_ext, I->ExtentMin, I->ExtentMin);
1575 max3f(max_ext, I->ExtentMax, I->ExtentMax);
1576 }
1577 }
1578 }
1579
1580 if(I->TTTFlag && I->ExtentFlag) {
1581 const float *ttt;
1582 double tttd[16];
1583 if(ObjectGetTTT(I, &ttt, -1)) {
1584 convertTTTfR44d(ttt, tttd);
1585 MatrixTransformExtentsR44d3f(tttd,
1586 I->ExtentMin, I->ExtentMax,
1587 I->ExtentMin, I->ExtentMax);
1588 }
1589 }
1590
1591 PRINTFD(I->G, FB_ObjectMap)
1592 " ObjectMapUpdateExtents-DEBUG: ExtentFlag %d\n", I->ExtentFlag ENDFD;
1593 }
1594
ObjectMapStateClamp(ObjectMapState * I,float clamp_floor,float clamp_ceiling)1595 void ObjectMapStateClamp(ObjectMapState * I, float clamp_floor, float clamp_ceiling)
1596 {
1597 int a, b, c;
1598 float *fp;
1599
1600 for(a = 0; a < I->FDim[0]; a++)
1601 for(b = 0; b < I->FDim[1]; b++)
1602 for(c = 0; c < I->FDim[2]; c++) {
1603 fp = F3Ptr(I->Field->data, a, b, c);
1604 if(*fp < clamp_floor)
1605 *fp = clamp_floor;
1606 else if(*fp > clamp_ceiling)
1607 *fp = clamp_ceiling;
1608 }
1609 }
1610
ObjectMapStateSetBorder(ObjectMapState * I,float level)1611 int ObjectMapStateSetBorder(ObjectMapState * I, float level)
1612 {
1613 int result = true;
1614 int a, b, c;
1615
1616 c = I->FDim[2] - 1;
1617 for(a = 0; a < I->FDim[0]; a++)
1618 for(b = 0; b < I->FDim[1]; b++) {
1619 F3(I->Field->data, a, b, 0) = level;
1620 F3(I->Field->data, a, b, c) = level;
1621 }
1622
1623 a = I->FDim[0] - 1;
1624 for(b = 0; b < I->FDim[1]; b++)
1625 for(c = 0; c < I->FDim[2]; c++) {
1626 F3(I->Field->data, 0, b, c) = level;
1627 F3(I->Field->data, a, b, c) = level;
1628 }
1629
1630 b = I->FDim[1] - 1;
1631 for(a = 0; a < I->FDim[0]; a++)
1632 for(c = 0; c < I->FDim[2]; c++) {
1633 F3(I->Field->data, a, 0, c) = level;
1634 F3(I->Field->data, a, b, c) = level;
1635 }
1636 return (result);
1637 }
1638
ObjectMapStatePurge(PyMOLGlobals * G,ObjectMapState * I)1639 void ObjectMapStatePurge(PyMOLGlobals * G, ObjectMapState * I)
1640 {
1641 ObjectStatePurge(I);
1642 I->Field = nullptr;
1643 I->Origin.clear();
1644 I->Dim.clear();
1645 I->Range.clear();
1646 I->Grid.clear();
1647 I->shaderCGO = nullptr;
1648
1649 I->Symmetry = nullptr;
1650
1651 I->Active = false;
1652 }
1653
update()1654 void ObjectMap::update()
1655 {
1656 auto I = this;
1657 if(!I->ExtentFlag) {
1658 ObjectMapUpdateExtents(I);
1659 if(I->ExtentFlag)
1660 SceneInvalidate(I->G);
1661 }
1662 }
1663
invalidate(int rep,int level,int state)1664 void ObjectMap::invalidate(int rep, int level, int state)
1665 {
1666 auto I = this;
1667 if(level >= cRepInvExtents) {
1668 I->ExtentFlag = false;
1669 }
1670 if((rep < 0) || (rep == cRepDot)) {
1671 int a;
1672 for(a = 0; a < I->State.size(); a++) {
1673 if(I->State[a].Active)
1674 I->State[a].have_range = false;
1675 I->State[a].shaderCGO = nullptr;
1676 }
1677 }
1678 SceneInvalidate(I->G);
1679 }
1680
1681 /* Has no prototype */
ObjectMapCGOGenerate(PyMOLGlobals * G,float * corner)1682 static CGO* ObjectMapCGOGenerate(PyMOLGlobals *G, float* corner)
1683 {
1684 int ok = true;
1685 CGO *convertCGO = NULL;
1686 CGO *shaderCGO = CGONewSized(G, 0);
1687 CGOBegin(shaderCGO, GL_LINES);
1688
1689 CGOVertexv(shaderCGO, corner + 3 * 0);
1690 CGOVertexv(shaderCGO, corner + 3 * 1);
1691
1692 CGOVertexv(shaderCGO, corner + 3 * 0);
1693 CGOVertexv(shaderCGO, corner + 3 * 2);
1694
1695 CGOVertexv(shaderCGO, corner + 3 * 2);
1696 CGOVertexv(shaderCGO, corner + 3 * 3);
1697
1698 CGOVertexv(shaderCGO, corner + 3 * 1);
1699 CGOVertexv(shaderCGO, corner + 3 * 3);
1700
1701 CGOVertexv(shaderCGO, corner + 3 * 0);
1702 CGOVertexv(shaderCGO, corner + 3 * 4);
1703
1704 CGOVertexv(shaderCGO, corner + 3 * 1);
1705 CGOVertexv(shaderCGO, corner + 3 * 5);
1706
1707 CGOVertexv(shaderCGO, corner + 3 * 2);
1708 CGOVertexv(shaderCGO, corner + 3 * 6);
1709
1710 CGOVertexv(shaderCGO, corner + 3 * 3);
1711 CGOVertexv(shaderCGO, corner + 3 * 7);
1712
1713 CGOVertexv(shaderCGO, corner + 3 * 4);
1714 CGOVertexv(shaderCGO, corner + 3 * 5);
1715
1716 CGOVertexv(shaderCGO, corner + 3 * 4);
1717 CGOVertexv(shaderCGO, corner + 3 * 6);
1718
1719 CGOVertexv(shaderCGO, corner + 3 * 6);
1720 CGOVertexv(shaderCGO, corner + 3 * 7);
1721
1722 CGOVertexv(shaderCGO, corner + 3 * 5);
1723 CGOVertexv(shaderCGO, corner + 3 * 7);
1724
1725 CGOEnd(shaderCGO);
1726
1727 CGOStop(shaderCGO);
1728
1729 convertCGO = CGOCombineBeginEnd(shaderCGO, 0);
1730 CHECKOK(ok, convertCGO);
1731 CGOFree(shaderCGO);
1732
1733 shaderCGO = convertCGO;
1734 if (ok)
1735 convertCGO = CGOOptimizeToVBONotIndexedWithReturnedData(shaderCGO, 0, 0, NULL);
1736 else
1737 return NULL;
1738 CHECKOK(ok, convertCGO);
1739 if (!ok)
1740 return NULL;
1741 CGOFree(shaderCGO);
1742 shaderCGO = convertCGO;
1743
1744 shaderCGO->use_shader = true;
1745
1746 return shaderCGO;
1747 }
1748
render(RenderInfo * info)1749 void ObjectMap::render(RenderInfo * info)
1750 {
1751 auto I = this;
1752 int state = info->state;
1753 CRay *ray = info->ray;
1754 auto pick = info->pick;
1755 int pass = info->pass;
1756 ObjectMapState *ms = NULL;
1757
1758 if(pass)
1759 return;
1760
1761 for(StateIterator iter(G, I->Setting, state, I->State.size());
1762 iter.next();) {
1763 state = iter.state;
1764 if(I->State[state].Active)
1765 ms = &I->State[state];
1766
1767 if(ms) {
1768 float *corner = ms->Corner;
1769 float tr_corner[24];
1770 ObjectPrepareContext(I, info);
1771
1772 if(!ms->Matrix.empty()) { /* transform the corners before drawing */
1773 int a;
1774 for(a = 0; a < 8; a++) {
1775 transform44d3f(ms->Matrix.data(), corner + 3 * a, tr_corner + 3 * a);
1776 }
1777 corner = tr_corner;
1778 }
1779
1780 if((I->visRep & cRepExtentBit)) {
1781 if(ray) {
1782 const float *vc;
1783 float radius = ray->PixelRadius / 1.4142F;
1784 vc = ColorGet(G, I->Color);
1785 ray->color3fv(vc);
1786 ray->sausage3fv(corner + 3 * 0, corner + 3 * 1, radius, vc, vc);
1787 ray->sausage3fv(corner + 3 * 0, corner + 3 * 2, radius, vc, vc);
1788 ray->sausage3fv(corner + 3 * 2, corner + 3 * 3, radius, vc, vc);
1789 ray->sausage3fv(corner + 3 * 1, corner + 3 * 3, radius, vc, vc);
1790 ray->sausage3fv(corner + 3 * 0, corner + 3 * 4, radius, vc, vc);
1791 ray->sausage3fv(corner + 3 * 1, corner + 3 * 5, radius, vc, vc);
1792 ray->sausage3fv(corner + 3 * 2, corner + 3 * 6, radius, vc, vc);
1793 ray->sausage3fv(corner + 3 * 3, corner + 3 * 7, radius, vc, vc);
1794 ray->sausage3fv(corner + 3 * 4, corner + 3 * 5, radius, vc, vc);
1795 ray->sausage3fv(corner + 3 * 4, corner + 3 * 6, radius, vc, vc);
1796 ray->sausage3fv(corner + 3 * 6, corner + 3 * 7, radius, vc, vc);
1797 ray->sausage3fv(corner + 3 * 5, corner + 3 * 7, radius, vc, vc);
1798 } else if(G->HaveGUI && G->ValidContext) {
1799 if(pick) {
1800 #ifndef PURE_OPENGL_ES_2
1801 } else if (!info->use_shaders) {
1802 // immediate
1803 ObjectUseColor(I);
1804 glDisable(GL_LIGHTING);
1805 glBegin(GL_LINES);
1806 glVertex3fv(corner + 3 * 0);
1807 glVertex3fv(corner + 3 * 1);
1808
1809 glVertex3fv(corner + 3 * 0);
1810 glVertex3fv(corner + 3 * 2);
1811
1812 glVertex3fv(corner + 3 * 2);
1813 glVertex3fv(corner + 3 * 3);
1814
1815 glVertex3fv(corner + 3 * 1);
1816 glVertex3fv(corner + 3 * 3);
1817
1818 glVertex3fv(corner + 3 * 0);
1819 glVertex3fv(corner + 3 * 4);
1820
1821 glVertex3fv(corner + 3 * 1);
1822 glVertex3fv(corner + 3 * 5);
1823
1824 glVertex3fv(corner + 3 * 2);
1825 glVertex3fv(corner + 3 * 6);
1826
1827 glVertex3fv(corner + 3 * 3);
1828 glVertex3fv(corner + 3 * 7);
1829
1830 glVertex3fv(corner + 3 * 4);
1831 glVertex3fv(corner + 3 * 5);
1832
1833 glVertex3fv(corner + 3 * 4);
1834 glVertex3fv(corner + 3 * 6);
1835
1836 glVertex3fv(corner + 3 * 6);
1837 glVertex3fv(corner + 3 * 7);
1838
1839 glVertex3fv(corner + 3 * 5);
1840 glVertex3fv(corner + 3 * 7);
1841
1842 glEnd();
1843 glEnable(GL_LIGHTING);
1844 } else {
1845 #endif
1846 // shader
1847 if (!ms->shaderCGO) {
1848 ms->shaderCGO.reset(ObjectMapCGOGenerate(G, corner));
1849 }
1850
1851 if (ms->shaderCGO) {
1852 CShaderPrg* shaderPrg = G->ShaderMgr->Enable_DefaultShader(info->pass);
1853 if (shaderPrg) {
1854 shaderPrg->SetLightingEnabled(0);
1855
1856 CGORenderGL(ms->shaderCGO.get(), ColorGet(G, I->Color),
1857 NULL, NULL, info, NULL);
1858 shaderPrg->Disable();
1859 }
1860 }
1861 }
1862 }
1863 }
1864
1865 if((I->visRep & cRepDotBit)) {
1866 if(!ms->have_range) {
1867 double sum = 0.0, sumsq = 0.0;
1868 CField *data = ms->Field->data.get();
1869 int cnt = data->dim[0] * data->dim[1] * data->dim[2];
1870 float *raw_data = (float *) data->data.data();
1871 int a;
1872 for(a = 0; a < cnt; a++) {
1873 double f_val = *(raw_data++);
1874 sum += f_val;
1875 sumsq += (f_val * f_val);
1876 }
1877 if(cnt) {
1878 float mean, stdev;
1879 mean = (float) (sum / cnt);
1880 stdev = (float) sqrt1d((sumsq - (sum * sum / cnt)) / (cnt));
1881 ms->high_cutoff = mean + stdev;
1882 ms->low_cutoff = mean - stdev;
1883 ms->have_range = true;
1884 }
1885 }
1886 if(ms->have_range && SettingGet_b(G, NULL, I->Setting, cSetting_dot_normals)) {
1887 IsofieldComputeGradients(G, ms->Field.get());
1888 }
1889 if(ms->have_range) {
1890 int a;
1891 CField *data = ms->Field->data.get();
1892 int cnt = data->dim[0] * data->dim[1] * data->dim[2];
1893 CField *points = ms->Field->points.get();
1894 CField *gradients = NULL;
1895
1896 if(SettingGet_b(G, NULL, I->Setting, cSetting_dot_normals)) {
1897 gradients = ms->Field->gradients.get();
1898 }
1899 if(data && points) {
1900 float *raw_data = (float *) data->data.data();
1901 float raw_point[3], *raw_point_ptr = (float *) points->data.data();
1902
1903 #define RAW_POINT_TRANSFORM(ptr, v3f) { \
1904 if(!ms->Matrix.empty()) \
1905 transform44d3f(ms->Matrix.data(), ptr, v3f); \
1906 else \
1907 copy3f(ptr, v3f); \
1908 ptr += 3; \
1909 }
1910
1911 float *raw_gradient = NULL;
1912 float high_cut = ms->high_cutoff, low_cut = ms->low_cutoff;
1913 float width =
1914 SettingGet_f(G, NULL, I->Setting, cSetting_dot_width);
1915
1916 if(ray) {
1917 float radius = ray->PixelRadius * width / 1.4142F;
1918 float vc[3];
1919 int color = I->Color;
1920 int ramped = ColorCheckRamped(G, I->Color);
1921
1922 {
1923 const float *tmp = ColorGet(G, I->Color);
1924 copy3f(tmp, vc);
1925 }
1926
1927 for(a = 0; a < cnt; a++) {
1928 float f_val = *(raw_data++);
1929 RAW_POINT_TRANSFORM(raw_point_ptr, raw_point);
1930 if((f_val >= high_cut) || (f_val <= low_cut)) {
1931 if(ramped) {
1932 ColorGetRamped(G, color, raw_point, vc, state);
1933 ray->color3fv(vc);
1934 }
1935 ray->sphere3fv(raw_point, radius);
1936 }
1937 }
1938 } else if(G->HaveGUI && G->ValidContext) {
1939 if(pick) {
1940 } else if (ALWAYS_IMMEDIATE_OR(!info->use_shaders)) {
1941 if(gradients) {
1942 raw_gradient = (float *) gradients->data.data();
1943 } else {
1944 glDisable(GL_LIGHTING);
1945 }
1946 {
1947 int ramped = ColorCheckRamped(G, I->Color);
1948 float vc[3];
1949 int color = I->Color;
1950 float gt[3];
1951
1952 glPointSize(width);
1953 glDisable(GL_POINT_SMOOTH);
1954 glBegin(GL_POINTS);
1955 ObjectUseColor(I);
1956 for(a = 0; a < cnt; a++) {
1957 float f_val = *(raw_data++);
1958 RAW_POINT_TRANSFORM(raw_point_ptr, raw_point);
1959 if(f_val >= high_cut) {
1960 if(raw_gradient) {
1961 normalize23f(raw_gradient, gt);
1962 invert3f(gt);
1963 glNormal3fv(gt);
1964 }
1965 if(ramped) {
1966 ColorGetRamped(G, color, raw_point, vc, state);
1967 glColor3fv(vc);
1968 }
1969 glVertex3fv(raw_point);
1970 } else if(f_val <= low_cut) {
1971 if(raw_gradient) {
1972 normalize23f(raw_gradient, gt);
1973 glNormal3fv(gt);
1974 }
1975 if(ramped) {
1976 ColorGetRamped(G, color, raw_point, vc, state);
1977 glColor3fv(vc);
1978 }
1979 glVertex3fv(raw_point);
1980 }
1981 if(raw_gradient)
1982 raw_gradient += 3;
1983 }
1984 glEnd();
1985 glEnable(GL_POINT_SMOOTH);
1986 }
1987 }
1988 }
1989 }
1990 }
1991 }
1992 }
1993 }
1994 }
1995
ObjectMapState(PyMOLGlobals * G)1996 ObjectMapState::ObjectMapState(PyMOLGlobals* G)
1997 : CObjectState(G)
1998 {
1999 auto I = this;
2000 ObjectMapStatePurge(G, I);
2001 ObjectStateInit(G, I);
2002 I->Symmetry.reset(new CSymmetry(G));
2003 I->Field = NULL;
2004 I->Origin.clear();
2005 I->Dim.clear();
2006 I->Range.clear();
2007 I->Grid.clear();
2008 I->MapSource = cMapSourceUndefined;
2009 I->have_range = false;
2010 }
2011
getNFrame() const2012 int ObjectMap::getNFrame() const
2013 {
2014 return State.size();
2015 }
2016
2017 /*========================================================================*/
_getObjectState(int state)2018 CObjectState* ObjectMap::_getObjectState(int state)
2019 {
2020 if (!State[state].Active)
2021 return nullptr;
2022 return &State[state];
2023 }
2024
2025 /*========================================================================*/
getSymmetry(int state) const2026 CSymmetry const* ObjectMap::getSymmetry(int state) const
2027 {
2028 auto* ms = ObjectMapGetState(const_cast<ObjectMap*>(this), state);
2029
2030 if (ms) {
2031 return ms->Symmetry.get();
2032 }
2033
2034 return nullptr;
2035 }
2036
2037 /*========================================================================*/
setSymmetry(CSymmetry const & symmetry,int state)2038 bool ObjectMap::setSymmetry(CSymmetry const& symmetry, int state)
2039 {
2040 bool success = false;
2041
2042 for (StateIterator iter(G, Setting, state, State.size()); iter.next();) {
2043 auto& oms = State[iter.state];
2044 if (oms.Active) {
2045 oms.Symmetry.reset(new CSymmetry(symmetry));
2046 success = true;
2047 }
2048 }
2049
2050 if (success) {
2051 ObjectMapRegeneratePoints(this);
2052 }
2053
2054 return success;
2055 }
2056
2057 /*========================================================================*/
ObjectMap(PyMOLGlobals * G)2058 ObjectMap::ObjectMap(PyMOLGlobals * G) : CObject(G)
2059 {
2060 auto I = this;
2061 I->type = cObjectMap;
2062
2063 I->visRep = cRepExtentBit;
2064 }
2065
2066
2067 /*========================================================================*/
ObjectMapNewStateFromDesc(PyMOLGlobals * G,ObjectMap * I,ObjectMapDesc * inp_md,int state,int quiet)2068 ObjectMapState *ObjectMapNewStateFromDesc(PyMOLGlobals * G, ObjectMap * I,
2069 ObjectMapDesc * inp_md, int state, int quiet)
2070 {
2071 int ok = true;
2072 float v[3];
2073 int a, b, c, d;
2074 float *fp;
2075 ObjectMapState *ms = NULL;
2076 ObjectMapDesc _md, *md;
2077 ms = ObjectMapStatePrime(I, state);
2078
2079 md = &_md;
2080 *(md) = *(inp_md);
2081
2082 if(I) {
2083 ms->Origin = std::vector<float>(3);
2084 ms->Range = std::vector<float>(3);
2085 ms->Grid = std::vector<float>(3);
2086 ms->MapSource = cMapSourceDesc;
2087 }
2088 switch (md->mode) {
2089 case cObjectMap_OrthoMinMaxGrid: /* Orthorhombic: min, max, spacing, centered over range */
2090
2091 subtract3f(md->MaxCorner, md->MinCorner, v);
2092 for(a = 0; a < 3; a++) {
2093 if(v[a] < 0.0)
2094 std::swap(md->MaxCorner[a], md->MinCorner[a]);
2095 };
2096 subtract3f(md->MaxCorner, md->MinCorner, v);
2097 for(a = 0; a < 3; a++) {
2098 md->Dim[a] = (int) (v[a] / md->Grid[a]);
2099 if(md->Dim[a] < 1)
2100 md->Dim[a] = 1;
2101 if((md->Dim[a] * md->Grid[a]) < v[a])
2102 md->Dim[a]++;
2103 }
2104
2105 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2106 " ObjectMap: Dim %d %d %d\n", md->Dim[0], md->Dim[1], md->Dim[2]
2107 ENDFB(I->G);
2108
2109 average3f(md->MaxCorner, md->MinCorner, v);
2110 for(a = 0; a < 3; a++) {
2111 md->MinCorner[a] = v[a] - 0.5F * md->Dim[a] * md->Grid[a];
2112 }
2113
2114 if(Feedback(I->G, FB_ObjectMap, FB_Blather)) {
2115 dump3f(md->MinCorner, " ObjectMap: MinCorner:");
2116 dump3f(md->MaxCorner, " ObjectMap: MaxCorner:");
2117 dump3f(md->Grid, " ObjectMap: Grid:");
2118 }
2119
2120 /* now populate the map data structure */
2121
2122 copy3f(md->MinCorner, ms->Origin.data());
2123 copy3f(md->Grid, ms->Grid.data());
2124 for(a = 0; a < 3; a++)
2125 ms->Range[a] = md->Grid[a] * (md->Dim[a] - 1);
2126
2127 /* these maps start at zero */
2128 for(a = 0; a < 3; a++) {
2129 ms->Min[a] = 0;
2130 ms->Max[a] = md->Dim[a] - 1;
2131 ms->Div[a] = md->Dim[a] - 1;
2132 }
2133
2134 /* define corners */
2135
2136 for(a = 0; a < 8; a++)
2137 copy3f(ms->Origin.data(), ms->Corner + 3 * a);
2138
2139 d = 0;
2140 for(c = 0; c < 2; c++) {
2141 {
2142 v[2] = (c ? ms->Range[2] : 0.0F);
2143 for(b = 0; b < 2; b++) {
2144 v[1] = (b ? ms->Range[1] : 0.0F);
2145 for(a = 0; a < 2; a++) {
2146 v[0] = (a ? ms->Range[0] : 0.0F);
2147 add3f(v, ms->Corner + 3 * d, ms->Corner + 3 * d);
2148 d++;
2149 }
2150 }
2151 }
2152 }
2153 for(a = 0; a < 3; a++)
2154 ms->FDim[a] = ms->Max[a] - ms->Min[a] + 1;
2155 ms->FDim[3] = 3;
2156
2157 ms->Field.reset(new Isofield(I->G, ms->FDim));
2158 if(!ms->Field)
2159 ok = false;
2160 else {
2161 for(a = 0; a < md->Dim[0]; a++) {
2162 v[0] = md->MinCorner[0] + a * md->Grid[0];
2163 for(b = 0; b < md->Dim[1]; b++) {
2164 v[1] = md->MinCorner[1] + b * md->Grid[1];
2165 for(c = 0; c < md->Dim[2]; c++) {
2166 v[2] = md->MinCorner[2] + c * md->Grid[2];
2167 fp = F4Ptr(ms->Field->points, a, b, c, 0);
2168 copy3f(v, fp);
2169 }
2170 }
2171 }
2172 }
2173 break;
2174 default:
2175 ok = false;
2176 }
2177 if(ok) {
2178 switch (md->init_mode) {
2179 case 0:
2180 for(a = 0; a < md->Dim[0]; a++) {
2181 for(b = 0; b < md->Dim[1]; b++) {
2182 for(c = 0; c < md->Dim[2]; c++) {
2183 F3(ms->Field->data, a, b, c) = 0.0F;
2184 }
2185 }
2186 }
2187 break;
2188 case 1:
2189 for(a = 0; a < md->Dim[0]; a++) {
2190 for(b = 0; b < md->Dim[1]; b++) {
2191 for(c = 0; c < md->Dim[2]; c++) {
2192 F3(ms->Field->data, a, b, c) = 1.0F;
2193 }
2194 }
2195 }
2196 break;
2197 case -2: /* for testing... */
2198 for(a = 0; a < md->Dim[0]; a++) {
2199 for(b = 0; b < md->Dim[1]; b++) {
2200 for(c = 0; c < md->Dim[2]; c++) {
2201 F3(ms->Field->data, a, b, c) = (float) sqrt1d(a * a + b * b + c * c);
2202 }
2203 }
2204 }
2205 break;
2206 }
2207 }
2208
2209 if(ok) {
2210 copy3f(ms->Origin.data(), ms->ExtentMin);
2211 copy3f(ms->Origin.data(), ms->ExtentMax);
2212 add3f(ms->Range.data(), ms->ExtentMax, ms->ExtentMax);
2213 ObjectMapUpdateExtents(I);
2214 }
2215 if(!ok) {
2216 ErrMessage(I->G, "ObjectMap", "Unable to create map");
2217 DeleteP(I);
2218 } else {
2219 if(!quiet) {
2220 PRINTFB(I->G, FB_ObjectMap, FB_Actions)
2221 " ObjectMap: Map created.\n" ENDFB(I->G);
2222 }
2223 }
2224
2225 return (ms);
2226 }
2227
2228
2229 /*========================================================================*/
ccp4_next_value(char ** pp,int mode)2230 static float ccp4_next_value(char ** pp, int mode) {
2231 char * p = *pp;
2232 switch(mode) {
2233 case 0:
2234 *pp += 1;
2235 return (float) *((int8_t *) p);
2236 case 1:
2237 *pp += 2;
2238 return (float) *((int16_t *) p);
2239 case 2:
2240 *pp += 4;
2241 return *((float *) p);
2242 }
2243 printf("ERROR unsupported mode\n");
2244 return 0.f;
2245 }
2246
2247 /* swaps n*width bytes in memory starting at p */
swap_endian(char * p,int n,int width)2248 static void swap_endian(char * p, int n, int width) {
2249 char tmp, *q, *pstop = p + (n - 1) * width + 1;
2250 int w2 = width/2, wm1 = width - 1;
2251 for(; p < pstop; p += w2) {
2252 for(q = p + wm1; p < q; p++, q--) {
2253 tmp = *p; *p = *q; *q = tmp;
2254 }
2255 }
2256 }
2257
ObjectMapCCP4StrToMap(ObjectMap * I,char * CCP4Str,int bytes,int state,int quiet,int format)2258 static int ObjectMapCCP4StrToMap(ObjectMap * I, char *CCP4Str, int bytes, int state,
2259 int quiet, int format)
2260 {
2261 auto G = I->G;
2262 char *p;
2263 int *i;
2264 size_t bytes_per_pt;
2265 char *q;
2266 float dens;
2267 int a, b, c, d, e;
2268 float v[3], vr[3], maxd, mind;
2269 int ok = true;
2270 int little_endian = 1, map_endian;
2271 /* CCP4 named from their docs */
2272 int nc, nr, ns;
2273 int map_mode;
2274 int ncstart, nrstart, nsstart;
2275 int nx, ny, nz;
2276 float xlen, ylen, zlen, alpha, beta, gamma;
2277 int ispg; // space group number
2278 int sym_skip;
2279 int mapc, mapr, maps;
2280 int cc[3];
2281 int n_pts;
2282 double sum, sumsq;
2283 float mean, stdev;
2284 int normalize;
2285 ObjectMapState *ms;
2286 int expectation;
2287
2288 switch (format) {
2289 case cLoadTypeCCP4Map:
2290 format = cLoadTypeCCP4Str;
2291 case cLoadTypeCCP4Str:
2292 case cLoadTypeCCP4Unspecified:
2293 case cLoadTypeMRC:
2294 break;
2295 default:
2296 ErrMessage(G, __func__, "wrong format");
2297 return false;
2298 }
2299
2300 if(bytes < 256 * sizeof(int)) {
2301 PRINTFB(I->G, FB_ObjectMap, FB_Errors)
2302 " ObjectMapCCP4: Map appears to be truncated -- aborting." ENDFB(I->G);
2303 return (0);
2304 }
2305
2306 /* state check */
2307 if(state < 0)
2308 state = I->State.size();
2309 /* alloc/init a new MapState */
2310 if(I->State.size() <= state) {
2311 VecCheckEmplace(I->State, state, I->G);
2312 }
2313 ms = &I->State[state];
2314
2315 normalize = SettingGetGlobal_b(I->G, cSetting_normalize_ccp4_maps);
2316
2317 p = CCP4Str;
2318 little_endian = *((char *) &little_endian);
2319 map_endian = (*p || *(p + 1)); // NOTE: this assumes 0x0 < NC < 0x10000
2320
2321 if(little_endian != map_endian) {
2322 if(!quiet) {
2323 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2324 " ObjectMapCCP4: Map appears to be reverse endian, swapping...\n" ENDFB(I->G);
2325 }
2326 swap_endian(p, 256, sizeof(int));
2327 }
2328
2329 i = (int *) p;
2330 nc = *(i++); /* columns */
2331 nr = *(i++); /* rows */
2332 ns = *(i++); /* sections */
2333 if(!quiet) {
2334 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2335 " ObjectMapCCP4: NC %d NR %d NS %d\n", nc, nr, ns ENDFB(I->G);
2336 }
2337 map_mode = *(i++); /* mode */
2338
2339 switch(map_mode) {
2340 case 0:
2341 bytes_per_pt = 1; // int8_t
2342 break;
2343 case 1:
2344 bytes_per_pt = 2; // int16_t
2345 break;
2346 case 2:
2347 bytes_per_pt = 4; // float
2348 break;
2349 default:
2350 PRINTFB(I->G, FB_ObjectMap, FB_Errors)
2351 "ObjectMapCCP4-ERR: Only map mode 0-2 currently supported (this map is mode %d)",
2352 map_mode ENDFB(I->G);
2353 return (0);
2354 }
2355
2356 if(!quiet) {
2357 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2358 " ObjectMapCCP4: Map is mode %d.\n", map_mode ENDFB(I->G);
2359 }
2360 ncstart = *(i++);
2361 nrstart = *(i++);
2362 nsstart = *(i++);
2363
2364 if(!quiet) {
2365 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2366 " ObjectMapCCP4: NCSTART %d NRSTART %d NSSTART %d\n",
2367 ncstart, nrstart, nsstart ENDFB(I->G);
2368 }
2369
2370 nx = *(i++);
2371 ny = *(i++);
2372 nz = *(i++);
2373
2374 if(!quiet) {
2375 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2376 " ObjectMapCCP4: NX %d NY %d NZ %d \n", nx, ny, nz ENDFB(I->G);
2377 }
2378
2379 xlen = *(float *) (i++);
2380 ylen = *(float *) (i++);
2381 zlen = *(float *) (i++);
2382
2383 if(!quiet) {
2384 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2385 " ObjectMapCCP4: X %8.3f Y %8.3f Z %8.3f \n", xlen, ylen, zlen ENDFB(I->G);
2386 }
2387
2388 alpha = *(float *) (i++);
2389 beta = *(float *) (i++);
2390 gamma = *(float *) (i++);
2391
2392 if(!quiet) {
2393 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2394 " ObjectMapCCP4: alpha %8.3f beta %8.3f gamma %8.3f \n",
2395 alpha, beta, gamma ENDFB(I->G);
2396 }
2397
2398 mapc = *(i++);
2399 mapr = *(i++);
2400 maps = *(i++);
2401
2402 if(!quiet) {
2403 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2404 " ObjectMapCCP4: MAPC %d MAPR %d MAPS %d \n", mapc, mapr, maps ENDFB(I->G);
2405 }
2406
2407 // AMIN, AMAX, AMEAN
2408 mind = *(float *) (i++);
2409 maxd = *(float *) (i++);
2410 mean = *(float *) (i++);
2411
2412 ispg = *(i++);
2413 sym_skip = *(i++);
2414
2415 // CCP4 Format:
2416 // 25-37 skew transformation
2417 // 38-52 future use "Storage space sometimes used by specific programs (default = 0)"
2418
2419 // MRC 2000 Format:
2420 // 25-49 "Extra, user defined storage space (default = 0)"
2421 // 50-52 ORIGIN
2422
2423 {
2424 int skew = *(i++);
2425
2426 if (format != cLoadTypeMRC && skew) {
2427 double matrix[16];
2428 auto i_float = (const float *) i;
2429
2430 // SKWMAT Skew matrix S (in order S11, S12, S13, S21 etc)
2431 copy33f44d(i_float, matrix);
2432 xx_matrix_invert(matrix, matrix, 4);
2433
2434 // SKWTRN Skew translation t
2435 matrix[3] = i_float[9];
2436 matrix[7] = i_float[10];
2437 matrix[11] = i_float[11];
2438
2439 // Xo(map) = S * (Xo(atoms) - t)
2440 ObjectStateSetMatrix(ms, matrix);
2441
2442 PRINTFB(I->G, FB_ObjectMap, FB_Details)
2443 " ObjectMapCCP4: Applied skew transformation\n"
2444 ENDFB(I->G);
2445 }
2446 }
2447
2448 // XORIGIN, YORIGIN, ZORIGIN (50 - 52)
2449 // TODO See "Origin Conventions" in http://situs.biomachina.org/fmap.pdf
2450 float * mrc2000origin = (float *)(i + 49 - 25);
2451 if (format != cLoadTypeCCP4Str && lengthsq3f(mrc2000origin) > R_SMALL4) {
2452 double matrix[] = {
2453 1., 0., 0., mrc2000origin[0],
2454 0., 1., 0., mrc2000origin[1],
2455 0., 0., 1., mrc2000origin[2],
2456 0., 0., 0., 1.};
2457
2458 ObjectStateSetMatrix(ms, matrix);
2459
2460 if (!quiet) {
2461 PRINTFB(I->G, FB_ObjectMap, FB_Details)
2462 " ObjectMapCCP4: Applied MRC 2000 ORIGIN %.2f %.2f %.2f\n",
2463 mrc2000origin[0], mrc2000origin[1], mrc2000origin[2]
2464 ENDFB(I->G);
2465
2466 if (ncstart || nrstart || nsstart) {
2467 PRINTFB(I->G, FB_ObjectMap, FB_Details)
2468 " ObjectMapCCP4: Ignoring N*START\n" ENDFB(I->G);
2469 }
2470 }
2471
2472 ncstart = 0;
2473 nrstart = 0;
2474 nsstart = 0;
2475 }
2476
2477 i += 54 - 25;
2478 stdev = *(float *) (i++);
2479
2480 if(!quiet) {
2481 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2482 " ObjectMapCCP4: AMIN %f AMAX %f AMEAN %f ARMS %f\n", mind, maxd, mean, stdev ENDFB(I->G);
2483 }
2484
2485 n_pts = nc * ns * nr;
2486
2487 /* at least one EM map encountered lacks NZ, so we'll try to guess it */
2488
2489 if((nx == ny) && (!nz)) {
2490 nz = nx;
2491
2492 if(!quiet) {
2493 PRINTFB(I->G, FB_ObjectMap, FB_Warnings)
2494 " ObjectMapCCP4: NZ value is zero, but NX = NY, so we'll guess NZ = NX = NY.\n"
2495 ENDFB(I->G);
2496 }
2497 }
2498
2499 expectation = sym_skip + sizeof(int) * 256 + bytes_per_pt * n_pts;
2500
2501 if(!quiet) {
2502 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2503 " ObjectMapCCP4: sym_skip %d bytes %d expectation %d\n",
2504 sym_skip, bytes, expectation ENDFB(I->G);
2505 }
2506
2507 if(bytes < expectation) {
2508 if(bytes == (expectation - sym_skip)) {
2509 /* accomodate bogus CCP4 map files with bad symmetry length information */
2510 if(!quiet) {
2511 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
2512 " ObjectMapCCP4: Map has invalid symmetry length -- working around.\n"
2513 ENDFB(I->G);
2514 }
2515
2516 expectation -= sym_skip;
2517 sym_skip = 0;
2518
2519 } else {
2520 PRINTFB(I->G, FB_ObjectMap, FB_Errors)
2521 " ObjectMapCCP4: Map appears to be truncated -- aborting.\n" ENDFB(I->G);
2522 return (0);
2523 }
2524 }
2525
2526 q = p + (sizeof(int) * 256) + sym_skip;
2527
2528 if(little_endian != map_endian && bytes_per_pt > 1) {
2529 swap_endian(q, n_pts, bytes_per_pt);
2530 }
2531
2532 // with normalize == 2, use mean and stdev from file header
2533 if(normalize == 1 && n_pts > 1) {
2534 c = n_pts;
2535 sum = 0.0;
2536 sumsq = 0.0;
2537 while(c--) {
2538 dens = ccp4_next_value(&q, map_mode);
2539 sumsq += dens * dens;
2540 sum += dens;
2541 }
2542 mean = (float) (sum / n_pts);
2543 stdev = (float) sqrt1d((sumsq - (sum * sum / n_pts)) / (n_pts - 1));
2544 if(stdev < 0.000001)
2545 stdev = 1.0;
2546 }
2547
2548 q = p + (sizeof(int) * 256) + sym_skip;
2549 mapc--; /* convert to C indexing... */
2550 mapr--;
2551 maps--;
2552
2553 ms->Div[0] = nx;
2554 ms->Div[1] = ny;
2555 ms->Div[2] = nz;
2556
2557 ms->FDim[mapc] = nc;
2558 ms->Min[mapc] = ncstart;
2559 ms->Max[mapc] = nc + ncstart - 1;
2560
2561 ms->FDim[mapr] = nr;
2562 ms->Min[mapr] = nrstart;
2563 ms->Max[mapr] = nr + nrstart - 1;
2564
2565 ms->FDim[maps] = ns;
2566 ms->Min[maps] = nsstart;
2567 ms->Max[maps] = ns + nsstart - 1;
2568
2569 if(!quiet) {
2570 if(Feedback(I->G, FB_ObjectMap, FB_Blather)) {
2571 dump3i(ms->Div, " ObjectMapCCP4: ms->Div");
2572 dump3i(ms->Min, " ObjectMapCCP4: ms->Min");
2573 dump3i(ms->Max, " ObjectMapCCP4: ms->Max");
2574 dump3i(ms->FDim, " ObjectMapCCP4: ms->FDim");
2575 }
2576 }
2577
2578 ms->Symmetry->Crystal.Dim[0] = xlen;
2579 ms->Symmetry->Crystal.Dim[1] = ylen;
2580 ms->Symmetry->Crystal.Dim[2] = zlen;
2581
2582 ms->Symmetry->Crystal.Angle[0] = alpha;
2583 ms->Symmetry->Crystal.Angle[1] = beta;
2584 ms->Symmetry->Crystal.Angle[2] = gamma;
2585
2586 if(ispg < n_space_group_numbers) {
2587 UtilNCopy(ms->Symmetry->SpaceGroup, space_group_numbers[ispg], WordLength);
2588 }
2589
2590 /* -- JV; for vol
2591 ms->Mean = mean;
2592 ms->SD = stdev;
2593 ms->MaxValue = maxd;
2594 ms->MinValue = mind;
2595 */
2596
2597 ms->FDim[3] = 3;
2598 if(!(ms->FDim[0] && ms->FDim[1] && ms->FDim[2]))
2599 ok = false;
2600 else {
2601 SymmetryUpdate(ms->Symmetry.get());
2602 /* CrystalDump(ms->Crystal); */
2603 ms->Field.reset(new Isofield(I->G, ms->FDim));
2604 ms->MapSource = cMapSourceCCP4;
2605 ms->Field->save_points = false;
2606
2607 maxd = -FLT_MAX;
2608 mind = FLT_MAX;
2609
2610 for(cc[maps] = 0; cc[maps] < ms->FDim[maps]; cc[maps]++) {
2611 v[maps] = (cc[maps] + ms->Min[maps]) / ((float) ms->Div[maps]);
2612
2613 for(cc[mapr] = 0; cc[mapr] < ms->FDim[mapr]; cc[mapr]++) {
2614 v[mapr] = (cc[mapr] + ms->Min[mapr]) / ((float) ms->Div[mapr]);
2615
2616 for(cc[mapc] = 0; cc[mapc] < ms->FDim[mapc]; cc[mapc]++) {
2617 v[mapc] = (cc[mapc] + ms->Min[mapc]) / ((float) ms->Div[mapc]);
2618
2619 dens = ccp4_next_value(&q, map_mode);
2620
2621 if(normalize)
2622 dens = (dens - mean) / stdev;
2623 F3(ms->Field->data, cc[0], cc[1], cc[2]) = dens;
2624 if(maxd < dens)
2625 maxd = dens;
2626 if(mind > dens)
2627 mind = dens;
2628 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
2629 for(e = 0; e < 3; e++)
2630 F4(ms->Field->points, cc[0], cc[1], cc[2], e) = vr[e];
2631 }
2632 }
2633 }
2634 }
2635 if(ok) {
2636 d = 0;
2637 for(c = 0; c < 2; c++) {
2638 v[2] = (c * (ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
2639 for(b = 0; b < 2; b++) {
2640 v[1] = (b * (ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
2641 for(a = 0; a < 2; a++) {
2642 v[0] = (a * (ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
2643 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
2644 copy3f(vr, ms->Corner + 3 * d);
2645 d++;
2646 }
2647 }
2648 }
2649 }
2650
2651 if(!quiet) {
2652 PRINTFB(I->G, FB_ObjectMap, FB_Details)
2653 " ObjectMapCCP4: Map Size %d x %d x %d\n", ms->FDim[0], ms->FDim[1], ms->FDim[2]
2654 ENDFB(I->G);
2655 }
2656
2657 if(!quiet) {
2658 if(n_pts > 1) {
2659 if(normalize) {
2660 PRINTFB(I->G, FB_ObjectMap, FB_Details)
2661 " ObjectMapCCP4: Normalizing with mean = %8.6f and stdev = %8.6f.\n",
2662 mean, stdev ENDFB(I->G);
2663 } else {
2664
2665 PRINTFB(I->G, FB_ObjectMap, FB_Details)
2666 " ObjectMapCCP4: Map will not be normalized.\n ObjectMapCCP4: Current mean = %8.6f and stdev = %8.6f.\n",
2667 mean, stdev ENDFB(I->G);
2668 }
2669 }
2670 }
2671
2672 if(ok) {
2673 v[2] = (ms->Min[2]) / ((float) ms->Div[2]);
2674 v[1] = (ms->Min[1]) / ((float) ms->Div[1]);
2675 v[0] = (ms->Min[0]) / ((float) ms->Div[0]);
2676
2677 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMin);
2678
2679 v[2] = ((ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
2680 v[1] = ((ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
2681 v[0] = ((ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
2682
2683 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMax);
2684 }
2685 #ifdef _UNDEFINED
2686 printf("%d %d %d %d %d %d %d %d %d\n",
2687 ms->Div[0],
2688 ms->Min[0],
2689 ms->Max[0],
2690 ms->Div[1], ms->Min[1], ms->Max[1], ms->Div[2], ms->Min[2], ms->Max[2]);
2691 printf("Okay? %d\n", ok);
2692 fflush(stdout);
2693 #endif
2694 if(!ok) {
2695 ErrMessage(I->G, "ObjectMap", "Error reading map");
2696 } else {
2697 ms->Active = true;
2698 ObjectMapUpdateExtents(I);
2699 if(!quiet) {
2700 PRINTFB(I->G, FB_ObjectMap, FB_Results)
2701 " ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(I->G);
2702 }
2703 }
2704
2705 return (ok);
2706 }
2707
2708
2709 /*========================================================================*/
getObjectMapState(PyMOLGlobals * G,const char * name,int state)2710 ObjectMapState * getObjectMapState(PyMOLGlobals * G, const char * name, int state) {
2711 auto* om = ExecutiveFindObject<ObjectMap>(G, name);
2712 return om ? om->getObjectMapState(state) : nullptr;
2713 }
2714
2715
2716 /*========================================================================*/
2717 /**
2718 * Export to CCP4 or MRC map format.
2719 *
2720 * CCP4:
2721 * - can have SKWMAT/SKWTRN, but almost no other progam supports them (Maestro does, to some extent)
2722 * - can't have ORIGIN
2723 *
2724 * MRC:
2725 * - can have ORIGIN, most programs support it
2726 * - can't have SKWMAT/SKWTRN
2727 * - can't have N*START (or shouldn't, because of weird ORIGIN conventions)
2728 * - can't have non-orthogonal axes (or shouldn't, because others might ignore it)
2729 */
ObjectMapStateToCCP4Str(const ObjectMapState * ms,int quiet,int format)2730 std::vector<char> ObjectMapStateToCCP4Str(const ObjectMapState * ms, int quiet, int format)
2731 {
2732 std::vector<char> buffer;
2733
2734 if (!ms || !ms->Active)
2735 return buffer; // empty
2736
2737 auto G = ms->G;
2738 auto field = ms->Field->data;
2739
2740 if (field->type != cFieldFloat ||
2741 field->base_size != 4) {
2742 PRINTFB(G, FB_ObjectMap, FB_Errors)
2743 " MapStateToCCP4-Error: Unsupported field type\n" ENDFB(G);
2744 return buffer; // empty
2745 }
2746
2747 switch (format) {
2748 case cLoadTypeCCP4Map:
2749 format = cLoadTypeCCP4Str;
2750 case cLoadTypeCCP4Str:
2751 case cLoadTypeCCP4Unspecified:
2752 case cLoadTypeMRC:
2753 break;
2754 default:
2755 ErrMessage(G, __func__, "wrong format");
2756 return {};
2757 }
2758
2759 buffer.resize(1024 + field->size(), 0);
2760 auto buffer_s = reinterpret_cast<char*>(&buffer.front());
2761 auto buffer_i = reinterpret_cast<int32_t*>(&buffer.front());
2762 auto buffer_f = reinterpret_cast<float*>(&buffer.front());
2763
2764 buffer_i[0] = ms->FDim[2]; // NC / NX
2765 buffer_i[1] = ms->FDim[1]; // NR / NY
2766 buffer_i[2] = ms->FDim[0]; // NS / NZ
2767 buffer_i[3] = 2; // MODE (2 = 32-bit reals)
2768 buffer_i[4] = ms->Min[2]; // NCSTART / NXSTART
2769 buffer_i[5] = ms->Min[1]; // NRSTART / NYSTART
2770 buffer_i[6] = ms->Min[0]; // NSSTART / NZSTART
2771 buffer_i[7] = ms->Div[0]; // NX / MX
2772 buffer_i[8] = ms->Div[1]; // NY / MY
2773 buffer_i[9] = ms->Div[2]; // NZ / MZ
2774
2775 if (ms->Div[0] == 0) {
2776 // fallback
2777 buffer_i[7] = ms->FDim[0] - 1; // NX / MX
2778 buffer_i[8] = ms->FDim[1] - 1; // NY / MY
2779 buffer_i[9] = ms->FDim[2] - 1; // NZ / MZ
2780 }
2781
2782 bool cell_fallback = true;
2783
2784 // Cell
2785 if (ms->Symmetry) {
2786 auto crystal = ms->Symmetry->Crystal;
2787 buffer_f[10] = crystal.Dim[0]; // X length / CELL A
2788 buffer_f[11] = crystal.Dim[1]; // Y length
2789 buffer_f[12] = crystal.Dim[2]; // Z length
2790 buffer_f[13] = crystal.Angle[0]; // Alpha / CELL B
2791 buffer_f[14] = crystal.Angle[1]; // Beta
2792 buffer_f[15] = crystal.Angle[2]; // Gamma
2793
2794 // check for 1x1x1 dummy cell
2795 cell_fallback = fabs(lengthsq3f(crystal.Dim) - 3.f) < R_SMALL4;
2796 }
2797
2798 if (cell_fallback) {
2799 // fallback: orthogonal extents box
2800 subtract3f(ms->ExtentMax, ms->ExtentMin, buffer_f + 10); // CELL A
2801 set3f(buffer_f + 13, 90.f, 90.f, 90.f); // CELL B
2802 }
2803
2804 buffer_i[16] = 3; // MAPC
2805 buffer_i[17] = 2; // MAPR
2806 buffer_i[18] = 1; // MAPS
2807
2808 // dummy values
2809 buffer_f[19] = -5.f; // AMIN
2810 buffer_f[20] = 5.f; // AMAX
2811 buffer_f[21] = 0.f; // AMEAN
2812
2813 // Space group number
2814 if (ms->Symmetry) {
2815 for (int ispg = 0; ispg < n_space_group_numbers; ++ispg) {
2816 if (strcmp(ms->Symmetry->SpaceGroup, space_group_numbers[ispg]) == 0) {
2817 buffer_i[22] = ispg; // ISPG
2818 break;
2819 }
2820 }
2821 }
2822
2823 buffer_i[23] = 0; // NSYMBT
2824
2825 float* b_skwmat = buffer_f + 25;
2826 float* b_skwtrn = buffer_f + 34;
2827 float* b_origin = buffer_f + 49;
2828
2829 // MRC typically only supports orthogonal axes
2830 const float _f3_90[3]{90.f, 90.f, 90.f};
2831 bool mrc_possible =
2832 format != cLoadTypeCCP4Str && equal3f(buffer_f + 13 /* CELL B */, _f3_90);
2833
2834 if (format == cLoadTypeMRC && !mrc_possible) {
2835 PRINTFB(G, FB_ObjectMap, FB_Warnings)
2836 " %s-Warning: MRC expects orthogonal axes\n", __func__ ENDFB(G);
2837 }
2838
2839 // skew transformation
2840 if (!ms->Matrix.empty()) {
2841 double m[16];
2842 copy44d(ms->Matrix.data(), m);
2843
2844 // Skew translation t
2845 set3f(buffer_f + 34, m[3], m[7], m[11]); // SKWTRN
2846
2847 // Skew matrix S (in order S11, S12, S13, S21 etc)
2848 m[3] = m[7] = m[11] = 0.;
2849 xx_matrix_invert(m, m, 4);
2850 copy44d33f(m, buffer_f + 25); // SKWMAT
2851
2852 if (mrc_possible) {
2853 mrc_possible = is_identityf(3, b_skwmat);
2854
2855 if (format == cLoadTypeMRC && !mrc_possible) {
2856 PRINTFB(G, FB_ObjectMap, FB_Warnings)
2857 " %s-Warning: MRC expects orthonormal map\n", __func__ ENDFB(G);
2858 }
2859 }
2860
2861 if (mrc_possible) {
2862 // translation only -> use origin instead of skew matrix
2863 copy3(b_skwtrn, b_origin); // MRC ORIGIN
2864 zero3(b_skwtrn);
2865 } else if (format != cLoadTypeMRC) {
2866 buffer_i[24] = 1; // LSKFLG
2867 }
2868 }
2869
2870 // origin (stored with skew transformation)
2871 if (!ms->Origin.empty() && lengthsq3f(ms->Origin.data()) > R_SMALL4) {
2872 if (!mrc_possible && !buffer_i[24] /* LSKFLG */) {
2873 identity33f(b_skwmat);
2874 buffer_i[24] = 1; // LSKFLG
2875 }
2876
2877 if (buffer_i[24] /* LSKFLG */) {
2878 add3f(ms->Origin.data(), b_skwtrn, b_skwtrn); // add to SKWTRN
2879 } else {
2880 add3f(ms->Origin.data(), b_origin, b_origin); // add to MRC ORIGIN
2881 }
2882 }
2883
2884 // Don't write ORIGIN and N*START. The convention is to ignore N*START if
2885 // ORIGIN != (0,0,0), likely because IMOD uses ORIGIN but ignores N*START.
2886 if (lengthsq3f(b_origin) > R_SMALL4 || format == cLoadTypeMRC) {
2887 // data origin in fractional coordinates
2888 float n_start_frac[3] = {
2889 float(buffer_i[6] /* NSSTART */) / (buffer_i[7] /* NX */),
2890 float(buffer_i[5] /* NRSTART */) / (buffer_i[8] /* NY */),
2891 float(buffer_i[4] /* NCSTART */) / (buffer_i[9] /* NZ */)};
2892
2893 // data origin in Angstrom
2894 float n_start_real[3];
2895 transform33f3f(ms->Symmetry->Crystal.FracToReal, n_start_frac,
2896 n_start_real);
2897
2898 // add to MRC origin
2899 add3f(n_start_real, b_origin, b_origin);
2900
2901 // get rid of N*START
2902 buffer_i[4] = 0;
2903 buffer_i[5] = 0;
2904 buffer_i[6] = 0;
2905 }
2906
2907 if (format == cLoadTypeCCP4Unspecified) {
2908 format = mrc_possible ? cLoadTypeMRC : cLoadTypeCCP4Str;
2909 }
2910
2911 // Character string 'MAP ' to identify file type
2912 memcpy(buffer_s + 52 * 4, "MAP ", 4); // MAP
2913
2914 // Machine stamp indicating the machine type which wrote file
2915 if (1 == *(int32_t*)"\1\0\0\0") {
2916 buffer_i[53] = 0x00004144; // MACHST (little endian)
2917 } else {
2918 buffer_i[53] = 0x11110000; // MACHST (big endian)
2919 }
2920
2921 buffer_f[54] = 1.f; // ARMS
2922 buffer_i[55] = 1; // NLABL
2923 sprintf(buffer_s + 56 * 4, "PyMOL %s format=%s", _PyMOL_VERSION,
2924 (format == cLoadTypeMRC) ? "MRC" : "CCP4"); // 57-256 LABEL(20,10)
2925
2926 // Map data array follows
2927 memcpy(buffer_s + 1024, field->data.data(), field->size());
2928
2929 return buffer;
2930 }
2931
2932
2933 /*========================================================================*/
ObjectMapPHIStrToMap(ObjectMap * I,char * PHIStr,int bytes,int state,int quiet)2934 static int ObjectMapPHIStrToMap(ObjectMap * I, char *PHIStr, int bytes, int state,
2935 int quiet)
2936 {
2937
2938 char *p;
2939 float dens, dens_rev;
2940 int a, b, c, d, e;
2941 float v[3], maxd, mind;
2942 int ok = true;
2943 int little_endian = 1;
2944 /* PHI named from their docs */
2945 int map_endian = 0;
2946 int map_dim;
2947 int map_bytes;
2948
2949 ObjectMapState *ms;
2950
2951 char cc[MAXLINELEN];
2952 char *rev;
2953
2954 little_endian = *((char *) &little_endian);
2955
2956 if(state < 0)
2957 state = I->State.size();
2958 if(I->State.size() <= state) {
2959 VecCheckEmplace(I->State, state, I->G);
2960 }
2961 ms = &I->State[state];
2962
2963 maxd = -FLT_MAX;
2964 mind = FLT_MAX;
2965 p = PHIStr;
2966
2967 if(*p) /* use FORMATTED IO record to determine map endiness */
2968 map_endian = 1;
2969 else
2970 map_endian = 0;
2971
2972 p += 4;
2973
2974 ParseNCopy(cc, p, 20);
2975 if(!quiet) {
2976 PRINTFB(I->G, FB_ObjectMap, FB_Details)
2977 " PHIStrToMap: %s\n", cc ENDFB(I->G);
2978 }
2979 p += 20;
2980 p += 4;
2981
2982 p += 4;
2983 ParseNCopy(cc, p, 10);
2984 if(!quiet) {
2985 PRINTFB(I->G, FB_ObjectMap, FB_Details)
2986 " PHIStrToMap: %s\n", cc ENDFB(I->G);
2987 }
2988 p += 10;
2989 ParseNCopy(cc, p, 60);
2990 if(!quiet) {
2991 PRINTFB(I->G, FB_ObjectMap, FB_Details)
2992 " PHIStrToMap: %s\n", cc ENDFB(I->G);
2993 }
2994 p += 60;
2995 p += 4;
2996
2997 rev = (char *) &dens_rev;
2998
2999 if(little_endian != map_endian) {
3000 rev[0] = p[3];
3001 rev[1] = p[2];
3002 rev[2] = p[1];
3003 rev[3] = p[0];
3004 } else {
3005 rev[0] = p[0]; /* gotta go char by char because of memory alignment issues ... */
3006 rev[1] = p[1];
3007 rev[2] = p[2];
3008 rev[3] = p[3];
3009 }
3010
3011 map_bytes = *((int *) rev);
3012
3013 map_dim = (int) (pow((map_bytes / 4.0), 1 / 3.0) + 0.5);
3014
3015 if((4 * map_dim * map_dim * map_dim) != map_bytes) /* consistency check */
3016 map_dim = 65;
3017
3018 if(!quiet) {
3019 PRINTFB(I->G, FB_ObjectMap, FB_Details)
3020 " PHIStrToMap: Map Size %d x %d x %d\n", map_dim, map_dim, map_dim ENDFB(I->G);
3021 }
3022 p += 4;
3023
3024 ms->FDim[0] = map_dim;
3025 ms->FDim[1] = map_dim;
3026 ms->FDim[2] = map_dim;
3027 ms->FDim[3] = 3;
3028
3029 ms->Div[0] = (map_dim - 1) / 2;
3030 ms->Div[1] = (map_dim - 1) / 2;
3031 ms->Div[2] = (map_dim - 1) / 2;
3032 ms->Min[0] = -ms->Div[0];
3033 ms->Min[1] = -ms->Div[1];
3034 ms->Min[2] = -ms->Div[2];
3035 ms->Max[0] = ms->Div[0];
3036 ms->Max[1] = ms->Div[1];
3037 ms->Max[2] = ms->Div[2];
3038
3039 ms->Field.reset(new Isofield(I->G, ms->FDim));
3040 ms->MapSource = cMapSourceGeneralPurpose;
3041 ms->Field->save_points = false;
3042
3043 for(c = 0; c < ms->FDim[2]; c++) { /* z y x ordering into c b a so that x = a, etc. */
3044 for(b = 0; b < ms->FDim[1]; b++) {
3045 for(a = 0; a < ms->FDim[0]; a++) {
3046
3047 if(little_endian != map_endian) {
3048 rev[0] = p[3];
3049 rev[1] = p[2];
3050 rev[2] = p[1];
3051 rev[3] = p[0];
3052 } else {
3053 rev[0] = p[0]; /* gotta go char by char because of memory alignment issues ... */
3054 rev[1] = p[1];
3055 rev[2] = p[2];
3056 rev[3] = p[3];
3057 }
3058 dens = *((float *) rev);
3059 F3(ms->Field->data, a, b, c) = dens;
3060 if(maxd < dens)
3061 maxd = dens;
3062 if(mind > dens)
3063 mind = dens;
3064 p += 4;
3065 }
3066 }
3067 }
3068 p += 4;
3069
3070 p += 4;
3071 ParseNCopy(cc, p, 16);
3072 if(!quiet) {
3073 PRINTFB(I->G, FB_ObjectMap, FB_Details)
3074 " PHIStrToMap: %s\n", cc ENDFB(I->G);
3075 }
3076 p += 16;
3077 p += 4;
3078
3079 ms->Grid = std::vector<float>(3);
3080 p += 4;
3081 if(little_endian != map_endian) {
3082 rev[0] = p[3];
3083 rev[1] = p[2];
3084 rev[2] = p[1];
3085 rev[3] = p[0];
3086 } else {
3087 rev[0] = p[0]; /* gotta go char by char because of memory alignment issues ... */
3088 rev[1] = p[1];
3089 rev[2] = p[2];
3090 rev[3] = p[3];
3091 }
3092 ms->Grid[0] = 1.0F / (*((float *) rev));
3093 ms->Grid[1] = ms->Grid[0];
3094 ms->Grid[2] = ms->Grid[0];
3095 p += 4;
3096
3097 ms->Origin = std::vector<float>(3);
3098 if(little_endian != map_endian) {
3099 rev[0] = p[3];
3100 rev[1] = p[2];
3101 rev[2] = p[1];
3102 rev[3] = p[0];
3103 } else {
3104 rev[0] = p[0]; /* gotta go char by char because of memory alignment issues ... */
3105 rev[1] = p[1];
3106 rev[2] = p[2];
3107 rev[3] = p[3];;
3108 }
3109 ms->Origin[0] = *((float *) rev);
3110 p += 4;
3111
3112 if(little_endian != map_endian) {
3113 rev[0] = p[3];
3114 rev[1] = p[2];
3115 rev[2] = p[1];
3116 rev[3] = p[0];
3117 } else {
3118 rev[0] = p[0]; /* gotta go char by char because of memory alignment issues ... */
3119 rev[1] = p[1];
3120 rev[2] = p[2];
3121 rev[3] = p[3];;
3122 }
3123 ms->Origin[1] = *((float *) rev);
3124 p += 4;
3125 if(little_endian != map_endian) {
3126 rev[0] = p[3];
3127 rev[1] = p[2];
3128 rev[2] = p[1];
3129 rev[3] = p[0];
3130 } else {
3131 rev[0] = p[0]; /* gotta go char by char because of memory alignment issues ... */
3132 rev[1] = p[1];
3133 rev[2] = p[2];
3134 rev[3] = p[3];;
3135 }
3136 ms->Origin[2] = *((float *) rev);
3137 p += 4;
3138
3139 p += 4;
3140
3141 if(ok) {
3142 for(e = 0; e < 3; e++) {
3143 ms->ExtentMin[e] = ms->Origin[e] + ms->Grid[e] * ms->Min[e];
3144 ms->ExtentMax[e] = ms->Origin[e] + ms->Grid[e] * ms->Max[e];
3145 }
3146 }
3147
3148 for(c = 0; c < ms->FDim[2]; c++) {
3149 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
3150 for(b = 0; b < ms->FDim[1]; b++) {
3151 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
3152 for(a = 0; a < ms->FDim[0]; a++) {
3153 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
3154 for(e = 0; e < 3; e++) {
3155 F4(ms->Field->points, a, b, c, e) = v[e];
3156 }
3157 }
3158 }
3159 }
3160
3161 d = 0;
3162 for(c = 0; c < ms->FDim[2]; c += ms->FDim[2] - 1) {
3163 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
3164
3165 for(b = 0; b < ms->FDim[1]; b += ms->FDim[1] - 1) {
3166 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
3167
3168 for(a = 0; a < ms->FDim[0]; a += ms->FDim[0] - 1) {
3169 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
3170 copy3f(v, ms->Corner + 3 * d);
3171 d++;
3172 }
3173 }
3174 }
3175
3176 /* interpolation test code
3177 {
3178 float test[3];
3179 float result;
3180 float cmp;
3181
3182 for(c=0;c<ms->FDim[0]-1;c++) {
3183 for(b=0;b<ms->FDim[1]-1;b++) {
3184 for(a=0;a<ms->FDim[2]-1;a++) {
3185 for(e=0;e<3;e++) {
3186 test[e] = (F4(ms->Field->points,a,b,c,e)+
3187 F4(ms->Field->points,a,b,c,e))/2.0;
3188 }
3189 ObjectMapStateInterpolate(ms,test,&result,1);
3190 cmp = (F3(ms->Field->data,a,b,c)+
3191 F3(ms->Field->data,a,b,c))/2.0;
3192 if(fabs(cmp-result)>0.001) {
3193 printf("%d %d %d\n",a,b,c);
3194 printf("%8.5f %8.5f\n",
3195 cmp,
3196 result);
3197 }
3198 }
3199 }
3200 }
3201 }
3202 */
3203
3204 if(!ok) {
3205 ErrMessage(I->G, "ObjectMap", "Error reading map");
3206 } else {
3207 ms->Active = true;
3208 ObjectMapUpdateExtents(I);
3209 if(!quiet) {
3210 PRINTFB(I->G, FB_ObjectMap, FB_Results)
3211 " ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(I->G);
3212 }
3213 }
3214 return (ok);
3215 }
3216
3217
3218 /*========================================================================*/
ObjectMapXPLORStrToMap(ObjectMap * I,char * XPLORStr,int state,int quiet)3219 static int ObjectMapXPLORStrToMap(ObjectMap * I, char *XPLORStr, int state, int quiet)
3220 {
3221
3222 char *p;
3223 int a, b, c, d, e;
3224 float v[3], vr[3], dens, maxd, mind;
3225 char cc[MAXLINELEN];
3226 int n;
3227 int ok = true;
3228 ObjectMapState *ms;
3229
3230 if(state < 0)
3231 state = I->State.size();
3232 if(I->State.size() <= state) {
3233 VecCheckEmplace(I->State, state, I->G);
3234 }
3235
3236 ms = &I->State[state];
3237
3238 maxd = -FLT_MAX;
3239 mind = FLT_MAX;
3240 p = XPLORStr;
3241
3242 while(*p) {
3243 p = ParseNCopy(cc, p, 8);
3244 if(!*cc)
3245 p = ParseNextLine(p);
3246 else if(sscanf(cc, "%i", &n) == 1) {
3247 p = ParseWordCopy(cc, p, MAXLINELEN);
3248 if(strstr(cc, "!NTITLE") || (!*cc)) {
3249 p = ParseNextLine(p);
3250 while(n--) {
3251 p = ParseNextLine(p);
3252 }
3253 } else if(strstr(cc, "REMARKS")) {
3254 p = ParseNextLine(p);
3255 } else {
3256 break;
3257 }
3258 }
3259 }
3260 if(*p) { /* n contains first dimension */
3261 ms->Div[0] = n;
3262 if(sscanf(cc, "%i", &ms->Min[0]) != 1)
3263 ok = false;
3264 p = ParseNCopy(cc, p, 8);
3265 if(sscanf(cc, "%i", &ms->Max[0]) != 1)
3266 ok = false;
3267 p = ParseNCopy(cc, p, 8);
3268 if(sscanf(cc, "%i", &ms->Div[1]) != 1)
3269 ok = false;
3270 p = ParseNCopy(cc, p, 8);
3271 if(sscanf(cc, "%i", &ms->Min[1]) != 1)
3272 ok = false;
3273 p = ParseNCopy(cc, p, 8);
3274 if(sscanf(cc, "%i", &ms->Max[1]) != 1)
3275 ok = false;
3276 p = ParseNCopy(cc, p, 8);
3277 if(sscanf(cc, "%i", &ms->Div[2]) != 1)
3278 ok = false;
3279 p = ParseNCopy(cc, p, 8);
3280 if(sscanf(cc, "%i", &ms->Min[2]) != 1)
3281 ok = false;
3282 p = ParseNCopy(cc, p, 8);
3283 if(sscanf(cc, "%i", &ms->Max[2]) != 1)
3284 ok = false;
3285 p = ParseNextLine(p);
3286 p = ParseNCopy(cc, p, 12);
3287 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Dim[0]) != 1)
3288 ok = false;
3289 p = ParseNCopy(cc, p, 12);
3290 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Dim[1]) != 1)
3291 ok = false;
3292 p = ParseNCopy(cc, p, 12);
3293 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Dim[2]) != 1)
3294 ok = false;
3295 p = ParseNCopy(cc, p, 12);
3296 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Angle[0]) != 1)
3297 ok = false;
3298 p = ParseNCopy(cc, p, 12);
3299 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Angle[1]) != 1)
3300 ok = false;
3301 p = ParseNCopy(cc, p, 12);
3302 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Angle[2]) != 1)
3303 ok = false;
3304 p = ParseNextLine(p);
3305 p = ParseNCopy(cc, p, 3);
3306 if(strcmp(cc, "ZYX"))
3307 ok = false;
3308 p = ParseNextLine(p);
3309
3310 } else {
3311 ok = false;
3312 }
3313 if(ok) {
3314
3315 ms->FDim[0] = ms->Max[0] - ms->Min[0] + 1;
3316 ms->FDim[1] = ms->Max[1] - ms->Min[1] + 1;
3317 ms->FDim[2] = ms->Max[2] - ms->Min[2] + 1;
3318 ms->FDim[3] = 3;
3319 if(!(ms->FDim[0] && ms->FDim[1] && ms->FDim[2]))
3320 ok = false;
3321 else {
3322 SymmetryUpdate(ms->Symmetry.get());
3323 ms->Field.reset(new Isofield(I->G, ms->FDim));
3324 ms->MapSource = cMapSourceCrystallographic;
3325 ms->Field->save_points = false;
3326 for(c = 0; c < ms->FDim[2]; c++) {
3327 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
3328 p = ParseNextLine(p);
3329 for(b = 0; b < ms->FDim[1]; b++) {
3330 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
3331 for(a = 0; a < ms->FDim[0]; a++) {
3332 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
3333 p = ParseNCopy(cc, p, 12);
3334 if(!cc[0]) {
3335 p = ParseNextLine(p);
3336 p = ParseNCopy(cc, p, 12);
3337 }
3338 if(sscanf(cc, "%f", &dens) != 1) {
3339 ok = false;
3340 } else {
3341 F3(ms->Field->data, a, b, c) = dens;
3342 if(maxd < dens)
3343 maxd = dens;
3344 if(mind > dens)
3345 mind = dens;
3346 }
3347 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
3348 for(e = 0; e < 3; e++) {
3349 F4(ms->Field->points, a, b, c, e) = vr[e];
3350 }
3351 }
3352 }
3353 p = ParseNextLine(p);
3354 }
3355 if(ok) {
3356 d = 0;
3357 for(c = 0; c < ms->FDim[2]; c += (ms->FDim[2] - 1)) {
3358 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
3359 for(b = 0; b < ms->FDim[1]; b += (ms->FDim[1] - 1)) {
3360 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
3361 for(a = 0; a < ms->FDim[0]; a += (ms->FDim[0] - 1)) {
3362 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
3363 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
3364 copy3f(vr, ms->Corner + 3 * d);
3365 d++;
3366 }
3367 }
3368 }
3369 }
3370 }
3371 }
3372
3373 if(ok) {
3374 v[2] = (ms->Min[2]) / ((float) ms->Div[2]);
3375 v[1] = (ms->Min[1]) / ((float) ms->Div[1]);
3376 v[0] = (ms->Min[0]) / ((float) ms->Div[0]);
3377
3378 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMin);
3379
3380 v[2] = ((ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
3381 v[1] = ((ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
3382 v[0] = ((ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
3383
3384 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMax);
3385
3386 }
3387 #ifdef _UNDEFINED
3388 printf("%d %d %d %d %d %d %d %d %d\n",
3389 ms->Div[0],
3390 ms->Min[0],
3391 ms->Max[0],
3392 ms->Div[1], ms->Min[1], ms->Max[1], ms->Div[2], ms->Min[2], ms->Max[2]);
3393 printf("Okay? %d\n", ok);
3394 fflush(stdout);
3395 #endif
3396 if(!ok) {
3397 ErrMessage(I->G, "ObjectMap", "Error reading map");
3398 } else {
3399 ms->Active = true;
3400 ObjectMapUpdateExtents(I);
3401 if(!quiet) {
3402 PRINTFB(I->G, FB_ObjectMap, FB_Results)
3403 " ObjectMap: Map read. Range = %5.3f to %5.3f\n", mind, maxd ENDFB(I->G);
3404 }
3405 }
3406
3407 return (ok);
3408 }
3409
3410
3411 /*========================================================================*/
ObjectMapFLDStrToMap(ObjectMap * I,char * PHIStr,int bytes,int state,int quiet)3412 static int ObjectMapFLDStrToMap(ObjectMap * I, char *PHIStr, int bytes, int state,
3413 int quiet)
3414 {
3415 char *p;
3416 float dens, dens_rev;
3417 int a, b, c, d, e;
3418 float v[3], maxd, mind;
3419 int ok = true;
3420 int little_endian = 1;
3421 /* PHI named from their docs */
3422 int map_endian = 1;
3423 char cc[MAXLINELEN];
3424 char *rev;
3425 int ndim = 0;
3426 int veclen = 0;
3427 int nspace = 0;
3428 ObjectMapState *ms;
3429 int got_ndim = false;
3430 int got_dim1 = false;
3431 int got_dim2 = false;
3432 int got_dim3 = false;
3433 int got_data = false;
3434 int got_veclen = false;
3435 int got_min_ext = false;
3436 int got_max_ext = false;
3437 int got_field = false;
3438 int got_nspace = false;
3439
3440 little_endian = *((char *) &little_endian);
3441 map_endian = little_endian;
3442
3443 if(state < 0)
3444 state = I->State.size();
3445 if(I->State.size() <= state) {
3446 VecCheckEmplace(I->State, state, I->G);
3447 }
3448 ms = &I->State[state];
3449
3450 maxd = -FLT_MAX;
3451 mind = FLT_MAX;
3452
3453 p = PHIStr;
3454
3455 while((*p) && (!(got_ndim &&
3456 got_dim1 &&
3457 got_dim2 &&
3458 got_dim3 &&
3459 got_veclen &&
3460 got_min_ext && got_max_ext && got_field && got_nspace && got_data))) {
3461
3462 if(!got_ndim) {
3463 ParseWordCopy(cc, p, 4);
3464 if(!strcmp(cc, "ndim")) {
3465 p = ParseSkipEquals(p);
3466 p = ParseWordCopy(cc, p, 50);
3467 if(sscanf(cc, "%d", &ndim) == 1)
3468 got_ndim = true;
3469 }
3470 }
3471
3472 if(!got_dim1) {
3473 ParseWordCopy(cc, p, 4);
3474 if(!strcmp(cc, "dim1")) {
3475 p = ParseSkipEquals(p);
3476 p = ParseWordCopy(cc, p, 50);
3477 if(sscanf(cc, "%d", &ms->FDim[0]) == 1)
3478 got_dim1 = true;
3479 }
3480 }
3481
3482 if(!got_dim2) {
3483 ParseWordCopy(cc, p, 4);
3484 if(!strcmp(cc, "dim2")) {
3485 p = ParseSkipEquals(p);
3486 p = ParseWordCopy(cc, p, 50);
3487 if(sscanf(cc, "%d", &ms->FDim[1]) == 1)
3488 got_dim2 = true;
3489 }
3490 }
3491
3492 if(!got_dim3) {
3493 ParseWordCopy(cc, p, 4);
3494 if(!strcmp(cc, "dim3")) {
3495 p = ParseSkipEquals(p);
3496 p = ParseWordCopy(cc, p, 50);
3497 if(sscanf(cc, "%d", &ms->FDim[2]) == 1)
3498 got_dim3 = true;
3499 }
3500 }
3501
3502 if(!got_nspace) {
3503 ParseWordCopy(cc, p, 6);
3504 if(!strcmp(cc, "nspace")) {
3505 p = ParseSkipEquals(p);
3506 p = ParseWordCopy(cc, p, 50);
3507 if(sscanf(cc, "%d", &nspace) == 1)
3508 got_nspace = true;
3509 }
3510 }
3511
3512 if(!got_veclen) {
3513 ParseWordCopy(cc, p, 6);
3514 if(!strcmp(cc, "veclen")) {
3515 p = ParseSkipEquals(p);
3516 p = ParseWordCopy(cc, p, 50);
3517 if(sscanf(cc, "%d", &veclen) == 1)
3518 got_veclen = true;
3519 }
3520 }
3521
3522 if(!got_data) {
3523 ParseWordCopy(cc, p, 4);
3524 if(!strcmp(cc, "data")) {
3525 p = ParseSkipEquals(p);
3526 p = ParseWordCopy(cc, p, 50);
3527 if(!strcmp(cc, "float"))
3528 got_data = true;
3529 }
3530 }
3531
3532 if(!got_min_ext) {
3533 ParseWordCopy(cc, p, 7);
3534 if(!strcmp(cc, "min_ext")) {
3535 p = ParseSkipEquals(p);
3536 p = ParseWordCopy(cc, p, 50);
3537 if(sscanf(cc, "%f", &ms->ExtentMin[0]) == 1) {
3538 p = ParseWordCopy(cc, p, 50);
3539 if(sscanf(cc, "%f", &ms->ExtentMin[1]) == 1) {
3540 p = ParseWordCopy(cc, p, 50);
3541 if(sscanf(cc, "%f", &ms->ExtentMin[2]) == 1) {
3542 got_min_ext = true;
3543 }
3544 }
3545 }
3546 }
3547 }
3548
3549 if(!got_max_ext) {
3550 ParseWordCopy(cc, p, 7);
3551 if(!strcmp(cc, "max_ext")) {
3552 p = ParseSkipEquals(p);
3553 p = ParseWordCopy(cc, p, 50);
3554 if(sscanf(cc, "%f", &ms->ExtentMax[0]) == 1) {
3555 p = ParseWordCopy(cc, p, 50);
3556 if(sscanf(cc, "%f", &ms->ExtentMax[1]) == 1) {
3557 p = ParseWordCopy(cc, p, 50);
3558 if(sscanf(cc, "%f", &ms->ExtentMax[2]) == 1) {
3559 got_max_ext = true;
3560 }
3561 }
3562 }
3563 }
3564 }
3565
3566 if(!got_field) {
3567 ParseWordCopy(cc, p, 5);
3568 if(!strcmp(cc, "field")) {
3569 p = ParseSkipEquals(p);
3570 p = ParseWordCopy(cc, p, 50);
3571 if(!strcmp(cc, "uniform"))
3572 got_field = true;
3573 }
3574 }
3575 p = ParseNextLine(p);
3576 }
3577
3578 if(got_ndim &&
3579 got_dim1 &&
3580 got_dim2 &&
3581 got_dim3 &&
3582 got_veclen && got_min_ext && got_max_ext && got_field && got_nspace && got_data) {
3583
3584 int pass = 0;
3585
3586 ms->Origin = std::vector<float>(3);
3587 ms->Range = std::vector<float>(3);
3588 ms->Grid = std::vector<float>(3);
3589
3590 copy3f(ms->ExtentMin, ms->Origin.data());
3591 subtract3f(ms->ExtentMax, ms->ExtentMin, ms->Range.data());
3592 ms->FDim[3] = 3;
3593
3594 PRINTFB(I->G, FB_ObjectMap, FB_Details)
3595 " FLDStrToMap: Map Size %d x %d x %d\n", ms->FDim[0], ms->FDim[1], ms->FDim[2]
3596 ENDFB(I->G);
3597
3598 for(a = 0; a < 3; a++) {
3599 ms->Min[a] = 0;
3600 ms->Max[a] = ms->FDim[a] - 1;
3601 ms->Div[a] = ms->FDim[a] - 1;
3602
3603 if(ms->FDim[a])
3604 ms->Grid[a] = ms->Range[a] / (ms->Max[a]);
3605 else
3606 ms->Grid[a] = 0.0F;
3607 }
3608
3609 ms->Field.reset(new Isofield(I->G, ms->FDim));
3610 ms->MapSource = cMapSourceFLD;
3611 ms->Field->save_points = false;
3612
3613 while(1) {
3614 maxd = -FLT_MAX;
3615 mind = FLT_MAX;
3616 p = PHIStr;
3617
3618 while(*p) { /* ^L^L sentinel */
3619 if((p[0] == 12) && (p[1] == 12)) {
3620 p += 2;
3621 break;
3622 }
3623 p++;
3624 }
3625
3626 rev = (char *) &dens_rev;
3627 for(c = 0; c < ms->FDim[2]; c++) { /* z y x ordering into c b a so that x = a, etc. */
3628 for(b = 0; b < ms->FDim[1]; b++) {
3629 for(a = 0; a < ms->FDim[0]; a++) {
3630
3631 if(little_endian != map_endian) {
3632 rev[0] = p[3];
3633 rev[1] = p[2];
3634 rev[2] = p[1];
3635 rev[3] = p[0];
3636 } else {
3637 rev[0] = p[0]; /* gotta go char by char because of memory alignment issues ... */
3638 rev[1] = p[1];
3639 rev[2] = p[2];
3640 rev[3] = p[3];
3641 }
3642 dens = *((float *) rev);
3643 F3(ms->Field->data, a, b, c) = dens;
3644 if(maxd < dens)
3645 maxd = dens;
3646 if(mind > dens)
3647 mind = dens;
3648 p += 4;
3649 }
3650 }
3651 }
3652
3653 /* There's no way to determine the original handedness of input
3654 field files. So instead, we simplymake an educated guess about
3655 whether we're byte-swapped based on the range of the density
3656 values obtained. */
3657
3658 if(((maxd / FLT_MAX) > 0.1F) && ((mind / (-FLT_MAX)) > 0.1F)) {
3659 if(pass == 0) {
3660 map_endian = (!map_endian); /* okay, try again swapped */
3661 } else if(pass == 1) {
3662 /* didn't help, so resort to original order */
3663 map_endian = (!map_endian);
3664 } else {
3665 break;
3666 }
3667 } else {
3668 break;
3669 }
3670 pass++;
3671 }
3672
3673 for(c = 0; c < ms->FDim[2]; c++) {
3674 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
3675 for(b = 0; b < ms->FDim[1]; b++) {
3676 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
3677 for(a = 0; a < ms->FDim[0]; a++) {
3678 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
3679 for(e = 0; e < 3; e++) {
3680 F4(ms->Field->points, a, b, c, e) = v[e];
3681 }
3682 }
3683 }
3684 }
3685
3686 d = 0;
3687 for(c = 0; c < ms->FDim[2]; c += ms->FDim[2] - 1) {
3688 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
3689 if(!c)
3690 v[2] = ms->ExtentMin[2];
3691 else
3692 v[2] = ms->ExtentMax[2];
3693
3694 for(b = 0; b < ms->FDim[1]; b += ms->FDim[1] - 1) {
3695 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
3696
3697 if(!b)
3698 v[1] = ms->ExtentMin[1];
3699 else
3700 v[1] = ms->ExtentMax[1];
3701
3702 for(a = 0; a < ms->FDim[0]; a += ms->FDim[0] - 1) {
3703 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
3704
3705 if(!a)
3706 v[0] = ms->ExtentMin[0];
3707 else
3708 v[0] = ms->ExtentMax[0];
3709
3710 copy3f(v, ms->Corner + 3 * d);
3711 d++;
3712 }
3713 }
3714 }
3715
3716 ms->Active = true;
3717 ObjectMapUpdateExtents(I);
3718 if(!quiet) {
3719 PRINTFB(I->G, FB_ObjectMap, FB_Results)
3720 " ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(I->G);
3721 }
3722 } else {
3723 PRINTFB(I->G, FB_ObjectMap, FB_Errors)
3724 " Error: unable to read FLD file.\n" ENDFB(I->G);
3725 /* printf(" got_ndim %d\n",got_ndim);
3726 printf(" got_dim1 %d\n",got_dim1);
3727 printf(" got_dim2 %d\n",got_dim2);
3728 printf(" got_dim3 %d\n",got_dim3);
3729 printf(" got_veclen %d\n",got_veclen);
3730 printf(" got_min_ext %d\n",got_min_ext);
3731 printf(" got_max_ext %d\n",got_max_ext);
3732 printf(" got_field %d\n",got_field);
3733 printf(" got_nspace %d\n",got_nspace);
3734 printf(" got_data %d\n",got_data);
3735 */
3736 }
3737
3738 return (ok);
3739
3740 }
3741
3742
3743 /*========================================================================*/
ObjectMapBRIXStrToMap(ObjectMap * I,char * BRIXStr,int bytes,int state,int quiet)3744 static int ObjectMapBRIXStrToMap(ObjectMap * I, char *BRIXStr, int bytes, int state,
3745 int quiet)
3746 {
3747 char *p, *pp;
3748 float dens;
3749 int a, b, c, d, e;
3750 float v[3], vr[3], maxd, mind;
3751 int ok = true;
3752 char cc[MAXLINELEN];
3753 ObjectMapState *ms;
3754 int got_origin = false;
3755 int got_extent = false;
3756 int got_grid = false;
3757 int got_cell = false;
3758 int got_prod = false;
3759 int got_plus = false;
3760 int got_sigma = false;
3761 int got_smiley = false;
3762 float prod = 1.0F;
3763 float plus = 0.0F;
3764 float sigma = 0.0F;
3765 float calc_sigma = 0.0F;
3766 float calc_mean = 0.0F;
3767 int normalize;
3768 int swap_bytes;
3769
3770 normalize = SettingGetGlobal_b(I->G, cSetting_normalize_o_maps);
3771 swap_bytes = SettingGetGlobal_b(I->G, cSetting_swap_dsn6_bytes);
3772 if(state < 0)
3773 state = I->State.size();
3774 if(I->State.size() <= state) {
3775 VecCheckEmplace(I->State, state, I->G);
3776 }
3777 ms = &I->State[state];
3778
3779 maxd = -FLT_MAX;
3780 mind = FLT_MAX;
3781
3782 p = BRIXStr;
3783
3784 ParseNCopy(cc, p, 3);
3785 got_smiley = (strcmp(cc, ":-)") == 0);
3786
3787 if(got_smiley) {
3788
3789 /* ASCII BRIX format header */
3790
3791 while((*p) && (!(got_origin &&
3792 got_extent &&
3793 got_grid && got_cell && got_prod && got_plus && got_sigma))) {
3794
3795 if(!got_origin) {
3796 pp = ParseWordCopy(cc, p, 6);
3797 if(WordMatch(I->G, "origin", cc, true) < 0) {
3798 p = ParseWordCopy(cc, pp, 50);
3799 if(sscanf(cc, "%d", &ms->Min[0]) == 1) {
3800 p = ParseWordCopy(cc, p, 50);
3801 if(sscanf(cc, "%d", &ms->Min[1]) == 1) {
3802 p = ParseWordCopy(cc, p, 50);
3803 if(sscanf(cc, "%d", &ms->Min[2]) == 1) {
3804 got_origin = true;
3805 }
3806 }
3807 }
3808 }
3809 }
3810
3811 if(!got_extent) {
3812 pp = ParseWordCopy(cc, p, 6);
3813 if(WordMatch(I->G, "extent", cc, true) < 0) {
3814 p = ParseWordCopy(cc, pp, 50);
3815 if(sscanf(cc, "%d", &ms->Max[0]) == 1) {
3816 p = ParseWordCopy(cc, p, 50);
3817 if(sscanf(cc, "%d", &ms->Max[1]) == 1) {
3818 p = ParseWordCopy(cc, p, 50);
3819 if(sscanf(cc, "%d", &ms->Max[2]) == 1) {
3820 got_extent = true;
3821 ms->Max[0] += ms->Min[0] - 1;
3822 ms->Max[1] += ms->Min[1] - 1;
3823 ms->Max[2] += ms->Min[2] - 1;
3824 }
3825 }
3826 }
3827 }
3828 }
3829
3830 if(!got_grid) {
3831 pp = ParseWordCopy(cc, p, 4);
3832 if(WordMatch(I->G, "grid", cc, true) < 0) {
3833 p = ParseWordCopy(cc, pp, 50);
3834 if(sscanf(cc, "%d", &ms->Div[0]) == 1) {
3835 p = ParseWordCopy(cc, p, 50);
3836 if(sscanf(cc, "%d", &ms->Div[1]) == 1) {
3837 p = ParseWordCopy(cc, p, 50);
3838 if(sscanf(cc, "%d", &ms->Div[2]) == 1) {
3839 got_grid = true;
3840 }
3841 }
3842 }
3843 }
3844 }
3845
3846 if(!got_cell) {
3847 pp = ParseWordCopy(cc, p, 4);
3848 if(WordMatch(I->G, "cell", cc, true) < 0) {
3849 p = ParseWordCopy(cc, pp, 50);
3850
3851 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Dim[0]) == 1) {
3852 p = ParseWordCopy(cc, p, 50);
3853 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Dim[1]) == 1) {
3854 p = ParseWordCopy(cc, p, 50);
3855 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Dim[2]) == 1) {
3856 p = ParseWordCopy(cc, p, 50);
3857 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Angle[0]) == 1) {
3858 p = ParseWordCopy(cc, p, 50);
3859 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Angle[1]) == 1) {
3860 p = ParseWordCopy(cc, p, 50);
3861 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Angle[2]) == 1) {
3862 got_cell = true;
3863 }
3864 }
3865 }
3866 }
3867 }
3868 }
3869 }
3870 }
3871
3872 if(!got_plus) {
3873 pp = ParseWordCopy(cc, p, 4);
3874 if(WordMatch(I->G, "plus", cc, true) < 0) {
3875 p = ParseWordCopy(cc, pp, 50);
3876 if(sscanf(cc, "%f", &plus) == 1) {
3877 got_plus = true;
3878 }
3879 }
3880 }
3881
3882 if(!got_prod) {
3883 pp = ParseWordCopy(cc, p, 4);
3884 if(WordMatch(I->G, "prod", cc, true) < 0) {
3885 p = ParseWordCopy(cc, pp, 50);
3886 if(sscanf(cc, "%f", &prod) == 1) {
3887 got_prod = true;
3888 }
3889 }
3890 }
3891
3892 if(!got_sigma) {
3893 pp = ParseWordCopy(cc, p, 5);
3894 if(WordMatch(I->G, "sigma", cc, true) < 0) {
3895 p = ParseWordCopy(cc, pp, 50);
3896 if(sscanf(cc, "%f", &sigma) == 1) {
3897 got_sigma = true;
3898 }
3899 }
3900 }
3901
3902 p++;
3903 }
3904 } else {
3905 /* Binary DSN6 format */
3906
3907 /* first, determine whether or not we need to byte-swap the header... */
3908
3909 int passed_endian_check = false;
3910 short int *shint_ptr;
3911 float scale_factor;
3912 char tmp_char;
3913 int a;
3914 int start_swap_at = 0;
3915 int stop_swap_at = bytes;
3916
3917 shint_ptr = (short int *) (p + 18 * 2);
3918 if(*shint_ptr == 100) {
3919 passed_endian_check = true;
3920 start_swap_at = 512; /* don't byte-swap header */
3921 }
3922
3923 if(!swap_bytes) {
3924 stop_swap_at = 512;
3925 }
3926 /* for DSN6, always swap the bricks */
3927
3928 p = BRIXStr + start_swap_at;
3929 for(a = start_swap_at; a < stop_swap_at; a += 2) {
3930 tmp_char = *p;
3931 *p = *(p + 1);
3932 *(p + 1) = tmp_char;
3933 p += 2;
3934 }
3935
3936 p = BRIXStr;
3937
3938 if(*shint_ptr == 100) {
3939 passed_endian_check = true;
3940 }
3941
3942 if(!passed_endian_check) {
3943 PRINTFB(I->G, FB_ObjectMap, FB_Errors)
3944 " Error: This looks like a DSN6 map file, but I can't match endianness.\n"
3945 ENDFB(I->G);
3946 } else {
3947 shint_ptr = (short int *) p;
3948
3949 ms->Min[0] = *(shint_ptr++);
3950 ms->Min[1] = *(shint_ptr++);
3951 ms->Min[2] = *(shint_ptr++);
3952 got_origin = true;
3953
3954 ms->Max[0] = *(shint_ptr++);
3955 ms->Max[1] = *(shint_ptr++);
3956 ms->Max[2] = *(shint_ptr++);
3957 got_extent = true;
3958 ms->Max[0] += ms->Min[0] - 1;
3959 ms->Max[1] += ms->Min[1] - 1;
3960 ms->Max[2] += ms->Min[2] - 1;
3961
3962 ms->Div[0] = *(shint_ptr++);
3963 ms->Div[1] = *(shint_ptr++);
3964 ms->Div[2] = *(shint_ptr++);
3965 got_grid = true;
3966
3967 ms->Symmetry->Crystal.Dim[0] = (float) (*(shint_ptr++));
3968 ms->Symmetry->Crystal.Dim[1] = (float) (*(shint_ptr++));
3969 ms->Symmetry->Crystal.Dim[2] = (float) (*(shint_ptr++));
3970 ms->Symmetry->Crystal.Angle[0] = (float) (*(shint_ptr++));
3971 ms->Symmetry->Crystal.Angle[1] = (float) (*(shint_ptr++));
3972 ms->Symmetry->Crystal.Angle[2] = (float) (*(shint_ptr++));
3973 got_cell = true;
3974
3975 prod = (float) (*(shint_ptr++)) / 100.0F;
3976 got_prod = true;
3977
3978 plus = (float) (*(shint_ptr++));
3979 got_plus = true;
3980
3981 scale_factor = (float) (*(shint_ptr++));
3982
3983 for(a = 0; a < 3; a++) {
3984 ms->Symmetry->Crystal.Dim[a] /= scale_factor;
3985 ms->Symmetry->Crystal.Angle[a] /= scale_factor;
3986 }
3987 }
3988 }
3989
3990 if(got_origin && got_extent && got_grid && got_cell && got_plus && got_prod) {
3991
3992 PRINTFB(I->G, FB_ObjectMap, FB_Blather)
3993 " BRIXStrToMap: Prod = %8.3f, Plus = %8.3f\n", prod, plus ENDFB(I->G);
3994
3995 ms->FDim[0] = ms->Max[0] - ms->Min[0] + 1;
3996 ms->FDim[1] = ms->Max[1] - ms->Min[1] + 1;
3997 ms->FDim[2] = ms->Max[2] - ms->Min[2] + 1;
3998 ms->FDim[3] = 3;
3999 if(!(ms->FDim[0] && ms->FDim[1] && ms->FDim[2]))
4000 ok = false;
4001 else {
4002 SymmetryUpdate(ms->Symmetry.get());
4003 ms->Field.reset(new Isofield(I->G, ms->FDim));
4004 ms->MapSource = cMapSourceBRIX;
4005 ms->Field->save_points = false;
4006
4007 {
4008 int block_size = 8;
4009 int xblocks = ((ms->FDim[0] - 1) / block_size) + 1;
4010 int yblocks = ((ms->FDim[1] - 1) / block_size) + 1;
4011 int zblocks = ((ms->FDim[2] - 1) / block_size) + 1;
4012 int cc, bb, aa, ia, ib, ic, xa, xb, xc;
4013 double sum = 0.0;
4014 double sumsq = 0.0;
4015 int n_pts = 0;
4016
4017 p = BRIXStr + 512;
4018 for(cc = 0; cc < zblocks; cc++) {
4019 for(bb = 0; bb < yblocks; bb++) {
4020 for(aa = 0; aa < xblocks; aa++) {
4021
4022 for(c = 0; c < block_size; c++) {
4023 xc = c + cc * block_size;
4024 ic = xc + ms->Min[2];
4025 v[2] = ic / ((float) ms->Div[2]);
4026
4027 for(b = 0; b < block_size; b++) {
4028 xb = b + bb * block_size;
4029 ib = xb + ms->Min[1];
4030 v[1] = ib / ((float) ms->Div[1]);
4031
4032 for(a = 0; a < block_size; a++) {
4033 xa = a + aa * block_size;
4034 ia = xa + ms->Min[0];
4035 v[0] = ia / ((float) ms->Div[0]);
4036
4037 dens = (((float) (*((unsigned char *) (p++)))) - plus) / prod;
4038 if((ia <= ms->Max[0]) && (ib <= ms->Max[1]) && (ic <= ms->Max[2])) {
4039 F3(ms->Field->data, xa, xb, xc) = dens;
4040 sumsq += dens * dens;
4041 sum += dens;
4042 n_pts++;
4043 if(maxd < dens)
4044 maxd = dens;
4045 if(mind > dens)
4046 mind = dens;
4047 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
4048 for(e = 0; e < 3; e++) {
4049 F4(ms->Field->points, xa, xb, xc, e) = vr[e];
4050 }
4051
4052 }
4053 }
4054 }
4055 }
4056 }
4057 }
4058 }
4059
4060 if(n_pts > 1) {
4061 calc_mean = (float) (sum / n_pts);
4062 calc_sigma = (float) sqrt1d((sumsq - (sum * sum / n_pts)) / (n_pts - 1));
4063 }
4064
4065 }
4066
4067 if(ok) {
4068 d = 0;
4069 for(c = 0; c < ms->FDim[2]; c += (ms->FDim[2] - 1)) {
4070 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
4071 for(b = 0; b < ms->FDim[1]; b += (ms->FDim[1] - 1)) {
4072 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
4073 for(a = 0; a < ms->FDim[0]; a += (ms->FDim[0] - 1)) {
4074 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
4075 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
4076 copy3f(vr, ms->Corner + 3 * d);
4077 d++;
4078 }
4079 }
4080 }
4081 }
4082
4083 if(ok && normalize) {
4084
4085 if(calc_sigma > R_SMALL8) {
4086 for(c = 0; c < ms->FDim[2]; c++) {
4087 for(b = 0; b < ms->FDim[1]; b++) {
4088 for(a = 0; a < ms->FDim[0]; a++) {
4089 F3(ms->Field->data, a, b, c) =
4090 (F3(ms->Field->data, a, b, c) - calc_mean) / calc_sigma;
4091 }
4092 }
4093 }
4094 }
4095 }
4096
4097 v[2] = (ms->Min[2]) / ((float) ms->Div[2]);
4098 v[1] = (ms->Min[1]) / ((float) ms->Div[1]);
4099 v[0] = (ms->Min[0]) / ((float) ms->Div[0]);
4100
4101 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMin);
4102
4103 v[2] = ((ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
4104 v[1] = ((ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
4105 v[0] = ((ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
4106
4107 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMax);
4108
4109 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4110 " BRIXStrToMap: Map Size %d x %d x %d\n", ms->FDim[0], ms->FDim[1], ms->FDim[2]
4111 ENDFB(I->G);
4112
4113 if(got_sigma) {
4114 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4115 " BRIXStrToMap: Reported Sigma = %8.3f\n", sigma ENDFB(I->G);
4116 }
4117
4118 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4119 " BRIXStrToMap: Range = %5.6f to %5.6f\n", mind, maxd ENDFB(I->G);
4120
4121 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4122 " BRIXStrToMap: Calculated Mean = %8.3f, Sigma = %8.3f\n", calc_mean, calc_sigma
4123 ENDFB(I->G);
4124
4125 if(normalize) {
4126 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4127 " BRIXStrToMap: Normalizing...\n" ENDFB(I->G);
4128 }
4129
4130 ms->Active = true;
4131 ObjectMapUpdateExtents(I);
4132
4133 }
4134 } else {
4135 PRINTFB(I->G, FB_ObjectMap, FB_Errors)
4136 " Error: unable to read BRIX/DSN6 file.\n" ENDFB(I->G);
4137 }
4138
4139 return (ok);
4140
4141 }
4142
4143
4144 /*========================================================================*/
ObjectMapGRDStrToMap(ObjectMap * I,char * GRDStr,int bytes,int state,int quiet)4145 static int ObjectMapGRDStrToMap(ObjectMap * I, char *GRDStr, int bytes, int state,
4146 int quiet)
4147 {
4148 /* NOTE: binary GRD reader not yet validated */
4149
4150 char *p;
4151 float dens;
4152 float *f = NULL;
4153 int a, b, c, d, e;
4154 float v[3], vr[3], maxd, mind;
4155 int ok = true;
4156 char cc[MAXLINELEN];
4157 ObjectMapState *ms;
4158 int got_cell = false;
4159 int got_brick = false;
4160 int fast_axis = 1; /* assumes fast X */
4161
4162 int got_ranges = false;
4163 int normalize = false;
4164 float mean, stdev;
4165 double sum = 0.0;
4166 double sumsq = 0.0;
4167 int n_pts = 0;
4168 int ascii = true;
4169 int little_endian = 1, map_endian = 0;
4170 union int_or_char_array {
4171 int block_len;
4172 char rev[4];
4173 };
4174 union int_or_char_array rev_union;
4175 rev_union.block_len = 0;
4176
4177 if(state < 0)
4178 state = I->State.size();
4179 if(I->State.size() <= state) {
4180 VecCheckEmplace(I->State, state, I->G);
4181 }
4182 ms = &I->State[state];
4183 normalize = SettingGetGlobal_b(I->G, cSetting_normalize_grd_maps);
4184 maxd = -FLT_MAX;
4185 mind = FLT_MAX;
4186
4187 p = GRDStr;
4188
4189 if(bytes > 4) {
4190 if((!p[0]) || (!p[1]) || (!p[2]) || (!p[3])) {
4191 /* if zero appears in the first four bytes, then this is a binary map */
4192 ascii = false;
4193
4194 little_endian = *((char *) &little_endian);
4195 map_endian = (*p || *(p + 1));
4196
4197 }
4198 }
4199 /* print map title */
4200
4201 if(ascii) {
4202 p = ParseNCopy(cc, p, 100);
4203 } else {
4204
4205
4206 /*
4207
4208 according to one site on the internet...GRD binary formats look like this:
4209
4210 write(14) title
4211 write(14) scale,oldmid
4212 write(14) ivary,nbyte,intdat,extent,extent
4213 write(14) extent,xang,yang,zang,xstart
4214 write(14) xend,ystart,yend,zstart,zend
4215 write(14) intx,inty,intz
4216
4217 do 100 i = 1,intz+1
4218 do 100 j = 1,inty+1
4219 100 write(14)(phimap(k,j,i),k=1,intx+1)
4220
4221 where: scale is the inverse grid spacing,
4222 oldmid are x,y,z coordinates of the grid center,
4223 intx,inty,intz = grid points per side - 1
4224 ivary = 0
4225 nbyte = 4
4226 intdat = 0
4227 xang,yang,zang = 90
4228 extent is the absolute value of the x,y, or z coordinate of the grid
4229 corners with the largest absolute value
4230 xstart,xend,ystart, yend,zstart,zend are the fractional limits (-1 to 1)
4231 of the molecule in each direction.
4232 title is a descriptive header for the molecule
4233
4234 However, that doesn't make sense...in the files I see, scale and oldmid don't seem to exist!
4235
4236 So here's another claim which I find more credible...
4237
4238 character*132 toplbl !ascii header
4239 integer*4 ivary !0 => x index varys most rapidly
4240 integer*4 nbyte !=4, # of bytes in data
4241 integer*4 inddat !=0, floating point data
4242 real*4 xang,yang,zang !=90,90,90 unit cell angles
4243 integer*4 intx,inty,intz !=igrid-1, # of intervals/grid side
4244 real*4 extent !maximum extent of grid
4245 real*4 xstart,xend !beginning, end of grid sides
4246 real*4 ystart,yend !in fractional
4247 real*4 zstart,zend !units of extent
4248 write(14)toplbl write(14)ivary, nbyte, intdat, extent, extent, extent, xang, yang, zang, xstart, xend, ystart, yend, zstart, zend, intx, inty, intz
4249 do k = 1,igrid
4250 do j = 1,igrid
4251 write(14)(phimap(i,j,k),i=1,igrid)
4252 end do
4253 end d
4254
4255 */
4256
4257 if(!quiet) {
4258 PRINTFB(I->G, FB_ObjectMap, FB_Warnings)
4259 " ObjectMapGRD-Warning: Binary GRD reader not yet validated.\n" ENDFB(I->G);
4260 }
4261
4262 if(little_endian != map_endian) {
4263 rev_union.rev[0] = p[3];
4264 rev_union.rev[1] = p[2];
4265 rev_union.rev[2] = p[1];
4266 rev_union.rev[3] = p[0];
4267 } else {
4268 rev_union.rev[0] = p[0]; /* gotta go char by char because of memory alignment issues ... */
4269 rev_union.rev[1] = p[1];
4270 rev_union.rev[2] = p[2];
4271 rev_union.rev[3] = p[3];
4272 }
4273 ParseNCopy(cc, p + 4, rev_union.block_len);
4274 /* now flip file to correct endianess */
4275 if(little_endian != map_endian) {
4276 char *c = p;
4277 int cnt = bytes >> 2;
4278 unsigned char c0, c1, c2, c3;
4279 while(cnt) {
4280 c0 = c[0];
4281 c1 = c[1];
4282 c2 = c[2];
4283 c3 = c[3];
4284 c[0] = c3;
4285 c[1] = c2;
4286 c[2] = c1;
4287 c[3] = c0;
4288 c += 4;
4289 cnt--;
4290 }
4291 }
4292 p += 4 + rev_union.block_len;
4293 f = (float *) p;
4294 }
4295
4296 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4297 " ObjectMap: %s\n", cc ENDFB(I->G);
4298
4299 if(ascii)
4300 p = ParseNextLine(p);
4301
4302 /* skip format -- we're reading float regardless... */
4303 if(ascii)
4304 p = ParseNextLine(p);
4305 else {
4306 {
4307 int block_len_check;
4308 int ivary;
4309 int nbyte;
4310 int intdat;
4311
4312 block_len_check = *((int *) (f++));
4313
4314 if(rev_union.block_len != block_len_check) {
4315 PRINTFB(I->G, FB_ObjectMap, FB_Warnings)
4316 " ObjectMapGRD-Warning: block length not matched -- not a true GRD binary?\n"
4317 ENDFB(I->G);
4318 }
4319
4320 rev_union.block_len = *((int *) (f++));
4321 ivary = *((int *) (f++));
4322 nbyte = *((int *) (f++));
4323 intdat = *((int *) (f++));
4324
4325 if((ivary) || (nbyte != 4) || (intdat)) {
4326 if(!quiet) {
4327 PRINTFB(I->G, FB_ObjectMap, FB_Warnings)
4328 " ObjectMapGRD-Warning: funky ivary, nbyte, intdat -- not a true GRD binary?\n"
4329 ENDFB(I->G);
4330 }
4331 }
4332 }
4333
4334 }
4335
4336 /* read unit cell */
4337
4338 if(ok) {
4339 if(ascii) {
4340 p = ParseWordCopy(cc, p, 50);
4341 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Dim[0]) == 1) {
4342 p = ParseWordCopy(cc, p, 50);
4343 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Dim[1]) == 1) {
4344 p = ParseWordCopy(cc, p, 50);
4345 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Dim[2]) == 1) {
4346 p = ParseWordCopy(cc, p, 50);
4347 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Angle[0]) == 1) {
4348 p = ParseWordCopy(cc, p, 50);
4349 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Angle[1]) == 1) {
4350 p = ParseWordCopy(cc, p, 50);
4351 if(sscanf(cc, "%f", &ms->Symmetry->Crystal.Angle[2]) == 1) {
4352 got_cell = true;
4353 }
4354 }
4355 }
4356 }
4357 }
4358 }
4359 ok = got_cell;
4360 } else {
4361 ms->Symmetry->Crystal.Dim[0] = *(f++); /* x-extent */
4362 ms->Symmetry->Crystal.Dim[1] = *(f++); /* y-extent */
4363 ms->Symmetry->Crystal.Dim[2] = *(f++); /* z-extent */
4364 ms->Symmetry->Crystal.Angle[0] = (*f++); /* xang */
4365 ms->Symmetry->Crystal.Angle[1] = (*f++); /* yang */
4366 ms->Symmetry->Crystal.Angle[2] = (*f++); /* zang */
4367 got_cell = 1;
4368 }
4369 }
4370
4371 if(ascii)
4372 p = ParseNextLine(p);
4373
4374 if(ascii) {
4375 /* read brick dimensions */
4376
4377 if(ok) {
4378 p = ParseWordCopy(cc, p, 50);
4379 if(sscanf(cc, "%d", &ms->FDim[0]) == 1) {
4380 p = ParseWordCopy(cc, p, 50);
4381 if(sscanf(cc, "%d", &ms->FDim[1]) == 1) {
4382 p = ParseWordCopy(cc, p, 50);
4383 if(sscanf(cc, "%d", &ms->FDim[2]) == 1) {
4384 got_brick = true;
4385 ms->FDim[0]++;
4386 ms->FDim[1]++;
4387 ms->FDim[2]++; /* stupid 0-based counters */
4388 }
4389 }
4390 }
4391 ok = got_brick;
4392 }
4393
4394 p = ParseNextLine(p);
4395
4396 /* read ranges */
4397
4398 if(ok) {
4399 p = ParseWordCopy(cc, p, 50);
4400 if(sscanf(cc, "%d", &fast_axis) == 1) {
4401 p = ParseWordCopy(cc, p, 50);
4402 if(sscanf(cc, "%d", &ms->Min[0]) == 1) {
4403 p = ParseWordCopy(cc, p, 50);
4404 if(sscanf(cc, "%d", &ms->Div[0]) == 1) {
4405 p = ParseWordCopy(cc, p, 50);
4406 if(sscanf(cc, "%d", &ms->Min[1]) == 1) {
4407 p = ParseWordCopy(cc, p, 50);
4408 if(sscanf(cc, "%d", &ms->Div[1]) == 1) {
4409 p = ParseWordCopy(cc, p, 50);
4410 if(sscanf(cc, "%d", &ms->Min[2]) == 1) {
4411 p = ParseWordCopy(cc, p, 50);
4412 if(sscanf(cc, "%d", &ms->Div[2]) == 1) {
4413 got_ranges = true;
4414 }
4415 }
4416 }
4417 }
4418 }
4419 }
4420 }
4421 ok = got_ranges;
4422 }
4423 } else {
4424 float fmin[3];
4425 float fmax[3];
4426
4427 fmin[0] = *(f++);
4428 fmax[0] = *(f++);
4429 fmin[1] = *(f++);
4430 fmax[1] = *(f++);
4431 fmin[2] = *(f++);
4432 fmax[2] = *(f++);
4433
4434 dump3f(fmin, "fmin");
4435 dump3f(fmax, "fmax");
4436 ms->FDim[0] = *((int *) f++) + 1;
4437 ms->FDim[1] = *((int *) f++) + 1;
4438 ms->FDim[2] = *((int *) f++) + 1;
4439
4440 ms->Div[0] = pymol_roundf((ms->FDim[0] - 1) / (fmax[0] - fmin[0]));
4441 ms->Div[1] = pymol_roundf((ms->FDim[1] - 1) / (fmax[1] - fmin[1]));
4442 ms->Div[2] = pymol_roundf((ms->FDim[2] - 1) / (fmax[2] - fmin[2]));
4443
4444 ms->Min[0] = pymol_roundf(ms->Div[0] * fmin[0]);
4445 ms->Min[1] = pymol_roundf(ms->Div[1] * fmin[1]);
4446 ms->Min[2] = pymol_roundf(ms->Div[2] * fmin[2]);
4447
4448 ms->Max[0] = ms->Min[0] + ms->FDim[0] - 1;
4449 ms->Max[1] = ms->Min[1] + ms->FDim[1] - 1;
4450 ms->Max[2] = ms->Min[2] + ms->FDim[2] - 1;
4451 got_ranges = true;
4452
4453 {
4454 int block_len_check;
4455
4456 block_len_check = *((int *) (f++));
4457 if(rev_union.block_len != block_len_check) {
4458 PRINTFB(I->G, FB_ObjectMap, FB_Warnings)
4459 " ObjectMapGRD-Warning: block length not matched -- not a true GRD binary?\n"
4460 ENDFB(I->G);
4461 }
4462 }
4463 }
4464
4465 if(ok) {
4466
4467 if(ascii) {
4468 ms->Div[0] = ms->Div[0] - ms->Min[0];
4469 ms->Div[1] = ms->Div[1] - ms->Min[1];
4470 ms->Div[2] = ms->Div[2] - ms->Min[2];
4471
4472 ms->Max[0] = ms->Min[0] + ms->FDim[0] - 1;
4473 ms->Max[1] = ms->Min[1] + ms->FDim[1] - 1;
4474 ms->Max[2] = ms->Min[2] + ms->FDim[2] - 1;
4475 }
4476
4477 ms->FDim[3] = 3;
4478
4479 if(Feedback(I->G, FB_ObjectMap, FB_Blather)) {
4480 dump3i(ms->Div, "ms->Div");
4481 dump3i(ms->Min, "ms->Min");
4482 dump3i(ms->Max, "ms->Max");
4483 dump3i(ms->FDim, "ms->FDim");
4484 }
4485
4486 SymmetryUpdate(ms->Symmetry.get());
4487 ms->Field.reset(new Isofield(I->G, ms->FDim));
4488 ms->MapSource = cMapSourceGRD;
4489 ms->Field->save_points = false;
4490
4491 switch (fast_axis) {
4492 case 3: /* Fast Y - BROKEN! */
4493 PRINTFB(I->G, FB_ObjectMap, FB_Warnings)
4494 " ObjectMapGRD-Warning: fast_axis %d unsupported!\n", fast_axis ENDFB(I->G);
4495 /* intentional fall though... */
4496 case 1: /* Fast X */
4497 default:
4498 for(c = 0; c < ms->FDim[2]; c++) {
4499 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
4500 for(b = 0; b < ms->FDim[1]; b++) {
4501 if(!ascii)
4502 f++; /* skip block delimiter */
4503 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
4504 for(a = 0; a < ms->FDim[0]; a++) {
4505 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
4506 if(ascii) {
4507 p = ParseNextLine(p);
4508 p = ParseNCopy(cc, p, 24);
4509 if(sscanf(cc, "%f", &dens) != 1) {
4510 ok = false;
4511 }
4512 } else {
4513 dens = *(f++);
4514 }
4515 if(ok) {
4516 sumsq += dens * dens;
4517 sum += dens;
4518 n_pts++;
4519 F3(ms->Field->data, a, b, c) = dens;
4520 if(maxd < dens)
4521 maxd = dens;
4522 if(mind > dens)
4523 mind = dens;
4524 }
4525 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
4526 for(e = 0; e < 3; e++) {
4527 F4(ms->Field->points, a, b, c, e) = vr[e];
4528 }
4529 }
4530 if(!ascii)
4531 f++; /* skip fortran block delimiter */
4532 }
4533 }
4534 break;
4535 }
4536 }
4537
4538 if(n_pts > 1) {
4539 mean = (float) (sum / n_pts);
4540 stdev = (float) sqrt1d((sumsq - (sum * sum / n_pts)) / (n_pts - 1));
4541
4542 if(normalize) {
4543 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4544 " ObjectMapGRDStrToMap: Normalizing: mean = %8.6f and stdev = %8.6f.\n",
4545 mean, stdev ENDFB(I->G);
4546 if(stdev < R_SMALL8)
4547 stdev = 1.0F;
4548 for(c = 0; c < ms->FDim[2]; c++)
4549 for(b = 0; b < ms->FDim[1]; b++) {
4550 for(a = 0; a < ms->FDim[0]; a++) {
4551 dens = F3(ms->Field->data, a, b, c);
4552 F3(ms->Field->data, a, b, c) = (dens - mean) / stdev;
4553 }
4554 }
4555 } else {
4556 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4557 " ObjectMapGRDStrToMap: Mean = %8.6f and stdev = %8.6f.\n",
4558 mean, stdev ENDFB(I->G);
4559 }
4560 }
4561
4562 if(ok) {
4563 d = 0;
4564 for(c = 0; c < ms->FDim[2]; c += (ms->FDim[2] - 1)) {
4565 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
4566 for(b = 0; b < ms->FDim[1]; b += (ms->FDim[1] - 1)) {
4567 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
4568 for(a = 0; a < ms->FDim[0]; a += (ms->FDim[0] - 1)) {
4569 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
4570 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
4571 copy3f(vr, ms->Corner + 3 * d);
4572 d++;
4573 }
4574 }
4575 }
4576
4577 v[2] = (ms->Min[2]) / ((float) ms->Div[2]);
4578 v[1] = (ms->Min[1]) / ((float) ms->Div[1]);
4579 v[0] = (ms->Min[0]) / ((float) ms->Div[0]);
4580
4581 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMin);
4582
4583 v[2] = ((ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
4584 v[1] = ((ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
4585 v[0] = ((ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
4586
4587 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMax);
4588
4589 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4590 " GRDXStrToMap: Map Size %d x %d x %d\n", ms->FDim[0], ms->FDim[1], ms->FDim[2]
4591 ENDFB(I->G);
4592
4593 ms->Active = true;
4594 ObjectMapUpdateExtents(I);
4595 if(!quiet) {
4596 PRINTFB(I->G, FB_ObjectMap, FB_Results)
4597 " ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(I->G);
4598 }
4599 } else {
4600 PRINTFB(I->G, FB_ObjectMap, FB_Errors)
4601 " Error: unable to read GRD file.\n" ENDFB(I->G);
4602 }
4603
4604 return (ok);
4605
4606 }
4607
4608
4609 /*========================================================================*/
ObjectMapReadXPLORStr(PyMOLGlobals * G,ObjectMap * I,char * XPLORStr,int state,int quiet)4610 static ObjectMap *ObjectMapReadXPLORStr(PyMOLGlobals * G, ObjectMap * I, char *XPLORStr,
4611 int state, int quiet)
4612 {
4613 int ok = true;
4614 int isNew = true;
4615
4616 if(!I)
4617 isNew = true;
4618 else
4619 isNew = false;
4620 if(ok) {
4621 if(isNew) {
4622 I = (ObjectMap *) new ObjectMap(G);
4623 isNew = true;
4624 } else {
4625 isNew = false;
4626 }
4627 ObjectMapXPLORStrToMap(I, XPLORStr, state, quiet);
4628
4629 SceneChanged(I->G);
4630 SceneCountFrames(I->G);
4631 }
4632 return (I);
4633 }
4634
4635
4636 /*========================================================================*/
ObjectMapReadCCP4Str(PyMOLGlobals * G,ObjectMap * I,char * XPLORStr,int bytes,int state,int quiet,int format)4637 static ObjectMap *ObjectMapReadCCP4Str(PyMOLGlobals * G, ObjectMap * I, char *XPLORStr,
4638 int bytes, int state, int quiet,
4639 int format)
4640 {
4641 int ok = true;
4642 int isNew = true;
4643
4644 if(!I)
4645 isNew = true;
4646 else
4647 isNew = false;
4648 if(ok) {
4649 if(isNew) {
4650 I = (ObjectMap *) new ObjectMap(G);
4651 isNew = true;
4652 } else {
4653 isNew = false;
4654 }
4655 ObjectMapCCP4StrToMap(I, XPLORStr, bytes, state, quiet, format);
4656 SceneChanged(G);
4657 SceneCountFrames(G);
4658 }
4659 return (I);
4660 }
4661
4662
4663 /*========================================================================*/
ObjectMapLoadCCP4(PyMOLGlobals * G,ObjectMap * obj,const char * fname,int state,int is_string,int bytes,int quiet,int format)4664 ObjectMap *ObjectMapLoadCCP4(PyMOLGlobals * G, ObjectMap * obj, const char *fname, int state,
4665 int is_string, int bytes, int quiet,
4666 int format)
4667 {
4668 ObjectMap *I = NULL;
4669 char *buffer;
4670 long size;
4671
4672 if(!is_string) {
4673 if (!quiet)
4674 PRINTFB(G, FB_ObjectMap, FB_Actions)
4675 " ObjectMapLoadCCP4File: Loading from '%s'.\n", fname ENDFB(G);
4676
4677 buffer = FileGetContents(fname, &size);
4678
4679 if(!buffer)
4680 ErrMessage(G, "ObjectMapLoadCCP4File", "Unable to open file!");
4681 } else {
4682 buffer = (char*) fname;
4683 size = (long) bytes;
4684 }
4685
4686 if (buffer) {
4687 I = ObjectMapReadCCP4Str(G, obj, buffer, size, state, quiet, format);
4688
4689 if(!is_string)
4690 mfree(buffer);
4691
4692 if(!quiet) {
4693 if(state < 0)
4694 state = I->State.size() - 1;
4695 if(state < I->State.size()) {
4696 ObjectMapState *ms;
4697 ms = &I->State[state];
4698 if(ms->Active) {
4699 CrystalDump(&ms->Symmetry->Crystal);
4700 }
4701 }
4702 }
4703 }
4704 return (I);
4705 }
4706
4707
4708 /*========================================================================*/
ObjectMapReadFLDStr(PyMOLGlobals * G,ObjectMap * I,char * MapStr,int bytes,int state,int quiet)4709 static ObjectMap *ObjectMapReadFLDStr(PyMOLGlobals * G, ObjectMap * I, char *MapStr,
4710 int bytes, int state, int quiet)
4711 {
4712 int ok = true;
4713 int isNew = true;
4714
4715 if(!I)
4716 isNew = true;
4717 else
4718 isNew = false;
4719 if(ok) {
4720 if(isNew) {
4721 I = (ObjectMap *) new ObjectMap(G);
4722 isNew = true;
4723 } else {
4724 isNew = false;
4725 }
4726 ObjectMapFLDStrToMap(I, MapStr, bytes, state, quiet);
4727 SceneChanged(G);
4728 SceneCountFrames(G);
4729 }
4730 return (I);
4731 }
4732
4733
4734 /*========================================================================*/
ObjectMapReadBRIXStr(PyMOLGlobals * G,ObjectMap * I,char * MapStr,int bytes,int state,int quiet)4735 static ObjectMap *ObjectMapReadBRIXStr(PyMOLGlobals * G, ObjectMap * I, char *MapStr,
4736 int bytes, int state, int quiet)
4737 {
4738 int ok = true;
4739 int isNew = true;
4740
4741 if(!I)
4742 isNew = true;
4743 else
4744 isNew = false;
4745 if(ok) {
4746 if(isNew) {
4747 I = (ObjectMap *) new ObjectMap(G);
4748 isNew = true;
4749 } else {
4750 isNew = false;
4751 }
4752 ObjectMapBRIXStrToMap(I, MapStr, bytes, state, quiet);
4753 SceneChanged(G);
4754 SceneCountFrames(G);
4755 }
4756 return (I);
4757 }
4758
4759
4760 /*========================================================================*/
ObjectMapReadGRDStr(PyMOLGlobals * G,ObjectMap * I,char * MapStr,int bytes,int state,int quiet)4761 static ObjectMap *ObjectMapReadGRDStr(PyMOLGlobals * G, ObjectMap * I, char *MapStr,
4762 int bytes, int state, int quiet)
4763 {
4764 int ok = true;
4765 int isNew = true;
4766
4767 if(!I)
4768 isNew = true;
4769 else
4770 isNew = false;
4771 if(ok) {
4772 if(isNew) {
4773 I = (ObjectMap *) new ObjectMap(G);
4774 isNew = true;
4775 } else {
4776 isNew = false;
4777 }
4778 ObjectMapGRDStrToMap(I, MapStr, bytes, state, quiet);
4779 SceneChanged(G);
4780 SceneCountFrames(G);
4781 }
4782 return (I);
4783 }
4784
4785
4786 /*========================================================================*/
ObjectMapReadPHIStr(PyMOLGlobals * G,ObjectMap * I,char * MapStr,int bytes,int state,int quiet)4787 static ObjectMap *ObjectMapReadPHIStr(PyMOLGlobals * G, ObjectMap * I, char *MapStr,
4788 int bytes, int state, int quiet)
4789 {
4790 int ok = true;
4791 int isNew = true;
4792
4793 if(!I)
4794 isNew = true;
4795 else
4796 isNew = false;
4797 if(ok) {
4798 if(isNew) {
4799 I = (ObjectMap *) new ObjectMap(G);
4800 isNew = true;
4801 } else {
4802 isNew = false;
4803 }
4804 ObjectMapPHIStrToMap(I, MapStr, bytes, state, quiet);
4805 SceneChanged(G);
4806 SceneCountFrames(G);
4807 }
4808 return (I);
4809 }
4810
4811
4812 /*========================================================================*/
ObjectMapLoadPHI(PyMOLGlobals * G,ObjectMap * obj,const char * fname,int state,int is_string,int bytes,int quiet)4813 ObjectMap *ObjectMapLoadPHI(PyMOLGlobals * G, ObjectMap * obj, const char *fname, int state,
4814 int is_string, int bytes, int quiet)
4815 {
4816
4817 ObjectMap *I = NULL;
4818 long size;
4819 char *buffer;
4820
4821 if(!is_string) {
4822 if (!quiet)
4823 PRINTFB(G, FB_ObjectMap, FB_Actions)
4824 " ObjectMapLoadPHIFile: Loading from '%s'.\n", fname ENDFB(G);
4825
4826 buffer = FileGetContents(fname, &size);
4827
4828 if(!buffer)
4829 ErrMessage(G, "ObjectMapLoadPHIFile", "Unable to open file!");
4830 } else {
4831 buffer = (char*) fname;
4832 size = (long) bytes;
4833 }
4834
4835 if(buffer) {
4836 I = ObjectMapReadPHIStr(G, obj, buffer, size, state, quiet);
4837
4838 if(!is_string)
4839 mfree(buffer);
4840
4841 }
4842 return (I);
4843
4844 }
4845
4846
4847 /*========================================================================*/
4848
is_number(char * p)4849 static int is_number(char *p)
4850 {
4851 int result = (*p != 0);
4852 char c;
4853 if(result)
4854 while((c = *(p++))) {
4855 if(!((c == '.') ||
4856 (c == '-') ||
4857 (c == '+') || (c == 'e') || (c == 'E') || ((c >= '0') && (c <= '9'))))
4858 result = false;
4859 break;
4860 }
4861 return result;
4862 }
4863
ObjectMapDXStrToMap(ObjectMap * I,char * DXStr,int bytes,int state,int quiet)4864 static int ObjectMapDXStrToMap(ObjectMap * I, char *DXStr, int bytes, int state,
4865 int quiet)
4866 {
4867
4868 int n_items = 0;
4869
4870 char *p, *pp;
4871 float dens;
4872 int a, b, c, d, e;
4873 float v[3], maxd, mind;
4874 int ok = true;
4875 /* DX named from their docs */
4876
4877 int stage = 0;
4878
4879 ObjectMapState *ms;
4880
4881 char cc[MAXLINELEN];
4882
4883 if(state < 0)
4884 state = I->State.size();
4885 if(I->State.size() <= state) {
4886 VecCheckEmplace(I->State, state, I->G);
4887 }
4888 ms = &I->State[state];
4889
4890 ms->Origin = std::vector<float>(3);
4891 ms->Grid = std::vector<float>(3);
4892
4893 maxd = -FLT_MAX;
4894 mind = FLT_MAX;
4895 p = DXStr;
4896
4897 /* get the dimensions */
4898
4899 ms->FDim[3] = 3;
4900
4901 while(ok && (*p) && (stage == 0)) {
4902 pp = p;
4903 p = ParseNCopy(cc, p, 35);
4904 if((strcmp(cc, "object 1 class gridpositions counts") == 0) || is_number(cc)) {
4905 if(is_number(cc))
4906 p = pp;
4907 p = ParseWordCopy(cc, p, 10);
4908 if(sscanf(cc, "%d", &ms->FDim[0]) == 1) {
4909 p = ParseWordCopy(cc, p, 10);
4910 if(sscanf(cc, "%d", &ms->FDim[1]) == 1) {
4911 p = ParseWordCopy(cc, p, 10);
4912 if(sscanf(cc, "%d", &ms->FDim[2]) == 1) {
4913 stage = 1;
4914 }
4915 }
4916 }
4917 }
4918 p = ParseNextLine(p);
4919 }
4920
4921 if(ok && (stage == 1)) {
4922 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4923 " DXStrToMap: Dimensions: %d %d %d\n", ms->FDim[0], ms->FDim[1], ms->FDim[2]
4924 ENDFB(I->G);
4925 }
4926
4927 /* get the origin */
4928
4929 while(ok && (*p) && (stage == 1)) {
4930 pp = p;
4931 p = ParseNCopy(cc, p, 6);
4932 if((strcmp(cc, "origin") == 0) || is_number(cc)) {
4933 if(is_number(cc))
4934 p = pp;
4935 p = ParseWordCopy(cc, p, 20);
4936 if(sscanf(cc, "%f", &ms->Origin[0]) == 1) {
4937 p = ParseWordCopy(cc, p, 20);
4938 if(sscanf(cc, "%f", &ms->Origin[1]) == 1) {
4939 p = ParseWordCopy(cc, p, 20);
4940 if(sscanf(cc, "%f", &ms->Origin[2]) == 1) {
4941 stage = 2;
4942 }
4943 }
4944 }
4945 }
4946 p = ParseNextLine(p);
4947 }
4948
4949 if(ok && (stage == 2)) {
4950 PRINTFB(I->G, FB_ObjectMap, FB_Details)
4951 " DXStrToMap: Origin %8.3f %8.3f %8.3f\n", ms->Origin[0], ms->Origin[1],
4952 ms->Origin[2]
4953 ENDFB(I->G);
4954 }
4955
4956 float delta[9];
4957 int delta_i = 0;
4958
4959 while(ok && (*p) && (stage == 2)) {
4960 pp = p;
4961 p = ParseNCopy(cc, p, 5);
4962
4963 if(strcmp(cc, "delta") != 0) {
4964 if(is_number(cc)) {
4965 p = pp;
4966 } else {
4967 p = ParseNextLine(p);
4968 continue;
4969 }
4970 }
4971
4972 if(3 != sscanf(p, " %f %f %f",
4973 delta + delta_i,
4974 delta + delta_i + 1,
4975 delta + delta_i + 2)) {
4976 // error
4977 break;
4978 }
4979
4980 p = ParseNextLine(p);
4981 delta_i += 3;
4982
4983 if (delta_i == 9) {
4984 stage = 3;
4985
4986 if (is_diagonalf(3, delta)) {
4987 ms->Grid[0] = delta[0];
4988 ms->Grid[1] = delta[4];
4989 ms->Grid[2] = delta[8];
4990 } else {
4991 if(ms->Matrix.empty())
4992 ms->Matrix = std::vector<double>(16);
4993
4994 copy33f44d(delta, ms->Matrix.data());
4995 ms->Matrix[3] = ms->Origin[0];
4996 ms->Matrix[7] = ms->Origin[1];
4997 ms->Matrix[11] = ms->Origin[2];
4998
4999 ones3f(ms->Grid.data());
5000 zero3f(ms->Origin.data());
5001 }
5002 }
5003 }
5004
5005 if(ok && (stage == 3)) {
5006 PRINTFB(I->G, FB_ObjectMap, FB_Details)
5007 " DXStrToMap: Grid %8.3f %8.3f %8.3f\n", ms->Grid[0], ms->Grid[1], ms->Grid[2]
5008 ENDFB(I->G);
5009 }
5010
5011 while(ok && (*p) && (stage == 3)) {
5012 p = ParseNCopy(cc, p, 6);
5013 if(strcmp(cc, "object") == 0) {
5014 p = ParseWordCopy(cc, p, 20);
5015 if (1 == sscanf(p, " class array type %*s rank %*s items %s", cc)) {
5016 if(sscanf(cc, "%d", &n_items) == 1) {
5017 if(n_items == ms->FDim[0] * ms->FDim[1] * ms->FDim[2])
5018 stage = 4;
5019 }
5020 }
5021 } else if(is_number(cc)) {
5022 n_items = ms->FDim[0] * ms->FDim[1] * ms->FDim[2];
5023 stage = 4;
5024 break;
5025 }
5026 p = ParseNextLine(p);
5027 }
5028
5029 if(stage == 4) {
5030
5031 if(ok && (stage == 4)) {
5032 PRINTFB(I->G, FB_ObjectMap, FB_Details)
5033 " DXStrToMap: %d data points.\n", n_items ENDFB(I->G);
5034 }
5035
5036 ms->Field.reset(new Isofield(I->G, ms->FDim));
5037 ms->MapSource = cMapSourceGeneralPurpose;
5038 ms->Field->save_points = false;
5039
5040 for(a = 0; a < 3; a++) {
5041 ms->Div[a] = ms->FDim[a] - 1;
5042 ms->Min[a] = 0;
5043 ms->Max[a] = ms->FDim[a] - 1;
5044 }
5045
5046 for(a = 0; a < ms->FDim[0]; a++) {
5047 for(b = 0; b < ms->FDim[1]; b++) {
5048 for(c = 0; c < ms->FDim[2]; c++) {
5049
5050 p = ParseWordCopy(cc, p, 20);
5051 if(!cc[0]) {
5052 p = ParseNextLine(p);
5053 p = ParseWordCopy(cc, p, 20);
5054 }
5055 if(sscanf(cc, "%f", &dens) == 1) {
5056 if(maxd < dens)
5057 maxd = dens;
5058 if(mind > dens)
5059 mind = dens;
5060 F3(ms->Field->data, a, b, c) = dens;
5061 } else {
5062 ok = false;
5063 }
5064 }
5065 }
5066
5067 }
5068
5069 for(e = 0; e < 3; e++) {
5070 ms->ExtentMin[e] = ms->Origin[e] + ms->Grid[e] * ms->Min[e];
5071 ms->ExtentMax[e] = ms->Origin[e] + ms->Grid[e] * ms->Max[e];
5072 }
5073
5074 for(c = 0; c < ms->FDim[2]; c++) {
5075 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
5076 for(b = 0; b < ms->FDim[1]; b++) {
5077 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
5078 for(a = 0; a < ms->FDim[0]; a++) {
5079 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
5080 for(e = 0; e < 3; e++) {
5081 F4(ms->Field->points, a, b, c, e) = v[e];
5082 }
5083 }
5084 }
5085 }
5086
5087 d = 0;
5088 for(c = 0; c < ms->FDim[2]; c += ms->FDim[2] - 1) {
5089 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
5090
5091 for(b = 0; b < ms->FDim[1]; b += ms->FDim[1] - 1) {
5092 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
5093
5094 for(a = 0; a < ms->FDim[0]; a += ms->FDim[0] - 1) {
5095 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
5096 copy3f(v, ms->Corner + 3 * d);
5097 d++;
5098 }
5099 }
5100 }
5101 if(ok)
5102 stage = 5;
5103 }
5104
5105 if(stage != 5)
5106 ok = false;
5107
5108 if(!ok) {
5109 ErrMessage(I->G, "ObjectMap", "Error reading map");
5110 } else {
5111 ms->Active = true;
5112 ObjectMapUpdateExtents(I);
5113 if(!quiet) {
5114 PRINTFB(I->G, FB_ObjectMap, FB_Results)
5115 " ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(I->G);
5116 }
5117 }
5118 return (ok);
5119 }
5120
5121
5122 /*========================================================================*/
ObjectMapReadDXStr(PyMOLGlobals * G,ObjectMap * I,char * MapStr,int bytes,int state,int quiet)5123 static ObjectMap *ObjectMapReadDXStr(PyMOLGlobals * G, ObjectMap * I,
5124 char *MapStr, int bytes, int state, int quiet)
5125 {
5126 int ok = true;
5127 int isNew = true;
5128
5129 if(!I)
5130 isNew = true;
5131 else
5132 isNew = false;
5133 if(ok) {
5134 if(isNew) {
5135 I = (ObjectMap *) new ObjectMap(G);
5136 isNew = true;
5137 } else {
5138 isNew = false;
5139 }
5140 ObjectMapDXStrToMap(I, MapStr, bytes, state, quiet);
5141 SceneChanged(G);
5142 SceneCountFrames(G);
5143 }
5144 return (I);
5145 }
5146
5147
5148 /*========================================================================*/
ObjectMapLoadDXFile(PyMOLGlobals * G,ObjectMap * obj,const char * fname,int state,int quiet)5149 ObjectMap *ObjectMapLoadDXFile(PyMOLGlobals * G, ObjectMap * obj, const char *fname, int state,
5150 int quiet)
5151 {
5152 ObjectMap *I = NULL;
5153 long size;
5154 char *buffer;
5155 float mat[9];
5156
5157 buffer = FileGetContents(fname, &size);
5158
5159 if(!buffer) {
5160 ErrMessage(G, "ObjectMapLoadDXFile", "Unable to open file!");
5161 PRINTFB(G, FB_ObjectMap, FB_Errors)
5162 "ObjectMapLoadDXFile: Does '%s' exist?\n", fname ENDFB(G);
5163 } else {
5164 if(Feedback(G, FB_ObjectMap, FB_Actions)) {
5165 printf(" ObjectMapLoadDXFile: Loading from '%s'.\n", fname);
5166 }
5167
5168 I = ObjectMapReadDXStr(G, obj, buffer, size, state, quiet);
5169
5170 mfree(buffer);
5171 if(state < 0)
5172 state = I->State.size() - 1;
5173 if(state < I->State.size()) {
5174 ObjectMapState *ms;
5175 ms = &I->State[state];
5176 if(ms->Active) {
5177 multiply33f33f(ms->Symmetry->Crystal.FracToReal, ms->Symmetry->Crystal.RealToFrac, mat);
5178 }
5179 }
5180 }
5181 return (I);
5182
5183 }
5184
5185
5186 /*========================================================================*/
ObjectMapACNTStrToMap(ObjectMap * I,char * ACNTStr,int bytes,int state,int quiet)5187 static int ObjectMapACNTStrToMap(ObjectMap * I, char *ACNTStr, int bytes, int state,
5188 int quiet)
5189 {
5190
5191 int n_items = 0;
5192
5193 char *p;
5194 float dens;
5195 int a, b, c, d, e;
5196 float v[3], maxd, mind;
5197 int ok = true;
5198
5199 int stage = 0;
5200
5201 ObjectMapState *ms;
5202
5203 char cc[MAXLINELEN];
5204
5205 if(state < 0)
5206 state = I->State.size();
5207 if(I->State.size() <= state) {
5208 VecCheckEmplace(I->State, state, I->G);
5209 }
5210
5211 ms = &I->State[state];
5212
5213 ms->Origin = std::vector<float>(3);
5214 ms->Grid = std::vector<float>(3);
5215
5216 maxd = -FLT_MAX;
5217 mind = FLT_MAX;
5218 p = ACNTStr;
5219
5220 /* skip header */
5221
5222 p = ParseNextLine(p);
5223
5224 /* get the origin, spacing, and dimensions */
5225
5226 ms->FDim[3] = 3;
5227
5228 /* read the map info */
5229
5230 {
5231 int i;
5232 for(i = 0; i < 3; i++) {
5233 int ii = (4 - i) % 3; /* y, x , z */
5234 p = ParseWordCopy(cc, p, MAXLINELEN);
5235 if(sscanf(cc, "%f", &ms->Origin[ii]) == 1) {
5236 p = ParseWordCopy(cc, p, MAXLINELEN);
5237 if(sscanf(cc, "%f", &ms->Grid[ii]) == 1) {
5238 p = ParseWordCopy(cc, p, MAXLINELEN);
5239 if(sscanf(cc, "%d", &ms->FDim[ii]) == 1) {
5240 p = ParseNextLine(p);
5241 stage++;
5242 }
5243 }
5244 }
5245 }
5246 }
5247
5248 /* skip the angles (ignored -- should probably check instead */
5249
5250 p = ParseNextLine(p);
5251
5252 /* read the data block */
5253
5254 if(ok && (stage == 3)) {
5255
5256 PRINTFB(I->G, FB_ObjectMap, FB_Details)
5257 " ACNTStrToMap: Dimensions: %d %d %d\n", ms->FDim[0], ms->FDim[1], ms->FDim[2]
5258 ENDFB(I->G);
5259 PRINTFB(I->G, FB_ObjectMap, FB_Details)
5260 " ACNTStrToMap: Origin %8.3f %8.3f %8.3f\n", ms->Origin[0], ms->Origin[1],
5261 ms->Origin[2]
5262 ENDFB(I->G);
5263 PRINTFB(I->G, FB_ObjectMap, FB_Details)
5264 " ACNTStrToMap: Grid %8.3f %8.3f %8.3f\n", ms->Grid[0], ms->Grid[1], ms->Grid[2]
5265 ENDFB(I->G);
5266
5267 n_items = ms->FDim[0] * ms->FDim[1] * ms->FDim[2];
5268
5269 if(ok && (stage == 1)) {
5270 PRINTFB(I->G, FB_ObjectMap, FB_Details)
5271 " ACNTStrToMap: %d data points.\n", n_items ENDFB(I->G);
5272 }
5273
5274 ms->Field.reset(new Isofield(I->G, ms->FDim));
5275 ms->MapSource = cMapSourceGeneralPurpose;
5276 ms->Field->save_points = false;
5277
5278 for(a = 0; a < 3; a++) {
5279 ms->Div[a] = ms->FDim[a] - 1;
5280 ms->Min[a] = 0;
5281 ms->Max[a] = ms->FDim[a] - 1;
5282 }
5283
5284 for(c = 0; c < ms->FDim[2]; c++) {
5285 for(a = 0; a < ms->FDim[0]; a++) {
5286 for(b = 0; b < ms->FDim[1]; b++) {
5287 p = ParseWordCopy(cc, p, MAXLINELEN);
5288 p = ParseNextLine(p);
5289 if(sscanf(cc, "%f", &dens) == 1) {
5290 if(maxd < dens)
5291 maxd = dens;
5292 if(mind > dens)
5293 mind = dens;
5294 F3(ms->Field->data, a, b, c) = dens;
5295 } else {
5296 ok = false;
5297 }
5298 }
5299 }
5300 }
5301
5302 for(e = 0; e < 3; e++) {
5303 ms->ExtentMin[e] = ms->Origin[e] + ms->Grid[e] * ms->Min[e];
5304 ms->ExtentMax[e] = ms->Origin[e] + ms->Grid[e] * ms->Max[e];
5305 }
5306
5307 for(c = 0; c < ms->FDim[2]; c++) {
5308 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
5309 for(b = 0; b < ms->FDim[1]; b++) {
5310 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
5311 for(a = 0; a < ms->FDim[0]; a++) {
5312 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
5313 for(e = 0; e < 3; e++) {
5314 F4(ms->Field->points, a, b, c, e) = v[e];
5315 }
5316 }
5317 }
5318 }
5319
5320 d = 0;
5321 for(c = 0; c < ms->FDim[2]; c += ms->FDim[2] - 1) {
5322 v[2] = ms->Origin[2] + ms->Grid[2] * (c + ms->Min[2]);
5323
5324 for(b = 0; b < ms->FDim[1]; b += ms->FDim[1] - 1) {
5325 v[1] = ms->Origin[1] + ms->Grid[1] * (b + ms->Min[1]);
5326
5327 for(a = 0; a < ms->FDim[0]; a += ms->FDim[0] - 1) {
5328 v[0] = ms->Origin[0] + ms->Grid[0] * (a + ms->Min[0]);
5329 copy3f(v, ms->Corner + 3 * d);
5330 d++;
5331 }
5332 }
5333 }
5334
5335 if(ok)
5336 stage = 4;
5337 }
5338
5339 if(stage != 4)
5340 ok = false;
5341
5342 if(!ok) {
5343 ErrMessage(I->G, "ObjectMap", "Error reading map");
5344 } else {
5345 ms->Active = true;
5346 ObjectMapUpdateExtents(I);
5347 if(!quiet) {
5348 PRINTFB(I->G, FB_ObjectMap, FB_Results)
5349 " ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(I->G);
5350 }
5351 }
5352 return (ok);
5353 }
5354
ObjectMapReadACNTStr(PyMOLGlobals * G,ObjectMap * I,char * MapStr,int bytes,int state,int quiet)5355 static ObjectMap *ObjectMapReadACNTStr(PyMOLGlobals * G, ObjectMap * I,
5356 char *MapStr, int bytes, int state, int quiet)
5357 {
5358 int ok = true;
5359 int isNew = true;
5360
5361 if(!I)
5362 isNew = true;
5363 else
5364 isNew = false;
5365 if(ok) {
5366 if(isNew) {
5367 I = (ObjectMap *) new ObjectMap(G);
5368 isNew = true;
5369 } else {
5370 isNew = false;
5371 }
5372 ObjectMapACNTStrToMap(I, MapStr, bytes, state, quiet);
5373 SceneChanged(G);
5374 SceneCountFrames(G);
5375 }
5376 return (I);
5377 }
5378
ObjectMapLoadACNTFile(PyMOLGlobals * G,ObjectMap * obj,const char * fname,int state,int quiet)5379 ObjectMap *ObjectMapLoadACNTFile(PyMOLGlobals * G, ObjectMap * obj, const char *fname,
5380 int state, int quiet)
5381 {
5382 ObjectMap *I = NULL;
5383 long size;
5384 char *buffer;
5385 float mat[9];
5386
5387 buffer = FileGetContents(fname, &size);
5388
5389 if(!buffer) {
5390 ErrMessage(G, "ObjectMapLoadACNTFile", "Unable to open file!");
5391 PRINTFB(G, FB_ObjectMap, FB_Errors)
5392 "ObjectMapLoadACNTFile: Does '%s' exist?\n", fname ENDFB(G);
5393 } else {
5394 if(Feedback(G, FB_ObjectMap, FB_Actions)) {
5395 printf(" ObjectMapLoadACNTFile: Loading from '%s'.\n", fname);
5396 }
5397
5398 I = ObjectMapReadACNTStr(G, obj, buffer, size, state, quiet);
5399
5400 mfree(buffer);
5401 if(state < 0)
5402 state = I->State.size() - 1;
5403 if(state < I->State.size()) {
5404 ObjectMapState *ms;
5405 ms = &I->State[state];
5406 if(ms->Active) {
5407 multiply33f33f(ms->Symmetry->Crystal.FracToReal, ms->Symmetry->Crystal.RealToFrac, mat);
5408 }
5409 }
5410 }
5411 return (I);
5412
5413 }
5414
5415
5416 /*========================================================================*/
ObjectMapLoadFLDFile(PyMOLGlobals * G,ObjectMap * obj,const char * fname,int state,int quiet)5417 ObjectMap *ObjectMapLoadFLDFile(PyMOLGlobals * G, ObjectMap * obj, const char *fname, int state,
5418 int quiet)
5419 {
5420 ObjectMap *I = NULL;
5421 long size;
5422 char *buffer;
5423 float mat[9];
5424
5425 buffer = FileGetContents(fname, &size);
5426
5427 if(!buffer) {
5428 ErrMessage(G, "ObjectMapLoadFLDFile", "Unable to open file!");
5429 } else {
5430 if(Feedback(G, FB_ObjectMap, FB_Actions)) {
5431 printf(" ObjectMapLoadFLDFile: Loading from '%s'.\n", fname);
5432 }
5433
5434 I = ObjectMapReadFLDStr(G, obj, buffer, size, state, quiet);
5435
5436 mfree(buffer);
5437 if(state < 0)
5438 state = I->State.size() - 1;
5439 if(state < I->State.size()) {
5440 ObjectMapState *ms;
5441 ms = &I->State[state];
5442 if(ms->Active) {
5443 multiply33f33f(ms->Symmetry->Crystal.FracToReal, ms->Symmetry->Crystal.RealToFrac, mat);
5444 }
5445 }
5446 }
5447 return (I);
5448
5449 }
5450
5451
5452 /*========================================================================*/
ObjectMapLoadBRIXFile(PyMOLGlobals * G,ObjectMap * obj,const char * fname,int state,int quiet)5453 ObjectMap *ObjectMapLoadBRIXFile(PyMOLGlobals * G, ObjectMap * obj, const char *fname,
5454 int state, int quiet)
5455 {
5456 ObjectMap *I = NULL;
5457 long size;
5458 char *buffer;
5459 float mat[9];
5460
5461 buffer = FileGetContents(fname, &size);
5462
5463 if(!buffer) {
5464 ErrMessage(G, "ObjectMapLoadBRIXFile", "Unable to open file!");
5465 } else {
5466 if(Feedback(G, FB_ObjectMap, FB_Actions)) {
5467 printf(" ObjectMapLoadBRIXFile: Loading from '%s'.\n", fname);
5468 }
5469
5470 I = ObjectMapReadBRIXStr(G, obj, buffer, size, state, quiet);
5471
5472 mfree(buffer);
5473 if(state < 0)
5474 state = I->State.size() - 1;
5475 if(state < I->State.size()) {
5476 ObjectMapState *ms;
5477 ms = &I->State[state];
5478 if(ms->Active) {
5479 CrystalDump(&ms->Symmetry->Crystal);
5480 multiply33f33f(ms->Symmetry->Crystal.FracToReal, ms->Symmetry->Crystal.RealToFrac, mat);
5481 }
5482 }
5483 }
5484 return (I);
5485
5486 }
5487
5488
5489 /*========================================================================*/
ObjectMapLoadGRDFile(PyMOLGlobals * G,ObjectMap * obj,const char * fname,int state,int quiet)5490 ObjectMap *ObjectMapLoadGRDFile(PyMOLGlobals * G, ObjectMap * obj, const char *fname, int state,
5491 int quiet)
5492 {
5493 ObjectMap *I = NULL;
5494 long size;
5495 char *buffer;
5496 float mat[9];
5497
5498 buffer = FileGetContents(fname, &size);
5499
5500 if(!buffer) {
5501 ErrMessage(G, "ObjectMapLoadGRDFile", "Unable to open file!");
5502 } else {
5503 if(Feedback(G, FB_ObjectMap, FB_Actions)) {
5504 printf(" ObjectMapLoadGRDFile: Loading from '%s'.\n", fname);
5505 }
5506
5507 I = ObjectMapReadGRDStr(G, obj, buffer, size, state, quiet);
5508
5509 mfree(buffer);
5510 if(state < 0)
5511 state = I->State.size() - 1;
5512 if(state < I->State.size()) {
5513 ObjectMapState *ms;
5514 ms = &I->State[state];
5515 if(ms->Active) {
5516 CrystalDump(&ms->Symmetry->Crystal);
5517 multiply33f33f(ms->Symmetry->Crystal.FracToReal, ms->Symmetry->Crystal.RealToFrac, mat);
5518 }
5519 }
5520 }
5521 return (I);
5522
5523 }
5524
5525
5526 /*========================================================================*/
ObjectMapLoadXPLOR(PyMOLGlobals * G,ObjectMap * obj,const char * fname,int state,int is_file,int quiet)5527 ObjectMap *ObjectMapLoadXPLOR(PyMOLGlobals * G, ObjectMap * obj, const char *fname,
5528 int state, int is_file, int quiet)
5529 {
5530 ObjectMap *I = NULL;
5531 long size;
5532 char *buffer;
5533
5534 if(is_file) {
5535 buffer = FileGetContents(fname, &size);
5536
5537 if(!buffer)
5538 ErrMessage(G, "ObjectMapLoadXPLOR", "Unable to open file!");
5539 } else {
5540 buffer = (char*) fname;
5541 }
5542
5543 if (buffer) {
5544 if((!quiet) && (Feedback(G, FB_ObjectMap, FB_Actions))) {
5545 if(is_file) {
5546 printf(" ObjectMapLoadXPLOR: Loading from '%s'.\n", fname);
5547 } else {
5548 printf(" ObjectMapLoadXPLOR: Loading...\n");
5549 }
5550 }
5551
5552 I = ObjectMapReadXPLORStr(G, obj, buffer, state, quiet);
5553
5554 if(is_file)
5555 mfree(buffer);
5556 if(!quiet) {
5557 if(Feedback(G, FB_ObjectMap, FB_Details)) {
5558 if(state < 0)
5559 state = I->State.size() - 1;
5560
5561 if(state < I->State.size()) {
5562 ObjectMapState *ms;
5563 ms = &I->State[state];
5564 if(ms->Active) {
5565 CrystalDump(&ms->Symmetry->Crystal);
5566 }
5567 }
5568 }
5569 }
5570 }
5571 return (I);
5572
5573 }
5574
5575
5576 /*========================================================================*/
ObjectMapSetBorder(ObjectMap * I,float level,int state)5577 int ObjectMapSetBorder(ObjectMap * I, float level, int state)
5578 {
5579 for (StateIterator iter(I, state); iter.next();) {
5580 auto* oms = &I->State[iter.state];
5581 if (oms->Active) {
5582 if (!ObjectMapStateSetBorder(oms, level))
5583 return false;
5584 }
5585 }
5586 return true;
5587 }
5588
5589
5590 /*========================================================================*/
5591 #ifndef _PYMOL_NOPY
ObjectMapNumPyArrayToMapState(PyMOLGlobals * G,ObjectMapState * ms,PyObject * ary,int quiet)5592 static int ObjectMapNumPyArrayToMapState(PyMOLGlobals * G, ObjectMapState * ms,
5593 PyObject * ary, int quiet)
5594 {
5595
5596 int a, b, c, d, e;
5597 float v[3], dens, maxd, mind;
5598 int ok = true;
5599 void * ptr;
5600
5601 #ifdef _PYMOL_NUMPY
5602 PyArrayObject * pao = (PyArrayObject *) ary;
5603 const int itemsize = PyArray_ITEMSIZE(pao);
5604 #endif
5605 maxd = -FLT_MAX;
5606 mind = FLT_MAX;
5607 if(ok) {
5608 ms->FDim[0] = ms->Dim[0];
5609 ms->FDim[1] = ms->Dim[1];
5610 ms->FDim[2] = ms->Dim[2];
5611 ms->FDim[3] = 3;
5612
5613 if(!(ms->FDim[0] && ms->FDim[1] && ms->FDim[2]))
5614 ok = false;
5615 else {
5616 ms->Field.reset(new Isofield(G, ms->FDim));
5617 for(c = 0; c < ms->FDim[2]; c++) {
5618 v[2] = ms->Origin[2] + ms->Grid[2] * c;
5619 for(b = 0; b < ms->FDim[1]; b++) {
5620 v[1] = ms->Origin[1] + ms->Grid[1] * b;
5621 for(a = 0; a < ms->FDim[0]; a++) {
5622 v[0] = ms->Origin[0] + ms->Grid[0] * a;
5623 #ifdef _PYMOL_NUMPY
5624 ptr = PyArray_GETPTR3(pao, a, b, c);
5625 switch(itemsize) {
5626 case sizeof(double):
5627 dens = (float) *((double*)ptr);
5628 break;
5629 case sizeof(float):
5630 dens = *((float*)ptr);
5631 break;
5632 default:
5633 dens = 0.0;
5634 printf("no itemsize match\n");
5635 }
5636 #else
5637 dens = 0.0;
5638 #endif
5639 F3(ms->Field->data, a, b, c) = dens;
5640 if(maxd < dens)
5641 maxd = dens;
5642 if(mind > dens)
5643 mind = dens;
5644 for(e = 0; e < 3; e++)
5645 F4(ms->Field->points, a, b, c, e) = v[e];
5646 }
5647 }
5648 }
5649 d = 0;
5650 for(c = 0; c < ms->FDim[2]; c += (ms->FDim[2] - 1)) {
5651 v[2] = ms->Origin[2] + ms->Grid[2] * c;
5652 for(b = 0; b < ms->FDim[1]; b += (ms->FDim[1] - 1)) {
5653 v[1] = ms->Origin[1] + ms->Grid[1] * b;
5654 for(a = 0; a < ms->FDim[0]; a += (ms->FDim[0] - 1)) {
5655 v[0] = ms->Origin[0] + ms->Grid[0] * a;
5656 copy3f(v, ms->Corner + 3 * d);
5657 d++;
5658 }
5659 }
5660 }
5661 }
5662 }
5663 if(ok) {
5664 copy3f(ms->Origin.data(), ms->ExtentMin);
5665 copy3f(ms->Origin.data(), ms->ExtentMax);
5666 add3f(ms->Range.data(), ms->ExtentMax, ms->ExtentMax);
5667 }
5668 if(!ok) {
5669 ErrMessage(G, "ObjectMap", "Error reading map");
5670 } else {
5671 ms->Active = true;
5672 if(!quiet) {
5673 PRINTFB(G, FB_ObjectMap, FB_Results)
5674 " ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(G);
5675 }
5676 }
5677 return (ok);
5678 }
5679 #endif
5680
5681 /*========================================================================*/
ObjectMapLoadChemPyBrick(PyMOLGlobals * G,ObjectMap * I,PyObject * Map,int state,int discrete,int quiet)5682 ObjectMap *ObjectMapLoadChemPyBrick(PyMOLGlobals * G, ObjectMap * I, PyObject * Map,
5683 int state, int discrete, int quiet)
5684 {
5685 #ifndef _PYMOL_NUMPY
5686 ErrMessage(G, "ObjectMap", "missing numpy support, required for chempy.brick");
5687 return NULL;
5688 #endif
5689
5690 #ifdef _PYMOL_NOPY
5691 return NULL;
5692 #else
5693
5694 int ok = true;
5695 int isNew = true;
5696 PyObject *tmp;
5697 ObjectMapState *ms;
5698
5699 if(!I)
5700 isNew = true;
5701 else
5702 isNew = false;
5703
5704 if(ok) {
5705
5706 if(isNew) {
5707 I = (ObjectMap *) new ObjectMap(G);
5708 isNew = true;
5709 } else {
5710 isNew = false;
5711 }
5712
5713 if(state < 0)
5714 state = I->State.size();
5715 if(I->State.size() <= state) {
5716 VecCheckEmplace(I->State, state, I->G);
5717 }
5718 ms = &I->State[state];
5719
5720 if(PyObject_HasAttrString(Map, "origin") &&
5721 PyObject_HasAttrString(Map, "dim") &&
5722 PyObject_HasAttrString(Map, "range") &&
5723 PyObject_HasAttrString(Map, "grid") && PyObject_HasAttrString(Map, "lvl")) {
5724 tmp = PyObject_GetAttrString(Map, "origin");
5725 if(tmp) {
5726 PConvFromPyObject(G, tmp, ms->Origin);
5727 Py_DECREF(tmp);
5728 } else
5729 ok = ErrMessage(G, "ObjectMap", "missing brick origin.");
5730 tmp = PyObject_GetAttrString(Map, "dim");
5731 if(tmp) {
5732 PConvFromPyObject(G, tmp, ms->Dim);
5733 Py_DECREF(tmp);
5734 } else
5735 ok = ErrMessage(G, "ObjectMap", "missing brick dimension.");
5736 tmp = PyObject_GetAttrString(Map, "range");
5737 if(tmp) {
5738 PConvFromPyObject(G, tmp, ms->Range);
5739 Py_DECREF(tmp);
5740 } else
5741 ok = ErrMessage(G, "ObjectMap", "missing brick range.");
5742 tmp = PyObject_GetAttrString(Map, "grid");
5743 if(tmp) {
5744 PConvFromPyObject(G, tmp, ms->Grid);
5745 Py_DECREF(tmp);
5746 } else
5747 ok = ErrMessage(G, "ObjectMap", "missing brick grid.");
5748 tmp = PyObject_GetAttrString(Map, "lvl");
5749 if(tmp) {
5750
5751 ObjectMapNumPyArrayToMapState(G, ms, tmp, quiet);
5752
5753 Py_DECREF(tmp);
5754 } else
5755 ok = ErrMessage(G, "ObjectMap", "missing brick density.");
5756
5757 } else
5758 ok = ErrMessage(G, "ObjectMap", "missing any brick attribute.");
5759
5760 SceneChanged(G);
5761 SceneCountFrames(G);
5762 if(ok) {
5763 int a;
5764 for(a = 0; a < 3; a++) {
5765 ms->Min[a] = 0;
5766 ms->Max[a] = ms->Dim[a] - 1;
5767 }
5768 ms->Active = true;
5769 ms->MapSource = cMapSourceChempyBrick;
5770 ObjectMapUpdateExtents(I);
5771
5772 }
5773 }
5774 return (I);
5775 #endif
5776 }
5777
5778
5779 /*========================================================================*/
ObjectMapLoadChemPyMap(PyMOLGlobals * G,ObjectMap * I,PyObject * Map,int state,int discrete,int quiet)5780 ObjectMap *ObjectMapLoadChemPyMap(PyMOLGlobals * G, ObjectMap * I, PyObject * Map,
5781 int state, int discrete, int quiet)
5782 {
5783 #ifdef _PYMOL_NOPY
5784 return NULL;
5785 #else
5786 int ok = true;
5787 int isNew = true;
5788 float *cobj;
5789 WordType format;
5790 float v[3], vr[3], dens, maxd, mind;
5791 int a, b, c, d, e;
5792 ObjectMapState *ms;
5793
5794 maxd = -FLT_MAX;
5795 mind = FLT_MAX;
5796
5797 if(!I)
5798 isNew = true;
5799 else
5800 isNew = false;
5801
5802 if(ok) {
5803
5804 if(isNew) {
5805 I = (ObjectMap *) new ObjectMap(G);
5806 isNew = true;
5807 } else {
5808 isNew = false;
5809 }
5810
5811 if(state < 0)
5812 state = I->State.size();
5813 if(I->State.size() <= state) {
5814 VecCheckEmplace(I->State, state, I->G);
5815 }
5816 ms = &I->State[state];
5817
5818 if(!PConvAttrToStrMaxLen(Map, "format", format, sizeof(WordType) - 1))
5819 ok = ErrMessage(G, "LoadChemPyMap", "bad 'format' parameter.");
5820 else if(!PConvAttrToFloatArrayInPlace(Map, "cell_dim", ms->Symmetry->Crystal.Dim, 3))
5821 ok = ErrMessage(G, "LoadChemPyMap", "bad 'cell_dim' parameter.");
5822 else if(!PConvAttrToFloatArrayInPlace(Map, "cell_ang", ms->Symmetry->Crystal.Angle, 3))
5823 ok = ErrMessage(G, "LoadChemPyMap", "bad 'cell_ang' parameter.");
5824 else if(!PConvAttrToIntArrayInPlace(Map, "cell_div", ms->Div, 3))
5825 ok = ErrMessage(G, "LoadChemPyMap", "bad 'cell_div' parameter.");
5826 else if(!PConvAttrToIntArrayInPlace(Map, "first", ms->Min, 3))
5827 ok = ErrMessage(G, "LoadChemPyMap", "bad 'first' parameter.");
5828 else if(!PConvAttrToIntArrayInPlace(Map, "last", ms->Max, 3))
5829 ok = ErrMessage(G, "LoadChemPyMap", "bad 'last' parameter.");
5830
5831 if(ok) {
5832 if(strcmp(format, "CObjectZYXfloat") == 0) {
5833 ok = PConvAttrToPtr(Map, "c_object", (void **) (void *) &cobj);
5834 if(!ok)
5835 ErrMessage(G, "LoadChemPyMap", "CObject unreadable.");
5836 } else {
5837 ok = ErrMessage(G, "LoadChemPyMap", "unsupported format.");
5838 }
5839 }
5840 /* good to go */
5841
5842 if(ok) {
5843 if(strcmp(format, "CObjectZYXfloat") == 0) {
5844
5845 ms->FDim[0] = ms->Max[0] - ms->Min[0] + 1;
5846 ms->FDim[1] = ms->Max[1] - ms->Min[1] + 1;
5847 ms->FDim[2] = ms->Max[2] - ms->Min[2] + 1;
5848 if(Feedback(G, FB_ObjectMap, FB_Actions)) {
5849 printf(" LoadChemPyMap: CObjectZYXdouble %dx%dx%d\n", ms->FDim[0], ms->FDim[1],
5850 ms->FDim[2]);
5851 }
5852 ms->FDim[3] = 3;
5853 if(!(ms->FDim[0] && ms->FDim[1] && ms->FDim[2]))
5854 ok = false;
5855 else {
5856 SymmetryUpdate(ms->Symmetry.get());
5857 ms->Field.reset(new Isofield(G, ms->FDim));
5858 for(c = 0; c < ms->FDim[2]; c++) {
5859 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
5860 for(b = 0; b < ms->FDim[1]; b++) {
5861 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
5862 for(a = 0; a < ms->FDim[0]; a++) {
5863 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
5864
5865 dens = *(cobj++);
5866
5867 F3(ms->Field->data, a, b, c) = dens;
5868 if(maxd < dens)
5869 maxd = dens;
5870 if(mind > dens)
5871 mind = dens;
5872 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
5873 for(e = 0; e < 3; e++)
5874 F4(ms->Field->points, a, b, c, e) = vr[e];
5875 }
5876 }
5877 }
5878
5879 if(ok) {
5880 d = 0;
5881 for(c = 0; c < ms->FDim[2]; c += (ms->FDim[2] - 1)) {
5882 v[2] = (c + ms->Min[2]) / ((float) ms->Div[2]);
5883 for(b = 0; b < ms->FDim[1]; b += (ms->FDim[1] - 1)) {
5884 v[1] = (b + ms->Min[1]) / ((float) ms->Div[1]);
5885 for(a = 0; a < ms->FDim[0]; a += (ms->FDim[0] - 1)) {
5886 v[0] = (a + ms->Min[0]) / ((float) ms->Div[0]);
5887 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, vr);
5888 copy3f(vr, ms->Corner + 3 * d);
5889 d++;
5890 }
5891 }
5892 }
5893 }
5894 }
5895 }
5896 }
5897
5898 if(ok) {
5899 CrystalDump(&ms->Symmetry->Crystal);
5900
5901 v[2] = (ms->Min[2]) / ((float) ms->Div[2]);
5902 v[1] = (ms->Min[1]) / ((float) ms->Div[1]);
5903 v[0] = (ms->Min[0]) / ((float) ms->Div[0]);
5904
5905 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMin);
5906
5907 v[2] = ((ms->FDim[2] - 1) + ms->Min[2]) / ((float) ms->Div[2]);
5908 v[1] = ((ms->FDim[1] - 1) + ms->Min[1]) / ((float) ms->Div[1]);
5909 v[0] = ((ms->FDim[0] - 1) + ms->Min[0]) / ((float) ms->Div[0]);
5910
5911 transform33f3f(ms->Symmetry->Crystal.FracToReal, v, ms->ExtentMax);
5912
5913 }
5914
5915 if(!ok) {
5916 ErrMessage(G, "ObjectMap", "Error reading map");
5917 } else {
5918 ms->Active = true;
5919 ObjectMapUpdateExtents(I);
5920 if(!quiet) {
5921 PRINTFB(I->G, FB_ObjectMap, FB_Results)
5922 " ObjectMap: Map read. Range: %5.3f to %5.3f\n", mind, maxd ENDFB(I->G);
5923 }
5924 }
5925
5926 if(ok) {
5927 SceneChanged(G);
5928 SceneCountFrames(G);
5929 }
5930 }
5931 return (I);
5932 #endif
5933 }
5934
ObjectMapDump(const ObjectMap * om,const char * fname,int state,int quiet)5935 void ObjectMapDump(const ObjectMap* om, const char* fname, int state, int quiet)
5936 {
5937 auto* oms = om->getObjectMapState(state);
5938
5939 if (!oms) {
5940 ErrMessage(om->G, __func__, "state out of range");
5941 return;
5942 }
5943
5944 auto file = fopen(fname, "wb");
5945 if (!file) {
5946 ErrMessage(om->G, __func__, "can't open file for writing");
5947 return;
5948 }
5949
5950 auto* field = oms->Field.get();
5951
5952 for (int xi = 0; xi < field->dimensions[0]; xi++) {
5953 for (int yi = 0; yi < field->dimensions[1]; yi++) {
5954 for (int zi = 0; zi < field->dimensions[2]; zi++) {
5955
5956 float x = Ffloat4(field->points, xi, yi, zi, 0);
5957 float y = Ffloat4(field->points, xi, yi, zi, 1);
5958 float z = Ffloat4(field->points, xi, yi, zi, 2);
5959
5960 switch (field->data->type) {
5961 case cFieldFloat: {
5962 float value = Ffloat3(field->data, xi, yi, zi);
5963 fprintf(file, "%10.4f%10.4f%10.4f%10.4f\n", x, y, z, value);
5964 break;
5965 }
5966 case cFieldInt: {
5967 int value = Fint3(field->data, xi, yi, zi);
5968 fprintf(file, "%10.4f%10.4f%10.4f%10d\n", x, y, z, value);
5969 break;
5970 }
5971 default:
5972 ErrMessage(om->G, __func__, "unknown field type");
5973 fclose(file);
5974 return;
5975 }
5976 }
5977 }
5978 }
5979 fclose(file);
5980 if (!quiet) {
5981 PRINTFB(om->G, FB_ObjectMap, FB_Actions)
5982 " ObjectMapDump: %s written to %s\n", om->Name, fname ENDFB(om->G);
5983 }
5984 }
5985
clone() const5986 CObject* ObjectMap::clone() const
5987 {
5988 return new ObjectMap(*this);
5989 }
5990
5991