ProSHADE  0.7.5.3 (FEB 2021)
Protein Shape Detection
ProSHADE_mapManip.cpp
Go to the documentation of this file.
1 
22 //==================================================== ProSHADE
23 #include "ProSHADE_mapManip.hpp"
24 
25 //==================================================== Define round for C++98
31 proshade_signed ProSHADE_internal_mapManip::myRound ( proshade_double x )
32 {
33 #if __cplusplus >= 201103L
34  return ( std::round ( x ) );
35 #else
36  return ( round ( x ) );
37 #endif
38 }
39 
57 void ProSHADE_internal_mapManip::determinePDBRanges ( gemmi::Structure pdbFile, proshade_single* xFrom, proshade_single* xTo, proshade_single* yFrom, proshade_single* yTo, proshade_single* zFrom, proshade_single* zTo, bool firstModel )
58 {
59  //================================================ Initialise structure crawl
60  bool firstAtom = true;
61 
62  //================================================ Use the first model, if it exists
63  if ( pdbFile.models.size() > 0 )
64  {
65  //============================================ For each model
66  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); sIt++ )
67  {
68  //======================================== Check if multiple models are allowed
69  if ( firstModel && ( sIt != 0 ) ) { break; }
70 
71  //======================================== Get model
72  gemmi::Model model = pdbFile.models.at(sIt);
73 
74  //======================================== For each chain
75  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model.chains.size() ); mIt++ )
76  {
77  //==================================== Get chain
78  gemmi::Chain chain = model.chains.at(mIt);
79 
80  //==================================== For each residue
81  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain.residues.size() ); rIt++ )
82  {
83  //================================ Get residue
84  gemmi::Residue residue = chain.residues.at(rIt);
85 
86  //================================ For each atom
87  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue.atoms.size() ); aIt++ )
88  {
89  //============================ Get atom
90  gemmi::Atom atom = residue.atoms.at(aIt);
91 
92  //============================ Ignore hydrogens, map computations ignore them anyway and inclusion here causes map - co-ordinate mismatches.
93  if ( atom.is_hydrogen() ) { continue; }
94 
95  //============================ Find the coordinate ranges
96  if ( firstAtom )
97  {
98  *xTo = static_cast<proshade_single> ( atom.pos.x );
99  *xFrom = static_cast<proshade_single> ( atom.pos.x );
100  *yTo = static_cast<proshade_single> ( atom.pos.y );
101  *yFrom = static_cast<proshade_single> ( atom.pos.y );
102  *zTo = static_cast<proshade_single> ( atom.pos.z );
103  *zFrom = static_cast<proshade_single> ( atom.pos.z );
104  firstAtom = false;
105  }
106  else
107  {
108  if ( static_cast<proshade_single> ( atom.pos.x ) > *xTo ) { *xTo = static_cast<proshade_single> ( atom.pos.x ); }
109  if ( static_cast<proshade_single> ( atom.pos.x ) < *xFrom ) { *xFrom = static_cast<proshade_single> ( atom.pos.x ); }
110  if ( static_cast<proshade_single> ( atom.pos.y ) > *yTo ) { *yTo = static_cast<proshade_single> ( atom.pos.y ); }
111  if ( static_cast<proshade_single> ( atom.pos.y ) < *yFrom ) { *yFrom = static_cast<proshade_single> ( atom.pos.y ); }
112  if ( static_cast<proshade_single> ( atom.pos.z ) > *zTo ) { *zTo = static_cast<proshade_single> ( atom.pos.z ); }
113  if ( static_cast<proshade_single> ( atom.pos.z ) < *zFrom ) { *zFrom = static_cast<proshade_single> ( atom.pos.z ); }
114  }
115  }
116  }
117  }
118  }
119  }
120  else
121  {
122  std::stringstream hlpSS;
123  hlpSS << "Found 0 models in input file " << pdbFile.name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
124  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
125  }
126 
127  //================================================ Done
128  return ;
129 
130 }
131 
142 void ProSHADE_internal_mapManip::findPDBCOMValues ( gemmi::Structure pdbFile, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, bool firstModel )
143 {
144  //================================================ Initialise structure crawl
145  proshade_double totAtoms = 0.0;
146  *xCom = 0.0;
147  *yCom = 0.0;
148  *zCom = 0.0;
149 
150  //================================================ Use the first model, if it exists
151  if ( pdbFile.models.size() > 0 )
152  {
153  //============================================ For each model
154  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); sIt++ )
155  {
156  //======================================== Get model
157  gemmi::Model model = pdbFile.models.at(sIt);
158 
159  //======================================== Check if multiple models are allowed
160  if ( firstModel && ( sIt != 0 ) ) { break; }
161 
162  //======================================== For each chain
163  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model.chains.size() ); mIt++ )
164  {
165  //==================================== Get chain
166  gemmi::Chain chain = model.chains.at(mIt);
167 
168  //==================================== For each residue
169  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain.residues.size() ); rIt++ )
170  {
171  //================================ Get residue
172  gemmi::Residue residue = chain.residues.at(rIt);
173 
174  //================================ For each atom
175  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue.atoms.size() ); aIt++ )
176  {
177  //============================ Get atom
178  gemmi::Atom atom = residue.atoms.at(aIt);
179 
180  //============================ Save the COM sums
181  *xCom += atom.pos.x * atom.element.weight();
182  *yCom += atom.pos.y * atom.element.weight();
183  *zCom += atom.pos.z * atom.element.weight();
184  totAtoms += atom.element.weight();
185  }
186  }
187  }
188  }
189  }
190  else
191  {
192  std::stringstream hlpSS;
193  hlpSS << "Found 0 models in input file " << pdbFile.name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
194  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
195  }
196 
197  //================================================ Normalise sums to COM
198  *xCom /= totAtoms;
199  *yCom /= totAtoms;
200  *zCom /= totAtoms;
201 
202  //================================================ Done
203  return ;
204 
205 }
206 
225 void ProSHADE_internal_mapManip::findMAPCOMValues ( proshade_double* map, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo )
226 {
227  //================================================ Initialise computation
228  proshade_double totDensity = 0.0;
229  *xCom = 0.0;
230  *yCom = 0.0;
231  *zCom = 0.0;
232  proshade_signed arrPos = 0;
233  proshade_double xSampRate = xAngs / static_cast<proshade_single> ( xTo - xFrom );
234  proshade_double ySampRate = yAngs / static_cast<proshade_single> ( yTo - yFrom );
235  proshade_double zSampRate = zAngs / static_cast<proshade_single> ( zTo - zFrom );
236 
237  //================================================ For each map point
238  for ( proshade_signed xIt = xFrom; xIt < xTo; xIt++ )
239  {
240  for ( proshade_signed yIt = yFrom; yIt < yTo; yIt++ )
241  {
242  for ( proshade_signed zIt = zFrom; zIt < zTo; zIt++ )
243  {
244  arrPos = (zIt-zFrom) + ( zTo - zFrom ) * ( (yIt-yFrom) + ( yTo - yFrom ) * (xIt-xFrom) );
245  totDensity += map[arrPos];
246  *xCom += static_cast<proshade_double> ( xIt * xSampRate ) * map[arrPos];
247  *yCom += static_cast<proshade_double> ( yIt * ySampRate ) * map[arrPos];
248  *zCom += static_cast<proshade_double> ( zIt * zSampRate ) * map[arrPos];
249  }
250  }
251  }
252 
253  //================================================ Normalise sums to COM
254  *xCom /= totDensity;
255  *yCom /= totDensity;
256  *zCom /= totDensity;
257 
258  //================================================ Done
259  return ;
260 
261 }
262 
277 void ProSHADE_internal_mapManip::rotatePDBCoordinates ( gemmi::Structure *pdbFile, proshade_double euA, proshade_double euB, proshade_double euG, proshade_double xCom,
278 proshade_double yCom, proshade_double zCom, bool firstModel )
279 {
280  //================================================ Convert Euler angles to rotation matrix
281  proshade_double *rotMat = new proshade_double[9];
282  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat, __FILE__, __LINE__, __func__ );
284 
285  //================================================ Initialise internal variables
286  proshade_single xTmp, yTmp, zTmp;
287 
288  //================================================ Use the first model, if it exists
289  if ( pdbFile->models.size() > 0 )
290  {
291  //============================================ For each model
292  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
293  {
294  //======================================== Get model
295  gemmi::Model *model = &pdbFile->models.at(sIt);
296 
297  //======================================== Check if multiple models are allowed
298  if ( firstModel && ( sIt != 0 ) ) { break; }
299 
300  //======================================== For each chain
301  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
302  {
303  //==================================== Get chain
304  gemmi::Chain *chain = &model->chains.at(mIt);
305 
306  //==================================== For each residue
307  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
308  {
309  //================================ Get residue
310  gemmi::Residue *residue = &chain->residues.at(rIt);
311 
312  //================================ For each atom
313  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
314  {
315  //============================ Get atom
316  gemmi::Atom *atom = &residue->atoms.at(aIt);
317 
318  //============================ Move to mid-point
319  xTmp = atom->pos.x - xCom;
320  yTmp = atom->pos.y - yCom;
321  zTmp = atom->pos.z - zCom;
322 
323  //============================ Rotate the atom position
324  atom->pos.x = ( xTmp * rotMat[0] ) + ( yTmp * rotMat[1] ) + ( zTmp * rotMat[2] );
325  atom->pos.y = ( xTmp * rotMat[3] ) + ( yTmp * rotMat[4] ) + ( zTmp * rotMat[5] );
326  atom->pos.z = ( xTmp * rotMat[6] ) + ( yTmp * rotMat[7] ) + ( zTmp * rotMat[8] );
327 
328  //============================ Move back
329  atom->pos.x = atom->pos.x + xCom;
330  atom->pos.y = atom->pos.y + yCom;
331  atom->pos.z = atom->pos.z + zCom;
332  }
333  }
334  }
335  }
336  }
337  else
338  {
339  std::stringstream hlpSS;
340  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
341  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
342  }
343 
344  //================================================ Release memory
345  delete[] rotMat;
346 
347  //================================================ Done
348  return ;
349 
350 }
351 
362 void ProSHADE_internal_mapManip::translatePDBCoordinates ( gemmi::Structure *pdbFile, proshade_double transX, proshade_double transY, proshade_double transZ, bool firstModel )
363 {
364  //================================================ Use the first model, if it exists
365  if ( pdbFile->models.size() > 0 )
366  {
367  //============================================ For each model
368  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
369  {
370  //======================================== Check if multiple models are allowed
371  if ( firstModel && ( sIt != 0 ) ) { break; }
372 
373  //======================================== Get model
374  gemmi::Model *model = &pdbFile->models.at(sIt);
375 
376  //======================================== For each chain
377  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
378  {
379  //==================================== Get chain
380  gemmi::Chain *chain = &model->chains.at(mIt);
381 
382  //==================================== For each residue
383  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
384  {
385  //================================ Get residue
386  gemmi::Residue *residue = &chain->residues.at(rIt);
387 
388  //================================ For each atom
389  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
390  {
391  //============================ Get atom
392  gemmi::Atom *atom = &residue->atoms.at(aIt);
393 
394  //============================ Translate
395  atom->pos.x += transX;
396  atom->pos.y += transY;
397  atom->pos.z += transZ;
398  }
399  }
400  }
401  }
402  }
403  else
404  {
405  std::stringstream hlpSS;
406  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
407  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
408  }
409 
410  //================================================ Done
411  return ;
412 
413 }
414 
425 void ProSHADE_internal_mapManip::changePDBBFactors ( gemmi::Structure *pdbFile, proshade_double newBFactorValue, bool firstModel )
426 {
427  //================================================ Use the first model, if it exists
428  if ( pdbFile->models.size() > 0 )
429  {
430  //============================================ For each model
431  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
432  {
433  //======================================== Check if multiple models are allowed
434  if ( firstModel && ( sIt != 0 ) ) { break; }
435 
436  //======================================== Get model
437  gemmi::Model *model = &pdbFile->models.at(sIt);
438 
439  //======================================== For each chain
440  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
441  {
442  //==================================== Get chain
443  gemmi::Chain *chain = &model->chains.at(mIt);
444 
445  //==================================== For each residue
446  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
447  {
448  //================================ Get residue
449  gemmi::Residue *residue = &chain->residues.at(rIt);
450 
451  //================================ For each atom
452  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
453  {
454  //============================ Get atom
455  gemmi::Atom *atom = &residue->atoms.at(aIt);
456 
457  //============================ Change the B-factors
458  atom->b_iso = newBFactorValue;
459  }
460  }
461  }
462  }
463  }
464  else
465  {
466  std::stringstream hlpSS;
467  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
468  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
469  }
470 
471  //================================================ Done
472  return ;
473 
474 }
475 
485 void ProSHADE_internal_mapManip::removeWaters ( gemmi::Structure *pdbFile, bool firstModel )
486 {
487  //================================================ Use the first model, if it exists
488  if ( pdbFile->models.size() > 0 )
489  {
490  //============================================ For each model
491  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
492  {
493  //======================================== Check if multiple models are allowed
494  if ( firstModel && ( sIt != 0 ) ) { break; }
495 
496  //======================================== Get model
497  gemmi::Model *model = &pdbFile->models.at(sIt);
498 
499  //======================================== For each chain
500  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
501  {
502  //==================================== Get chain
503  gemmi::Chain *chain = &model->chains.at(mIt);
504 
505  //==================================== Initialise del vector
506  std::vector< proshade_unsign > delVec;
507 
508  //==================================== For each residue
509  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
510  {
511  //================================ Get residue
512  gemmi::Residue *residue = &chain->residues.at(rIt);
513 
514  //================================ If residue is water
515  if ( residue->is_water() )
516  {
518  }
519  }
520 
521  //==================================== Delete from end to avoid indexing issues
522  std::sort ( delVec.begin(), delVec.end(), std::greater<int>() );
523  for ( proshade_unsign vecIt = 0; vecIt < static_cast<proshade_unsign> ( delVec.size() ); vecIt++ )
524  {
525  chain->residues.erase ( chain->residues.begin() + delVec.at(vecIt) );
526  }
527  }
528  }
529  }
530  else
531  {
532  std::stringstream hlpSS;
533  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
534  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
535  }
536 
537  //================================================ Done
538  return ;
539 
540 }
541 
554 void ProSHADE_internal_mapManip::movePDBForMapCalc ( gemmi::Structure *pdbFile, proshade_single xMov, proshade_single yMov, proshade_single zMov, bool firstModel )
555 {
556  //================================================ Use the first model, if it exists
557  if ( pdbFile->models.size() > 0 )
558  {
559  //============================================ For each model
560  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
561  {
562  //======================================== Check if multiple models are allowed
563  if ( firstModel && ( sIt != 0 ) ) { break; }
564 
565  //======================================== Get model
566  gemmi::Model *model = &pdbFile->models.at(sIt);
567 
568  //======================================== For each chain
569  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
570  {
571  //==================================== Get chain
572  gemmi::Chain *chain = &model->chains.at(mIt);
573 
574  //==================================== For each residue
575  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
576  {
577  //================================ Get residue
578  gemmi::Residue *residue = &chain->residues.at(rIt);
579 
580  //================================ For each atom
581  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
582  {
583  //============================ Get atom
584  gemmi::Atom *atom = &residue->atoms.at(aIt);
585 
586  //============================ Move the atoms
587  atom->pos = gemmi::Position ( atom->pos.x + xMov, atom->pos.y + yMov, atom->pos.z + zMov );
588  }
589  }
590  }
591 
592  }
593  }
594  else
595  {
596  std::stringstream hlpSS;
597  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
598  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
599  }
600 
601  //================================================ Done
602  return ;
603 }
604 
625 void ProSHADE_internal_mapManip::generateMapFromPDB ( gemmi::Structure pdbFile, proshade_double*& map, proshade_single requestedResolution, proshade_single xCell, proshade_single yCell, proshade_single zCell, proshade_signed* xTo, proshade_signed* yTo, proshade_signed* zTo, bool forceP1, bool firstModel )
626 {
627  //================================================ Set cell dimensions from the increased ranges (we need to add some space) and re-calculate cell properties
628  if ( forceP1 ) { pdbFile.cell = gemmi::UnitCell(); }
629  pdbFile.cell.a = xCell;
630  pdbFile.cell.b = yCell;
631  pdbFile.cell.c = zCell;
632  pdbFile.cell.calculate_properties ( );
633 
634  //================================================ Get elements in Gemmi format
635  std::string totElString;
636  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); mIt++ )
637  {
638  //============================================ Check if multiple models are allowed
639  if ( firstModel && ( mIt != 0 ) )
640  {
641  std::stringstream hlpSS;
642  hlpSS << "!!! ProSHADE WARNING !!! Found multiple models (" << pdbFile.models.size() << ") in input file " << pdbFile.name << ", while the settings state that only the first PDB file model should be used. If all models should be used, please supply ProSHADE with the \"-x\" option.";
643  ProSHADE_internal_messages::printWarningMessage ( 0, hlpSS.str(), "WP00055" );
644  break;
645  }
646 
647  std::string hlpStr = pdbFile.models[mIt].present_elements ( ).to_string<char,std::char_traits<char>,std::allocator<char> >();
648  totElString = totElString + hlpStr;
649  }
650  std::bitset<(size_t)gemmi::El::END> present_elems ( totElString );
651 
652  //================================================ Sanity checks
653  if ( present_elems[static_cast<int> ( gemmi::El::X )] )
654  {
655  throw ProSHADE_exception ( "Found unknown element in input file.", "EP00051", __FILE__, __LINE__, __func__, "Gemmi library does not recognise some of the elements in\n : the co-ordinate file. Please check the file for not being\n : corrupted and containing standard elements." );
656  }
657 
658  for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( present_elems.size() ); elIt++ )
659  {
660  if ( present_elems[elIt] && !gemmi::IT92<double>::has ( static_cast<gemmi::El> ( elIt ) ) )
661  {
662  std::stringstream hlpSS;
663  hlpSS << "Missing form factor for element " << element_name ( static_cast<gemmi::El> ( elIt ) );
664  throw ProSHADE_exception ( hlpSS.str().c_str(), "EP00052", __FILE__, __LINE__, __func__, "Gemmi library does not have a form factor value for this\n : reported element. Please report this to the author." );
665  }
666  }
667 
668  //================================================ Compute the f's
669  double wavelength = 0.1;
670  double energy = gemmi::hc() / wavelength;
671 
672  //================================================ Create the density calculator object and fill it in
673  gemmi::DensityCalculator<gemmi::IT92<double>, float> dencalc;
674 
675  dencalc.d_min = static_cast<double> ( requestedResolution );
676  for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( present_elems.size() ); elIt++ ) { if ( present_elems[elIt] ) { dencalc.addends.set ( static_cast<gemmi::El> ( elIt ), gemmi::cromer_libermann ( elIt, energy, nullptr ) ); } }
677  dencalc.set_grid_cell_and_spacegroup ( pdbFile );
678 
679  //================================================ Force P1 spacegroup
680  if ( forceP1 ) { dencalc.grid.spacegroup = &gemmi::get_spacegroup_p1(); }
681 
682  //================================================ Compute the theoretical map for each model
683  dencalc.grid.data.clear ( );
684  dencalc.grid.set_size_from_spacing ( dencalc.d_min / ( 2.0 * dencalc.rate), true );
685  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); mIt++ )
686  {
687  if ( firstModel && ( mIt != 0 ) ) { break; }
688  dencalc.add_model_density_to_grid ( pdbFile.models[mIt] );
689  dencalc.grid.symmetrize ( [](float a, float b) { return a + b; } );
690  }
691 
692  //================================================ Get the map
693  const gemmi::Grid<float>& grid = dencalc.grid;
694 
695  //================================================ Save the map dimensions
696  *xTo = grid.nu;
697  *yTo = grid.nv;
698  *zTo = grid.nw;
699 
700  //================================================ Copy the gemmi::Grid to my map format
701  map = new proshade_double [(*xTo) * (*yTo) * (*zTo)];
702  ProSHADE_internal_misc::checkMemoryAllocation ( map, __FILE__, __LINE__, __func__ );
703 
704  proshade_signed arrPos = 0;
705  for ( proshade_signed uIt = 0; uIt < (*xTo); uIt++ )
706  {
707  for ( proshade_signed vIt = 0; vIt < (*yTo); vIt++ )
708  {
709  for ( proshade_signed wIt = 0; wIt < (*zTo); wIt++ )
710  {
711  arrPos = wIt + (*zTo) * ( vIt + (*yTo) * uIt );
712  map[arrPos] = grid.get_value_q( uIt, vIt, wIt );
713  }
714  }
715  }
716 
717  //================================================ Done
718  return ;
719 
720 }
721 
744 void ProSHADE_internal_mapManip::moveMapByIndices ( proshade_single* xMov, proshade_single* yMov, proshade_single* zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed* xFrom, proshade_signed* xTo, proshade_signed* yFrom, proshade_signed* yTo, proshade_signed* zFrom, proshade_signed* zTo, proshade_signed* xOrigin, proshade_signed* yOrigin, proshade_signed* zOrigin )
745 {
746  //================================================ Compute movement in indices
747  proshade_signed xIndMove = static_cast<proshade_signed> ( std::floor ( -(*xMov) / ( xAngs / (*xTo - *xFrom + 1) ) ) );
748  proshade_signed yIndMove = static_cast<proshade_signed> ( std::floor ( -(*yMov) / ( yAngs / (*yTo - *yFrom + 1) ) ) );
749  proshade_signed zIndMove = static_cast<proshade_signed> ( std::floor ( -(*zMov) / ( zAngs / (*zTo - *zFrom + 1) ) ) );
750 
751  //================================================ Set the movs to the remainder
752  *xMov = -( *xMov ) - ( xIndMove * ( xAngs / (*xTo - *xFrom + 1) ) );
753  *yMov = -( *yMov ) - ( yIndMove * ( yAngs / (*yTo - *yFrom + 1) ) );
754  *zMov = -( *zMov ) - ( zIndMove * ( zAngs / (*zTo - *zFrom + 1) ) );
755 
756  //================================================ Move indices by as much
757  *xFrom += xIndMove;
758  *xTo += xIndMove;
759  *yFrom += yIndMove;
760  *yTo += yIndMove;
761  *zFrom += zIndMove;
762  *zTo += zIndMove;
763 
764  //================================================ And set origin to reflect the changes
765  *xOrigin = *xFrom;
766  *yOrigin = *yFrom;
767  *zOrigin = *zFrom;
768 
769  //================================================ Done
770  return ;
771 
772 }
773 
792 void ProSHADE_internal_mapManip::moveMapByFourier ( proshade_double*& map, proshade_single xMov, proshade_single yMov, proshade_single zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim )
793 {
794  //================================================ Local variables initialisation
795  proshade_signed arrayPos = 0;
796  proshade_signed h, k, l;
797  proshade_double real = 0.0;
798  proshade_double imag = 0.0;
799  proshade_double trCoeffReal, trCoeffImag;
800  proshade_double normFactor = static_cast<proshade_double> ( xDim * yDim * zDim );
801  proshade_double exponent = 0.0;
802  proshade_double hlpArrReal;
803  proshade_double hlpArrImag;
804 
805  //================================================ Create Fourier map variable
806  fftw_complex *fCoeffs = new fftw_complex [xDim * yDim * zDim];
807  fftw_complex *translatedMap = new fftw_complex [xDim * yDim * zDim];
808 
809  //================================================ Check memory allocation
810  ProSHADE_internal_misc::checkMemoryAllocation ( fCoeffs, __FILE__, __LINE__, __func__ );
811  ProSHADE_internal_misc::checkMemoryAllocation ( translatedMap, __FILE__, __LINE__, __func__ );
812 
813  //================================================ Create plans
814  fftw_plan planForwardFourier = fftw_plan_dft_3d ( xDim, yDim, zDim, translatedMap, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
815  fftw_plan planBackwardFourier = fftw_plan_dft_3d ( xDim, yDim, zDim, fCoeffs, translatedMap, FFTW_BACKWARD, FFTW_ESTIMATE );
816 
817  //================================================ Copy map to complex format
818  for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
819  {
820  for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
821  {
822  for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
823  {
824  arrayPos = wIt + zDim * ( vIt + yDim * uIt );
825 
826  if ( map[arrayPos] == map[arrayPos] ) { translatedMap[arrayPos][0] = map[arrayPos]; }
827  else { translatedMap[arrayPos][0] = 0.0; }
828  translatedMap[arrayPos][1] = 0.0;
829  }
830  }
831  }
832 
833  //================================================ Compute Forward Fourier
834  fftw_execute ( planForwardFourier );
835 
836  //================================================ Add the required shift
837  for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
838  {
839  for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
840  {
841  for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
842  {
843  //==================================== Var init
844  arrayPos = wIt + zDim * ( vIt + yDim * uIt );
845  real = fCoeffs[arrayPos][0];
846  imag = fCoeffs[arrayPos][1];
847 
848  //==================================== Change the B-factors, if required
849  if ( uIt > static_cast<proshade_unsign> ( (xDim+1) / 2) ) { h = uIt - static_cast <proshade_signed> ( xDim ); } else { h = uIt; }
850  if ( vIt > static_cast<proshade_unsign> ( (yDim+1) / 2) ) { k = vIt - static_cast <proshade_signed> ( yDim ); } else { k = vIt; }
851  if ( wIt > static_cast<proshade_unsign> ( (zDim+1) / 2) ) { l = wIt - static_cast <proshade_signed> ( zDim ); } else { l = wIt; }
852 
853  //==================================== Get translation coefficient change
854  exponent = ( ( ( static_cast <proshade_double> ( h ) / static_cast <proshade_double> ( xAngs ) ) * (-xMov) ) +
855  ( ( static_cast <proshade_double> ( k ) / static_cast <proshade_double> ( yAngs ) ) * (-yMov) ) +
856  ( ( static_cast <proshade_double> ( l ) / static_cast <proshade_double> ( zAngs ) ) * (-zMov) ) ) * 2.0 * M_PI;
857 
858  trCoeffReal = cos ( exponent );
859  trCoeffImag = sin ( exponent );
860  ProSHADE_internal_maths::complexMultiplication ( &real, &imag, &trCoeffReal, &trCoeffImag, &hlpArrReal, &hlpArrImag );
861 
862  //==================================== Save the mask data
863  fCoeffs[arrayPos][0] = hlpArrReal / normFactor;
864  fCoeffs[arrayPos][1] = hlpArrImag / normFactor;
865  }
866  }
867  }
868 
869  //================================================ Compute inverse Fourier
870  fftw_execute ( planBackwardFourier );
871 
872  //================================================ Copy back to map
873  for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
874  {
875  for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
876  {
877  for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
878  {
879  arrayPos = wIt + zDim * ( vIt + yDim * uIt );
880  map[arrayPos] = translatedMap[arrayPos][0];
881  }
882  }
883  }
884 
885  //================================================ Release memory
886  fftw_destroy_plan ( planForwardFourier );
887  fftw_destroy_plan ( planBackwardFourier );
888  delete[] fCoeffs;
889  delete[] translatedMap;
890 
891  //================================================ Done
892  return ;
893 
894 }
895 
912 void ProSHADE_internal_mapManip::blurSharpenMap ( proshade_double*& map, proshade_double*& blurredMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single blurringFactor )
913 {
914  //================================================ Set local variables
915  proshade_signed xDim = xDimS;
916  proshade_signed yDim = yDimS;
917  proshade_signed zDim = zDimS;
918  proshade_double real, imag, S, mag, phase;
919  proshade_signed h, k, l;
920  proshade_unsign arrayPos = 0;
921  proshade_double normFactor = static_cast<proshade_double> ( xDim * yDim * zDim );
922 
923  //================================================ Copy map for processing
924  fftw_complex* mapCoeffs = new fftw_complex[xDim * yDim * zDim];
925  fftw_complex* mapMask = new fftw_complex[xDim * yDim * zDim];
926 
927  //================================================ Check memory allocation
928  ProSHADE_internal_misc::checkMemoryAllocation ( mapCoeffs, __FILE__, __LINE__, __func__ );
929  ProSHADE_internal_misc::checkMemoryAllocation ( mapMask, __FILE__, __LINE__, __func__ );
930 
931  //================================================ Copy data to mask
932  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> (xDim * yDim * zDim); iter++ )
933  {
934  mapMask[iter][0] = map[iter];
935  mapMask[iter][1] = 0.0;
936  }
937 
938  //================================================ Prepare FFTW plans
939  fftw_plan forward = fftw_plan_dft_3d ( xDim, yDim, zDim, mapMask, mapCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
940  fftw_plan inverse = fftw_plan_dft_3d ( xDim, yDim, zDim, mapCoeffs, mapMask, FFTW_BACKWARD, FFTW_ESTIMATE );
941 
942  //================================================ Run forward Fourier
943  fftw_execute ( forward );
944 
945  //================================================ Blur the coeffs
946  for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
947  {
948  for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
949  {
950  for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
951  {
952  //==================================== Var init
953  arrayPos = wIt + zDim * ( vIt + yDim * uIt );
954  real = mapCoeffs[arrayPos][0];
955  imag = mapCoeffs[arrayPos][1];
956 
957  //==================================== Change the B-factors, if required
958  if ( uIt > static_cast<proshade_unsign> ( (xDim+1) / 2) ) { h = uIt - static_cast <proshade_signed> ( xDim ); } else { h = uIt; }
959  if ( vIt > static_cast<proshade_unsign> ( (yDim+1) / 2) ) { k = vIt - static_cast <proshade_signed> ( yDim ); } else { k = vIt; }
960  if ( wIt > static_cast<proshade_unsign> ( (zDim+1) / 2) ) { l = wIt - static_cast <proshade_signed> ( zDim ); } else { l = wIt; }
961 
962  //====================================Get magnitude and phase with mask parameters
963  S = ( pow( static_cast<proshade_double> ( h ) / xAngs, 2.0 ) +
964  pow( static_cast<proshade_double> ( k ) / yAngs, 2.0 ) +
965  pow( static_cast<proshade_double> ( l ) / zAngs, 2.0 ) );
966  mag = std::sqrt ( (real*real) + (imag*imag) ) * std::exp ( - ( ( blurringFactor * S ) / 4.0 ) );
967  phase = std::atan2 ( imag, real );
968 
969  //==================================== Save the mask data
970  mapCoeffs[arrayPos][0] = ( mag * cos(phase) ) / normFactor;
971  mapCoeffs[arrayPos][1] = ( mag * sin(phase) ) / normFactor;
972  }
973  }
974  }
975 
976  //================================================ Run inverse Fourier
977  fftw_execute ( inverse );
978 
979  //================================================ Save the results
980  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> (xDim * yDim * zDim); iter++ )
981  {
982  blurredMap[iter] = mapMask[iter][0];
983  }
984 
985  //================================================ Release memory
986  delete[] mapMask;
987  delete[] mapCoeffs;
988 
989  //================================================ Delete FFTW plans
990  fftw_destroy_plan ( forward );
991  fftw_destroy_plan ( inverse );
992 
993  //================================================ Done
994  return ;
995 
996 }
997 
1012 void ProSHADE_internal_mapManip::getMaskFromBlurr ( proshade_double*& blurMap, proshade_double*& outMap, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single noIQRs )
1013 {
1014  //================================================ Initialise vector for map data
1015  std::vector<proshade_double> mapVals ( xDim * yDim * zDim, 0.0 );
1016 
1017  //================================================ Save map values in vector
1018  for ( proshade_unsign iter = 0; iter < ( xDim * yDim * zDim ); iter++ )
1019  {
1020  mapVals.at(iter) = blurMap[iter];
1021  }
1022 
1023  //================================================ Find median and IQRs
1024  proshade_double* medAndIQR = new proshade_double[2];
1025  ProSHADE_internal_maths::vectorMedianAndIQR ( &mapVals, medAndIQR );
1026 
1027  //================================================ Find the threshold
1028  proshade_double maskThreshold = medAndIQR[0] + ( medAndIQR[1] * static_cast<proshade_double> ( noIQRs ) );
1029 
1030  //================================================ Apply threshold
1031  for ( proshade_unsign iter = 0; iter < ( xDim * yDim * zDim ); iter++ )
1032  {
1033  if ( blurMap[iter] < maskThreshold )
1034  {
1035  outMap[iter] = 0.0;
1036  blurMap[iter] = 0.0;
1037  }
1038  }
1039 
1040  //================================================ Release vector values
1041  mapVals.clear ( );
1042 
1043  //================================================ Release memory
1044  delete[] medAndIQR;
1045 
1046  //================================================ Done
1047  return ;
1048 
1049 }
1050 
1065 void ProSHADE_internal_mapManip::getNonZeroBounds ( proshade_double* map, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed*& ret )
1066 {
1067  //================================================ Initialise local variables
1068  proshade_unsign arrayPos = 0;
1069 
1070  //================================================ Initialise result variable
1071  ret[0] = xDim;
1072  ret[1] = 0;
1073  ret[2] = yDim;
1074  ret[3] = 0;
1075  ret[4] = zDim;
1076  ret[5] = 0;
1077 
1078  //================================================ Iterate through map and check bounds
1079  for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1080  {
1081  for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1082  {
1083  for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1084  {
1085  //==================================== Var init
1086  arrayPos = zIt + zDim * ( yIt + yDim * xIt );
1087 
1088  //==================================== Check bounds
1089  if ( map[arrayPos] > 0.001 )
1090  {
1091  if ( xIt < ret[0] ) { ret[0] = xIt; }
1092  if ( xIt > ret[1] ) { ret[1] = xIt; }
1093  if ( yIt < ret[2] ) { ret[2] = yIt; }
1094  if ( yIt > ret[3] ) { ret[3] = yIt; }
1095  if ( zIt < ret[4] ) { ret[4] = zIt; }
1096  if ( zIt > ret[5] ) { ret[5] = zIt; }
1097  }
1098  }
1099  }
1100  }
1101 
1102  //================================================ Done
1103  return ;
1104 
1105 }
1106 
1122 void ProSHADE_internal_mapManip::addExtraBoundSpace ( proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed*& bounds, proshade_single extraSpace )
1123 {
1124  //================================================ Convert angstrom distance to indices
1125  proshade_signed xExtraInds = ProSHADE_internal_mapManip::myRound ( extraSpace / ( xAngs / static_cast<proshade_single> ( xDim ) ) );
1126  proshade_signed yExtraInds = ProSHADE_internal_mapManip::myRound ( extraSpace / ( yAngs / static_cast<proshade_single> ( yDim ) ) );
1127  proshade_signed zExtraInds = ProSHADE_internal_mapManip::myRound ( extraSpace / ( zAngs / static_cast<proshade_single> ( zDim ) ) );
1128 
1129  //=============================================== Changed bounds even if exceeding physical map - this will be deal with in the map creation part
1130  bounds[0] = bounds[0] - xExtraInds;
1131  bounds[1] = bounds[1] + xExtraInds;
1132  bounds[2] = bounds[2] - yExtraInds;
1133  bounds[3] = bounds[3] + yExtraInds;
1134  bounds[4] = bounds[4] - zExtraInds;
1135  bounds[5] = bounds[5] + zExtraInds;
1136 
1137  //================================================ Done
1138  return ;
1139 
1140 }
1141 
1158 void ProSHADE_internal_mapManip::reSampleMapToResolutionTrilinear ( proshade_double*& map, proshade_single resolution, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single*& corrs )
1159 {
1160  //================================================ Sanity check - the resolution needs to be set
1161  if ( resolution <= 0.0 )
1162  {
1163  throw ProSHADE_exception ( "Requested resolution not set for map re-sampling.", "EM00015", __FILE__, __LINE__, __func__, "There is no resolution value set, but map re-sampling to\n : this unset resolution value is required. This error\n : occurs when a task with no resolution requirement is\n : requested on a map data and the map resolution change is\n : set to \'on\'. Either supply a resolution value, or do not\n : re-sample the map." );
1164  }
1165 
1166  //================================================ Initialise local variables
1167  proshade_signed xDim = static_cast<proshade_signed> ( xDimS );
1168  proshade_signed yDim = static_cast<proshade_signed> ( yDimS );
1169  proshade_signed zDim = static_cast<proshade_signed> ( zDimS );
1170  proshade_single oldXSample = ( xAngs / static_cast<proshade_single> ( xDim ) );
1171  proshade_single oldYSample = ( yAngs / static_cast<proshade_single> ( yDim ) );
1172  proshade_single oldZSample = ( zAngs / static_cast<proshade_single> ( zDim ) );
1173  proshade_single newXSample = resolution / 2.0;
1174  proshade_single newYSample = resolution / 2.0;
1175  proshade_single newZSample = resolution / 2.0;
1176 
1177  //================================================ Compute required grid size
1178  proshade_signed newXDim = static_cast<proshade_signed> ( std::ceil ( xAngs / newXSample ) );
1179  proshade_signed newYDim = static_cast<proshade_signed> ( std::ceil ( yAngs / newYSample ) );
1180  proshade_signed newZDim = static_cast<proshade_signed> ( std::ceil ( zAngs / newZSample ) );
1181 
1182  //================================================ Create a new map variable
1183  proshade_double* newMap = new proshade_double [newXDim * newYDim * newZDim];
1184 
1185  //================================================ For each new map point
1186  proshade_signed xBottom = 0, xTop, yBottom = 0, yTop, zBottom = 0, zTop, oldMapIndex, newMapIndex;
1187  std::vector<proshade_double> c000 = std::vector<proshade_double> ( 4, 0.0 );
1188  std::vector<proshade_double> c001 = std::vector<proshade_double> ( 4, 0.0 );
1189  std::vector<proshade_double> c010 = std::vector<proshade_double> ( 4, 0.0 );
1190  std::vector<proshade_double> c011 = std::vector<proshade_double> ( 4, 0.0 );
1191  std::vector<proshade_double> c100 = std::vector<proshade_double> ( 4, 0.0 );
1192  std::vector<proshade_double> c101 = std::vector<proshade_double> ( 4, 0.0 );
1193  std::vector<proshade_double> c110 = std::vector<proshade_double> ( 4, 0.0 );
1194  std::vector<proshade_double> c111 = std::vector<proshade_double> ( 4, 0.0 );
1195  std::vector<proshade_double> c00 = std::vector<proshade_double> ( 4, 0.0 );
1196  std::vector<proshade_double> c01 = std::vector<proshade_double> ( 4, 0.0 );
1197  std::vector<proshade_double> c10 = std::vector<proshade_double> ( 4, 0.0 );
1198  std::vector<proshade_double> c11 = std::vector<proshade_double> ( 4, 0.0 );
1199  std::vector<proshade_double> c0 = std::vector<proshade_double> ( 4, 0.0 );
1200  std::vector<proshade_double> c1 = std::vector<proshade_double> ( 4, 0.0 );
1201  proshade_double xRelative, yRelative, zRelative;
1202 
1203  for ( proshade_signed xIt = 0; xIt < newXDim; xIt++ )
1204  {
1205  for ( proshade_signed yIt = 0; yIt < newYDim; yIt++ )
1206  {
1207  for ( proshade_signed zIt = 0; zIt < newZDim; zIt++ )
1208  {
1209  //==================================== Get this point's index
1210  newMapIndex = zIt + newZDim * ( yIt + newYDim * xIt );
1211 
1212  //==================================== Find this points bottom and top positions in the old map (including periodicity)
1213  for ( proshade_unsign ox = 0; ox < (xDimS-1); ox++ ) { if ( ( ( xIt * newXSample ) >= ( ox * oldXSample ) ) && ( ( xIt * newXSample ) <= ( (ox+1) * oldXSample ) ) ) { xBottom = static_cast<proshade_signed> ( ox ); break; } }
1214  for ( proshade_unsign oy = 0; oy < (yDimS-1); oy++ ) { if ( ( ( yIt * newYSample ) >= ( oy * oldYSample ) ) && ( ( yIt * newYSample ) <= ( (oy+1) * oldYSample ) ) ) { yBottom = static_cast<proshade_signed> ( oy ); break; } }
1215  for ( proshade_unsign oz = 0; oz < (zDimS-1); oz++ ) { if ( ( ( zIt * newZSample ) >= ( oz * oldZSample ) ) && ( ( zIt * newZSample ) <= ( (oz+1) * oldZSample ) ) ) { zBottom = static_cast<proshade_signed> ( oz ); break; } }
1216  xTop = xBottom + 1;
1217  yTop = yBottom + 1;
1218  zTop = zBottom + 1;
1219 
1220  //==================================== Find the surrounding point's values from the original map
1221  oldMapIndex = zBottom + zDimS * ( yBottom + yDimS * xBottom );
1222  c000.at(0) = static_cast<proshade_double> ( xBottom * oldXSample );
1223  c000.at(1) = static_cast<proshade_double> ( yBottom * oldYSample );
1224  c000.at(2) = static_cast<proshade_double> ( zBottom * oldZSample );
1225  c000.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1226 
1227  oldMapIndex = zTop + zDimS * ( yBottom + yDimS * xBottom );
1228  c001.at(0) = static_cast<proshade_double> ( xBottom * oldXSample );
1229  c001.at(1) = static_cast<proshade_double> ( yBottom * oldYSample );
1230  c001.at(2) = static_cast<proshade_double> ( zTop * oldZSample );
1231  c001.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1232 
1233  oldMapIndex = zBottom + zDimS * ( yTop + yDimS * xBottom );
1234  c010.at(0) = static_cast<proshade_double> ( xBottom * oldXSample );
1235  c010.at(1) = static_cast<proshade_double> ( yTop * oldYSample );
1236  c010.at(2) = static_cast<proshade_double> ( zBottom * oldZSample );
1237  c010.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1238 
1239  oldMapIndex = zTop + zDimS * ( yTop + yDimS * xBottom );
1240  c011.at(0) = static_cast<proshade_double> ( xBottom * oldXSample );
1241  c011.at(1) = static_cast<proshade_double> ( yTop * oldYSample );
1242  c011.at(2) = static_cast<proshade_double> ( zTop * oldZSample );
1243  c011.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1244 
1245  oldMapIndex = zBottom + zDimS * ( yBottom + yDimS * xTop );
1246  c100.at(0) = static_cast<proshade_double> ( xTop * oldXSample );
1247  c100.at(1) = static_cast<proshade_double> ( yBottom * oldYSample );
1248  c100.at(2) = static_cast<proshade_double> ( zBottom * oldZSample );
1249  c100.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1250 
1251  oldMapIndex = zTop + zDimS * ( yBottom + yDimS * xTop );
1252  c101.at(0) = static_cast<proshade_double> ( xTop * oldXSample );
1253  c101.at(1) = static_cast<proshade_double> ( yBottom * oldYSample );
1254  c101.at(2) = static_cast<proshade_double> ( zTop * oldZSample );
1255  c101.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1256 
1257  oldMapIndex = zBottom + zDimS * ( yTop + yDimS * xTop );
1258  c110.at(0) = static_cast<proshade_double> ( xTop * oldXSample );
1259  c110.at(1) = static_cast<proshade_double> ( yTop * oldYSample );
1260  c110.at(2) = static_cast<proshade_double> ( zBottom * oldZSample );
1261  c110.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1262 
1263  oldMapIndex = zTop + zDimS * ( yTop + yDimS * xTop );
1264  c111.at(0) = static_cast<proshade_double> ( xTop * oldXSample );
1265  c111.at(1) = static_cast<proshade_double> ( yTop * oldYSample );
1266  c111.at(2) = static_cast<proshade_double> ( zTop * oldZSample );
1267  c111.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1268 
1269  //==================================== Interpolate to the new grid along X
1270  xRelative = ( ( xIt * newXSample ) - ( xBottom * oldXSample ) ) / ( ( xTop * oldXSample ) - ( xBottom * oldXSample ) );
1271 
1272  //==================================== Interpolate for the less less point
1273  c00.at(0) = ( newXSample * xRelative ) + c000.at(0);
1274  c00.at(1) = c000.at(1);
1275  c00.at(2) = c000.at(2);
1276  c00.at(3) = ( c000.at(3) * ( 1.0 - xRelative ) ) + ( c100.at(3) * xRelative );
1277 
1278  //==================================== Interpolate for the less more point
1279  c01.at(0) = ( newXSample * xRelative ) + c001.at(0);
1280  c01.at(1) = c001.at(1);
1281  c01.at(2) = c001.at(2);
1282  c01.at(3) = ( c001.at(3) * ( 1.0 - xRelative ) ) + ( c101.at(3) * xRelative );
1283 
1284  //==================================== Interpolate for the more less point
1285  c10.at(0) = ( newXSample * xRelative ) + c010.at(0);
1286  c10.at(1) = c010.at(1);
1287  c10.at(2) = c010.at(2);
1288  c10.at(3) = ( c010.at(3) * ( 1.0 - xRelative ) ) + ( c110.at(3) * xRelative );
1289 
1290  //==================================== Interpolate for the more more point
1291  c11.at(0) = ( newXSample * xRelative ) + c011.at(0);
1292  c11.at(1) = c011.at(1);
1293  c11.at(2) = c011.at(2);
1294  c11.at(3) = ( c011.at(3) * ( 1.0 - xRelative ) ) + ( c111.at(3) * xRelative );
1295 
1296  //==================================== Interpolate to the new grid along Y
1297  yRelative = ( ( yIt * newYSample ) - ( yBottom * oldYSample ) ) / ( ( yTop * oldYSample ) - ( yBottom * oldYSample ) );
1298 
1299  //==================================== Interpolate for the less point
1300  c0.at(0) = c00.at(0);
1301  c0.at(1) = ( newYSample * yRelative ) + c00.at(1);
1302  c0.at(2) = c00.at(2);
1303  c0.at(3) = ( c00.at(3) * ( 1.0 - yRelative ) ) + ( c10.at(3) * yRelative );
1304 
1305  //==================================== Interpolate for the more point
1306  c1.at(0) = c01.at(0);
1307  c1.at(1) = ( newYSample * yRelative ) + c01.at(1);
1308  c1.at(2) = c01.at(2);
1309  c1.at(3) = ( c01.at(3) * ( 1.0 - yRelative ) ) + ( c11.at(3) * yRelative );
1310 
1311  //==================================== Interpolate to the new grid along Z
1312  zRelative = ( ( zIt * newZSample ) - ( zBottom * oldZSample ) ) / ( ( zTop * oldZSample ) - ( zBottom * oldZSample ) );
1313  newMap[newMapIndex] = ( c0.at(3) * ( 1.0 - zRelative ) ) + ( c1.at(3) * zRelative );
1314  }
1315  }
1316  }
1317 
1318  //================================================ Delete old map and allocate new memory
1319  delete[] map;
1320  map = new proshade_double [newXDim * newYDim * newZDim];
1321 
1322  //================================================ Copy map
1323  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ )
1324  {
1325  map[iter] = newMap[iter];
1326  }
1327 
1328  //================================================ Release memory
1329  delete[] newMap;
1330 
1331  //================================================ Define change in indices and return it
1332  corrs[0] = newXDim - xDim;
1333  corrs[1] = newYDim - yDim;
1334  corrs[2] = newZDim - zDim;
1335  corrs[3] = newXDim * newXSample;
1336  corrs[4] = newYDim * newYSample;
1337  corrs[5] = newZDim * newZSample;
1338 
1339  //======================================== Done
1340  return ;
1341 
1342 }
1343 
1359 void ProSHADE_internal_mapManip::reSampleMapToResolutionFourier ( proshade_double*& map, proshade_single resolution, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single*& corrs )
1360 {
1361  //================================================ Sanity check - the resolution needs to be set
1362  if ( resolution <= 0.0 )
1363  {
1364  throw ProSHADE_exception ( "Requested resolution not set for map re-sampling.", "EM00015", __FILE__, __LINE__, __func__, "There is no resolution value set, but map re-sampling to\n : this unset resolution value is required. This error\n : occurs when a task with no resolution requirement is\n : requested on a map data and the map resolution change is\n : set to \'on\'. Either supply a resolution value, or do not\n : re-sample the map." );
1365  }
1366 
1367  //================================================ Initialise variables
1368  proshade_unsign newXDim = static_cast<proshade_unsign> ( std::ceil ( xAngs / ( resolution / 2.0 ) ) );
1369  proshade_unsign newYDim = static_cast<proshade_unsign> ( std::ceil ( yAngs / ( resolution / 2.0 ) ) );
1370  proshade_unsign newZDim = static_cast<proshade_unsign> ( std::ceil ( zAngs / ( resolution / 2.0 ) ) );
1371 
1372  if ( newXDim % 2 != 0 ) { newXDim += 1; }
1373  if ( newYDim % 2 != 0 ) { newYDim += 1; }
1374  if ( newZDim % 2 != 0 ) { newZDim += 1; }
1375 
1376  proshade_signed preXChange, preYChange, preZChange;
1377  if ( ( xDimS % 2 ) == 0 ) { preXChange = std::ceil ( ( static_cast<proshade_signed> ( xDimS ) - static_cast<proshade_signed> ( newXDim ) ) / 2.0 ); }
1378  else { preXChange = std::floor ( ( static_cast<proshade_signed> ( xDimS ) - static_cast<proshade_signed> ( newXDim ) ) / 2.0 ); }
1379  if ( ( yDimS % 2 ) == 0 ) { preYChange = std::ceil ( ( static_cast<proshade_signed> ( yDimS ) - static_cast<proshade_signed> ( newYDim ) ) / 2.0 ); }
1380  else { preYChange = std::floor ( ( static_cast<proshade_signed> ( yDimS ) - static_cast<proshade_signed> ( newYDim ) ) / 2.0 ); }
1381  if ( ( zDimS % 2 ) == 0 ) { preZChange = std::ceil ( ( static_cast<proshade_signed> ( zDimS ) - static_cast<proshade_signed> ( newZDim ) ) / 2.0 ); }
1382  else { preZChange = std::floor ( ( static_cast<proshade_signed> ( zDimS ) - static_cast<proshade_signed> ( newZDim ) ) / 2.0 ); }
1383 
1384  proshade_signed postXChange = static_cast<proshade_signed> ( xDimS ) - ( preXChange + static_cast<proshade_signed> ( newXDim ) );
1385  proshade_signed postYChange = static_cast<proshade_signed> ( yDimS ) - ( preYChange + static_cast<proshade_signed> ( newYDim ) );
1386  proshade_signed postZChange = static_cast<proshade_signed> ( zDimS ) - ( preZChange + static_cast<proshade_signed> ( newZDim ) );
1387 
1388  proshade_signed origSizeArr = 0, newSizeArr = 0;
1389  proshade_double normFactor = static_cast<proshade_double> ( xDimS * yDimS * zDimS );
1390 
1391  //================================================ Manage memory
1392  fftw_complex *origMap, *fCoeffs, *newFCoeffs, *newMap;
1393  fftw_plan planForwardFourier, planBackwardRescaledFourier;
1394  allocateResolutionFourierMemory ( origMap, fCoeffs, newFCoeffs, newMap, planForwardFourier, planBackwardRescaledFourier,
1395  xDimS, yDimS, zDimS, newXDim, newYDim, newZDim );
1396 
1397  //================================================ Fill maps with data and zeroes
1398  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDimS * yDimS * zDimS ); iter++ ) { origMap[iter][0] = map[iter]; origMap[iter][1] = 0.0; }
1399  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ ) { newFCoeffs[iter][0] = 0.0; newFCoeffs[iter][1] = 0.0; }
1400 
1401  //================================================ Get the Fourier coeffs
1402  fftw_execute ( planForwardFourier );
1403 
1404  //================================================ Change the order of Fourier coefficients
1405  changeFourierOrder ( fCoeffs, xDimS, yDimS, zDimS, true );
1406 
1407  //================================================ Re-sample the coefficients by removing high frequencies or adding these with 0 values
1408  for ( proshade_unsign xIt = 0; xIt < newXDim; xIt++ )
1409  {
1410  for ( proshade_unsign yIt = 0; yIt < newYDim; yIt++ )
1411  {
1412  for ( proshade_unsign zIt = 0; zIt < newZDim; zIt++ )
1413  {
1414  //==================================== Find the array positions
1415  origSizeArr = (zIt + preZChange) + zDimS * ( (yIt + preYChange) + yDimS * (xIt + preXChange) );
1416  newSizeArr = zIt + newZDim * ( yIt + newYDim * xIt );
1417 
1418  //==================================== If original coefficient for this new coefficient position exists, copy
1419  if ( ( ( -1 < static_cast<proshade_signed> ( xIt + preXChange ) ) && ( -1 < static_cast<proshade_signed> ( yIt + preYChange ) ) && ( -1 < static_cast<proshade_signed> ( zIt + preZChange ) ) ) &&
1420  ( ( xIt < static_cast<proshade_unsign> ( newXDim + postXChange ) ) && ( yIt < static_cast<proshade_unsign> ( newYDim + postYChange ) ) && ( zIt < static_cast<proshade_unsign> ( newZDim + postZChange ) ) ) )
1421  {
1422  //================================ Copy the Fourier coeff
1423  newFCoeffs[newSizeArr][0] = fCoeffs[origSizeArr][0] / normFactor;
1424  newFCoeffs[newSizeArr][1] = fCoeffs[origSizeArr][1] / normFactor;
1425  }
1426  }
1427  }
1428  }
1429 
1430  //================================================ Change the order of the re-sampled Fourier coefficients
1431  changeFourierOrder ( newFCoeffs, newXDim, newYDim, newZDim, false );
1432 
1433  //================================================ Get the new map from the re-sized Fourier coefficients
1434  fftw_execute ( planBackwardRescaledFourier );
1435 
1436  //================================================ Delete the old map and create a new, re-sized one. Then copy the new map values into this new map memory.
1437  delete map;
1438  map = new proshade_double [newXDim * newYDim * newZDim];
1439  ProSHADE_internal_misc::checkMemoryAllocation ( map, __FILE__, __LINE__, __func__ );
1440  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ ) { map[iter] = newMap[iter][0]; }
1441 
1442  //================================================ Release memory
1443  releaseResolutionFourierMemory ( origMap, fCoeffs, newFCoeffs, newMap, planForwardFourier, planBackwardRescaledFourier );
1444 
1445  //================================================ Define change in indices and return it
1446  corrs[0] = static_cast<proshade_signed> ( newXDim ) - static_cast<proshade_signed> ( xDimS );
1447  corrs[1] = static_cast<proshade_signed> ( newYDim ) - static_cast<proshade_signed> ( yDimS );
1448  corrs[2] = static_cast<proshade_signed> ( newZDim ) - static_cast<proshade_signed> ( zDimS );
1449  corrs[3] = static_cast<proshade_signed> ( newXDim ) * ( resolution / 2.0 );
1450  corrs[4] = static_cast<proshade_signed> ( newYDim ) * ( resolution / 2.0 );
1451  corrs[5] = static_cast<proshade_signed> ( newZDim ) * ( resolution / 2.0 );
1452 
1453  //======================================== Done
1454  return ;
1455 
1456 }
1457 
1476 void ProSHADE_internal_mapManip::allocateResolutionFourierMemory ( fftw_complex*& origMap, fftw_complex*& fCoeffs, fftw_complex*& newFCoeffs, fftw_complex*& newMap, fftw_plan& planForwardFourier, fftw_plan& planBackwardRescaledFourier, proshade_unsign xDimOld, proshade_unsign yDimOld, proshade_unsign zDimOld, proshade_unsign xDimNew, proshade_unsign yDimNew, proshade_unsign zDimNew )
1477 {
1478  //================================================ Initialise memory
1479  origMap = new fftw_complex [xDimOld * yDimOld * zDimOld];
1480  fCoeffs = new fftw_complex [xDimOld * yDimOld * zDimOld];
1481  newFCoeffs = new fftw_complex [xDimNew * yDimNew * zDimNew];
1482  newMap = new fftw_complex [xDimNew * yDimNew * zDimNew];
1483 
1484  //================================================ Check memory allocation
1485  ProSHADE_internal_misc::checkMemoryAllocation ( origMap, __FILE__, __LINE__, __func__ );
1486  ProSHADE_internal_misc::checkMemoryAllocation ( fCoeffs, __FILE__, __LINE__, __func__ );
1487  ProSHADE_internal_misc::checkMemoryAllocation ( newFCoeffs, __FILE__, __LINE__, __func__ );
1488  ProSHADE_internal_misc::checkMemoryAllocation ( newMap, __FILE__, __LINE__, __func__ );
1489 
1490  //================================================ Create plans
1491  planForwardFourier = fftw_plan_dft_3d ( xDimOld, yDimOld, zDimOld, origMap, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
1492  planBackwardRescaledFourier = fftw_plan_dft_3d ( xDimNew, yDimNew, zDimNew, newFCoeffs, newMap, FFTW_BACKWARD, FFTW_ESTIMATE );
1493 
1494  //================================================ Done
1495  return ;
1496 
1497 }
1498 
1510 void ProSHADE_internal_mapManip::releaseResolutionFourierMemory ( fftw_complex*& origMap, fftw_complex*& fCoeffs, fftw_complex*& newFCoeffs, fftw_complex*& newMap, fftw_plan& planForwardFourier, fftw_plan& planBackwardRescaledFourier )
1511 {
1512  //================================================ Delete the FFTW plans
1513  fftw_destroy_plan ( planForwardFourier );
1514  fftw_destroy_plan ( planBackwardRescaledFourier );
1515 
1516  //================================================ Delete the complex arrays
1517  delete[] origMap;
1518  delete[] fCoeffs;
1519  delete[] newFCoeffs;
1520  delete[] newMap;
1521 
1522  //================================================ Done
1523  return ;
1524 
1525 }
1526 
1539 void ProSHADE_internal_mapManip::changeFourierOrder ( fftw_complex*& fCoeffs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, bool negativeFirst )
1540 {
1541  //================================================ Initialise local variables
1542  proshade_signed h = 0, k = 0, l = 0, origSizeArr = 0, newSizeArr = 0;
1543  proshade_signed xSeq1FreqStart, ySeq1FreqStart, zSeq1FreqStart, xSeq2FreqStart, ySeq2FreqStart, zSeq2FreqStart;
1544 
1545  //================================================ Find the positive and negative indices cot-offs
1546  if ( negativeFirst )
1547  {
1548  if ( ( xDim % 2 ) == 0 ) { xSeq1FreqStart = xDim / 2; xSeq2FreqStart = xDim / 2; } else { xSeq1FreqStart = (xDim / 2) + 1; xSeq2FreqStart = xDim / 2; }
1549  if ( ( yDim % 2 ) == 0 ) { ySeq1FreqStart = yDim / 2; ySeq2FreqStart = yDim / 2; } else { ySeq1FreqStart = (yDim / 2) + 1; ySeq2FreqStart = yDim / 2; }
1550  if ( ( zDim % 2 ) == 0 ) { zSeq1FreqStart = zDim / 2; zSeq2FreqStart = zDim / 2; } else { zSeq1FreqStart = (zDim / 2) + 1; zSeq2FreqStart = zDim / 2; }
1551  }
1552  else
1553  {
1554  if ( ( xDim % 2 ) == 0 ) { xSeq1FreqStart = xDim / 2; xSeq2FreqStart = xDim / 2; } else { xSeq1FreqStart = (xDim / 2); xSeq2FreqStart = xDim / 2 + 1; }
1555  if ( ( yDim % 2 ) == 0 ) { ySeq1FreqStart = yDim / 2; ySeq2FreqStart = yDim / 2; } else { ySeq1FreqStart = (yDim / 2); ySeq2FreqStart = yDim / 2 + 1; }
1556  if ( ( zDim % 2 ) == 0 ) { zSeq1FreqStart = zDim / 2; zSeq2FreqStart = zDim / 2; } else { zSeq1FreqStart = (zDim / 2); zSeq2FreqStart = zDim / 2 + 1; }
1557  }
1558 
1559  //================================================ Allocate helper array memory
1560  fftw_complex *hlpFCoeffs = new fftw_complex [xDim * yDim * zDim];
1561  ProSHADE_internal_misc::checkMemoryAllocation ( hlpFCoeffs, __FILE__, __LINE__, __func__ );
1562 
1563  //================================================ Change the coefficients order
1564  for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1565  {
1566  //============================================ Find x frequency
1567  if ( xIt < xSeq1FreqStart ) { h = xIt + xSeq2FreqStart; } else { h = xIt - xSeq1FreqStart; }
1568  for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1569  {
1570  //======================================== Find y frequency
1571  if ( yIt < ySeq1FreqStart ) { k = yIt + ySeq2FreqStart; } else { k = yIt - ySeq1FreqStart; }
1572 
1573  for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1574  {
1575  //==================================== Find z frequency
1576  if ( zIt < zSeq1FreqStart ) { l = zIt + zSeq2FreqStart; } else { l = zIt - zSeq1FreqStart; }
1577 
1578  //==================================== Find array positions
1579  newSizeArr = l + zDim * ( k + yDim * h );
1580  origSizeArr = zIt + zDim * ( yIt + yDim * xIt );
1581 
1582  //==================================== Copy vals
1583  hlpFCoeffs[newSizeArr][0] = fCoeffs[origSizeArr][0];
1584  hlpFCoeffs[newSizeArr][1] = fCoeffs[origSizeArr][1];
1585  }
1586  }
1587  }
1588 
1589  //================================================ Copy the helper array to the input Fourier coefficients array
1590  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { fCoeffs[iter][0] = hlpFCoeffs[iter][0]; fCoeffs[iter][1] = hlpFCoeffs[iter][1]; }
1591 
1592  //================================================ Release helper array memory
1593  delete[] hlpFCoeffs;
1594 
1595  //================================================ Done
1596  return ;
1597 
1598 }
1599 
1611 void ProSHADE_internal_mapManip::removeMapPhase ( fftw_complex*& mapCoeffs, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim )
1612 {
1613  //================================================ Set local variables
1614  proshade_double real, imag, mag, phase;
1615  proshade_unsign arrayPos = 0;
1616  proshade_double normFactor = static_cast<proshade_double> ( xDim * yDim * zDim );
1617 
1618  //================================================ Iterate through the map
1619  for ( proshade_unsign uIt = 0; uIt < xDim; uIt++ )
1620  {
1621  for ( proshade_unsign vIt = 0; vIt < yDim; vIt++ )
1622  {
1623  for ( proshade_unsign wIt = 0; wIt < zDim; wIt++ )
1624  {
1625  //==================================== Var init
1626  arrayPos = wIt + zDim * ( vIt + yDim * uIt );
1627  real = mapCoeffs[arrayPos][0];
1628  imag = mapCoeffs[arrayPos][1];
1629 
1630  //==================================== Get magnitude and phase with mask parameters
1631  mag = std::sqrt ( (real*real) + (imag*imag) );;
1632  phase = 0.0; // This would be std::atan2 ( imag, real ); - but here we remove the phase.
1633 
1634  //==================================== Save the phaseless data
1635  mapCoeffs[arrayPos][0] = ( mag * cos(phase) ) / normFactor;
1636  mapCoeffs[arrayPos][1] = ( mag * sin(phase) ) / normFactor;
1637  }
1638  }
1639  }
1640 
1641  //================================================ Done
1642  return ;
1643 
1644 }
1645 
1660 void ProSHADE_internal_mapManip::getFakeHalfMap ( proshade_double*& map, proshade_double*& fakeHalfMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_signed fakeMapKernel )
1661 {
1662  //================================================ Set local variables
1663  proshade_signed xDim = xDimS;
1664  proshade_signed yDim = yDimS;
1665  proshade_signed zDim = zDimS;
1666  proshade_signed currentPos, neighArrPos, neighXPos, neighYPos, neighZPos;
1667  proshade_double neighSum;
1668  proshade_double neighCount = pow ( ( ( fakeMapKernel * 2 ) + 1 ), 3.0 ) - 1.0;
1669 
1670  //================================================ Blur the coeffs
1671  for ( proshade_signed uIt = 0; uIt < xDim; uIt++ )
1672  {
1673  for ( proshade_signed vIt = 0; vIt < yDim; vIt++ )
1674  {
1675  for ( proshade_signed wIt = 0; wIt < zDim; wIt++ )
1676  {
1677  //==================================== Var init
1678  currentPos = wIt + zDim * ( vIt + yDim * uIt );
1679  neighSum = 0.0;
1680 
1681  //==================================== Average neighbours
1682  for ( proshade_signed xCh = -fakeMapKernel; xCh <= +fakeMapKernel; xCh++ )
1683  {
1684  for ( proshade_signed yCh = -fakeMapKernel; yCh <= +fakeMapKernel; yCh++ )
1685  {
1686  for ( proshade_signed zCh = -fakeMapKernel; zCh <= +fakeMapKernel; zCh++ )
1687  {
1688  if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) { continue; }
1689 
1690  //======================== Find the nieghbour peak indices (with periodicity)
1691  neighXPos = uIt + xCh; if ( neighXPos >= xDim ) { neighXPos -= xDim; }; if ( neighXPos < 0 ) { neighXPos += xDim; }
1692  neighYPos = vIt + yCh; if ( neighYPos >= yDim ) { neighYPos -= yDim; }; if ( neighYPos < 0 ) { neighYPos += yDim; }
1693  neighZPos = wIt + zCh; if ( neighZPos >= zDim ) { neighZPos -= zDim; }; if ( neighZPos < 0 ) { neighZPos += zDim; }
1694  neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1695 
1696  //======================== Add to average
1697  neighSum += map[neighArrPos];
1698  }
1699  }
1700  }
1701 
1702  //==================================== Save the average to "fake" half-map
1703  fakeHalfMap[currentPos] = neighSum / neighCount;
1704  }
1705  }
1706  }
1707 
1708  //================================================ Done
1709  return ;
1710 
1711 }
1712 
1725 void ProSHADE_internal_mapManip::getCorrelationMapMask ( proshade_double*& map, proshade_double*& fakeHalfMap, proshade_double*& correlationMask, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_signed corrMaskKernel )
1726 {
1727  //================================================ Set local variables
1728  proshade_signed xDim = xDimS, yDim = yDimS, zDim = zDimS, currentPos, neighArrPos, neighXPos, neighYPos, neighZPos, corrIter;
1729  proshade_unsign noCorrVals = static_cast<proshade_unsign> ( pow ( ( ( corrMaskKernel * 2 ) + 1 ), 3 ) );
1730 
1731  //================================================ Alocate memory
1732  proshade_double *origMap = new proshade_double [noCorrVals];
1733  proshade_double *fakeHM = new proshade_double [noCorrVals];
1734 
1735  //================================================ Check memory allocation
1736  ProSHADE_internal_misc::checkMemoryAllocation ( origMap, __FILE__, __LINE__, __func__ );
1737  ProSHADE_internal_misc::checkMemoryAllocation ( fakeHM, __FILE__, __LINE__, __func__ );
1738 
1739  //================================================ Blur the coeffs
1740  for ( proshade_signed uIt = 0; uIt < xDim; uIt++ )
1741  {
1742  for ( proshade_signed vIt = 0; vIt < yDim; vIt++ )
1743  {
1744  for ( proshade_signed wIt = 0; wIt < zDim; wIt++ )
1745  {
1746  //==================================== Var init
1747  currentPos = wIt + zDim * ( vIt + yDim * uIt );
1748  corrIter = 0;
1749 
1750  //==================================== Average neighbours
1751  for ( proshade_signed xCh = -corrMaskKernel; xCh <= +corrMaskKernel; xCh++ )
1752  {
1753  for ( proshade_signed yCh = -corrMaskKernel; yCh <= +corrMaskKernel; yCh++ )
1754  {
1755  for ( proshade_signed zCh = -corrMaskKernel; zCh <= +corrMaskKernel; zCh++ )
1756  {
1757  //======================== Find the nieghbour peak indices (with periodicity)
1758  neighXPos = uIt + xCh; if ( neighXPos >= xDim ) { neighXPos -= xDim; }; if ( neighXPos < 0 ) { neighXPos += xDim; }
1759  neighYPos = vIt + yCh; if ( neighYPos >= yDim ) { neighYPos -= yDim; }; if ( neighYPos < 0 ) { neighYPos += yDim; }
1760  neighZPos = wIt + zCh; if ( neighZPos >= zDim ) { neighZPos -= zDim; }; if ( neighZPos < 0 ) { neighZPos += zDim; }
1761  neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1762 
1763  //======================== Add to correct arrays
1764  origMap[corrIter] = map[neighArrPos];
1765  fakeHM[corrIter] = fakeHalfMap[neighArrPos];
1766  corrIter += 1;
1767  }
1768  }
1769  }
1770 
1771  //==================================== Save the correlation comparison result
1772  correlationMask[currentPos] = ProSHADE_internal_maths::pearsonCorrCoeff ( origMap, fakeHM, noCorrVals );
1773  }
1774  }
1775  }
1776 
1777  //================================================ Done
1778  return ;
1779 
1780 }
1781 
1796 proshade_single ProSHADE_internal_mapManip::getIndicesFromAngstroms ( proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single dist )
1797 {
1798  //================================================ Compute
1799  proshade_unsign ret = ProSHADE_internal_mapManip::myRound ( std::max ( dist / ( xAngs / static_cast<proshade_single> ( xDim ) ),
1800  std::max ( dist / ( yAngs / static_cast<proshade_single> ( yDim ) ),
1801  dist / ( zAngs / static_cast<proshade_single> ( zDim ) ) ) ) );
1802 
1803  //================================================ Done
1804  return ( ret );
1805 
1806 }
1807 
1819 void ProSHADE_internal_mapManip::connectMaskBlobs ( proshade_double*& mask, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single maskThres )
1820 {
1821  //================================================ Initialise variables
1822  proshade_double* hlpMap = new proshade_double[xDim * yDim * zDim];
1823  proshade_signed addSurroundingPoints = std::max ( 3L, static_cast<proshade_signed> ( std::ceil ( getIndicesFromAngstroms( xDim, yDim, zDim, xAngs, yAngs, zAngs, std::max( xAngs, std::max( yAngs, zAngs ) ) * 0.1 ) ) ) );
1824  proshade_signed currPos, neighXPos, neighYPos, neighZPos, neighArrPos;
1825 
1826  //================================================ Check memory allocation
1827  ProSHADE_internal_misc::checkMemoryAllocation ( hlpMap, __FILE__, __LINE__, __func__ );
1828 
1829  //================================================ Copy the mask map
1830  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { hlpMap[iter] = mask[iter]; }
1831 
1832  //================================================ Repeat as many times as needed
1833  for ( proshade_signed it = 0; it < addSurroundingPoints; it++ )
1834  {
1835  //============================================ For each x, y and z
1836  for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1837  {
1838  for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1839  {
1840  for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1841  {
1842  //================================ Current position
1843  currPos = zIt + zDim * ( yIt + yDim * xIt );
1844 
1845  //================================ If zero, ignore
1846  if ( hlpMap[currPos] < maskThres ) { continue; }
1847 
1848  //================================ Check neighbours
1849  for ( proshade_signed xCh = -1; xCh <= +1; xCh++ )
1850  {
1851  for ( proshade_signed yCh = -1; yCh <= +1; yCh++ )
1852  {
1853  for ( proshade_signed zCh = -1; zCh <= +1; zCh++ )
1854  {
1855  if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) { continue; }
1856 
1857  //==================== Find the nieghbour indices (without periodicity)
1858  neighXPos = xIt + xCh; if ( neighXPos < 0 ) { continue; } if ( neighXPos >= xDim ) { continue; }
1859  neighYPos = yIt + yCh; if ( neighYPos < 0 ) { continue; } if ( neighYPos >= yDim ) { continue; }
1860  neighZPos = zIt + zCh; if ( neighZPos < 0 ) { continue; } if ( neighZPos >= zDim ) { continue; }
1861  neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1862 
1863  //==================== Add to mask if this point is below it (as it is a neighbour to a point which is part of the mask)
1864  if ( hlpMap[neighArrPos] < maskThres ) { mask[neighArrPos] = maskThres; }
1865  }
1866  }
1867  }
1868  }
1869  }
1870  }
1871 
1872  //============================================ Now copy the updated mask to the temporary helper, unless last iteration
1873  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { hlpMap[iter] = mask[iter]; }
1874  }
1875 
1876  //================================================ Release memory
1877  delete[] hlpMap;
1878 
1879  //================================================ Done
1880  return ;
1881 
1882 }
1883 
1897 void ProSHADE_internal_mapManip::beautifyBoundaries ( proshade_signed*& bounds, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_signed boundsDiffThres, proshade_signed verbose )
1898 {
1899  //================================================ If extra bounds space pushed the bounds over the physical map, freely add up to 10 indices over the current value
1900  while ( bounds[1] >= static_cast<proshade_signed> ( xDim ) ) { xDim += 10; }
1901  while ( bounds[3] >= static_cast<proshade_signed> ( yDim ) ) { yDim += 10; }
1902  while ( bounds[5] >= static_cast<proshade_signed> ( zDim ) ) { zDim += 10; }
1903 
1904  //================================================ Find if better bouds exist in terms of prime numbers
1905  proshade_signed addToX = betterClosePrimeFactors ( bounds[1] - bounds[0] + 1, xDim );
1906  proshade_signed addToY = betterClosePrimeFactors ( bounds[3] - bounds[2] + 1, yDim );
1907  proshade_signed addToZ = betterClosePrimeFactors ( bounds[5] - bounds[4] + 1, zDim );
1908 
1909  //================================================ Figure if similar sizes should not be forced to be the same
1910  proshade_signed XtoY = std::abs ( addToX - addToY );
1911  proshade_signed XtoZ = std::abs ( addToX - addToZ );
1912  proshade_signed YtoZ = std::abs ( addToY - addToZ );
1913 
1914  if ( ( ( XtoY < boundsDiffThres ) && ( XtoZ < boundsDiffThres ) ) ||
1915  ( ( XtoY < boundsDiffThres ) && ( YtoZ < boundsDiffThres ) ) ||
1916  ( ( XtoZ < boundsDiffThres ) && ( YtoZ < boundsDiffThres ) ) )
1917  {
1918  //============================================ Dealing with chain ( a <= b <= c type ) where at least two dista are smaller than threshold
1919  proshade_signed maxSize = std::max ( addToX, std::max ( addToY, addToZ ) );
1920  addToX = maxSize;
1921  addToY = maxSize;
1922  addToZ = maxSize;
1923  }
1924  else
1925  {
1926  //============================================ Only two may be similar enough
1927  if ( XtoY <= boundsDiffThres )
1928  {
1929  proshade_signed maxSize = std::max ( addToX, addToY );
1930  addToX = maxSize;
1931  addToY = maxSize;
1932  }
1933  if ( XtoZ <= boundsDiffThres )
1934  {
1935  proshade_signed maxSize = std::max ( addToX, addToZ );
1936  addToX = maxSize;
1937  addToZ = maxSize;
1938  }
1939  if ( YtoZ <= boundsDiffThres )
1940  {
1941  proshade_signed maxSize = std::max ( addToY, addToZ );
1942  addToY = maxSize;
1943  addToZ = maxSize;
1944  }
1945  }
1946 
1947  //================================================ Add the extra space appropriately
1948  distributeSpaceToBoundaries ( bounds[0], bounds[1], ( bounds[1] - bounds[0] + 1 ), addToX );
1949  distributeSpaceToBoundaries ( bounds[2], bounds[3], ( bounds[3] - bounds[2] + 1 ), addToY );
1950  distributeSpaceToBoundaries ( bounds[4], bounds[5], ( bounds[5] - bounds[4] + 1 ), addToZ );
1951 
1952  //================================================ Done
1953  return ;
1954 
1955 }
1956 
1967 proshade_signed ProSHADE_internal_mapManip::betterClosePrimeFactors ( proshade_signed fromRange, proshade_signed toRange )
1968 {
1969  //================================================ Initialise variables
1970  proshade_signed ret = fromRange;
1971  std::vector < proshade_signed > posibles, hlp;
1972  proshade_signed sum;
1973 
1974  //================================================ Find the sums of prime factors for each number in the whole range
1975  for ( proshade_signed iter = fromRange; iter < toRange; iter++ )
1976  {
1977  sum = 0;
1979  for ( proshade_unsign i = 0; i < static_cast<proshade_unsign> ( hlp.size() ); i++ ) { sum += hlp.at(i); }
1980  hlp.clear ( );
1981  ProSHADE_internal_misc::addToSignedVector ( &posibles, sum );
1982  }
1983 
1984  //================================================ Find a better number
1985  for ( proshade_signed iter = fromRange; iter < toRange; iter++ )
1986  {
1987  //============================================ Ignore odd numbers
1988  if ( iter %2 != 0 ) { continue; }
1989 
1990  //============================================ Better number needs to have lower prime factor sum, but also needs not to be too far
1991  if ( posibles.at(iter-fromRange) < ( posibles.at(ret-fromRange) - ( iter - ret ) ) ) { ret = iter; }
1992  }
1993 
1994  //================================================ In the case of large prime number, just add one for even number
1995  if ( ( ret % 2 != 0 ) && ( ret < ( toRange - 1 ) ) ) { ret += 1; }
1996 
1997  //================================================ Done
1998  return ( ret );
1999 
2000 }
2001 
2013 void ProSHADE_internal_mapManip::distributeSpaceToBoundaries ( proshade_signed& minBound, proshade_signed& maxBound, proshade_signed oldBoundRange, proshade_signed newBoundRange )
2014 {
2015  //================================================ Only bother when sometings is to be added
2016  if ( newBoundRange > oldBoundRange )
2017  {
2018  //============================================ How much are we adding?
2019  proshade_signed distributeThis = newBoundRange - oldBoundRange;
2020 
2021  while ( distributeThis != 0 )
2022  {
2023  minBound -= 1;
2024  distributeThis -= 1;
2025 
2026  if ( distributeThis != 0 )
2027  {
2028  maxBound += 1;
2029  distributeThis -= 1;
2030  }
2031  }
2032  }
2033 
2034  //================================================ Done
2035  return ;
2036 
2037 }
2038 
2063 void ProSHADE_internal_mapManip::copyMapByBounds ( proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo, proshade_signed origXFrom, proshade_signed origYFrom, proshade_signed origZFrom, proshade_signed yDimIndices, proshade_signed zDimIndices, proshade_signed origXDimIndices, proshade_signed origYDimIndices, proshade_signed origZDimIndices, proshade_double*& newMap, proshade_double* origMap )
2064 {
2065  //================================================ Initialise variables
2066  proshade_signed newMapIndex, oldMapIndex, oldX, oldY, oldZ, newX, newY, newZ;
2067 
2068  //=============================================== For all values in the new map
2069  for ( proshade_signed xIt = xFrom; xIt <= xTo; xIt++ )
2070  {
2071  //============================================ Index position init
2072  newX = ( xIt - xFrom );
2073  oldX = ( newX + ( xFrom - origXFrom ) );
2074 
2075  for ( proshade_signed yIt = yFrom; yIt <= yTo; yIt++ )
2076  {
2077  //======================================== Index position init
2078  newY = ( yIt - yFrom );
2079  oldY = ( newY + ( yFrom - origYFrom ) );
2080 
2081  for ( proshade_signed zIt = zFrom; zIt <= zTo; zIt++ )
2082  {
2083  //==================================== Index position init
2084  newZ = ( zIt - zFrom );
2085  oldZ = ( newZ + ( zFrom - origZFrom ) );
2086 
2087  //==================================== Index arrays
2088  newMapIndex = newZ + zDimIndices * ( newY + yDimIndices * newX );
2089  oldMapIndex = oldZ + origZDimIndices * ( oldY + origYDimIndices * oldX );
2090 
2091  //============================ Check if we are adding new point
2092  if ( ( ( oldX < 0 ) || ( oldX >= origXDimIndices ) ) ||
2093  ( ( oldY < 0 ) || ( oldY >= origYDimIndices ) ) ||
2094  ( ( oldZ < 0 ) || ( oldZ >= origZDimIndices ) ) )
2095  {
2096  //================================ Yes, this is a new point, no known value for it
2097  newMap[newMapIndex] = 0.0;
2098  continue;
2099  }
2100 
2101  //==================================== If we got here, this should be within the old map, so just copy value
2102  newMap[newMapIndex] = origMap[oldMapIndex];
2103  }
2104  }
2105  }
2106 
2107  //================================================ Done
2108  return ;
2109 
2110 }
ProSHADE_internal_mapManip::addExtraBoundSpace
void addExtraBoundSpace(proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed *&bounds, proshade_single extraSpace)
This function takes a set of bounds and adds a given number of Angstroms to them.
Definition: ProSHADE_mapManip.cpp:1122
ProSHADE_internal_mapManip::determinePDBRanges
void determinePDBRanges(gemmi::Structure pdbFile, proshade_single *xFrom, proshade_single *xTo, proshade_single *yFrom, proshade_single *yTo, proshade_single *zFrom, proshade_single *zTo, bool firstModel)
Function for finding the PDB file ranges.
Definition: ProSHADE_mapManip.cpp:57
ProSHADE_internal_mapManip::changePDBBFactors
void changePDBBFactors(gemmi::Structure *pdbFile, proshade_double newBFactorValue, bool firstModel)
Function for changing the PDB B-factor values to a specific single value.
Definition: ProSHADE_mapManip.cpp:425
ProSHADE_exception
This class is the representation of ProSHADE exception.
Definition: ProSHADE_exceptions.hpp:37
ProSHADE_internal_mapManip::distributeSpaceToBoundaries
void distributeSpaceToBoundaries(proshade_signed &minBound, proshade_signed &maxBound, proshade_signed oldBoundRange, proshade_signed newBoundRange)
Function for adding space to boundaries within the map confines.
Definition: ProSHADE_mapManip.cpp:2013
ProSHADE_internal_mapManip::getMaskFromBlurr
void getMaskFromBlurr(proshade_double *&blurMap, proshade_double *&outMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single noIQRs)
Function for computing mask from blurred map.
Definition: ProSHADE_mapManip.cpp:1012
ProSHADE_mapManip.hpp
This header file declares the ProSHADE_internal_mapManip namespace, which groups functions for intern...
ProSHADE_internal_mapManip::findMAPCOMValues
void findMAPCOMValues(proshade_double *map, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo)
This function finds the Centre of Mass for a map.
Definition: ProSHADE_mapManip.cpp:225
ProSHADE_internal_mapManip::blurSharpenMap
void blurSharpenMap(proshade_double *&map, proshade_double *&maskedMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single blurringFactor)
Function for blurring/sharpening maps.
Definition: ProSHADE_mapManip.cpp:912
ProSHADE_internal_mapManip::moveMapByFourier
void moveMapByFourier(proshade_double *&map, proshade_single xMov, proshade_single yMov, proshade_single zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim)
Function for moving map back to original PDB location by using Fourier transformation.
Definition: ProSHADE_mapManip.cpp:792
ProSHADE_internal_mapManip::movePDBForMapCalc
void movePDBForMapCalc(gemmi::Structure *pdbFile, proshade_single xMov, proshade_single yMov, proshade_single zMov, bool firstModel)
Function for moving co-ordinate atoms to better suit theoretical map computation.
Definition: ProSHADE_mapManip.cpp:554
ProSHADE_internal_maths::getRotationMatrixFromEulerZXZAngles
void getRotationMatrixFromEulerZXZAngles(proshade_double eulerAlpha, proshade_double eulerBeta, proshade_double eulerGamma, proshade_double *matrix)
Function to find the rotation matrix from Euler angles (ZXZ convention).
Definition: ProSHADE_maths.cpp:1005
ProSHADE_internal_maths::complexMultiplication
void complexMultiplication(proshade_double *r1, proshade_double *i1, proshade_double *r2, proshade_double *i2, proshade_double *retReal, proshade_double *retImag)
Function to multiply two complex numbers.
Definition: ProSHADE_maths.cpp:38
ProSHADE_internal_messages::printWarningMessage
void printWarningMessage(proshade_signed verbose, std::string message, std::string warnCode)
General stderr message printing (used for warnings).
Definition: ProSHADE_messages.cpp:101
ProSHADE_internal_mapManip::reSampleMapToResolutionTrilinear
void reSampleMapToResolutionTrilinear(proshade_double *&map, proshade_single resolution, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single *&corrs)
This function re-samples a map to conform to given resolution using tri-linear interpolation.
Definition: ProSHADE_mapManip.cpp:1158
ProSHADE_internal_mapManip::beautifyBoundaries
void beautifyBoundaries(proshade_signed *&bounds, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_signed boundsDiffThres, proshade_signed verbose)
Function for modifying boundaries to a mathematically more pleasant values.
Definition: ProSHADE_mapManip.cpp:1897
ProSHADE_internal_mapManip::getCorrelationMapMask
void getCorrelationMapMask(proshade_double *&map, proshade_double *&fakeHalfMap, proshade_double *&correlationMask, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_signed corrMaskKernel)
Function for creating the correlation mask.
Definition: ProSHADE_mapManip.cpp:1725
ProSHADE_internal_mapManip::changeFourierOrder
void changeFourierOrder(fftw_complex *&fCoeffs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, bool negativeFirst)
This function changes the order of Fourier coefficients in a 3D array between positive first (default...
Definition: ProSHADE_mapManip.cpp:1539
ProSHADE_internal_mapManip::releaseResolutionFourierMemory
void releaseResolutionFourierMemory(fftw_complex *&origMap, fftw_complex *&fCoeffs, fftw_complex *&newFCoeffs, fftw_complex *&newMap, fftw_plan &planForwardFourier, fftw_plan &planBackwardRescaledFourier)
This function releases the memory required by the Fourier resampling.
Definition: ProSHADE_mapManip.cpp:1510
ProSHADE_internal_mapManip::removeWaters
void removeWaters(gemmi::Structure *pdbFile, bool firstModel)
This function removed all waters from PDB input file.
Definition: ProSHADE_mapManip.cpp:485
ProSHADE_internal_maths::pearsonCorrCoeff
proshade_double pearsonCorrCoeff(proshade_double *valSet1, proshade_double *valSet2, proshade_unsign length)
Function for computing the Pearson's correlation coefficient.
Definition: ProSHADE_maths.cpp:244
ProSHADE_internal_mapManip::findPDBCOMValues
void findPDBCOMValues(gemmi::Structure pdbFile, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, bool firstModel)
This function finds the Centre of Mass for the co-ordinate file.
Definition: ProSHADE_mapManip.cpp:142
ProSHADE_internal_misc::addToSignedVector
void addToSignedVector(std::vector< proshade_signed > *vecToAddTo, proshade_signed elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:121
ProSHADE_internal_mapManip::getIndicesFromAngstroms
proshade_single getIndicesFromAngstroms(proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single dist)
This function converts distance in Angstroms to distance in map indices.
Definition: ProSHADE_mapManip.cpp:1796
ProSHADE_internal_mapManip::getNonZeroBounds
void getNonZeroBounds(proshade_double *map, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed *&ret)
Function for finding the map boundaries enclosing positive only values.
Definition: ProSHADE_mapManip.cpp:1065
ProSHADE_internal_mapManip::copyMapByBounds
void copyMapByBounds(proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo, proshade_signed origXFrom, proshade_signed origYFrom, proshade_signed origZFrom, proshade_signed yDimIndices, proshade_signed zDimIndices, proshade_signed origXDimIndices, proshade_signed origYDimIndices, proshade_signed origZDimIndices, proshade_double *&newMap, proshade_double *origMap)
This function copies an old map to a new map with different boundaries.
Definition: ProSHADE_mapManip.cpp:2063
ProSHADE_internal_mapManip::removeMapPhase
void removeMapPhase(fftw_complex *&mapCoeffs, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim)
This function removes the phase from reciprocal (frequency) map.
Definition: ProSHADE_mapManip.cpp:1611
ProSHADE_internal_mapManip::allocateResolutionFourierMemory
void allocateResolutionFourierMemory(fftw_complex *&origMap, fftw_complex *&fCoeffs, fftw_complex *&newFCoeffs, fftw_complex *&newMap, fftw_plan &planForwardFourier, fftw_plan &planBackwardRescaledFourier, proshade_unsign xDimOld, proshade_unsign yDimOld, proshade_unsign zDimOld, proshade_unsign xDimNew, proshade_unsign yDimNew, proshade_unsign zDimNew)
This function allocates and checks the allocatio of the memory required by the Fourier resampling.
Definition: ProSHADE_mapManip.cpp:1476
ProSHADE_internal_maths::primeFactorsDecomp
std::vector< proshade_signed > primeFactorsDecomp(proshade_signed number)
Function to find prime factors of an integer.
Definition: ProSHADE_maths.cpp:1639
ProSHADE_internal_mapManip::betterClosePrimeFactors
proshade_signed betterClosePrimeFactors(proshade_signed fromRange, proshade_signed toRange)
Function for finding close numbers with better prime factors.
Definition: ProSHADE_mapManip.cpp:1967
ProSHADE_internal_mapManip::myRound
proshade_signed myRound(proshade_double x)
Calls the appropriate version of round function depending on compiler version.
Definition: ProSHADE_mapManip.cpp:31
ProSHADE_internal_mapManip::rotatePDBCoordinates
void rotatePDBCoordinates(gemmi::Structure *pdbFile, proshade_double euA, proshade_double euB, proshade_double euG, proshade_double xCom, proshade_double yCom, proshade_double zCom, bool firstModel)
Function for rotating the PDB file co-ordinates by Euler angles.
Definition: ProSHADE_mapManip.cpp:277
ProSHADE_internal_mapManip::generateMapFromPDB
void generateMapFromPDB(gemmi::Structure pdbFile, proshade_double *&map, proshade_single requestedResolution, proshade_single xCell, proshade_single yCell, proshade_single zCell, proshade_signed *xTo, proshade_signed *yTo, proshade_signed *zTo, bool forceP1, bool firstModel)
This function generates a theoretical map from co-ordinate input files.
Definition: ProSHADE_mapManip.cpp:625
ProSHADE_internal_misc::checkMemoryAllocation
void checkMemoryAllocation(chVar checkVar, std::string fileP, unsigned int lineP, std::string funcP, std::string infoP="This error may occurs when ProSHADE requests memory to be\n : allocated to it and this operation fails. This could\n : happen when not enough memory is available, either due to\n : other processes using a lot of memory, or when the machine\n : does not have sufficient memory available. Re-run to see\n : if this problem persists.")
Checks if memory was allocated properly.
Definition: ProSHADE_misc.hpp:65
ProSHADE_internal_mapManip::translatePDBCoordinates
void translatePDBCoordinates(gemmi::Structure *pdbFile, proshade_double transX, proshade_double transY, proshade_double transZ, bool firstModel)
Function for translating the PDB file co-ordinates by given distances in Angstroms.
Definition: ProSHADE_mapManip.cpp:362
ProSHADE_internal_maths::vectorMedianAndIQR
void vectorMedianAndIQR(std::vector< proshade_double > *vec, proshade_double *&ret)
Function to get vector median and inter-quartile range.
Definition: ProSHADE_maths.cpp:147
ProSHADE_internal_mapManip::getFakeHalfMap
void getFakeHalfMap(proshade_double *&map, proshade_double *&fakeHalfMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_signed fakeMapKernel)
Function for creating "fake" half-maps.
Definition: ProSHADE_mapManip.cpp:1660
ProSHADE_internal_misc::addToUnsignVector
void addToUnsignVector(std::vector< proshade_unsign > *vecToAddTo, proshade_unsign elementToAdd)
Adds the element to the vector.
Definition: ProSHADE_misc.cpp:99
ProSHADE_internal_mapManip::reSampleMapToResolutionFourier
void reSampleMapToResolutionFourier(proshade_double *&map, proshade_single resolution, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single *&corrs)
This function re-samples a map to conform to given resolution using Fourier.
Definition: ProSHADE_mapManip.cpp:1359
ProSHADE_internal_mapManip::connectMaskBlobs
void connectMaskBlobs(proshade_double *&mask, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single maskThres)
This function connects blobs in mask.
Definition: ProSHADE_mapManip.cpp:1819
ProSHADE_internal_mapManip::moveMapByIndices
void moveMapByIndices(proshade_single *xMov, proshade_single *yMov, proshade_single *zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed *xFrom, proshade_signed *xTo, proshade_signed *yFrom, proshade_signed *yTo, proshade_signed *zFrom, proshade_signed *zTo, proshade_signed *xOrigin, proshade_signed *yOrigin, proshade_signed *zOrigin)
Function for moving map back to original PDB location by changing the indices.
Definition: ProSHADE_mapManip.cpp:744