ProSHADE  0.7.6.1 (AUG 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 ( static_cast< proshade_signed > ( std::round ( x ) ) );
35 #else
36  return ( static_cast< proshade_signed > ( round ( x ) ) );
37 #endif
38 }
39 
40 //==================================================== Define round for C++98
46 proshade_signed ProSHADE_internal_mapManip::myRound ( proshade_single x )
47 {
48 #if __cplusplus >= 201103L
49  return ( static_cast< proshade_signed > ( std::round ( x ) ) );
50 #else
51  return ( static_cast< proshade_signed > ( round ( x ) ) );
52 #endif
53 }
54 
72 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 )
73 {
74  //================================================ Initialise structure crawl
75  bool firstAtom = true;
76 
77  //================================================ Use the first model, if it exists
78  if ( pdbFile.models.size() > 0 )
79  {
80  //============================================ For each model
81  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); sIt++ )
82  {
83  //======================================== Check if multiple models are allowed
84  if ( firstModel && ( sIt != 0 ) ) { break; }
85 
86  //======================================== Get model
87  gemmi::Model model = pdbFile.models.at(sIt);
88 
89  //======================================== For each chain
90  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model.chains.size() ); mIt++ )
91  {
92  //==================================== Get chain
93  gemmi::Chain chain = model.chains.at(mIt);
94 
95  //==================================== For each residue
96  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain.residues.size() ); rIt++ )
97  {
98  //================================ Get residue
99  gemmi::Residue residue = chain.residues.at(rIt);
100 
101  //================================ For each atom
102  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue.atoms.size() ); aIt++ )
103  {
104  //============================ Get atom
105  gemmi::Atom atom = residue.atoms.at(aIt);
106 
107  //============================ Ignore hydrogens, map computations ignore them anyway and inclusion here causes map - co-ordinate mismatches.
108  if ( atom.is_hydrogen() ) { continue; }
109 
110  //============================ Find the coordinate ranges
111  if ( firstAtom )
112  {
113  *xTo = static_cast<proshade_single> ( atom.pos.x );
114  *xFrom = static_cast<proshade_single> ( atom.pos.x );
115  *yTo = static_cast<proshade_single> ( atom.pos.y );
116  *yFrom = static_cast<proshade_single> ( atom.pos.y );
117  *zTo = static_cast<proshade_single> ( atom.pos.z );
118  *zFrom = static_cast<proshade_single> ( atom.pos.z );
119  firstAtom = false;
120  }
121  else
122  {
123  if ( static_cast<proshade_single> ( atom.pos.x ) > *xTo ) { *xTo = static_cast<proshade_single> ( atom.pos.x ); }
124  if ( static_cast<proshade_single> ( atom.pos.x ) < *xFrom ) { *xFrom = static_cast<proshade_single> ( atom.pos.x ); }
125  if ( static_cast<proshade_single> ( atom.pos.y ) > *yTo ) { *yTo = static_cast<proshade_single> ( atom.pos.y ); }
126  if ( static_cast<proshade_single> ( atom.pos.y ) < *yFrom ) { *yFrom = static_cast<proshade_single> ( atom.pos.y ); }
127  if ( static_cast<proshade_single> ( atom.pos.z ) > *zTo ) { *zTo = static_cast<proshade_single> ( atom.pos.z ); }
128  if ( static_cast<proshade_single> ( atom.pos.z ) < *zFrom ) { *zFrom = static_cast<proshade_single> ( atom.pos.z ); }
129  }
130  }
131  }
132  }
133  }
134  }
135  else
136  {
137  std::stringstream hlpSS;
138  hlpSS << "Found 0 models in input file " << pdbFile.name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
139  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
140  }
141 
142  //================================================ Done
143  return ;
144 
145 }
146 
157 void ProSHADE_internal_mapManip::findPDBCOMValues ( gemmi::Structure pdbFile, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, bool firstModel )
158 {
159  //================================================ Initialise structure crawl
160  proshade_double totAtoms = 0.0;
161  *xCom = 0.0;
162  *yCom = 0.0;
163  *zCom = 0.0;
164 
165  //================================================ Use the first model, if it exists
166  if ( pdbFile.models.size() > 0 )
167  {
168  //============================================ For each model
169  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); sIt++ )
170  {
171  //======================================== Get model
172  gemmi::Model model = pdbFile.models.at(sIt);
173 
174  //======================================== Check if multiple models are allowed
175  if ( firstModel && ( sIt != 0 ) ) { break; }
176 
177  //======================================== For each chain
178  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model.chains.size() ); mIt++ )
179  {
180  //==================================== Get chain
181  gemmi::Chain chain = model.chains.at(mIt);
182 
183  //==================================== For each residue
184  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain.residues.size() ); rIt++ )
185  {
186  //================================ Get residue
187  gemmi::Residue residue = chain.residues.at(rIt);
188 
189  //================================ For each atom
190  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue.atoms.size() ); aIt++ )
191  {
192  //============================ Get atom
193  gemmi::Atom atom = residue.atoms.at(aIt);
194 
195  //============================ Save the COM sums
196  *xCom += atom.pos.x * atom.element.weight();
197  *yCom += atom.pos.y * atom.element.weight();
198  *zCom += atom.pos.z * atom.element.weight();
199  totAtoms += atom.element.weight();
200  }
201  }
202  }
203  }
204  }
205  else
206  {
207  std::stringstream hlpSS;
208  hlpSS << "Found 0 models in input file " << pdbFile.name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
209  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
210  }
211 
212  //================================================ Normalise sums to COM
213  *xCom /= totAtoms;
214  *yCom /= totAtoms;
215  *zCom /= totAtoms;
216 
217  //================================================ Done
218  return ;
219 
220 }
221 
240 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 )
241 {
242  //================================================ Initialise computation
243  proshade_double totDensity = 0.0;
244  *xCom = 0.0;
245  *yCom = 0.0;
246  *zCom = 0.0;
247  proshade_signed arrPos = 0;
248  proshade_single xSampRate = xAngs / static_cast< proshade_single > ( xTo - xFrom );
249  proshade_single ySampRate = yAngs / static_cast< proshade_single > ( yTo - yFrom );
250  proshade_single zSampRate = zAngs / static_cast< proshade_single > ( zTo - zFrom );
251 
252  //================================================ For each map point
253  for ( proshade_signed xIt = xFrom; xIt <= xTo; xIt++ )
254  {
255  for ( proshade_signed yIt = yFrom; yIt <= yTo; yIt++ )
256  {
257  for ( proshade_signed zIt = zFrom; zIt <= zTo; zIt++ )
258  {
259  arrPos = (zIt-zFrom) + ( zTo - zFrom + 1 ) * ( ( yIt - yFrom ) + ( yTo - yFrom + 1 ) * ( xIt - xFrom ) );
260 
261  if ( map[arrPos] > 0.0 )
262  {
263  totDensity += map[arrPos];
264  *xCom += static_cast<proshade_double> ( static_cast< proshade_single > ( xIt ) * xSampRate ) * map[arrPos];
265  *yCom += static_cast<proshade_double> ( static_cast< proshade_single > ( yIt ) * ySampRate ) * map[arrPos];
266  *zCom += static_cast<proshade_double> ( static_cast< proshade_single > ( zIt ) * zSampRate ) * map[arrPos];
267  }
268  }
269  }
270  }
271 
272  //================================================ Normalise sums to COM
273  *xCom /= totDensity;
274  *yCom /= totDensity;
275  *zCom /= totDensity;
276 
277  //================================================ Done
278  return ;
279 
280 }
281 
296 void ProSHADE_internal_mapManip::rotatePDBCoordinates ( gemmi::Structure *pdbFile, proshade_double euA, proshade_double euB, proshade_double euG, proshade_double xCom,
297 proshade_double yCom, proshade_double zCom, bool firstModel )
298 {
299  //================================================ Convert Euler angles to rotation matrix
300  proshade_double *rotMat = new proshade_double[9];
301  ProSHADE_internal_misc::checkMemoryAllocation ( rotMat, __FILE__, __LINE__, __func__ );
303 
304  //================================================ Initialise internal variables
305  proshade_double xTmp, yTmp, zTmp;
306 
307  //================================================ Use the first model, if it exists
308  if ( pdbFile->models.size() > 0 )
309  {
310  //============================================ For each model
311  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
312  {
313  //======================================== Get model
314  gemmi::Model *model = &pdbFile->models.at(sIt);
315 
316  //======================================== Check if multiple models are allowed
317  if ( firstModel && ( sIt != 0 ) ) { break; }
318 
319  //======================================== For each chain
320  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
321  {
322  //==================================== Get chain
323  gemmi::Chain *chain = &model->chains.at(mIt);
324 
325  //==================================== For each residue
326  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
327  {
328  //================================ Get residue
329  gemmi::Residue *residue = &chain->residues.at(rIt);
330 
331  //================================ For each atom
332  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
333  {
334  //============================ Get atom
335  gemmi::Atom *atom = &residue->atoms.at(aIt);
336 
337  //============================ Move to mid-point
338  xTmp = static_cast< proshade_double > ( atom->pos.x - xCom );
339  yTmp = static_cast< proshade_double > ( atom->pos.y - yCom );
340  zTmp = static_cast< proshade_double > ( atom->pos.z - zCom );
341 
342  //============================ Rotate the atom position
343  atom->pos.x = ( xTmp * rotMat[0] ) + ( yTmp * rotMat[1] ) + ( zTmp * rotMat[2] );
344  atom->pos.y = ( xTmp * rotMat[3] ) + ( yTmp * rotMat[4] ) + ( zTmp * rotMat[5] );
345  atom->pos.z = ( xTmp * rotMat[6] ) + ( yTmp * rotMat[7] ) + ( zTmp * rotMat[8] );
346 
347  //============================ Move back
348  atom->pos.x = atom->pos.x + xCom;
349  atom->pos.y = atom->pos.y + yCom;
350  atom->pos.z = atom->pos.z + zCom;
351  }
352  }
353  }
354  }
355  }
356  else
357  {
358  std::stringstream hlpSS;
359  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
360  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
361  }
362 
363  //================================================ Release memory
364  delete[] rotMat;
365 
366  //================================================ Done
367  return ;
368 
369 }
370 
381 void ProSHADE_internal_mapManip::translatePDBCoordinates ( gemmi::Structure *pdbFile, proshade_double transX, proshade_double transY, proshade_double transZ, bool firstModel )
382 {
383  //================================================ Use the first model, if it exists
384  if ( pdbFile->models.size() > 0 )
385  {
386  //============================================ For each model
387  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
388  {
389  //======================================== Check if multiple models are allowed
390  if ( firstModel && ( sIt != 0 ) ) { break; }
391 
392  //======================================== Get model
393  gemmi::Model *model = &pdbFile->models.at(sIt);
394 
395  //======================================== For each chain
396  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
397  {
398  //==================================== Get chain
399  gemmi::Chain *chain = &model->chains.at(mIt);
400 
401  //==================================== For each residue
402  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
403  {
404  //================================ Get residue
405  gemmi::Residue *residue = &chain->residues.at(rIt);
406 
407  //================================ For each atom
408  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
409  {
410  //============================ Get atom
411  gemmi::Atom *atom = &residue->atoms.at(aIt);
412 
413  //============================ Translate
414  atom->pos.x += transX;
415  atom->pos.y += transY;
416  atom->pos.z += transZ;
417  }
418  }
419  }
420  }
421  }
422  else
423  {
424  std::stringstream hlpSS;
425  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
426  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
427  }
428 
429  //================================================ Done
430  return ;
431 
432 }
433 
444 void ProSHADE_internal_mapManip::changePDBBFactors ( gemmi::Structure *pdbFile, proshade_double newBFactorValue, bool firstModel )
445 {
446  //================================================ Use the first model, if it exists
447  if ( pdbFile->models.size() > 0 )
448  {
449  //============================================ For each model
450  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
451  {
452  //======================================== Check if multiple models are allowed
453  if ( firstModel && ( sIt != 0 ) ) { break; }
454 
455  //======================================== Get model
456  gemmi::Model *model = &pdbFile->models.at(sIt);
457 
458  //======================================== For each chain
459  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
460  {
461  //==================================== Get chain
462  gemmi::Chain *chain = &model->chains.at(mIt);
463 
464  //==================================== For each residue
465  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
466  {
467  //================================ Get residue
468  gemmi::Residue *residue = &chain->residues.at(rIt);
469 
470  //================================ For each atom
471  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
472  {
473  //============================ Get atom
474  gemmi::Atom *atom = &residue->atoms.at(aIt);
475 
476  //============================ Change the B-factors
477  atom->b_iso = static_cast< float > ( newBFactorValue );
478  }
479  }
480  }
481  }
482  }
483  else
484  {
485  std::stringstream hlpSS;
486  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
487  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
488  }
489 
490  //================================================ Done
491  return ;
492 
493 }
494 
504 void ProSHADE_internal_mapManip::removeWaters ( gemmi::Structure *pdbFile, bool firstModel )
505 {
506  //================================================ Use the first model, if it exists
507  if ( pdbFile->models.size() > 0 )
508  {
509  //============================================ For each model
510  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
511  {
512  //======================================== Check if multiple models are allowed
513  if ( firstModel && ( sIt != 0 ) ) { break; }
514 
515  //======================================== Get model
516  gemmi::Model *model = &pdbFile->models.at(sIt);
517 
518  //======================================== For each chain
519  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
520  {
521  //==================================== Get chain
522  gemmi::Chain *chain = &model->chains.at(mIt);
523 
524  //==================================== Initialise del vector
525  std::vector< proshade_unsign > delVec;
526 
527  //==================================== For each residue
528  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
529  {
530  //================================ Get residue
531  gemmi::Residue *residue = &chain->residues.at(rIt);
532 
533  //================================ If residue is water
534  if ( residue->is_water() )
535  {
537  }
538  }
539 
540  //==================================== Delete from end to avoid indexing issues
541  std::sort ( delVec.begin(), delVec.end(), std::greater<int>() );
542  for ( proshade_unsign vecIt = 0; vecIt < static_cast<proshade_unsign> ( delVec.size() ); vecIt++ )
543  {
544  chain->residues.erase ( chain->residues.begin() + static_cast< long int > ( delVec.at(vecIt) ) );
545  }
546  }
547  }
548  }
549  else
550  {
551  std::stringstream hlpSS;
552  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
553  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
554  }
555 
556  //================================================ Done
557  return ;
558 
559 }
560 
573 void ProSHADE_internal_mapManip::movePDBForMapCalc ( gemmi::Structure *pdbFile, proshade_single xMov, proshade_single yMov, proshade_single zMov, bool firstModel )
574 {
575  //================================================ Use the first model, if it exists
576  if ( pdbFile->models.size() > 0 )
577  {
578  //============================================ For each model
579  for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
580  {
581  //======================================== Check if multiple models are allowed
582  if ( firstModel && ( sIt != 0 ) ) { break; }
583 
584  //======================================== Get model
585  gemmi::Model *model = &pdbFile->models.at(sIt);
586 
587  //======================================== For each chain
588  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
589  {
590  //==================================== Get chain
591  gemmi::Chain *chain = &model->chains.at(mIt);
592 
593  //==================================== For each residue
594  for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
595  {
596  //================================ Get residue
597  gemmi::Residue *residue = &chain->residues.at(rIt);
598 
599  //================================ For each atom
600  for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
601  {
602  //============================ Get atom
603  gemmi::Atom *atom = &residue->atoms.at(aIt);
604 
605  //============================ Move the atoms
606  atom->pos = gemmi::Position ( atom->pos.x + static_cast< proshade_double > ( xMov ), atom->pos.y + static_cast< proshade_double > ( yMov ), atom->pos.z + static_cast< proshade_double > ( zMov ) );
607  }
608  }
609  }
610 
611  }
612  }
613  else
614  {
615  std::stringstream hlpSS;
616  hlpSS << "Found 0 models in input file " << pdbFile->name << ".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
617  throw ProSHADE_exception ( "Found no model in co-ordinate file.", "EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
618  }
619 
620  //================================================ Done
621  return ;
622 }
623 
644 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 )
645 {
646  //================================================ Set cell dimensions from the increased ranges (we need to add some space) and re-calculate cell properties
647  if ( forceP1 ) { pdbFile.cell = gemmi::UnitCell(); }
648  pdbFile.cell.a = static_cast< proshade_double > ( xCell );
649  pdbFile.cell.b = static_cast< proshade_double > ( yCell );
650  pdbFile.cell.c = static_cast< proshade_double > ( zCell );
651  pdbFile.cell.calculate_properties ( );
652 
653  //================================================ Get elements in Gemmi format
654  std::string totElString;
655  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); mIt++ )
656  {
657  //============================================ Check if multiple models are allowed
658  if ( firstModel && ( mIt != 0 ) )
659  {
660  std::stringstream hlpSS;
661  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.";
662  ProSHADE_internal_messages::printWarningMessage ( 0, hlpSS.str(), "WP00055" );
663  break;
664  }
665 
666  std::string hlpStr = pdbFile.models[mIt].present_elements ( ).to_string<char,std::char_traits<char>,std::allocator<char> >();
667  totElString = totElString + hlpStr;
668  }
669  std::bitset< static_cast< size_t > ( gemmi::El::END )> present_elems ( totElString );
670 
671  //================================================ Sanity checks
672  if ( present_elems[static_cast<int> ( gemmi::El::X )] )
673  {
674  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." );
675  }
676 
677  for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( present_elems.size() ); elIt++ )
678  {
679  if ( present_elems[elIt] && !gemmi::IT92<double>::has ( static_cast<gemmi::El> ( elIt ) ) )
680  {
681  std::stringstream hlpSS;
682  hlpSS << "Missing form factor for element " << element_name ( static_cast<gemmi::El> ( elIt ) );
683  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." );
684  }
685  }
686 
687  //================================================ Compute the f's
688  double wavelength = 0.1;
689  double energy = gemmi::hc() / wavelength;
690 
691  //================================================ Create the density calculator object and fill it in
692  gemmi::DensityCalculator<gemmi::IT92<double>, float> dencalc;
693 
694  dencalc.d_min = static_cast< double > ( requestedResolution );
695  for ( size_t elIt = 0; elIt < present_elems.size(); elIt++ ) { if ( present_elems[elIt] ) { dencalc.addends.set ( static_cast< gemmi::El > ( elIt ), static_cast< float > ( gemmi::cromer_liberman ( static_cast< int > ( elIt ), energy, nullptr ) ) ); } }
696  dencalc.set_grid_cell_and_spacegroup ( pdbFile );
697 
698  //================================================ Force P1 spacegroup
699  if ( forceP1 ) { dencalc.grid.spacegroup = &gemmi::get_spacegroup_p1(); }
700 
701  //================================================ Compute the theoretical map for each model
702  dencalc.grid.data.clear ( );
703  dencalc.grid.set_size_from_spacing ( dencalc.d_min / ( 2.0 * dencalc.rate), true );
704  for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); mIt++ )
705  {
706  if ( firstModel && ( mIt != 0 ) ) { break; }
707  dencalc.add_model_density_to_grid ( pdbFile.models[mIt] );
708  dencalc.grid.symmetrize ( [](float a, float b) { return a + b; } );
709  }
710 
711  //================================================ Get the map
712  const gemmi::Grid<float>& grid = dencalc.grid;
713 
714  //================================================ Save the map dimensions
715  *xTo = grid.nu;
716  *yTo = grid.nv;
717  *zTo = grid.nw;
718 
719  //================================================ Copy the gemmi::Grid to my map format
720  map = new proshade_double [(*xTo) * (*yTo) * (*zTo)];
721  ProSHADE_internal_misc::checkMemoryAllocation ( map, __FILE__, __LINE__, __func__ );
722 
723  proshade_signed arrPos = 0;
724  for ( proshade_signed uIt = 0; uIt < (*xTo); uIt++ )
725  {
726  for ( proshade_signed vIt = 0; vIt < (*yTo); vIt++ )
727  {
728  for ( proshade_signed wIt = 0; wIt < (*zTo); wIt++ )
729  {
730  arrPos = wIt + (*zTo) * ( vIt + (*yTo) * uIt );
731  map[arrPos] = static_cast< proshade_double > ( grid.get_value_q( static_cast< int > ( uIt ), static_cast< int > ( vIt ), static_cast< int > ( wIt ) ) );
732  }
733  }
734  }
735 
736  //================================================ Done
737  return ;
738 
739 }
740 
763 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 )
764 {
765  //================================================ Compute movement in indices
766  proshade_single xIndMove = std::floor ( -(*xMov) / ( xAngs / ( static_cast< proshade_single > ( *xTo ) - static_cast< proshade_single > ( *xFrom ) ) ) );
767  proshade_single yIndMove = std::floor ( -(*yMov) / ( yAngs / ( static_cast< proshade_single > ( *yTo ) - static_cast< proshade_single > ( *yFrom ) ) ) );
768  proshade_single zIndMove = std::floor ( -(*zMov) / ( zAngs / ( static_cast< proshade_single > ( *zTo ) - static_cast< proshade_single > ( *zFrom ) ) ) );
769 
770  //================================================ Set the movs to the remainder
771  *xMov = -( *xMov ) - ( xIndMove * ( xAngs / ( static_cast< proshade_single > ( *xTo ) - static_cast< proshade_single > ( *xFrom ) ) ) );
772  *yMov = -( *yMov ) - ( yIndMove * ( yAngs / ( static_cast< proshade_single > ( *yTo ) - static_cast< proshade_single > ( *yFrom ) ) ) );
773  *zMov = -( *zMov ) - ( zIndMove * ( zAngs / ( static_cast< proshade_single > ( *zTo ) - static_cast< proshade_single > ( *zFrom ) ) ) );
774 
775  //================================================ Move indices by as much
776  *xFrom += static_cast< proshade_signed > ( xIndMove );
777  *xTo += static_cast< proshade_signed > ( xIndMove );
778  *yFrom += static_cast< proshade_signed > ( yIndMove );
779  *yTo += static_cast< proshade_signed > ( yIndMove );
780  *zFrom += static_cast< proshade_signed > ( zIndMove );
781  *zTo += static_cast< proshade_signed > ( zIndMove );
782 
783  //================================================ And set origin to reflect the changes
784  *xOrigin = *xFrom;
785  *yOrigin = *yFrom;
786  *zOrigin = *zFrom;
787 
788  //================================================ Done
789  return ;
790 
791 }
792 
811 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 )
812 {
813  //================================================ Local variables initialisation
814  proshade_unsign arrayPos = 0;
815  proshade_signed h, k, l;
816  proshade_double real = 0.0;
817  proshade_double imag = 0.0;
818  proshade_double trCoeffReal, trCoeffImag;
819  proshade_double normFactor = static_cast< proshade_double > ( xDim * yDim * zDim );
820  proshade_double exponent = 0.0;
821  proshade_double hlpArrReal;
822  proshade_double hlpArrImag;
823 
824  //================================================ Create Fourier map variable
825  fftw_complex *fCoeffs = new fftw_complex [xDim * yDim * zDim];
826  fftw_complex *translatedMap = new fftw_complex [xDim * yDim * zDim];
827 
828  //================================================ Check memory allocation
829  ProSHADE_internal_misc::checkMemoryAllocation ( fCoeffs, __FILE__, __LINE__, __func__ );
830  ProSHADE_internal_misc::checkMemoryAllocation ( translatedMap, __FILE__, __LINE__, __func__ );
831 
832  //================================================ Create plans
833  fftw_plan planForwardFourier = fftw_plan_dft_3d ( static_cast< int > ( xDim ), static_cast< int > ( yDim ), static_cast< int > ( zDim ), translatedMap, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
834  fftw_plan planBackwardFourier = fftw_plan_dft_3d ( static_cast< int > ( xDim ), static_cast< int > ( yDim ), static_cast< int > ( zDim ), fCoeffs, translatedMap, FFTW_BACKWARD, FFTW_ESTIMATE );
835 
836  //================================================ Copy map to complex format
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  arrayPos = wIt + static_cast< proshade_unsign > ( zDim ) * ( vIt + static_cast< proshade_unsign > ( yDim ) * uIt );
844 
845  const FloatingPoint< proshade_double > lhs ( map[arrayPos] ), rhs ( map[arrayPos] );
846  if ( lhs.AlmostEquals ( rhs ) ) { translatedMap[arrayPos][0] = map[arrayPos]; }
847  else { translatedMap[arrayPos][0] = 0.0; }
848  translatedMap[arrayPos][1] = 0.0;
849  }
850  }
851  }
852 
853  //================================================ Compute Forward Fourier
854  fftw_execute ( planForwardFourier );
855 
856  //================================================ Add the required shift
857  for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
858  {
859  for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
860  {
861  for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
862  {
863  //==================================== Var init
864  arrayPos = wIt + static_cast< proshade_unsign > ( zDim ) * ( vIt + static_cast< proshade_unsign > ( yDim ) * uIt );
865  real = fCoeffs[arrayPos][0];
866  imag = fCoeffs[arrayPos][1];
867 
868  //==================================== Change the B-factors, if required
869  if ( uIt > static_cast< proshade_unsign > ( (xDim+1) / 2) ) { h = static_cast < proshade_signed > ( uIt ) - xDim; } else { h = static_cast < proshade_signed > ( uIt ); }
870  if ( vIt > static_cast< proshade_unsign > ( (yDim+1) / 2) ) { k = static_cast < proshade_signed > ( vIt ) - yDim; } else { k = static_cast < proshade_signed > ( vIt ); }
871  if ( wIt > static_cast< proshade_unsign > ( (zDim+1) / 2) ) { l = static_cast < proshade_signed > ( wIt ) - zDim; } else { l = static_cast < proshade_signed > ( wIt ); }
872 
873  //==================================== Get translation coefficient change
874  exponent = ( ( ( static_cast <proshade_double> ( h ) / static_cast <proshade_double> ( xAngs ) ) * static_cast< proshade_double > ( -xMov ) ) +
875  ( ( static_cast <proshade_double> ( k ) / static_cast <proshade_double> ( yAngs ) ) * static_cast< proshade_double > ( -yMov ) ) +
876  ( ( static_cast <proshade_double> ( l ) / static_cast <proshade_double> ( zAngs ) ) * static_cast< proshade_double > ( -zMov ) ) ) * 2.0 * M_PI;
877 
878  trCoeffReal = cos ( exponent );
879  trCoeffImag = sin ( exponent );
880  ProSHADE_internal_maths::complexMultiplication ( &real, &imag, &trCoeffReal, &trCoeffImag, &hlpArrReal, &hlpArrImag );
881 
882  //==================================== Save the mask data
883  fCoeffs[arrayPos][0] = hlpArrReal / normFactor;
884  fCoeffs[arrayPos][1] = hlpArrImag / normFactor;
885  }
886  }
887  }
888 
889  //================================================ Compute inverse Fourier
890  fftw_execute ( planBackwardFourier );
891 
892  //================================================ Copy back to map
893  for ( proshade_unsign uIt = 0; uIt < static_cast< proshade_unsign > ( xDim ); uIt++ )
894  {
895  for ( proshade_unsign vIt = 0; vIt < static_cast< proshade_unsign > ( yDim ); vIt++ )
896  {
897  for ( proshade_unsign wIt = 0; wIt < static_cast< proshade_unsign > ( zDim ); wIt++ )
898  {
899  arrayPos = wIt + static_cast< proshade_unsign > ( zDim ) * ( vIt + static_cast< proshade_unsign > ( yDim ) * uIt );
900  map[arrayPos] = translatedMap[arrayPos][0];
901  }
902  }
903  }
904 
905  //================================================ Release memory
906  fftw_destroy_plan ( planForwardFourier );
907  fftw_destroy_plan ( planBackwardFourier );
908  delete[] fCoeffs;
909  delete[] translatedMap;
910 
911  //================================================ Done
912  return ;
913 
914 }
915 
932 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 )
933 {
934  //================================================ Set local variables
935  proshade_signed xDim = static_cast< proshade_signed > ( xDimS );
936  proshade_signed yDim = static_cast< proshade_signed > ( yDimS );
937  proshade_signed zDim = static_cast< proshade_signed > ( zDimS );
938  proshade_double real, imag, S, mag, phase;
939  proshade_signed h, k, l;
940  proshade_unsign arrayPos = 0;
941  proshade_double normFactor = static_cast<proshade_double> ( xDim * yDim * zDim );
942 
943  //================================================ Copy map for processing
944  fftw_complex* mapCoeffs = new fftw_complex[xDim * yDim * zDim];
945  fftw_complex* mapMask = new fftw_complex[xDim * yDim * zDim];
946 
947  //================================================ Check memory allocation
948  ProSHADE_internal_misc::checkMemoryAllocation ( mapCoeffs, __FILE__, __LINE__, __func__ );
949  ProSHADE_internal_misc::checkMemoryAllocation ( mapMask, __FILE__, __LINE__, __func__ );
950 
951  //================================================ Copy data to mask
952  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> (xDim * yDim * zDim); iter++ )
953  {
954  mapMask[iter][0] = map[iter];
955  mapMask[iter][1] = 0.0;
956  }
957 
958  //================================================ Prepare FFTW plans
959  fftw_plan forward = fftw_plan_dft_3d ( static_cast< int > ( xDim ), static_cast< int > ( yDim ), static_cast< int > ( zDim ), mapMask, mapCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
960  fftw_plan inverse = fftw_plan_dft_3d ( static_cast< int > ( xDim ), static_cast< int > ( yDim ), static_cast< int > ( zDim ), mapCoeffs, mapMask, FFTW_BACKWARD, FFTW_ESTIMATE );
961 
962  //================================================ Run forward Fourier
963  fftw_execute ( forward );
964 
965  //================================================ Blur the coeffs
966  for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
967  {
968  for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
969  {
970  for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
971  {
972  //==================================== Var init
973  arrayPos = wIt + static_cast< proshade_unsign > ( zDim ) * ( vIt + static_cast< proshade_unsign > ( yDim ) * uIt );
974  real = mapCoeffs[arrayPos][0];
975  imag = mapCoeffs[arrayPos][1];
976 
977  //==================================== Convert to HKL
978  if ( uIt > static_cast< proshade_unsign > ( (xDim+1) / 2) ) { h = static_cast < proshade_signed > ( uIt ) - xDim; } else { h = static_cast < proshade_signed > ( uIt ); }
979  if ( vIt > static_cast< proshade_unsign > ( (yDim+1) / 2) ) { k = static_cast < proshade_signed > ( vIt ) - yDim; } else { k = static_cast < proshade_signed > ( vIt ); }
980  if ( wIt > static_cast< proshade_unsign > ( (zDim+1) / 2) ) { l = static_cast < proshade_signed > ( wIt ) - zDim; } else { l = static_cast < proshade_signed > ( wIt ); }
981 
982  //====================================Get magnitude and phase with mask parameters
983  S = ( pow( static_cast< proshade_double > ( h ) / static_cast< proshade_double > ( xAngs ), 2.0 ) +
984  pow( static_cast< proshade_double > ( k ) / static_cast< proshade_double > ( yAngs ), 2.0 ) +
985  pow( static_cast< proshade_double > ( l ) / static_cast< proshade_double > ( zAngs ), 2.0 ) );
986  mag = std::sqrt ( (real*real) + (imag*imag) ) * std::exp ( - ( ( static_cast< proshade_double > ( blurringFactor ) * S ) / 4.0 ) );
987  phase = std::atan2 ( imag, real );
988 
989  //==================================== Save the mask data
990  mapCoeffs[arrayPos][0] = ( mag * cos(phase) ) / normFactor;
991  mapCoeffs[arrayPos][1] = ( mag * sin(phase) ) / normFactor;
992  }
993  }
994  }
995 
996  //================================================ Run inverse Fourier
997  fftw_execute ( inverse );
998 
999  //================================================ Save the results
1000  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> (xDim * yDim * zDim); iter++ )
1001  {
1002  blurredMap[iter] = mapMask[iter][0];
1003  }
1004 
1005  //================================================ Release memory
1006  delete[] mapMask;
1007  delete[] mapCoeffs;
1008 
1009  //================================================ Delete FFTW plans
1010  fftw_destroy_plan ( forward );
1011  fftw_destroy_plan ( inverse );
1012 
1013  //================================================ Done
1014  return ;
1015 
1016 }
1017 
1032 void ProSHADE_internal_mapManip::getMaskFromBlurr ( proshade_double*& blurMap, proshade_double*& outMap, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_single noIQRs )
1033 {
1034  //================================================ Initialise vector for map data
1035  std::vector<proshade_double> mapVals ( xDim * yDim * zDim, 0.0 );
1036 
1037  //================================================ Save map values in vector
1038  for ( proshade_unsign iter = 0; iter < ( xDim * yDim * zDim ); iter++ )
1039  {
1040  mapVals.at(iter) = blurMap[iter];
1041  }
1042 
1043  //================================================ Find median and IQRs
1044  proshade_double* medAndIQR = new proshade_double[2];
1045  ProSHADE_internal_maths::vectorMedianAndIQR ( &mapVals, medAndIQR );
1046 
1047  //================================================ Find the threshold
1048  proshade_double maskThreshold = medAndIQR[0] + ( medAndIQR[1] * static_cast<proshade_double> ( noIQRs ) );
1049 
1050  //================================================ Apply threshold
1051  for ( proshade_unsign iter = 0; iter < ( xDim * yDim * zDim ); iter++ )
1052  {
1053  if ( blurMap[iter] < maskThreshold )
1054  {
1055  outMap[iter] = 0.0;
1056  blurMap[iter] = 0.0;
1057  }
1058  }
1059 
1060  //================================================ Release vector values
1061  mapVals.clear ( );
1062 
1063  //================================================ Release memory
1064  delete[] medAndIQR;
1065 
1066  //================================================ Done
1067  return ;
1068 
1069 }
1070 
1082 void ProSHADE_internal_mapManip::getNonZeroBounds ( proshade_double* map, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_signed*& ret )
1083 {
1084  //================================================ Initialise local variables
1085  proshade_signed arrayPos = 0;
1086 
1087  //================================================ Initialise result variable
1088  ret[0] = xDim;
1089  ret[1] = 0;
1090  ret[2] = yDim;
1091  ret[3] = 0;
1092  ret[4] = zDim;
1093  ret[5] = 0;
1094 
1095  //================================================ Iterate through map and check bounds
1096  for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1097  {
1098  for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1099  {
1100  for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1101  {
1102  //==================================== Var init
1103  arrayPos = zIt + zDim * ( yIt + yDim * xIt );
1104 
1105  //==================================== Check bounds
1106  if ( map[arrayPos] > 0.001 )
1107  {
1108  if ( xIt < ret[0] ) { ret[0] = xIt; }
1109  if ( xIt > ret[1] ) { ret[1] = xIt; }
1110  if ( yIt < ret[2] ) { ret[2] = yIt; }
1111  if ( yIt > ret[3] ) { ret[3] = yIt; }
1112  if ( zIt < ret[4] ) { ret[4] = zIt; }
1113  if ( zIt > ret[5] ) { ret[5] = zIt; }
1114  }
1115  }
1116  }
1117  }
1118 
1119  //================================================ Done
1120  return ;
1121 
1122 }
1123 
1139 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 )
1140 {
1141  //================================================ Convert angstrom distance to indices
1142  proshade_signed xExtraInds = ProSHADE_internal_mapManip::myRound ( extraSpace / ( xAngs / static_cast<proshade_single> ( xDim ) ) );
1143  proshade_signed yExtraInds = ProSHADE_internal_mapManip::myRound ( extraSpace / ( yAngs / static_cast<proshade_single> ( yDim ) ) );
1144  proshade_signed zExtraInds = ProSHADE_internal_mapManip::myRound ( extraSpace / ( zAngs / static_cast<proshade_single> ( zDim ) ) );
1145 
1146  //=============================================== Changed bounds even if exceeding physical map - this will be deal with in the map creation part
1147  bounds[0] = bounds[0] - xExtraInds;
1148  bounds[1] = bounds[1] + xExtraInds;
1149  bounds[2] = bounds[2] - yExtraInds;
1150  bounds[3] = bounds[3] + yExtraInds;
1151  bounds[4] = bounds[4] - zExtraInds;
1152  bounds[5] = bounds[5] + zExtraInds;
1153 
1154  //================================================ Done
1155  return ;
1156 
1157 }
1158 
1175 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 )
1176 {
1177  //================================================ Sanity check - the resolution needs to be set
1178  if ( resolution <= 0.0f )
1179  {
1180  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." );
1181  }
1182 
1183  //================================================ Initialise local variables
1184  proshade_signed xDim = static_cast<proshade_signed> ( xDimS );
1185  proshade_signed yDim = static_cast<proshade_signed> ( yDimS );
1186  proshade_signed zDim = static_cast<proshade_signed> ( zDimS );
1187  proshade_single oldXSample = ( xAngs / static_cast<proshade_single> ( xDim ) );
1188  proshade_single oldYSample = ( yAngs / static_cast<proshade_single> ( yDim ) );
1189  proshade_single oldZSample = ( zAngs / static_cast<proshade_single> ( zDim ) );
1190  proshade_single newXSample = static_cast< proshade_single > ( resolution / 2.0f );
1191  proshade_single newYSample = static_cast< proshade_single > ( resolution / 2.0f );
1192  proshade_single newZSample = static_cast< proshade_single > ( resolution / 2.0f );
1193 
1194  //================================================ Compute required grid size
1195  proshade_signed newXDim = static_cast<proshade_signed> ( std::ceil ( xAngs / newXSample ) );
1196  proshade_signed newYDim = static_cast<proshade_signed> ( std::ceil ( yAngs / newYSample ) );
1197  proshade_signed newZDim = static_cast<proshade_signed> ( std::ceil ( zAngs / newZSample ) );
1198 
1199  //================================================ Create a new map variable
1200  proshade_double* newMap = new proshade_double [newXDim * newYDim * newZDim];
1201 
1202  //================================================ For each new map point
1203  proshade_signed xBottom = 0, xTop, yBottom = 0, yTop, zBottom = 0, zTop, oldMapIndex, newMapIndex;
1204  std::vector<proshade_double> c000 = std::vector<proshade_double> ( 4, 0.0 );
1205  std::vector<proshade_double> c001 = std::vector<proshade_double> ( 4, 0.0 );
1206  std::vector<proshade_double> c010 = std::vector<proshade_double> ( 4, 0.0 );
1207  std::vector<proshade_double> c011 = std::vector<proshade_double> ( 4, 0.0 );
1208  std::vector<proshade_double> c100 = std::vector<proshade_double> ( 4, 0.0 );
1209  std::vector<proshade_double> c101 = std::vector<proshade_double> ( 4, 0.0 );
1210  std::vector<proshade_double> c110 = std::vector<proshade_double> ( 4, 0.0 );
1211  std::vector<proshade_double> c111 = std::vector<proshade_double> ( 4, 0.0 );
1212  std::vector<proshade_double> c00 = std::vector<proshade_double> ( 4, 0.0 );
1213  std::vector<proshade_double> c01 = std::vector<proshade_double> ( 4, 0.0 );
1214  std::vector<proshade_double> c10 = std::vector<proshade_double> ( 4, 0.0 );
1215  std::vector<proshade_double> c11 = std::vector<proshade_double> ( 4, 0.0 );
1216  std::vector<proshade_double> c0 = std::vector<proshade_double> ( 4, 0.0 );
1217  std::vector<proshade_double> c1 = std::vector<proshade_double> ( 4, 0.0 );
1218  proshade_double xRelative, yRelative, zRelative;
1219 
1220  for ( proshade_signed xIt = 0; xIt < newXDim; xIt++ )
1221  {
1222  for ( proshade_signed yIt = 0; yIt < newYDim; yIt++ )
1223  {
1224  for ( proshade_signed zIt = 0; zIt < newZDim; zIt++ )
1225  {
1226  //==================================== Get this point's index
1227  newMapIndex = zIt + newZDim * ( yIt + newYDim * xIt );
1228 
1229  //==================================== Find this points bottom and top positions in the old map (including periodicity)
1230  for ( proshade_signed ox = 0; ox < ( static_cast< proshade_signed > ( xDimS ) - 1 ); ox++ ) { if ( ( ( static_cast< proshade_single > ( xIt ) * newXSample ) >= ( static_cast< proshade_single > ( ox ) * oldXSample ) ) && ( ( static_cast< proshade_single > ( xIt ) * newXSample ) <= ( ( static_cast< proshade_single > ( ox ) + 1 ) * oldXSample ) ) ) { xBottom = ox; break; } }
1231  for ( proshade_signed oy = 0; oy < ( static_cast< proshade_signed > ( yDimS ) - 1 ); oy++ ) { if ( ( ( static_cast< proshade_single > ( yIt ) * newYSample ) >= ( static_cast< proshade_single > ( oy ) * oldYSample ) ) && ( ( static_cast< proshade_single > ( yIt ) * newYSample ) <= ( ( static_cast< proshade_single > ( oy ) + 1 ) * oldYSample ) ) ) { yBottom = oy; break; } }
1232  for ( proshade_signed oz = 0; oz < ( static_cast< proshade_signed > ( zDimS ) - 1 ); oz++ ) { if ( ( ( static_cast< proshade_single > ( zIt ) * newZSample ) >= ( static_cast< proshade_single > ( oz ) * oldZSample ) ) && ( ( static_cast< proshade_single > ( zIt ) * newZSample ) <= ( ( static_cast< proshade_single > ( oz ) + 1 ) * oldZSample ) ) ) { zBottom = oz; break; } }
1233  xTop = xBottom + 1;
1234  yTop = yBottom + 1;
1235  zTop = zBottom + 1;
1236 
1237  //==================================== Find the surrounding point's values from the original map
1238  oldMapIndex = zBottom + static_cast< proshade_signed > ( zDimS ) * ( yBottom + static_cast< proshade_signed > ( yDimS ) * xBottom );
1239  c000.at(0) = static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample );
1240  c000.at(1) = static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample );
1241  c000.at(2) = static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample );
1242  c000.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1243 
1244  oldMapIndex = zTop + static_cast< proshade_signed > ( zDimS ) * ( yBottom + static_cast< proshade_signed > ( yDimS ) * xBottom );
1245  c001.at(0) = static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample );
1246  c001.at(1) = static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample );
1247  c001.at(2) = static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample );
1248  c001.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1249 
1250  oldMapIndex = zBottom + static_cast< proshade_signed > ( zDimS ) * ( yTop + static_cast< proshade_signed > ( yDimS ) * xBottom );
1251  c010.at(0) = static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample );
1252  c010.at(1) = static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample );
1253  c010.at(2) = static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample );
1254  c010.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1255 
1256  oldMapIndex = zTop + static_cast< proshade_signed > ( zDimS ) * ( yTop + static_cast< proshade_signed > ( yDimS ) * xBottom );
1257  c011.at(0) = static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample );
1258  c011.at(1) = static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample );
1259  c011.at(2) = static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample );
1260  c011.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1261 
1262  oldMapIndex = zBottom + static_cast< proshade_signed > ( zDimS ) * ( yBottom + static_cast< proshade_signed > ( yDimS ) * xTop );
1263  c100.at(0) = static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample );
1264  c100.at(1) = static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample );
1265  c100.at(2) = static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample );
1266  c100.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1267 
1268  oldMapIndex = zTop + static_cast< proshade_signed > ( zDimS ) * ( yBottom + static_cast< proshade_signed > ( yDimS ) * xTop );
1269  c101.at(0) = static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample );
1270  c101.at(1) = static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample );
1271  c101.at(2) = static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample );
1272  c101.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1273 
1274  oldMapIndex = zBottom + static_cast< proshade_signed > ( zDimS ) * ( yTop + static_cast< proshade_signed > ( yDimS ) * xTop );
1275  c110.at(0) = static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample );
1276  c110.at(1) = static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample );
1277  c110.at(2) = static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample );
1278  c110.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1279 
1280  oldMapIndex = zTop + static_cast< proshade_signed > ( zDimS ) * ( yTop + static_cast< proshade_signed > ( yDimS ) * xTop );
1281  c111.at(0) = static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample );
1282  c111.at(1) = static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample );
1283  c111.at(2) = static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample );
1284  c111.at(3) = static_cast<proshade_double> ( map[oldMapIndex] );
1285 
1286  //==================================== Interpolate to the new grid along X
1287  xRelative = ( ( static_cast<proshade_double> ( xIt ) * static_cast<proshade_double> ( newXSample ) ) - ( static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample ) ) ) / ( ( static_cast<proshade_double> ( xTop ) * static_cast<proshade_double> ( oldXSample ) ) - ( static_cast<proshade_double> ( xBottom ) * static_cast<proshade_double> ( oldXSample ) ) );
1288 
1289  //==================================== Interpolate for the less less point
1290  c00.at(0) = ( static_cast< proshade_double > ( newXSample ) * xRelative ) + c000.at(0);
1291  c00.at(1) = c000.at(1);
1292  c00.at(2) = c000.at(2);
1293  c00.at(3) = ( c000.at(3) * ( 1.0 - xRelative ) ) + ( c100.at(3) * xRelative );
1294 
1295  //==================================== Interpolate for the less more point
1296  c01.at(0) = ( static_cast< proshade_double > ( newXSample ) * xRelative ) + c001.at(0);
1297  c01.at(1) = c001.at(1);
1298  c01.at(2) = c001.at(2);
1299  c01.at(3) = ( c001.at(3) * ( 1.0 - xRelative ) ) + ( c101.at(3) * xRelative );
1300 
1301  //==================================== Interpolate for the more less point
1302  c10.at(0) = ( static_cast< proshade_double > ( newXSample ) * xRelative ) + c010.at(0);
1303  c10.at(1) = c010.at(1);
1304  c10.at(2) = c010.at(2);
1305  c10.at(3) = ( c010.at(3) * ( 1.0 - xRelative ) ) + ( c110.at(3) * xRelative );
1306 
1307  //==================================== Interpolate for the more more point
1308  c11.at(0) = ( static_cast< proshade_double > ( newXSample ) * xRelative ) + c011.at(0);
1309  c11.at(1) = c011.at(1);
1310  c11.at(2) = c011.at(2);
1311  c11.at(3) = ( c011.at(3) * ( 1.0 - xRelative ) ) + ( c111.at(3) * xRelative );
1312 
1313  //==================================== Interpolate to the new grid along Y
1314  yRelative = ( ( static_cast<proshade_double> ( yIt ) * static_cast<proshade_double> ( newYSample ) ) - ( static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample ) ) ) / ( ( static_cast<proshade_double> ( yTop ) * static_cast<proshade_double> ( oldYSample ) ) - ( static_cast<proshade_double> ( yBottom ) * static_cast<proshade_double> ( oldYSample ) ) );
1315 
1316  //==================================== Interpolate for the less point
1317  c0.at(0) = c00.at(0);
1318  c0.at(1) = ( static_cast< proshade_double > ( newYSample ) * yRelative ) + c00.at(1);
1319  c0.at(2) = c00.at(2);
1320  c0.at(3) = ( c00.at(3) * ( 1.0 - yRelative ) ) + ( c10.at(3) * yRelative );
1321 
1322  //==================================== Interpolate for the more point
1323  c1.at(0) = c01.at(0);
1324  c1.at(1) = ( static_cast< proshade_double > ( newYSample ) * yRelative ) + c01.at(1);
1325  c1.at(2) = c01.at(2);
1326  c1.at(3) = ( c01.at(3) * ( 1.0 - yRelative ) ) + ( c11.at(3) * yRelative );
1327 
1328  //==================================== Interpolate to the new grid along Z
1329  zRelative = ( ( static_cast<proshade_double> ( zIt ) * static_cast< proshade_double > ( newZSample ) ) - ( static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample ) ) ) / static_cast< proshade_double > ( ( static_cast<proshade_double> ( zTop ) * static_cast<proshade_double> ( oldZSample ) ) - ( static_cast<proshade_double> ( zBottom ) * static_cast<proshade_double> ( oldZSample ) ) );
1330  newMap[newMapIndex] = ( c0.at(3) * ( 1.0 - zRelative ) ) + ( c1.at(3) * zRelative );
1331  }
1332  }
1333  }
1334 
1335  //================================================ Delete old map and allocate new memory
1336  delete[] map;
1337  map = new proshade_double [newXDim * newYDim * newZDim];
1338 
1339  //================================================ Copy map
1340  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ )
1341  {
1342  map[iter] = newMap[iter];
1343  }
1344 
1345  //================================================ Release memory
1346  delete[] newMap;
1347 
1348  //================================================ Define change in indices and return it
1349  corrs[0] = static_cast< proshade_single > ( newXDim - xDim );
1350  corrs[1] = static_cast< proshade_single > ( newYDim - yDim );
1351  corrs[2] = static_cast< proshade_single > ( newZDim - zDim );
1352  corrs[3] = static_cast< proshade_single > ( newXDim ) * static_cast< proshade_single > ( newXSample );
1353  corrs[4] = static_cast< proshade_single > ( newYDim ) * static_cast< proshade_single > ( newYSample );
1354  corrs[5] = static_cast< proshade_single > ( newZDim ) * static_cast< proshade_single > ( newZSample );
1355 
1356  //======================================== Done
1357  return ;
1358 
1359 }
1360 
1376 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 )
1377 {
1378  //================================================ Sanity check - the resolution needs to be set
1379  if ( resolution <= 0.0f )
1380  {
1381  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." );
1382  }
1383 
1384  //================================================ Initialise variables
1385  proshade_unsign newXDim = static_cast<proshade_unsign> ( ProSHADE_internal_mapManip::myRound ( xAngs / ( resolution / 2.0f ) ) );
1386  proshade_unsign newYDim = static_cast<proshade_unsign> ( ProSHADE_internal_mapManip::myRound ( yAngs / ( resolution / 2.0f ) ) );
1387  proshade_unsign newZDim = static_cast<proshade_unsign> ( ProSHADE_internal_mapManip::myRound ( zAngs / ( resolution / 2.0f ) ) );
1388 
1389  if ( newXDim % 2 != 0 ) { newXDim += 1; }
1390  if ( newYDim % 2 != 0 ) { newYDim += 1; }
1391  if ( newZDim % 2 != 0 ) { newZDim += 1; }
1392 
1393  proshade_signed preXChange, preYChange, preZChange;
1394  if ( ( xDimS % 2 ) == 0 ) { preXChange = static_cast< proshade_signed > ( std::ceil ( ( static_cast<proshade_signed> ( xDimS ) - static_cast<proshade_signed> ( newXDim ) ) / 2 ) ); }
1395  else { preXChange = static_cast< proshade_signed > ( std::floor ( ( static_cast<proshade_signed> ( xDimS ) - static_cast<proshade_signed> ( newXDim ) ) / 2 ) ); }
1396  if ( ( yDimS % 2 ) == 0 ) { preYChange = static_cast< proshade_signed > ( std::ceil ( ( static_cast<proshade_signed> ( yDimS ) - static_cast<proshade_signed> ( newYDim ) ) / 2 ) ); }
1397  else { preYChange = static_cast< proshade_signed > ( std::floor ( ( static_cast<proshade_signed> ( yDimS ) - static_cast<proshade_signed> ( newYDim ) ) / 2 ) ); }
1398  if ( ( zDimS % 2 ) == 0 ) { preZChange = static_cast< proshade_signed > ( std::ceil ( ( static_cast<proshade_signed> ( zDimS ) - static_cast<proshade_signed> ( newZDim ) ) / 2 ) ); }
1399  else { preZChange = static_cast< proshade_signed > ( std::floor ( ( static_cast<proshade_signed> ( zDimS ) - static_cast<proshade_signed> ( newZDim ) ) / 2 ) ); }
1400 
1401  proshade_signed postXChange = static_cast<proshade_signed> ( xDimS ) - ( preXChange + static_cast<proshade_signed> ( newXDim ) );
1402  proshade_signed postYChange = static_cast<proshade_signed> ( yDimS ) - ( preYChange + static_cast<proshade_signed> ( newYDim ) );
1403  proshade_signed postZChange = static_cast<proshade_signed> ( zDimS ) - ( preZChange + static_cast<proshade_signed> ( newZDim ) );
1404 
1405  proshade_unsign origSizeArr = 0, newSizeArr = 0;
1406  proshade_double normFactor = static_cast<proshade_double> ( xDimS * yDimS * zDimS );
1407 
1408  //================================================ Manage memory
1409  fftw_complex *origMap, *fCoeffs, *newFCoeffs, *newMap;
1410  fftw_plan planForwardFourier, planBackwardRescaledFourier;
1411  allocateResolutionFourierMemory ( origMap, fCoeffs, newFCoeffs, newMap, planForwardFourier, planBackwardRescaledFourier,
1412  xDimS, yDimS, zDimS, newXDim, newYDim, newZDim );
1413 
1414  //================================================ Fill maps with data and zeroes
1415  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDimS * yDimS * zDimS ); iter++ ) { origMap[iter][0] = map[iter]; origMap[iter][1] = 0.0; }
1416  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ ) { newFCoeffs[iter][0] = 0.0; newFCoeffs[iter][1] = 0.0; }
1417 
1418  //================================================ Get the Fourier coeffs
1419  fftw_execute ( planForwardFourier );
1420 
1421  //================================================ Change the order of Fourier coefficients
1422  changeFourierOrder ( fCoeffs, static_cast< proshade_signed > ( xDimS ), static_cast< proshade_signed > ( yDimS ), static_cast< proshade_signed > ( zDimS ), true );
1423 
1424  //================================================ Re-sample the coefficients by removing high frequencies or adding these with 0 values
1425  for ( proshade_unsign xIt = 0; xIt < newXDim; xIt++ )
1426  {
1427  for ( proshade_unsign yIt = 0; yIt < newYDim; yIt++ )
1428  {
1429  for ( proshade_unsign zIt = 0; zIt < newZDim; zIt++ )
1430  {
1431  //==================================== Find the array positions
1432  origSizeArr = ( ( zIt + static_cast< proshade_unsign > ( preZChange ) ) + zDimS *
1433  ( ( yIt + static_cast< proshade_unsign > ( preYChange ) ) + yDimS *
1434  ( xIt + static_cast< proshade_unsign > ( preXChange ) ) ) );
1435  newSizeArr = zIt + newZDim * ( yIt + newYDim * xIt );
1436 
1437  //==================================== If original coefficient for this new coefficient position exists, copy
1438  if ( ( ( -1 < static_cast< proshade_signed > ( xIt ) + preXChange ) && ( -1 < static_cast<proshade_signed> ( yIt ) + preYChange ) && ( -1 < static_cast<proshade_signed> ( zIt ) + preZChange ) ) &&
1439  ( ( xIt < newXDim + static_cast<proshade_unsign> ( postXChange ) ) && ( yIt < newYDim + static_cast<proshade_unsign> ( postYChange ) ) && ( zIt < newZDim + static_cast<proshade_unsign> ( postZChange ) ) ) )
1440  {
1441  //================================ Copy the Fourier coeff
1442  newFCoeffs[newSizeArr][0] = fCoeffs[origSizeArr][0] / normFactor;
1443  newFCoeffs[newSizeArr][1] = fCoeffs[origSizeArr][1] / normFactor;
1444  }
1445  }
1446  }
1447  }
1448 
1449  //================================================ Change the order of the re-sampled Fourier coefficients
1450  changeFourierOrder ( newFCoeffs, static_cast< proshade_signed > ( newXDim ), static_cast< proshade_signed > ( newYDim ), static_cast< proshade_signed > ( newZDim ), false );
1451 
1452  //================================================ Get the new map from the re-sized Fourier coefficients
1453  fftw_execute ( planBackwardRescaledFourier );
1454 
1455  //================================================ Delete the old map and create a new, re-sized one. Then copy the new map values into this new map memory.
1456  delete map;
1457  map = new proshade_double [newXDim * newYDim * newZDim];
1458  ProSHADE_internal_misc::checkMemoryAllocation ( map, __FILE__, __LINE__, __func__ );
1459  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ ) { map[iter] = newMap[iter][0]; }
1460 
1461  //================================================ Release memory
1462  releaseResolutionFourierMemory ( origMap, fCoeffs, newFCoeffs, newMap, planForwardFourier, planBackwardRescaledFourier );
1463 
1464  //================================================ Define change in indices and return it
1465  corrs[0] = static_cast< proshade_single > ( newXDim ) - static_cast< proshade_single > ( xDimS );
1466  corrs[1] = static_cast< proshade_single > ( newYDim ) - static_cast< proshade_single > ( yDimS );
1467  corrs[2] = static_cast< proshade_single > ( newZDim ) - static_cast< proshade_single > ( zDimS );
1468  corrs[3] = static_cast< proshade_single > ( newXDim ) * static_cast< proshade_single > ( resolution / 2.0f );
1469  corrs[4] = static_cast< proshade_single > ( newYDim ) * static_cast< proshade_single > ( resolution / 2.0f );
1470  corrs[5] = static_cast< proshade_single > ( newZDim ) * static_cast< proshade_single > ( resolution / 2.0f );
1471 
1472 
1473  //======================================== Done
1474  return ;
1475 
1476 }
1477 
1496 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 )
1497 {
1498  //================================================ Initialise memory
1499  origMap = new fftw_complex [xDimOld * yDimOld * zDimOld];
1500  fCoeffs = new fftw_complex [xDimOld * yDimOld * zDimOld];
1501  newFCoeffs = new fftw_complex [xDimNew * yDimNew * zDimNew];
1502  newMap = new fftw_complex [xDimNew * yDimNew * zDimNew];
1503 
1504  //================================================ Check memory allocation
1505  ProSHADE_internal_misc::checkMemoryAllocation ( origMap, __FILE__, __LINE__, __func__ );
1506  ProSHADE_internal_misc::checkMemoryAllocation ( fCoeffs, __FILE__, __LINE__, __func__ );
1507  ProSHADE_internal_misc::checkMemoryAllocation ( newFCoeffs, __FILE__, __LINE__, __func__ );
1508  ProSHADE_internal_misc::checkMemoryAllocation ( newMap, __FILE__, __LINE__, __func__ );
1509 
1510  //================================================ Create plans
1511  planForwardFourier = fftw_plan_dft_3d ( static_cast< int > ( xDimOld ), static_cast< int > ( yDimOld ), static_cast< int > ( zDimOld ), origMap, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
1512  planBackwardRescaledFourier = fftw_plan_dft_3d ( static_cast< int > ( xDimNew ), static_cast< int > ( yDimNew ), static_cast< int > ( zDimNew ), newFCoeffs, newMap, FFTW_BACKWARD, FFTW_ESTIMATE );
1513 
1514  //================================================ Done
1515  return ;
1516 
1517 }
1518 
1530 void ProSHADE_internal_mapManip::releaseResolutionFourierMemory ( fftw_complex*& origMap, fftw_complex*& fCoeffs, fftw_complex*& newFCoeffs, fftw_complex*& newMap, fftw_plan& planForwardFourier, fftw_plan& planBackwardRescaledFourier )
1531 {
1532  //================================================ Delete the FFTW plans
1533  fftw_destroy_plan ( planForwardFourier );
1534  fftw_destroy_plan ( planBackwardRescaledFourier );
1535 
1536  //================================================ Delete the complex arrays
1537  delete[] origMap;
1538  delete[] fCoeffs;
1539  delete[] newFCoeffs;
1540  delete[] newMap;
1541 
1542  //================================================ Done
1543  return ;
1544 
1545 }
1546 
1559 void ProSHADE_internal_mapManip::changeFourierOrder ( fftw_complex*& fCoeffs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, bool negativeFirst )
1560 {
1561  //================================================ Initialise local variables
1562  proshade_signed h = 0, k = 0, l = 0, origSizeArr = 0, newSizeArr = 0;
1563  proshade_signed xSeq1FreqStart, ySeq1FreqStart, zSeq1FreqStart, xSeq2FreqStart, ySeq2FreqStart, zSeq2FreqStart;
1564 
1565  //================================================ Find the positive and negative indices cot-offs
1566  if ( negativeFirst )
1567  {
1568  if ( ( xDim % 2 ) == 0 ) { xSeq1FreqStart = xDim / 2; xSeq2FreqStart = xDim / 2; } else { xSeq1FreqStart = (xDim / 2) + 1; xSeq2FreqStart = xDim / 2; }
1569  if ( ( yDim % 2 ) == 0 ) { ySeq1FreqStart = yDim / 2; ySeq2FreqStart = yDim / 2; } else { ySeq1FreqStart = (yDim / 2) + 1; ySeq2FreqStart = yDim / 2; }
1570  if ( ( zDim % 2 ) == 0 ) { zSeq1FreqStart = zDim / 2; zSeq2FreqStart = zDim / 2; } else { zSeq1FreqStart = (zDim / 2) + 1; zSeq2FreqStart = zDim / 2; }
1571  }
1572  else
1573  {
1574  if ( ( xDim % 2 ) == 0 ) { xSeq1FreqStart = xDim / 2; xSeq2FreqStart = xDim / 2; } else { xSeq1FreqStart = (xDim / 2); xSeq2FreqStart = xDim / 2 + 1; }
1575  if ( ( yDim % 2 ) == 0 ) { ySeq1FreqStart = yDim / 2; ySeq2FreqStart = yDim / 2; } else { ySeq1FreqStart = (yDim / 2); ySeq2FreqStart = yDim / 2 + 1; }
1576  if ( ( zDim % 2 ) == 0 ) { zSeq1FreqStart = zDim / 2; zSeq2FreqStart = zDim / 2; } else { zSeq1FreqStart = (zDim / 2); zSeq2FreqStart = zDim / 2 + 1; }
1577  }
1578 
1579  //================================================ Allocate helper array memory
1580  fftw_complex *hlpFCoeffs = new fftw_complex [xDim * yDim * zDim];
1581  ProSHADE_internal_misc::checkMemoryAllocation ( hlpFCoeffs, __FILE__, __LINE__, __func__ );
1582 
1583  //================================================ Change the coefficients order
1584  for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1585  {
1586  //============================================ Find x frequency
1587  if ( xIt < xSeq1FreqStart ) { h = xIt + xSeq2FreqStart; } else { h = xIt - xSeq1FreqStart; }
1588  for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1589  {
1590  //======================================== Find y frequency
1591  if ( yIt < ySeq1FreqStart ) { k = yIt + ySeq2FreqStart; } else { k = yIt - ySeq1FreqStart; }
1592 
1593  for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1594  {
1595  //==================================== Find z frequency
1596  if ( zIt < zSeq1FreqStart ) { l = zIt + zSeq2FreqStart; } else { l = zIt - zSeq1FreqStart; }
1597 
1598  //==================================== Find array positions
1599  newSizeArr = l + zDim * ( k + yDim * h );
1600  origSizeArr = zIt + zDim * ( yIt + yDim * xIt );
1601 
1602  //==================================== Copy vals
1603  hlpFCoeffs[newSizeArr][0] = fCoeffs[origSizeArr][0];
1604  hlpFCoeffs[newSizeArr][1] = fCoeffs[origSizeArr][1];
1605  }
1606  }
1607  }
1608 
1609  //================================================ Copy the helper array to the input Fourier coefficients array
1610  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]; }
1611 
1612  //================================================ Release helper array memory
1613  delete[] hlpFCoeffs;
1614 
1615  //================================================ Done
1616  return ;
1617 
1618 }
1619 
1631 void ProSHADE_internal_mapManip::removeMapPhase ( fftw_complex*& mapCoeffs, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim )
1632 {
1633  //================================================ Set local variables
1634  proshade_double real, imag, mag, phase;
1635  proshade_unsign arrayPos = 0;
1636  proshade_double normFactor = static_cast<proshade_double> ( xDim * yDim * zDim );
1637 
1638  //================================================ Iterate through the map
1639  for ( proshade_unsign uIt = 0; uIt < xDim; uIt++ )
1640  {
1641  for ( proshade_unsign vIt = 0; vIt < yDim; vIt++ )
1642  {
1643  for ( proshade_unsign wIt = 0; wIt < zDim; wIt++ )
1644  {
1645  //==================================== Var init
1646  arrayPos = wIt + zDim * ( vIt + yDim * uIt );
1647  real = mapCoeffs[arrayPos][0];
1648  imag = mapCoeffs[arrayPos][1];
1649 
1650  //==================================== Get magnitude and phase with mask parameters
1651  mag = std::sqrt ( (real*real) + (imag*imag) );;
1652  phase = 0.0; // This would be std::atan2 ( imag, real ); - but here we remove the phase.
1653 
1654  //==================================== Save the phaseless data
1655  mapCoeffs[arrayPos][0] = ( mag * cos(phase) ) / normFactor;
1656  mapCoeffs[arrayPos][1] = ( mag * sin(phase) ) / normFactor;
1657  }
1658  }
1659  }
1660 
1661  //================================================ Done
1662  return ;
1663 
1664 }
1665 
1680 void ProSHADE_internal_mapManip::getFakeHalfMap ( proshade_double*& map, proshade_double*& fakeHalfMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_signed fakeMapKernel )
1681 {
1682  //================================================ Set local variables
1683  proshade_signed xDim = static_cast< proshade_signed > ( xDimS );
1684  proshade_signed yDim = static_cast< proshade_signed > ( yDimS );
1685  proshade_signed zDim = static_cast< proshade_signed > ( zDimS );
1686  proshade_signed currentPos, neighArrPos, neighXPos, neighYPos, neighZPos;
1687  proshade_double neighSum;
1688  proshade_double neighCount = pow ( ( ( fakeMapKernel * 2 ) + 1 ), 3.0 ) - 1.0;
1689 
1690  //================================================ Blur the coeffs
1691  for ( proshade_signed uIt = 0; uIt < xDim; uIt++ )
1692  {
1693  for ( proshade_signed vIt = 0; vIt < yDim; vIt++ )
1694  {
1695  for ( proshade_signed wIt = 0; wIt < zDim; wIt++ )
1696  {
1697  //==================================== Var init
1698  currentPos = wIt + zDim * ( vIt + yDim * uIt );
1699  neighSum = 0.0;
1700 
1701  //==================================== Average neighbours
1702  for ( proshade_signed xCh = -fakeMapKernel; xCh <= +fakeMapKernel; xCh++ )
1703  {
1704  for ( proshade_signed yCh = -fakeMapKernel; yCh <= +fakeMapKernel; yCh++ )
1705  {
1706  for ( proshade_signed zCh = -fakeMapKernel; zCh <= +fakeMapKernel; zCh++ )
1707  {
1708  if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) { continue; }
1709 
1710  //======================== Find the nieghbour peak indices (with periodicity)
1711  neighXPos = uIt + xCh; if ( neighXPos >= xDim ) { neighXPos -= xDim; }; if ( neighXPos < 0 ) { neighXPos += xDim; }
1712  neighYPos = vIt + yCh; if ( neighYPos >= yDim ) { neighYPos -= yDim; }; if ( neighYPos < 0 ) { neighYPos += yDim; }
1713  neighZPos = wIt + zCh; if ( neighZPos >= zDim ) { neighZPos -= zDim; }; if ( neighZPos < 0 ) { neighZPos += zDim; }
1714  neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1715 
1716  //======================== Add to average
1717  neighSum += map[neighArrPos];
1718  }
1719  }
1720  }
1721 
1722  //==================================== Save the average to "fake" half-map
1723  fakeHalfMap[currentPos] = neighSum / neighCount;
1724  }
1725  }
1726  }
1727 
1728  //================================================ Done
1729  return ;
1730 
1731 }
1732 
1745 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 )
1746 {
1747  //================================================ Set local variables
1748  proshade_signed xDim = static_cast< proshade_signed > ( xDimS ), yDim = static_cast< proshade_signed > ( yDimS ), zDim = static_cast< proshade_signed > ( zDimS ), currentPos, neighArrPos, neighXPos, neighYPos, neighZPos, corrIter;
1749  proshade_unsign noCorrVals = static_cast<proshade_unsign> ( pow ( ( ( corrMaskKernel * 2 ) + 1 ), 3 ) );
1750 
1751  //================================================ Alocate memory
1752  proshade_double *origMap = new proshade_double [noCorrVals];
1753  proshade_double *fakeHM = new proshade_double [noCorrVals];
1754 
1755  //================================================ Check memory allocation
1756  ProSHADE_internal_misc::checkMemoryAllocation ( origMap, __FILE__, __LINE__, __func__ );
1757  ProSHADE_internal_misc::checkMemoryAllocation ( fakeHM, __FILE__, __LINE__, __func__ );
1758 
1759  //================================================ Blur the coeffs
1760  for ( proshade_signed uIt = 0; uIt < xDim; uIt++ )
1761  {
1762  for ( proshade_signed vIt = 0; vIt < yDim; vIt++ )
1763  {
1764  for ( proshade_signed wIt = 0; wIt < zDim; wIt++ )
1765  {
1766  //==================================== Var init
1767  currentPos = wIt + zDim * ( vIt + yDim * uIt );
1768  corrIter = 0;
1769 
1770  //==================================== Average neighbours
1771  for ( proshade_signed xCh = -corrMaskKernel; xCh <= +corrMaskKernel; xCh++ )
1772  {
1773  for ( proshade_signed yCh = -corrMaskKernel; yCh <= +corrMaskKernel; yCh++ )
1774  {
1775  for ( proshade_signed zCh = -corrMaskKernel; zCh <= +corrMaskKernel; zCh++ )
1776  {
1777  //======================== Find the nieghbour peak indices (with periodicity)
1778  neighXPos = uIt + xCh; if ( neighXPos >= xDim ) { neighXPos -= xDim; }; if ( neighXPos < 0 ) { neighXPos += xDim; }
1779  neighYPos = vIt + yCh; if ( neighYPos >= yDim ) { neighYPos -= yDim; }; if ( neighYPos < 0 ) { neighYPos += yDim; }
1780  neighZPos = wIt + zCh; if ( neighZPos >= zDim ) { neighZPos -= zDim; }; if ( neighZPos < 0 ) { neighZPos += zDim; }
1781  neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1782 
1783  //======================== Add to correct arrays
1784  origMap[corrIter] = map[neighArrPos];
1785  fakeHM[corrIter] = fakeHalfMap[neighArrPos];
1786  corrIter += 1;
1787  }
1788  }
1789  }
1790 
1791  //==================================== Save the correlation comparison result
1792  correlationMask[currentPos] = ProSHADE_internal_maths::pearsonCorrCoeff ( origMap, fakeHM, noCorrVals );
1793  }
1794  }
1795  }
1796 
1797  //================================================ Done
1798  return ;
1799 
1800 }
1801 
1816 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 )
1817 {
1818  //================================================ Compute
1819  proshade_single ret = static_cast< proshade_single > ( ProSHADE_internal_mapManip::myRound ( std::max ( dist / ( xAngs / static_cast<proshade_single> ( xDim ) ),
1820  std::max ( dist / ( yAngs / static_cast<proshade_single> ( yDim ) ),
1821  dist / ( zAngs / static_cast<proshade_single> ( zDim ) ) ) ) ) );
1822 
1823  //================================================ Done
1824  return ( ret );
1825 
1826 }
1827 
1839 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 )
1840 {
1841  //================================================ Initialise variables
1842  proshade_double* hlpMap = new proshade_double[xDim * yDim * zDim];
1843  proshade_signed addSurroundingPoints = static_cast< proshade_signed > ( std::max ( 3L, static_cast<proshade_signed> ( std::ceil ( getIndicesFromAngstroms( static_cast< proshade_unsign > ( xDim ), static_cast< proshade_unsign > ( yDim ), static_cast< proshade_unsign > ( zDim ), xAngs, yAngs, zAngs, static_cast< proshade_single > ( std::max( xAngs, std::max( yAngs, zAngs ) ) * 0.1f ) ) ) ) ) );
1844  proshade_signed currPos, neighXPos, neighYPos, neighZPos, neighArrPos;
1845 
1846  //================================================ Check memory allocation
1847  ProSHADE_internal_misc::checkMemoryAllocation ( hlpMap, __FILE__, __LINE__, __func__ );
1848 
1849  //================================================ Copy the mask map
1850  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { hlpMap[iter] = mask[iter]; }
1851 
1852  //================================================ Repeat as many times as needed
1853  for ( proshade_signed it = 0; it < addSurroundingPoints; it++ )
1854  {
1855  //============================================ For each x, y and z
1856  for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1857  {
1858  for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1859  {
1860  for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1861  {
1862  //================================ Current position
1863  currPos = zIt + zDim * ( yIt + yDim * xIt );
1864 
1865  //================================ If zero, ignore
1866  if ( hlpMap[currPos] < static_cast< proshade_double > ( maskThres ) ) { continue; }
1867 
1868  //================================ Check neighbours
1869  for ( proshade_signed xCh = -1; xCh <= +1; xCh++ )
1870  {
1871  for ( proshade_signed yCh = -1; yCh <= +1; yCh++ )
1872  {
1873  for ( proshade_signed zCh = -1; zCh <= +1; zCh++ )
1874  {
1875  if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) { continue; }
1876 
1877  //==================== Find the nieghbour indices (without periodicity)
1878  neighXPos = xIt + xCh; if ( neighXPos < 0 ) { continue; } if ( neighXPos >= xDim ) { continue; }
1879  neighYPos = yIt + yCh; if ( neighYPos < 0 ) { continue; } if ( neighYPos >= yDim ) { continue; }
1880  neighZPos = zIt + zCh; if ( neighZPos < 0 ) { continue; } if ( neighZPos >= zDim ) { continue; }
1881  neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1882 
1883  //==================== Add to mask if this point is below it (as it is a neighbour to a point which is part of the mask)
1884  if ( hlpMap[neighArrPos] < static_cast< proshade_double > ( maskThres ) ) { mask[neighArrPos] = static_cast< proshade_double > ( maskThres ); }
1885  }
1886  }
1887  }
1888  }
1889  }
1890  }
1891 
1892  //============================================ Now copy the updated mask to the temporary helper, unless last iteration
1893  for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { hlpMap[iter] = mask[iter]; }
1894  }
1895 
1896  //================================================ Release memory
1897  delete[] hlpMap;
1898 
1899  //================================================ Done
1900  return ;
1901 
1902 }
1903 
1916 void ProSHADE_internal_mapManip::beautifyBoundaries ( proshade_signed*& bounds, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_signed boundsDiffThres )
1917 {
1918  //================================================ If extra bounds space pushed the bounds over the physical map, freely add up to 10 indices over the current value
1919  while ( bounds[1] >= static_cast<proshade_signed> ( xDim ) ) { xDim += 10; }
1920  while ( bounds[3] >= static_cast<proshade_signed> ( yDim ) ) { yDim += 10; }
1921  while ( bounds[5] >= static_cast<proshade_signed> ( zDim ) ) { zDim += 10; }
1922 
1923  //================================================ Find if better bouds exist in terms of prime numbers
1924  proshade_signed addToX = betterClosePrimeFactors ( bounds[1] - bounds[0] + 1, static_cast< proshade_signed > ( xDim ) );
1925  proshade_signed addToY = betterClosePrimeFactors ( bounds[3] - bounds[2] + 1, static_cast< proshade_signed > ( yDim ) );
1926  proshade_signed addToZ = betterClosePrimeFactors ( bounds[5] - bounds[4] + 1, static_cast< proshade_signed > ( zDim ) );
1927 
1928  //================================================ Figure if similar sizes should not be forced to be the same
1929  proshade_signed XtoY = std::abs ( addToX - addToY );
1930  proshade_signed XtoZ = std::abs ( addToX - addToZ );
1931  proshade_signed YtoZ = std::abs ( addToY - addToZ );
1932 
1933  if ( ( ( XtoY < boundsDiffThres ) && ( XtoZ < boundsDiffThres ) ) ||
1934  ( ( XtoY < boundsDiffThres ) && ( YtoZ < boundsDiffThres ) ) ||
1935  ( ( XtoZ < boundsDiffThres ) && ( YtoZ < boundsDiffThres ) ) )
1936  {
1937  //============================================ Dealing with chain ( a <= b <= c type ) where at least two dista are smaller than threshold
1938  proshade_signed maxSize = std::max ( addToX, std::max ( addToY, addToZ ) );
1939  addToX = maxSize;
1940  addToY = maxSize;
1941  addToZ = maxSize;
1942  }
1943  else
1944  {
1945  //============================================ Only two may be similar enough
1946  if ( XtoY <= boundsDiffThres )
1947  {
1948  proshade_signed maxSize = std::max ( addToX, addToY );
1949  addToX = maxSize;
1950  addToY = maxSize;
1951  }
1952  if ( XtoZ <= boundsDiffThres )
1953  {
1954  proshade_signed maxSize = std::max ( addToX, addToZ );
1955  addToX = maxSize;
1956  addToZ = maxSize;
1957  }
1958  if ( YtoZ <= boundsDiffThres )
1959  {
1960  proshade_signed maxSize = std::max ( addToY, addToZ );
1961  addToY = maxSize;
1962  addToZ = maxSize;
1963  }
1964  }
1965 
1966  //================================================ Add the extra space appropriately
1967  distributeSpaceToBoundaries ( bounds[0], bounds[1], ( bounds[1] - bounds[0] + 1 ), addToX );
1968  distributeSpaceToBoundaries ( bounds[2], bounds[3], ( bounds[3] - bounds[2] + 1 ), addToY );
1969  distributeSpaceToBoundaries ( bounds[4], bounds[5], ( bounds[5] - bounds[4] + 1 ), addToZ );
1970 
1971  //================================================ Done
1972  return ;
1973 
1974 }
1975 
1986 proshade_signed ProSHADE_internal_mapManip::betterClosePrimeFactors ( proshade_signed fromRange, proshade_signed toRange )
1987 {
1988  //================================================ Initialise variables
1989  proshade_signed ret = fromRange;
1990  std::vector < proshade_signed > posibles, hlp;
1991  proshade_signed sum;
1992 
1993  //================================================ Find the sums of prime factors for each number in the whole range
1994  for ( proshade_signed iter = fromRange; iter < toRange; iter++ )
1995  {
1996  sum = 0;
1998  for ( proshade_unsign i = 0; i < static_cast<proshade_unsign> ( hlp.size() ); i++ ) { sum += hlp.at(i); }
1999  hlp.clear ( );
2000  ProSHADE_internal_misc::addToSignedVector ( &posibles, sum );
2001  }
2002 
2003  //================================================ Find a better number
2004  for ( proshade_signed iter = fromRange; iter < toRange; iter++ )
2005  {
2006  //============================================ Ignore odd numbers
2007  if ( iter %2 != 0 ) { continue; }
2008 
2009  //============================================ Better number needs to have lower prime factor sum, but also needs not to be too far
2010  if ( posibles.at( static_cast< size_t > ( iter - fromRange ) ) < ( posibles.at( static_cast< size_t > ( ret - fromRange ) ) - ( iter - ret ) ) ) { ret = iter; }
2011  }
2012 
2013  //================================================ In the case of large prime number, just add one for even number
2014  if ( ( ret % 2 != 0 ) && ( ret < ( toRange - 1 ) ) ) { ret += 1; }
2015 
2016  //================================================ Done
2017  return ( ret );
2018 
2019 }
2020 
2032 void ProSHADE_internal_mapManip::distributeSpaceToBoundaries ( proshade_signed& minBound, proshade_signed& maxBound, proshade_signed oldBoundRange, proshade_signed newBoundRange )
2033 {
2034  //================================================ Only bother when sometings is to be added
2035  if ( newBoundRange > oldBoundRange )
2036  {
2037  //============================================ How much are we adding?
2038  proshade_signed distributeThis = newBoundRange - oldBoundRange;
2039 
2040  while ( distributeThis != 0 )
2041  {
2042  minBound -= 1;
2043  distributeThis -= 1;
2044 
2045  if ( distributeThis != 0 )
2046  {
2047  maxBound += 1;
2048  distributeThis -= 1;
2049  }
2050  }
2051  }
2052 
2053  //================================================ Done
2054  return ;
2055 
2056 }
2057 
2082 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_unsign yDimIndices, proshade_unsign zDimIndices, proshade_unsign origXDimIndices, proshade_unsign origYDimIndices, proshade_unsign origZDimIndices, proshade_double*& newMap, proshade_double* origMap )
2083 {
2084  //================================================ Initialise variables
2085  proshade_signed newMapIndex, oldMapIndex, oldX, oldY, oldZ, newX, newY, newZ;
2086 
2087  //=============================================== For all values in the new map
2088  for ( proshade_signed xIt = xFrom; xIt <= xTo; xIt++ )
2089  {
2090  //============================================ Index position init
2091  newX = ( xIt - xFrom );
2092  oldX = ( newX + ( xFrom - origXFrom ) );
2093 
2094  for ( proshade_signed yIt = yFrom; yIt <= yTo; yIt++ )
2095  {
2096  //======================================== Index position init
2097  newY = ( yIt - yFrom );
2098  oldY = ( newY + ( yFrom - origYFrom ) );
2099 
2100  for ( proshade_signed zIt = zFrom; zIt <= zTo; zIt++ )
2101  {
2102  //==================================== Index position init
2103  newZ = ( zIt - zFrom );
2104  oldZ = ( newZ + ( zFrom - origZFrom ) );
2105 
2106  //==================================== Index arrays
2107  newMapIndex = newZ + static_cast< proshade_signed > ( zDimIndices ) * ( newY + static_cast< proshade_signed > ( yDimIndices ) * newX );
2108  oldMapIndex = oldZ + static_cast< proshade_signed > ( origZDimIndices ) * ( oldY + static_cast< proshade_signed > ( origYDimIndices ) * oldX );
2109 
2110  //============================ Check if we are adding new point
2111  if ( ( ( oldX < 0 ) || ( oldX >= static_cast< proshade_signed > ( origXDimIndices ) ) ) ||
2112  ( ( oldY < 0 ) || ( oldY >= static_cast< proshade_signed > ( origYDimIndices ) ) ) ||
2113  ( ( oldZ < 0 ) || ( oldZ >= static_cast< proshade_signed > ( origZDimIndices ) ) ) )
2114  {
2115  //================================ Yes, this is a new point, no known value for it
2116  newMap[newMapIndex] = 0.0;
2117  continue;
2118  }
2119 
2120  //==================================== If we got here, this should be within the old map, so just copy value
2121  newMap[newMapIndex] = origMap[oldMapIndex];
2122  }
2123  }
2124  }
2125 
2126  //================================================ Done
2127  return ;
2128 
2129 }
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:1139
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:72
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:444
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:2032
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:1032
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:240
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:932
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:811
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:573
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:1007
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:1175
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:1745
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:1559
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:1530
ProSHADE_internal_mapManip::removeWaters
void removeWaters(gemmi::Structure *pdbFile, bool firstModel)
This function removed all waters from PDB input file.
Definition: ProSHADE_mapManip.cpp:504
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:246
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:157
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:1816
ProSHADE_internal_mapManip::getNonZeroBounds
void getNonZeroBounds(proshade_double *map, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_signed *&ret)
Function for finding the map boundaries enclosing positive only values.
Definition: ProSHADE_mapManip.cpp:1082
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:1631
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:1496
ProSHADE_internal_maths::primeFactorsDecomp
std::vector< proshade_signed > primeFactorsDecomp(proshade_signed number)
Function to find prime factors of an integer.
Definition: ProSHADE_maths.cpp:1707
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:1986
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::beautifyBoundaries
void beautifyBoundaries(proshade_signed *&bounds, proshade_unsign xDim, proshade_unsign yDim, proshade_unsign zDim, proshade_signed boundsDiffThres)
Function for modifying boundaries to a mathematically more pleasant values.
Definition: ProSHADE_mapManip.cpp:1916
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:296
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:644
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:67
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:381
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:149
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:1680
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:1376
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_unsign yDimIndices, proshade_unsign zDimIndices, proshade_unsign origXDimIndices, proshade_unsign origYDimIndices, proshade_unsign origZDimIndices, proshade_double *&newMap, proshade_double *origMap)
This function copies an old map to a new map with different boundaries.
Definition: ProSHADE_mapManip.cpp:2082
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:1839
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:763