33 #if __cplusplus >= 201103L
34 return ( std::round ( x ) );
36 return ( round ( x ) );
60 bool firstAtom =
true;
63 if ( pdbFile.models.size() > 0 )
66 for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); sIt++ )
69 if ( firstModel && ( sIt != 0 ) ) {
break; }
72 gemmi::Model model = pdbFile.models.at(sIt);
75 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model.chains.size() ); mIt++ )
78 gemmi::Chain chain = model.chains.at(mIt);
81 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain.residues.size() ); rIt++ )
84 gemmi::Residue residue = chain.residues.at(rIt);
87 for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue.atoms.size() ); aIt++ )
90 gemmi::Atom atom = residue.atoms.at(aIt);
93 if ( atom.is_hydrogen() ) {
continue; }
98 *xTo =
static_cast<proshade_single
> ( atom.pos.x );
99 *xFrom =
static_cast<proshade_single
> ( atom.pos.x );
100 *yTo =
static_cast<proshade_single
> ( atom.pos.y );
101 *yFrom =
static_cast<proshade_single
> ( atom.pos.y );
102 *zTo =
static_cast<proshade_single
> ( atom.pos.z );
103 *zFrom =
static_cast<proshade_single
> ( atom.pos.z );
108 if (
static_cast<proshade_single
> ( atom.pos.x ) > *xTo ) { *xTo =
static_cast<proshade_single
> ( atom.pos.x ); }
109 if (
static_cast<proshade_single
> ( atom.pos.x ) < *xFrom ) { *xFrom =
static_cast<proshade_single
> ( atom.pos.x ); }
110 if (
static_cast<proshade_single
> ( atom.pos.y ) > *yTo ) { *yTo =
static_cast<proshade_single
> ( atom.pos.y ); }
111 if (
static_cast<proshade_single
> ( atom.pos.y ) < *yFrom ) { *yFrom =
static_cast<proshade_single
> ( atom.pos.y ); }
112 if (
static_cast<proshade_single
> ( atom.pos.z ) > *zTo ) { *zTo =
static_cast<proshade_single
> ( atom.pos.z ); }
113 if (
static_cast<proshade_single
> ( atom.pos.z ) < *zFrom ) { *zFrom =
static_cast<proshade_single
> ( atom.pos.z ); }
122 std::stringstream hlpSS;
123 hlpSS <<
"Found 0 models in input file " << pdbFile.name <<
".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
124 throw ProSHADE_exception (
"Found no model in co-ordinate file.",
"EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
145 proshade_double totAtoms = 0.0;
151 if ( pdbFile.models.size() > 0 )
154 for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); sIt++ )
157 gemmi::Model model = pdbFile.models.at(sIt);
160 if ( firstModel && ( sIt != 0 ) ) {
break; }
163 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model.chains.size() ); mIt++ )
166 gemmi::Chain chain = model.chains.at(mIt);
169 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain.residues.size() ); rIt++ )
172 gemmi::Residue residue = chain.residues.at(rIt);
175 for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue.atoms.size() ); aIt++ )
178 gemmi::Atom atom = residue.atoms.at(aIt);
181 *xCom += atom.pos.x * atom.element.weight();
182 *yCom += atom.pos.y * atom.element.weight();
183 *zCom += atom.pos.z * atom.element.weight();
184 totAtoms += atom.element.weight();
192 std::stringstream hlpSS;
193 hlpSS <<
"Found 0 models in input file " << pdbFile.name <<
".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
194 throw ProSHADE_exception (
"Found no model in co-ordinate file.",
"EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
225 void ProSHADE_internal_mapManip::findMAPCOMValues ( proshade_double* map, proshade_double *xCom, proshade_double *yCom, proshade_double *zCom, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo )
228 proshade_double totDensity = 0.0;
232 proshade_signed arrPos = 0;
233 proshade_double xSampRate = xAngs /
static_cast<proshade_single
> ( xTo - xFrom );
234 proshade_double ySampRate = yAngs /
static_cast<proshade_single
> ( yTo - yFrom );
235 proshade_double zSampRate = zAngs /
static_cast<proshade_single
> ( zTo - zFrom );
238 for ( proshade_signed xIt = xFrom; xIt < xTo; xIt++ )
240 for ( proshade_signed yIt = yFrom; yIt < yTo; yIt++ )
242 for ( proshade_signed zIt = zFrom; zIt < zTo; zIt++ )
244 arrPos = (zIt-zFrom) + ( zTo - zFrom ) * ( (yIt-yFrom) + ( yTo - yFrom ) * (xIt-xFrom) );
245 totDensity += map[arrPos];
246 *xCom +=
static_cast<proshade_double
> ( xIt * xSampRate ) * map[arrPos];
247 *yCom +=
static_cast<proshade_double
> ( yIt * ySampRate ) * map[arrPos];
248 *zCom +=
static_cast<proshade_double
> ( zIt * zSampRate ) * map[arrPos];
278 proshade_double yCom, proshade_double zCom,
bool firstModel )
281 proshade_double *rotMat =
new proshade_double[9];
286 proshade_single xTmp, yTmp, zTmp;
289 if ( pdbFile->models.size() > 0 )
292 for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
295 gemmi::Model *model = &pdbFile->models.at(sIt);
298 if ( firstModel && ( sIt != 0 ) ) {
break; }
301 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
304 gemmi::Chain *chain = &model->chains.at(mIt);
307 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
310 gemmi::Residue *residue = &chain->residues.at(rIt);
313 for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
316 gemmi::Atom *atom = &residue->atoms.at(aIt);
319 xTmp = atom->pos.x - xCom;
320 yTmp = atom->pos.y - yCom;
321 zTmp = atom->pos.z - zCom;
324 atom->pos.x = ( xTmp * rotMat[0] ) + ( yTmp * rotMat[1] ) + ( zTmp * rotMat[2] );
325 atom->pos.y = ( xTmp * rotMat[3] ) + ( yTmp * rotMat[4] ) + ( zTmp * rotMat[5] );
326 atom->pos.z = ( xTmp * rotMat[6] ) + ( yTmp * rotMat[7] ) + ( zTmp * rotMat[8] );
329 atom->pos.x = atom->pos.x + xCom;
330 atom->pos.y = atom->pos.y + yCom;
331 atom->pos.z = atom->pos.z + zCom;
339 std::stringstream hlpSS;
340 hlpSS <<
"Found 0 models in input file " << pdbFile->name <<
".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
341 throw ProSHADE_exception (
"Found no model in co-ordinate file.",
"EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
365 if ( pdbFile->models.size() > 0 )
368 for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
371 if ( firstModel && ( sIt != 0 ) ) {
break; }
374 gemmi::Model *model = &pdbFile->models.at(sIt);
377 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
380 gemmi::Chain *chain = &model->chains.at(mIt);
383 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
386 gemmi::Residue *residue = &chain->residues.at(rIt);
389 for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
392 gemmi::Atom *atom = &residue->atoms.at(aIt);
395 atom->pos.x += transX;
396 atom->pos.y += transY;
397 atom->pos.z += transZ;
405 std::stringstream hlpSS;
406 hlpSS <<
"Found 0 models in input file " << pdbFile->name <<
".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
407 throw ProSHADE_exception (
"Found no model in co-ordinate file.",
"EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
428 if ( pdbFile->models.size() > 0 )
431 for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
434 if ( firstModel && ( sIt != 0 ) ) {
break; }
437 gemmi::Model *model = &pdbFile->models.at(sIt);
440 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
443 gemmi::Chain *chain = &model->chains.at(mIt);
446 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
449 gemmi::Residue *residue = &chain->residues.at(rIt);
452 for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
455 gemmi::Atom *atom = &residue->atoms.at(aIt);
458 atom->b_iso = newBFactorValue;
466 std::stringstream hlpSS;
467 hlpSS <<
"Found 0 models in input file " << pdbFile->name <<
".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
468 throw ProSHADE_exception (
"Found no model in co-ordinate file.",
"EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
488 if ( pdbFile->models.size() > 0 )
491 for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
494 if ( firstModel && ( sIt != 0 ) ) {
break; }
497 gemmi::Model *model = &pdbFile->models.at(sIt);
500 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
503 gemmi::Chain *chain = &model->chains.at(mIt);
506 std::vector< proshade_unsign > delVec;
509 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
512 gemmi::Residue *residue = &chain->residues.at(rIt);
515 if ( residue->is_water() )
522 std::sort ( delVec.begin(), delVec.end(), std::greater<int>() );
523 for ( proshade_unsign vecIt = 0; vecIt < static_cast<proshade_unsign> ( delVec.size() ); vecIt++ )
525 chain->residues.erase ( chain->residues.begin() + delVec.at(vecIt) );
532 std::stringstream hlpSS;
533 hlpSS <<
"Found 0 models in input file " << pdbFile->name <<
".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
534 throw ProSHADE_exception (
"Found no model in co-ordinate file.",
"EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
557 if ( pdbFile->models.size() > 0 )
560 for ( proshade_unsign sIt = 0; sIt < static_cast<proshade_unsign> ( pdbFile->models.size() ); sIt++ )
563 if ( firstModel && ( sIt != 0 ) ) {
break; }
566 gemmi::Model *model = &pdbFile->models.at(sIt);
569 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( model->chains.size() ); mIt++ )
572 gemmi::Chain *chain = &model->chains.at(mIt);
575 for ( proshade_unsign rIt = 0; rIt < static_cast<proshade_unsign> ( chain->residues.size() ); rIt++ )
578 gemmi::Residue *residue = &chain->residues.at(rIt);
581 for ( proshade_unsign aIt = 0; aIt < static_cast<proshade_unsign> ( residue->atoms.size() ); aIt++ )
584 gemmi::Atom *atom = &residue->atoms.at(aIt);
587 atom->pos = gemmi::Position ( atom->pos.x + xMov, atom->pos.y + yMov, atom->pos.z + zMov );
596 std::stringstream hlpSS;
597 hlpSS <<
"Found 0 models in input file " << pdbFile->name <<
".\n : This suggests that the input co-ordinate file is\n : corrupted or mis-formatted.";
598 throw ProSHADE_exception (
"Found no model in co-ordinate file.",
"EP00050", __FILE__, __LINE__, __func__, hlpSS.str() );
625 void ProSHADE_internal_mapManip::generateMapFromPDB ( gemmi::Structure pdbFile, proshade_double*& map, proshade_single requestedResolution, proshade_single xCell, proshade_single yCell, proshade_single zCell, proshade_signed* xTo, proshade_signed* yTo, proshade_signed* zTo,
bool forceP1,
bool firstModel )
628 if ( forceP1 ) { pdbFile.cell = gemmi::UnitCell(); }
629 pdbFile.cell.a = xCell;
630 pdbFile.cell.b = yCell;
631 pdbFile.cell.c = zCell;
632 pdbFile.cell.calculate_properties ( );
635 std::string totElString;
636 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); mIt++ )
639 if ( firstModel && ( mIt != 0 ) )
641 std::stringstream hlpSS;
642 hlpSS <<
"!!! ProSHADE WARNING !!! Found multiple models (" << pdbFile.models.size() <<
") in input file " << pdbFile.name <<
", while the settings state that only the first PDB file model should be used. If all models should be used, please supply ProSHADE with the \"-x\" option.";
647 std::string hlpStr = pdbFile.models[mIt].present_elements ( ).to_string<char,std::char_traits<char>,std::allocator<char> >();
648 totElString = totElString + hlpStr;
650 std::bitset<(size_t)gemmi::El::END> present_elems ( totElString );
653 if ( present_elems[
static_cast<int> ( gemmi::El::X )] )
655 throw ProSHADE_exception (
"Found unknown element in input file.",
"EP00051", __FILE__, __LINE__, __func__,
"Gemmi library does not recognise some of the elements in\n : the co-ordinate file. Please check the file for not being\n : corrupted and containing standard elements." );
658 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( present_elems.size() ); elIt++ )
660 if ( present_elems[elIt] && !gemmi::IT92<double>::has (
static_cast<gemmi::El
> ( elIt ) ) )
662 std::stringstream hlpSS;
663 hlpSS <<
"Missing form factor for element " << element_name (
static_cast<gemmi::El
> ( elIt ) );
664 throw ProSHADE_exception ( hlpSS.str().c_str(),
"EP00052", __FILE__, __LINE__, __func__,
"Gemmi library does not have a form factor value for this\n : reported element. Please report this to the author." );
669 double wavelength = 0.1;
670 double energy = gemmi::hc() / wavelength;
673 gemmi::DensityCalculator<gemmi::IT92<double>,
float> dencalc;
675 dencalc.d_min =
static_cast<double> ( requestedResolution );
676 for ( proshade_unsign elIt = 0; elIt < static_cast<proshade_unsign> ( present_elems.size() ); elIt++ ) {
if ( present_elems[elIt] ) { dencalc.addends.set (
static_cast<gemmi::El
> ( elIt ), gemmi::cromer_libermann ( elIt, energy,
nullptr ) ); } }
677 dencalc.set_grid_cell_and_spacegroup ( pdbFile );
680 if ( forceP1 ) { dencalc.grid.spacegroup = &gemmi::get_spacegroup_p1(); }
683 dencalc.grid.data.clear ( );
684 dencalc.grid.set_size_from_spacing ( dencalc.d_min / ( 2.0 * dencalc.rate),
true );
685 for ( proshade_unsign mIt = 0; mIt < static_cast<proshade_unsign> ( pdbFile.models.size() ); mIt++ )
687 if ( firstModel && ( mIt != 0 ) ) {
break; }
688 dencalc.add_model_density_to_grid ( pdbFile.models[mIt] );
689 dencalc.grid.symmetrize ( [](
float a,
float b) {
return a + b; } );
693 const gemmi::Grid<float>& grid = dencalc.grid;
701 map =
new proshade_double [(*xTo) * (*yTo) * (*zTo)];
704 proshade_signed arrPos = 0;
705 for ( proshade_signed uIt = 0; uIt < (*xTo); uIt++ )
707 for ( proshade_signed vIt = 0; vIt < (*yTo); vIt++ )
709 for ( proshade_signed wIt = 0; wIt < (*zTo); wIt++ )
711 arrPos = wIt + (*zTo) * ( vIt + (*yTo) * uIt );
712 map[arrPos] = grid.get_value_q( uIt, vIt, wIt );
744 void ProSHADE_internal_mapManip::moveMapByIndices ( proshade_single* xMov, proshade_single* yMov, proshade_single* zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed* xFrom, proshade_signed* xTo, proshade_signed* yFrom, proshade_signed* yTo, proshade_signed* zFrom, proshade_signed* zTo, proshade_signed* xOrigin, proshade_signed* yOrigin, proshade_signed* zOrigin )
747 proshade_signed xIndMove =
static_cast<proshade_signed
> ( std::floor ( -(*xMov) / ( xAngs / (*xTo - *xFrom) ) ) );
748 proshade_signed yIndMove =
static_cast<proshade_signed
> ( std::floor ( -(*yMov) / ( yAngs / (*yTo - *yFrom) ) ) );
749 proshade_signed zIndMove =
static_cast<proshade_signed
> ( std::floor ( -(*zMov) / ( zAngs / (*zTo - *zFrom) ) ) );
752 *xMov = -( *xMov ) - ( xIndMove * ( xAngs / (*xTo - *xFrom) ) );
753 *yMov = -( *yMov ) - ( yIndMove * ( yAngs / (*yTo - *yFrom) ) );
754 *zMov = -( *zMov ) - ( zIndMove * ( zAngs / (*zTo - *zFrom) ) );
792 void ProSHADE_internal_mapManip::moveMapByFourier ( proshade_double*& map, proshade_single xMov, proshade_single yMov, proshade_single zMov, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim )
795 proshade_signed arrayPos = 0;
796 proshade_signed h, k, l;
797 proshade_double real = 0.0;
798 proshade_double imag = 0.0;
799 proshade_double trCoeffReal, trCoeffImag;
800 proshade_double normFactor =
static_cast<proshade_double
> ( xDim * yDim * zDim );
801 proshade_double exponent = 0.0;
802 proshade_double hlpArrReal;
803 proshade_double hlpArrImag;
806 fftw_complex *fCoeffs =
new fftw_complex [xDim * yDim * zDim];
807 fftw_complex *translatedMap =
new fftw_complex [xDim * yDim * zDim];
814 fftw_plan planForwardFourier = fftw_plan_dft_3d ( xDim, yDim, zDim, translatedMap, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
815 fftw_plan planBackwardFourier = fftw_plan_dft_3d ( xDim, yDim, zDim, fCoeffs, translatedMap, FFTW_BACKWARD, FFTW_ESTIMATE );
818 for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
820 for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
822 for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
824 arrayPos = wIt + zDim * ( vIt + yDim * uIt );
826 if ( map[arrayPos] == map[arrayPos] ) { translatedMap[arrayPos][0] = map[arrayPos]; }
827 else { translatedMap[arrayPos][0] = 0.0; }
828 translatedMap[arrayPos][1] = 0.0;
834 fftw_execute ( planForwardFourier );
837 for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
839 for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
841 for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
844 arrayPos = wIt + zDim * ( vIt + yDim * uIt );
845 real = fCoeffs[arrayPos][0];
846 imag = fCoeffs[arrayPos][1];
849 if ( uIt >
static_cast<proshade_unsign
> ( (xDim+1) / 2) ) { h = uIt -
static_cast <proshade_signed
> ( xDim ); }
else { h = uIt; }
850 if ( vIt >
static_cast<proshade_unsign
> ( (yDim+1) / 2) ) { k = vIt -
static_cast <proshade_signed
> ( yDim ); }
else { k = vIt; }
851 if ( wIt >
static_cast<proshade_unsign
> ( (zDim+1) / 2) ) { l = wIt -
static_cast <proshade_signed
> ( zDim ); }
else { l = wIt; }
854 exponent = ( ( (
static_cast <proshade_double
> ( h ) /
static_cast <proshade_double
> ( xAngs ) ) * (-xMov) ) +
855 ( (
static_cast <proshade_double
> ( k ) /
static_cast <proshade_double
> ( yAngs ) ) * (-yMov) ) +
856 ( (
static_cast <proshade_double
> ( l ) /
static_cast <proshade_double
> ( zAngs ) ) * (-zMov) ) ) * 2.0 * M_PI;
858 trCoeffReal = cos ( exponent );
859 trCoeffImag = sin ( exponent );
863 fCoeffs[arrayPos][0] = hlpArrReal / normFactor;
864 fCoeffs[arrayPos][1] = hlpArrImag / normFactor;
870 fftw_execute ( planBackwardFourier );
873 for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
875 for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
877 for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
879 arrayPos = wIt + zDim * ( vIt + yDim * uIt );
880 map[arrayPos] = translatedMap[arrayPos][0];
886 fftw_destroy_plan ( planForwardFourier );
887 fftw_destroy_plan ( planBackwardFourier );
889 delete[] translatedMap;
912 void ProSHADE_internal_mapManip::blurSharpenMap ( proshade_double*& map, proshade_double*& blurredMap, proshade_unsign xDimS, proshade_unsign yDimS, proshade_unsign zDimS, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single blurringFactor )
915 proshade_signed xDim = xDimS;
916 proshade_signed yDim = yDimS;
917 proshade_signed zDim = zDimS;
918 proshade_double real, imag, S, mag, phase;
919 proshade_signed h, k, l;
920 proshade_unsign arrayPos = 0;
921 proshade_double normFactor =
static_cast<proshade_double
> ( xDim * yDim * zDim );
924 fftw_complex* mapCoeffs =
new fftw_complex[xDim * yDim * zDim];
925 fftw_complex* mapMask =
new fftw_complex[xDim * yDim * zDim];
932 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> (xDim * yDim * zDim); iter++ )
934 mapMask[iter][0] = map[iter];
935 mapMask[iter][1] = 0.0;
939 fftw_plan forward = fftw_plan_dft_3d ( xDim, yDim, zDim, mapMask, mapCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
940 fftw_plan inverse = fftw_plan_dft_3d ( xDim, yDim, zDim, mapCoeffs, mapMask, FFTW_BACKWARD, FFTW_ESTIMATE );
943 fftw_execute ( forward );
946 for ( proshade_unsign uIt = 0; uIt < static_cast<proshade_unsign> ( xDim ); uIt++ )
948 for ( proshade_unsign vIt = 0; vIt < static_cast<proshade_unsign> ( yDim ); vIt++ )
950 for ( proshade_unsign wIt = 0; wIt < static_cast<proshade_unsign> ( zDim ); wIt++ )
953 arrayPos = wIt + zDim * ( vIt + yDim * uIt );
954 real = mapCoeffs[arrayPos][0];
955 imag = mapCoeffs[arrayPos][1];
958 if ( uIt >
static_cast<proshade_unsign
> ( (xDim+1) / 2) ) { h = uIt -
static_cast <proshade_signed
> ( xDim ); }
else { h = uIt; }
959 if ( vIt >
static_cast<proshade_unsign
> ( (yDim+1) / 2) ) { k = vIt -
static_cast <proshade_signed
> ( yDim ); }
else { k = vIt; }
960 if ( wIt >
static_cast<proshade_unsign
> ( (zDim+1) / 2) ) { l = wIt -
static_cast <proshade_signed
> ( zDim ); }
else { l = wIt; }
963 S = ( pow(
static_cast<proshade_double
> ( h ) / xAngs, 2.0 ) +
964 pow(
static_cast<proshade_double
> ( k ) / yAngs, 2.0 ) +
965 pow(
static_cast<proshade_double
> ( l ) / zAngs, 2.0 ) );
966 mag = std::sqrt ( (real*real) + (imag*imag) ) * std::exp ( - ( ( blurringFactor * S ) / 4.0 ) );
967 phase = std::atan2 ( imag, real );
970 mapCoeffs[arrayPos][0] = ( mag * cos(phase) ) / normFactor;
971 mapCoeffs[arrayPos][1] = ( mag * sin(phase) ) / normFactor;
977 fftw_execute ( inverse );
980 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> (xDim * yDim * zDim); iter++ )
982 blurredMap[iter] = mapMask[iter][0];
990 fftw_destroy_plan ( forward );
991 fftw_destroy_plan ( inverse );
1015 std::vector<proshade_double> mapVals ( xDim * yDim * zDim, 0.0 );
1018 for ( proshade_unsign iter = 0; iter < ( xDim * yDim * zDim ); iter++ )
1020 mapVals.at(iter) = blurMap[iter];
1024 proshade_double* medAndIQR =
new proshade_double[2];
1028 proshade_double maskThreshold = medAndIQR[0] + ( medAndIQR[1] *
static_cast<proshade_double
> ( noIQRs ) );
1031 for ( proshade_unsign iter = 0; iter < ( xDim * yDim * zDim ); iter++ )
1033 if ( blurMap[iter] < maskThreshold )
1036 blurMap[iter] = 0.0;
1068 proshade_unsign arrayPos = 0;
1079 for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1081 for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1083 for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1086 arrayPos = zIt + zDim * ( yIt + yDim * xIt );
1089 if ( map[arrayPos] > 0.001 )
1091 if ( xIt < ret[0] ) { ret[0] = xIt; }
1092 if ( xIt > ret[1] ) { ret[1] = xIt; }
1093 if ( yIt < ret[2] ) { ret[2] = yIt; }
1094 if ( yIt > ret[3] ) { ret[3] = yIt; }
1095 if ( zIt < ret[4] ) { ret[4] = zIt; }
1096 if ( zIt > ret[5] ) { ret[5] = zIt; }
1130 bounds[0] = bounds[0] - xExtraInds;
1131 bounds[1] = bounds[1] + xExtraInds;
1132 bounds[2] = bounds[2] - yExtraInds;
1133 bounds[3] = bounds[3] + yExtraInds;
1134 bounds[4] = bounds[4] - zExtraInds;
1135 bounds[5] = bounds[5] + zExtraInds;
1161 if ( resolution <= 0.0 )
1163 throw ProSHADE_exception (
"Requested resolution not set for map re-sampling.",
"EM00015", __FILE__, __LINE__, __func__,
"There is no resolution value set, but map re-sampling to\n : this unset resolution value is required. This error\n : occurs when a task with no resolution requirement is\n : requested on a map data and the map resolution change is\n : set to \'on\'. Either supply a resolution value, or do not\n : re-sample the map." );
1167 proshade_signed xDim =
static_cast<proshade_signed
> ( xDimS );
1168 proshade_signed yDim =
static_cast<proshade_signed
> ( yDimS );
1169 proshade_signed zDim =
static_cast<proshade_signed
> ( zDimS );
1170 proshade_single oldXSample = ( xAngs /
static_cast<proshade_single
> ( xDim ) );
1171 proshade_single oldYSample = ( yAngs /
static_cast<proshade_single
> ( yDim ) );
1172 proshade_single oldZSample = ( zAngs /
static_cast<proshade_single
> ( zDim ) );
1173 proshade_single newXSample = resolution / 2.0;
1174 proshade_single newYSample = resolution / 2.0;
1175 proshade_single newZSample = resolution / 2.0;
1178 proshade_signed newXDim =
static_cast<proshade_signed
> ( std::ceil ( xAngs / newXSample ) );
1179 proshade_signed newYDim =
static_cast<proshade_signed
> ( std::ceil ( yAngs / newYSample ) );
1180 proshade_signed newZDim =
static_cast<proshade_signed
> ( std::ceil ( zAngs / newZSample ) );
1183 proshade_double* newMap =
new proshade_double [newXDim * newYDim * newZDim];
1186 proshade_signed xBottom = 0, xTop, yBottom = 0, yTop, zBottom = 0, zTop, oldMapIndex, newMapIndex;
1187 std::vector<proshade_double> c000 = std::vector<proshade_double> ( 4, 0.0 );
1188 std::vector<proshade_double> c001 = std::vector<proshade_double> ( 4, 0.0 );
1189 std::vector<proshade_double> c010 = std::vector<proshade_double> ( 4, 0.0 );
1190 std::vector<proshade_double> c011 = std::vector<proshade_double> ( 4, 0.0 );
1191 std::vector<proshade_double> c100 = std::vector<proshade_double> ( 4, 0.0 );
1192 std::vector<proshade_double> c101 = std::vector<proshade_double> ( 4, 0.0 );
1193 std::vector<proshade_double> c110 = std::vector<proshade_double> ( 4, 0.0 );
1194 std::vector<proshade_double> c111 = std::vector<proshade_double> ( 4, 0.0 );
1195 std::vector<proshade_double> c00 = std::vector<proshade_double> ( 4, 0.0 );
1196 std::vector<proshade_double> c01 = std::vector<proshade_double> ( 4, 0.0 );
1197 std::vector<proshade_double> c10 = std::vector<proshade_double> ( 4, 0.0 );
1198 std::vector<proshade_double> c11 = std::vector<proshade_double> ( 4, 0.0 );
1199 std::vector<proshade_double> c0 = std::vector<proshade_double> ( 4, 0.0 );
1200 std::vector<proshade_double> c1 = std::vector<proshade_double> ( 4, 0.0 );
1201 proshade_double xRelative, yRelative, zRelative;
1203 for ( proshade_signed xIt = 0; xIt < newXDim; xIt++ )
1205 for ( proshade_signed yIt = 0; yIt < newYDim; yIt++ )
1207 for ( proshade_signed zIt = 0; zIt < newZDim; zIt++ )
1210 newMapIndex = zIt + newZDim * ( yIt + newYDim * xIt );
1213 for ( proshade_unsign ox = 0; ox < (xDimS-1); ox++ ) {
if ( ( ( xIt * newXSample ) >= ( ox * oldXSample ) ) && ( ( xIt * newXSample ) <= ( (ox+1) * oldXSample ) ) ) { xBottom =
static_cast<proshade_signed
> ( ox );
break; } }
1214 for ( proshade_unsign oy = 0; oy < (yDimS-1); oy++ ) {
if ( ( ( yIt * newYSample ) >= ( oy * oldYSample ) ) && ( ( yIt * newYSample ) <= ( (oy+1) * oldYSample ) ) ) { yBottom =
static_cast<proshade_signed
> ( oy );
break; } }
1215 for ( proshade_unsign oz = 0; oz < (zDimS-1); oz++ ) {
if ( ( ( zIt * newZSample ) >= ( oz * oldZSample ) ) && ( ( zIt * newZSample ) <= ( (oz+1) * oldZSample ) ) ) { zBottom =
static_cast<proshade_signed
> ( oz );
break; } }
1221 oldMapIndex = zBottom + zDimS * ( yBottom + yDimS * xBottom );
1222 c000.at(0) =
static_cast<proshade_double
> ( xBottom * oldXSample );
1223 c000.at(1) =
static_cast<proshade_double
> ( yBottom * oldYSample );
1224 c000.at(2) =
static_cast<proshade_double
> ( zBottom * oldZSample );
1225 c000.at(3) =
static_cast<proshade_double
> ( map[oldMapIndex] );
1227 oldMapIndex = zTop + zDimS * ( yBottom + yDimS * xBottom );
1228 c001.at(0) =
static_cast<proshade_double
> ( xBottom * oldXSample );
1229 c001.at(1) =
static_cast<proshade_double
> ( yBottom * oldYSample );
1230 c001.at(2) =
static_cast<proshade_double
> ( zTop * oldZSample );
1231 c001.at(3) =
static_cast<proshade_double
> ( map[oldMapIndex] );
1233 oldMapIndex = zBottom + zDimS * ( yTop + yDimS * xBottom );
1234 c010.at(0) =
static_cast<proshade_double
> ( xBottom * oldXSample );
1235 c010.at(1) =
static_cast<proshade_double
> ( yTop * oldYSample );
1236 c010.at(2) =
static_cast<proshade_double
> ( zBottom * oldZSample );
1237 c010.at(3) =
static_cast<proshade_double
> ( map[oldMapIndex] );
1239 oldMapIndex = zTop + zDimS * ( yTop + yDimS * xBottom );
1240 c011.at(0) =
static_cast<proshade_double
> ( xBottom * oldXSample );
1241 c011.at(1) =
static_cast<proshade_double
> ( yTop * oldYSample );
1242 c011.at(2) =
static_cast<proshade_double
> ( zTop * oldZSample );
1243 c011.at(3) =
static_cast<proshade_double
> ( map[oldMapIndex] );
1245 oldMapIndex = zBottom + zDimS * ( yBottom + yDimS * xTop );
1246 c100.at(0) =
static_cast<proshade_double
> ( xTop * oldXSample );
1247 c100.at(1) =
static_cast<proshade_double
> ( yBottom * oldYSample );
1248 c100.at(2) =
static_cast<proshade_double
> ( zBottom * oldZSample );
1249 c100.at(3) =
static_cast<proshade_double
> ( map[oldMapIndex] );
1251 oldMapIndex = zTop + zDimS * ( yBottom + yDimS * xTop );
1252 c101.at(0) =
static_cast<proshade_double
> ( xTop * oldXSample );
1253 c101.at(1) =
static_cast<proshade_double
> ( yBottom * oldYSample );
1254 c101.at(2) =
static_cast<proshade_double
> ( zTop * oldZSample );
1255 c101.at(3) =
static_cast<proshade_double
> ( map[oldMapIndex] );
1257 oldMapIndex = zBottom + zDimS * ( yTop + yDimS * xTop );
1258 c110.at(0) =
static_cast<proshade_double
> ( xTop * oldXSample );
1259 c110.at(1) =
static_cast<proshade_double
> ( yTop * oldYSample );
1260 c110.at(2) =
static_cast<proshade_double
> ( zBottom * oldZSample );
1261 c110.at(3) =
static_cast<proshade_double
> ( map[oldMapIndex] );
1263 oldMapIndex = zTop + zDimS * ( yTop + yDimS * xTop );
1264 c111.at(0) =
static_cast<proshade_double
> ( xTop * oldXSample );
1265 c111.at(1) =
static_cast<proshade_double
> ( yTop * oldYSample );
1266 c111.at(2) =
static_cast<proshade_double
> ( zTop * oldZSample );
1267 c111.at(3) =
static_cast<proshade_double
> ( map[oldMapIndex] );
1270 xRelative = ( ( xIt * newXSample ) - ( xBottom * oldXSample ) ) / ( ( xTop * oldXSample ) - ( xBottom * oldXSample ) );
1273 c00.at(0) = ( newXSample * xRelative ) + c000.at(0);
1274 c00.at(1) = c000.at(1);
1275 c00.at(2) = c000.at(2);
1276 c00.at(3) = ( c000.at(3) * ( 1.0 - xRelative ) ) + ( c100.at(3) * xRelative );
1279 c01.at(0) = ( newXSample * xRelative ) + c001.at(0);
1280 c01.at(1) = c001.at(1);
1281 c01.at(2) = c001.at(2);
1282 c01.at(3) = ( c001.at(3) * ( 1.0 - xRelative ) ) + ( c101.at(3) * xRelative );
1285 c10.at(0) = ( newXSample * xRelative ) + c010.at(0);
1286 c10.at(1) = c010.at(1);
1287 c10.at(2) = c010.at(2);
1288 c10.at(3) = ( c010.at(3) * ( 1.0 - xRelative ) ) + ( c110.at(3) * xRelative );
1291 c11.at(0) = ( newXSample * xRelative ) + c011.at(0);
1292 c11.at(1) = c011.at(1);
1293 c11.at(2) = c011.at(2);
1294 c11.at(3) = ( c011.at(3) * ( 1.0 - xRelative ) ) + ( c111.at(3) * xRelative );
1297 yRelative = ( ( yIt * newYSample ) - ( yBottom * oldYSample ) ) / ( ( yTop * oldYSample ) - ( yBottom * oldYSample ) );
1300 c0.at(0) = c00.at(0);
1301 c0.at(1) = ( newYSample * yRelative ) + c00.at(1);
1302 c0.at(2) = c00.at(2);
1303 c0.at(3) = ( c00.at(3) * ( 1.0 - yRelative ) ) + ( c10.at(3) * yRelative );
1306 c1.at(0) = c01.at(0);
1307 c1.at(1) = ( newYSample * yRelative ) + c01.at(1);
1308 c1.at(2) = c01.at(2);
1309 c1.at(3) = ( c01.at(3) * ( 1.0 - yRelative ) ) + ( c11.at(3) * yRelative );
1312 zRelative = ( ( zIt * newZSample ) - ( zBottom * oldZSample ) ) / ( ( zTop * oldZSample ) - ( zBottom * oldZSample ) );
1313 newMap[newMapIndex] = ( c0.at(3) * ( 1.0 - zRelative ) ) + ( c1.at(3) * zRelative );
1320 map =
new proshade_double [newXDim * newYDim * newZDim];
1323 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ )
1325 map[iter] = newMap[iter];
1332 corrs[0] = newXDim - xDim;
1333 corrs[1] = newYDim - yDim;
1334 corrs[2] = newZDim - zDim;
1335 corrs[3] = newXDim * newXSample;
1336 corrs[4] = newYDim * newYSample;
1337 corrs[5] = newZDim * newZSample;
1362 if ( resolution <= 0.0 )
1364 throw ProSHADE_exception (
"Requested resolution not set for map re-sampling.",
"EM00015", __FILE__, __LINE__, __func__,
"There is no resolution value set, but map re-sampling to\n : this unset resolution value is required. This error\n : occurs when a task with no resolution requirement is\n : requested on a map data and the map resolution change is\n : set to \'on\'. Either supply a resolution value, or do not\n : re-sample the map." );
1372 if ( newXDim % 2 != 0 ) { newXDim += 1; }
1373 if ( newYDim % 2 != 0 ) { newYDim += 1; }
1374 if ( newZDim % 2 != 0 ) { newZDim += 1; }
1376 proshade_signed preXChange, preYChange, preZChange;
1377 if ( ( xDimS % 2 ) == 0 ) { preXChange = std::ceil ( (
static_cast<proshade_signed
> ( xDimS ) -
static_cast<proshade_signed
> ( newXDim ) ) / 2.0 ); }
1378 else { preXChange = std::floor ( (
static_cast<proshade_signed
> ( xDimS ) -
static_cast<proshade_signed
> ( newXDim ) ) / 2.0 ); }
1379 if ( ( yDimS % 2 ) == 0 ) { preYChange = std::ceil ( (
static_cast<proshade_signed
> ( yDimS ) -
static_cast<proshade_signed
> ( newYDim ) ) / 2.0 ); }
1380 else { preYChange = std::floor ( (
static_cast<proshade_signed
> ( yDimS ) -
static_cast<proshade_signed
> ( newYDim ) ) / 2.0 ); }
1381 if ( ( zDimS % 2 ) == 0 ) { preZChange = std::ceil ( (
static_cast<proshade_signed
> ( zDimS ) -
static_cast<proshade_signed
> ( newZDim ) ) / 2.0 ); }
1382 else { preZChange = std::floor ( (
static_cast<proshade_signed
> ( zDimS ) -
static_cast<proshade_signed
> ( newZDim ) ) / 2.0 ); }
1384 proshade_signed postXChange =
static_cast<proshade_signed
> ( xDimS ) - ( preXChange +
static_cast<proshade_signed
> ( newXDim ) );
1385 proshade_signed postYChange =
static_cast<proshade_signed
> ( yDimS ) - ( preYChange +
static_cast<proshade_signed
> ( newYDim ) );
1386 proshade_signed postZChange =
static_cast<proshade_signed
> ( zDimS ) - ( preZChange +
static_cast<proshade_signed
> ( newZDim ) );
1388 proshade_signed origSizeArr = 0, newSizeArr = 0;
1389 proshade_double normFactor =
static_cast<proshade_double
> ( xDimS * yDimS * zDimS );
1392 fftw_complex *origMap, *fCoeffs, *newFCoeffs, *newMap;
1393 fftw_plan planForwardFourier, planBackwardRescaledFourier;
1395 xDimS, yDimS, zDimS, newXDim, newYDim, newZDim );
1398 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDimS * yDimS * zDimS ); iter++ ) { origMap[iter][0] = map[iter]; origMap[iter][1] = 0.0; }
1399 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ ) { newFCoeffs[iter][0] = 0.0; newFCoeffs[iter][1] = 0.0; }
1402 fftw_execute ( planForwardFourier );
1408 for ( proshade_unsign xIt = 0; xIt < newXDim; xIt++ )
1410 for ( proshade_unsign yIt = 0; yIt < newYDim; yIt++ )
1412 for ( proshade_unsign zIt = 0; zIt < newZDim; zIt++ )
1415 origSizeArr = (zIt + preZChange) + zDimS * ( (yIt + preYChange) + yDimS * (xIt + preXChange) );
1416 newSizeArr = zIt + newZDim * ( yIt + newYDim * xIt );
1419 if ( ( ( -1 <
static_cast<proshade_signed
> ( xIt + preXChange ) ) && ( -1 <
static_cast<proshade_signed
> ( yIt + preYChange ) ) && ( -1 <
static_cast<proshade_signed
> ( zIt + preZChange ) ) ) &&
1420 ( ( xIt < static_cast<proshade_unsign> ( newXDim + postXChange ) ) && ( yIt <
static_cast<proshade_unsign
> ( newYDim + postYChange ) ) && ( zIt <
static_cast<proshade_unsign
> ( newZDim + postZChange ) ) ) )
1423 newFCoeffs[newSizeArr][0] = fCoeffs[origSizeArr][0] / normFactor;
1424 newFCoeffs[newSizeArr][1] = fCoeffs[origSizeArr][1] / normFactor;
1434 fftw_execute ( planBackwardRescaledFourier );
1438 map =
new proshade_double [newXDim * newYDim * newZDim];
1440 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( newXDim * newYDim * newZDim ); iter++ ) { map[iter] = newMap[iter][0]; }
1446 corrs[0] =
static_cast<proshade_signed
> ( newXDim ) -
static_cast<proshade_signed
> ( xDimS );
1447 corrs[1] =
static_cast<proshade_signed
> ( newYDim ) -
static_cast<proshade_signed
> ( yDimS );
1448 corrs[2] =
static_cast<proshade_signed
> ( newZDim ) -
static_cast<proshade_signed
> ( zDimS );
1449 corrs[3] =
static_cast<proshade_signed
> ( newXDim ) * ( resolution / 2.0 );
1450 corrs[4] =
static_cast<proshade_signed
> ( newYDim ) * ( resolution / 2.0 );
1451 corrs[5] =
static_cast<proshade_signed
> ( newZDim ) * ( resolution / 2.0 );
1476 void ProSHADE_internal_mapManip::allocateResolutionFourierMemory ( fftw_complex*& origMap, fftw_complex*& fCoeffs, fftw_complex*& newFCoeffs, fftw_complex*& newMap, fftw_plan& planForwardFourier, fftw_plan& planBackwardRescaledFourier, proshade_unsign xDimOld, proshade_unsign yDimOld, proshade_unsign zDimOld, proshade_unsign xDimNew, proshade_unsign yDimNew, proshade_unsign zDimNew )
1479 origMap =
new fftw_complex [xDimOld * yDimOld * zDimOld];
1480 fCoeffs =
new fftw_complex [xDimOld * yDimOld * zDimOld];
1481 newFCoeffs =
new fftw_complex [xDimNew * yDimNew * zDimNew];
1482 newMap =
new fftw_complex [xDimNew * yDimNew * zDimNew];
1491 planForwardFourier = fftw_plan_dft_3d ( xDimOld, yDimOld, zDimOld, origMap, fCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
1492 planBackwardRescaledFourier = fftw_plan_dft_3d ( xDimNew, yDimNew, zDimNew, newFCoeffs, newMap, FFTW_BACKWARD, FFTW_ESTIMATE );
1513 fftw_destroy_plan ( planForwardFourier );
1514 fftw_destroy_plan ( planBackwardRescaledFourier );
1519 delete[] newFCoeffs;
1542 proshade_signed h = 0, k = 0, l = 0, origSizeArr = 0, newSizeArr = 0;
1543 proshade_signed xSeq1FreqStart, ySeq1FreqStart, zSeq1FreqStart, xSeq2FreqStart, ySeq2FreqStart, zSeq2FreqStart;
1546 if ( negativeFirst )
1548 if ( ( xDim % 2 ) == 0 ) { xSeq1FreqStart = xDim / 2; xSeq2FreqStart = xDim / 2; }
else { xSeq1FreqStart = (xDim / 2) + 1; xSeq2FreqStart = xDim / 2; }
1549 if ( ( yDim % 2 ) == 0 ) { ySeq1FreqStart = yDim / 2; ySeq2FreqStart = yDim / 2; }
else { ySeq1FreqStart = (yDim / 2) + 1; ySeq2FreqStart = yDim / 2; }
1550 if ( ( zDim % 2 ) == 0 ) { zSeq1FreqStart = zDim / 2; zSeq2FreqStart = zDim / 2; }
else { zSeq1FreqStart = (zDim / 2) + 1; zSeq2FreqStart = zDim / 2; }
1554 if ( ( xDim % 2 ) == 0 ) { xSeq1FreqStart = xDim / 2; xSeq2FreqStart = xDim / 2; }
else { xSeq1FreqStart = (xDim / 2); xSeq2FreqStart = xDim / 2 + 1; }
1555 if ( ( yDim % 2 ) == 0 ) { ySeq1FreqStart = yDim / 2; ySeq2FreqStart = yDim / 2; }
else { ySeq1FreqStart = (yDim / 2); ySeq2FreqStart = yDim / 2 + 1; }
1556 if ( ( zDim % 2 ) == 0 ) { zSeq1FreqStart = zDim / 2; zSeq2FreqStart = zDim / 2; }
else { zSeq1FreqStart = (zDim / 2); zSeq2FreqStart = zDim / 2 + 1; }
1560 fftw_complex *hlpFCoeffs =
new fftw_complex [xDim * yDim * zDim];
1564 for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1567 if ( xIt < xSeq1FreqStart ) { h = xIt + xSeq2FreqStart; }
else { h = xIt - xSeq1FreqStart; }
1568 for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1571 if ( yIt < ySeq1FreqStart ) { k = yIt + ySeq2FreqStart; }
else { k = yIt - ySeq1FreqStart; }
1573 for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1576 if ( zIt < zSeq1FreqStart ) { l = zIt + zSeq2FreqStart; }
else { l = zIt - zSeq1FreqStart; }
1579 newSizeArr = l + zDim * ( k + yDim * h );
1580 origSizeArr = zIt + zDim * ( yIt + yDim * xIt );
1583 hlpFCoeffs[newSizeArr][0] = fCoeffs[origSizeArr][0];
1584 hlpFCoeffs[newSizeArr][1] = fCoeffs[origSizeArr][1];
1590 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { fCoeffs[iter][0] = hlpFCoeffs[iter][0]; fCoeffs[iter][1] = hlpFCoeffs[iter][1]; }
1593 delete[] hlpFCoeffs;
1614 proshade_double real, imag, mag, phase;
1615 proshade_unsign arrayPos = 0;
1616 proshade_double normFactor =
static_cast<proshade_double
> ( xDim * yDim * zDim );
1619 for ( proshade_unsign uIt = 0; uIt < xDim; uIt++ )
1621 for ( proshade_unsign vIt = 0; vIt < yDim; vIt++ )
1623 for ( proshade_unsign wIt = 0; wIt < zDim; wIt++ )
1626 arrayPos = wIt + zDim * ( vIt + yDim * uIt );
1627 real = mapCoeffs[arrayPos][0];
1628 imag = mapCoeffs[arrayPos][1];
1631 mag = std::sqrt ( (real*real) + (imag*imag) );;
1635 mapCoeffs[arrayPos][0] = ( mag * cos(phase) ) / normFactor;
1636 mapCoeffs[arrayPos][1] = ( mag * sin(phase) ) / normFactor;
1663 proshade_signed xDim = xDimS;
1664 proshade_signed yDim = yDimS;
1665 proshade_signed zDim = zDimS;
1666 proshade_signed currentPos, neighArrPos, neighXPos, neighYPos, neighZPos;
1667 proshade_double neighSum;
1668 proshade_double neighCount = pow ( ( ( fakeMapKernel * 2 ) + 1 ), 3.0 ) - 1.0;
1671 for ( proshade_signed uIt = 0; uIt < xDim; uIt++ )
1673 for ( proshade_signed vIt = 0; vIt < yDim; vIt++ )
1675 for ( proshade_signed wIt = 0; wIt < zDim; wIt++ )
1678 currentPos = wIt + zDim * ( vIt + yDim * uIt );
1682 for ( proshade_signed xCh = -fakeMapKernel; xCh <= +fakeMapKernel; xCh++ )
1684 for ( proshade_signed yCh = -fakeMapKernel; yCh <= +fakeMapKernel; yCh++ )
1686 for ( proshade_signed zCh = -fakeMapKernel; zCh <= +fakeMapKernel; zCh++ )
1688 if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) {
continue; }
1691 neighXPos = uIt + xCh;
if ( neighXPos >= xDim ) { neighXPos -= xDim; };
if ( neighXPos < 0 ) { neighXPos += xDim; }
1692 neighYPos = vIt + yCh;
if ( neighYPos >= yDim ) { neighYPos -= yDim; };
if ( neighYPos < 0 ) { neighYPos += yDim; }
1693 neighZPos = wIt + zCh;
if ( neighZPos >= zDim ) { neighZPos -= zDim; };
if ( neighZPos < 0 ) { neighZPos += zDim; }
1694 neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1697 neighSum += map[neighArrPos];
1703 fakeHalfMap[currentPos] = neighSum / neighCount;
1728 proshade_signed xDim = xDimS, yDim = yDimS, zDim = zDimS, currentPos, neighArrPos, neighXPos, neighYPos, neighZPos, corrIter;
1729 proshade_unsign noCorrVals =
static_cast<proshade_unsign
> ( pow ( ( ( corrMaskKernel * 2 ) + 1 ), 3 ) );
1732 proshade_double *origMap =
new proshade_double [noCorrVals];
1733 proshade_double *fakeHM =
new proshade_double [noCorrVals];
1740 for ( proshade_signed uIt = 0; uIt < xDim; uIt++ )
1742 for ( proshade_signed vIt = 0; vIt < yDim; vIt++ )
1744 for ( proshade_signed wIt = 0; wIt < zDim; wIt++ )
1747 currentPos = wIt + zDim * ( vIt + yDim * uIt );
1751 for ( proshade_signed xCh = -corrMaskKernel; xCh <= +corrMaskKernel; xCh++ )
1753 for ( proshade_signed yCh = -corrMaskKernel; yCh <= +corrMaskKernel; yCh++ )
1755 for ( proshade_signed zCh = -corrMaskKernel; zCh <= +corrMaskKernel; zCh++ )
1758 neighXPos = uIt + xCh;
if ( neighXPos >= xDim ) { neighXPos -= xDim; };
if ( neighXPos < 0 ) { neighXPos += xDim; }
1759 neighYPos = vIt + yCh;
if ( neighYPos >= yDim ) { neighYPos -= yDim; };
if ( neighYPos < 0 ) { neighYPos += yDim; }
1760 neighZPos = wIt + zCh;
if ( neighZPos >= zDim ) { neighZPos -= zDim; };
if ( neighZPos < 0 ) { neighZPos += zDim; }
1761 neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1764 origMap[corrIter] = map[neighArrPos];
1765 fakeHM[corrIter] = fakeHalfMap[neighArrPos];
1800 std::max ( dist / ( yAngs /
static_cast<proshade_single
> ( yDim ) ),
1801 dist / ( zAngs /
static_cast<proshade_single
> ( zDim ) ) ) ) );
1819 void ProSHADE_internal_mapManip::connectMaskBlobs ( proshade_double*& mask, proshade_signed xDim, proshade_signed yDim, proshade_signed zDim, proshade_single xAngs, proshade_single yAngs, proshade_single zAngs, proshade_single maskThres )
1822 proshade_double* hlpMap =
new proshade_double[xDim * yDim * zDim];
1823 proshade_signed addSurroundingPoints = std::max ( 3L,
static_cast<proshade_signed
> ( std::ceil (
getIndicesFromAngstroms( xDim, yDim, zDim, xAngs, yAngs, zAngs, std::max( xAngs, std::max( yAngs, zAngs ) ) * 0.1 ) ) ) );
1824 proshade_signed currPos, neighXPos, neighYPos, neighZPos, neighArrPos;
1830 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { hlpMap[iter] = mask[iter]; }
1833 for ( proshade_signed it = 0; it < addSurroundingPoints; it++ )
1836 for ( proshade_signed xIt = 0; xIt < xDim; xIt++ )
1838 for ( proshade_signed yIt = 0; yIt < yDim; yIt++ )
1840 for ( proshade_signed zIt = 0; zIt < zDim; zIt++ )
1843 currPos = zIt + zDim * ( yIt + yDim * xIt );
1846 if ( hlpMap[currPos] < maskThres ) {
continue; }
1849 for ( proshade_signed xCh = -1; xCh <= +1; xCh++ )
1851 for ( proshade_signed yCh = -1; yCh <= +1; yCh++ )
1853 for ( proshade_signed zCh = -1; zCh <= +1; zCh++ )
1855 if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) {
continue; }
1858 neighXPos = xIt + xCh;
if ( neighXPos < 0 ) {
continue; }
if ( neighXPos >= xDim ) {
continue; }
1859 neighYPos = yIt + yCh;
if ( neighYPos < 0 ) {
continue; }
if ( neighYPos >= yDim ) {
continue; }
1860 neighZPos = zIt + zCh;
if ( neighZPos < 0 ) {
continue; }
if ( neighZPos >= zDim ) {
continue; }
1861 neighArrPos = neighZPos + zDim * ( neighYPos + yDim * neighXPos );
1864 if ( hlpMap[neighArrPos] < maskThres ) { mask[neighArrPos] = maskThres; }
1873 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( xDim * yDim * zDim ); iter++ ) { hlpMap[iter] = mask[iter]; }
1900 while ( bounds[1] >=
static_cast<proshade_signed
> ( xDim ) ) { xDim += 10; }
1901 while ( bounds[3] >=
static_cast<proshade_signed
> ( yDim ) ) { yDim += 10; }
1902 while ( bounds[5] >=
static_cast<proshade_signed
> ( zDim ) ) { zDim += 10; }
1910 proshade_signed XtoY = std::abs ( addToX - addToY );
1911 proshade_signed XtoZ = std::abs ( addToX - addToZ );
1912 proshade_signed YtoZ = std::abs ( addToY - addToZ );
1914 if ( ( ( XtoY < boundsDiffThres ) && ( XtoZ < boundsDiffThres ) ) ||
1915 ( ( XtoY < boundsDiffThres ) && ( YtoZ < boundsDiffThres ) ) ||
1916 ( ( XtoZ < boundsDiffThres ) && ( YtoZ < boundsDiffThres ) ) )
1919 proshade_signed maxSize = std::max ( addToX, std::max ( addToY, addToZ ) );
1927 if ( XtoY <= boundsDiffThres )
1929 proshade_signed maxSize = std::max ( addToX, addToY );
1933 if ( XtoZ <= boundsDiffThres )
1935 proshade_signed maxSize = std::max ( addToX, addToZ );
1939 if ( YtoZ <= boundsDiffThres )
1941 proshade_signed maxSize = std::max ( addToY, addToZ );
1970 proshade_signed ret = fromRange;
1971 std::vector < proshade_signed > posibles, hlp;
1972 proshade_signed sum;
1975 for ( proshade_signed iter = fromRange; iter < toRange; iter++ )
1979 for ( proshade_unsign i = 0; i < static_cast<proshade_unsign> ( hlp.size() ); i++ ) { sum += hlp.at(i); }
1985 for ( proshade_signed iter = fromRange; iter < toRange; iter++ )
1988 if ( iter %2 != 0 ) {
continue; }
1991 if ( posibles.at(iter-fromRange) < ( posibles.at(ret-fromRange) - ( iter - ret ) ) ) { ret = iter; }
1995 if ( ( ret % 2 != 0 ) && ( ret < ( toRange - 1 ) ) ) { ret += 1; }
2016 if ( newBoundRange > oldBoundRange )
2019 proshade_signed distributeThis = newBoundRange - oldBoundRange;
2021 while ( distributeThis != 0 )
2024 distributeThis -= 1;
2026 if ( distributeThis != 0 )
2029 distributeThis -= 1;
2063 void ProSHADE_internal_mapManip::copyMapByBounds ( proshade_signed xFrom, proshade_signed xTo, proshade_signed yFrom, proshade_signed yTo, proshade_signed zFrom, proshade_signed zTo, proshade_signed origXFrom, proshade_signed origYFrom, proshade_signed origZFrom, proshade_signed yDimIndices, proshade_signed zDimIndices, proshade_signed origXDimIndices, proshade_signed origYDimIndices, proshade_signed origZDimIndices, proshade_double*& newMap, proshade_double* origMap )
2066 proshade_signed newMapIndex, oldMapIndex, oldX, oldY, oldZ, newX, newY, newZ;
2069 for ( proshade_signed xIt = xFrom; xIt <= xTo; xIt++ )
2072 newX = ( xIt - xFrom );
2073 oldX = ( newX + ( xFrom - origXFrom ) );
2075 for ( proshade_signed yIt = yFrom; yIt <= yTo; yIt++ )
2078 newY = ( yIt - yFrom );
2079 oldY = ( newY + ( yFrom - origYFrom ) );
2081 for ( proshade_signed zIt = zFrom; zIt <= zTo; zIt++ )
2084 newZ = ( zIt - zFrom );
2085 oldZ = ( newZ + ( zFrom - origZFrom ) );
2088 newMapIndex = newZ + zDimIndices * ( newY + yDimIndices * newX );
2089 oldMapIndex = oldZ + origZDimIndices * ( oldY + origYDimIndices * oldX );
2092 if ( ( ( oldX < 0 ) || ( oldX >= origXDimIndices ) ) ||
2093 ( ( oldY < 0 ) || ( oldY >= origYDimIndices ) ) ||
2094 ( ( oldZ < 0 ) || ( oldZ >= origZDimIndices ) ) )
2097 newMap[newMapIndex] = 0.0;
2102 newMap[newMapIndex] = origMap[oldMapIndex];