26 #include <wrap_fftw.h> 27 #include <makeweights.h> 28 #include <s2_primitive.h> 29 #include <s2_cospmls.h> 30 #include <s2_legendreTransforms.h> 31 #include <s2_semi_fly.h> 32 #include <rotate_so3_utils.h> 33 #include <utils_so3.h> 34 #include <soft_fftw.h> 35 #include <rotate_so3_fftw.h> 41 #include <rvapi_interface.h> 61 if ( !this->_trSigmaPreComputed )
63 std::cerr <<
"!!! ProSHADE ERROR !!! Finding the optimal overlay requires the data from pre-computation of the TrSigma descriptor. This is not really logical for the user, but it saves a lot of time in case both overlay and TrSigma are to be computed. Therefore, you need to call the precomputeTrSigmaDescriptor function. Sorry about the inconvenience." << std::endl;
67 std::stringstream hlpSS;
68 hlpSS <<
"<font color=\"red\">" <<
"Finding the optimal overlay requires the data from pre-computation of the integral over the spheres (E matrices), which has not been done. This looks like an internal bug, please report this case." <<
"</font>";
69 rvapi_set_text ( hlpSS.str().c_str(),
83 this->_so3Coeffs =
reinterpret_cast<fftw_complex*
> ( fftw_malloc (
sizeof ( fftw_complex ) * ( ( 4 * pow( static_cast<double> ( this->_bandwidthLimit ),3.0 ) - static_cast<double> ( this->_bandwidthLimit ) ) / 3.0 ) ) );
86 memset ( this->_so3Coeffs,
88 sizeof ( fftw_complex ) * totalCoeffs_so3 ( this->_bandwidthLimit ) );
92 unsigned int indexO = 0;
93 double signValue = 1.0;
94 bool toBeIgnored =
false;
95 unsigned int bandUsageIndex = 0;
98 for (
unsigned int bandIter = 0; bandIter < this->_bandwidthLimit; bandIter++ )
101 if ( !this->_keepOrRemove ) {
if ( ( bandIter % 2 ) != 0 ) {
continue; } }
105 for (
unsigned int igIt = 0; igIt < static_cast<unsigned int> ( this->_lsToIgnore.size() ); igIt++ ) {
if ( bandIter == static_cast<unsigned int> ( this->_lsToIgnore.at(igIt) ) ) { toBeIgnored =
true; } }
106 if ( toBeIgnored ) {
continue; }
109 wigNorm = 2.0 * M_PI * sqrt ( 2.0 / (2.0 * bandIter + 1.0 ) ) ;
112 for (
unsigned int ordIter = 0; ordIter < static_cast<unsigned int> ( ( bandIter * 2 ) + 1 ); ordIter++ )
115 if ( ( ordIter-bandIter + bandIter ) % 2 ) { signValue = -1.0 ; }
116 else { signValue = 1.0 ; }
119 for (
unsigned int ordIter2 = 0; ordIter2 < static_cast<unsigned int> ( ( bandIter * 2 ) + 1 ); ordIter2++ )
122 indexO = so3CoefLoc ( ordIter-bandIter,
125 this->_bandwidthLimit );
128 this->_so3Coeffs[indexO][0] = this->_trSigmaEMatrix[bandUsageIndex][ordIter][ordIter2][0] * wigNorm * signValue;
129 this->_so3Coeffs[indexO][1] = this->_trSigmaEMatrix[bandUsageIndex][ordIter][ordIter2][1] * wigNorm * signValue;
139 this->_so3InvCoeffs =
reinterpret_cast<fftw_complex*
> ( fftw_malloc (
sizeof ( fftw_complex ) * 8 * pow( static_cast<double> ( this->_bandwidthLimit ), 3.0 ) ) );
140 this->_so3Workspace1 =
reinterpret_cast<fftw_complex*
> ( fftw_malloc (
sizeof ( fftw_complex ) * 8 * pow( static_cast<double> ( this->_bandwidthLimit ), 3.0 ) ) );
141 this->_so3Workspace2 =
reinterpret_cast<fftw_complex*
> ( fftw_malloc (
sizeof ( fftw_complex ) * 14 * pow( static_cast<double> ( this->_bandwidthLimit ), 2.0 ) + 48 * this->_bandwidthLimit ) );
142 this->_so3Workspace3 =
new double [
static_cast<unsigned int> ( 24 * this->_bandwidthLimit + 2 * pow( static_cast<double> ( this->_bandwidthLimit ), 2.0 ) )];
145 fftw_plan inverseSO3;
147 int howmany = 4 * this->_bandwidthLimit * this->_bandwidthLimit;
148 int idist = 2 * this->_bandwidthLimit;
149 int odist = 2 * this->_bandwidthLimit;
152 int inembed[2], onembed[2];
153 inembed[0] = 2 * this->_bandwidthLimit;
154 inembed[1] = 4 * this->_bandwidthLimit * this->_bandwidthLimit;
155 onembed[0] = 2 * this->_bandwidthLimit;
156 onembed[1] = 4 * this->_bandwidthLimit * this->_bandwidthLimit;
163 na[1] = 2 * this->_bandwidthLimit;;
165 inverseSO3 = fftw_plan_many_dft ( rank,
168 this->_so3Workspace1,
180 Inverse_SO3_Naive_fftw ( this->_bandwidthLimit,
183 this->_so3Workspace1,
184 this->_so3Workspace2,
185 this->_so3Workspace3,
190 fftw_free ( this->_so3Coeffs );
191 fftw_free ( this->_so3Workspace1 );
192 fftw_free ( this->_so3Workspace2 );
193 delete this->_so3Workspace3;
195 this->_so3Coeffs =
nullptr;
196 this->_so3Workspace1 =
nullptr;
197 this->_so3Workspace2 =
nullptr;
198 this->_so3Workspace3 =
nullptr;
200 fftw_destroy_plan ( inverseSO3 );
203 this->_so3InvMapComputed =
true;
220 if ( !this->_trSigmaPreComputed )
222 std::cerr <<
"!!! ProSHADE ERROR !!! Finding the optimal overlay requires the data from pre-computation of the TrSigma descriptor. This is not really logical for the user, but it saves a lot of time in case both overlay and TrSigma are to be computed. Therefore, you need to call the precomputeTrSigmaDescriptor function. Sorry about the inconvenience." << std::endl;
226 std::stringstream hlpSS;
227 hlpSS <<
"<font color=\"red\">" <<
"Finding the optimal overlay requires the data from pre-computation of the integral over the spheres (E matrices), which has not been done. This looks like an internal bug, please report this case." <<
"</font>";
228 rvapi_set_text ( hlpSS.str().c_str(),
242 this->_so3InverseCoeffs = std::vector<fftw_complex*> ( this->all->size() );
243 std::vector<fftw_complex*> _so3Coeffs = std::vector<fftw_complex*> ( this->all->size() );
246 double wigNorm = 0.0;
247 unsigned int indexO = 0;
248 double signValue = 1.0;
249 bool toBeIgnored =
false;
250 unsigned int bandUsageIndex = 0;
251 unsigned int fileUsageIndex = 0;
252 unsigned int localComparisonBandLim = 0;
255 for (
unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt++ )
258 if ( this->_enLevelsDoNotFollow.at(strIt) == 1 ) {
continue; }
259 if ( this->_trSigmaDoNotFollow.at(strIt) == 1 ) { fileUsageIndex += 1;
continue; }
262 localComparisonBandLim = std::min( this->one->_bandwidthLimit, this->all->at(strIt)->_bandwidthLimit );
265 _so3Coeffs.at(strIt) =
reinterpret_cast<fftw_complex*
> ( fftw_malloc (
sizeof ( fftw_complex ) * ( ( 4 * pow( static_cast<double> ( localComparisonBandLim ), 3.0 ) - static_cast<double> ( localComparisonBandLim ) ) / 3.0 ) ) );
268 memset ( _so3Coeffs.at(strIt),
270 sizeof ( fftw_complex ) * totalCoeffs_so3 ( localComparisonBandLim ) );
273 for (
unsigned int bandIter = 0; bandIter < localComparisonBandLim; bandIter++ )
276 if ( !this->_keepOrRemove ) {
if ( ( bandIter % 2 ) != 0 ) {
continue; } }
280 for (
unsigned int igIt = 0; igIt < static_cast<unsigned int> ( this->_lsToIgnore.size() ); igIt++ ) {
if ( bandIter == static_cast<unsigned int> ( this->_lsToIgnore.at(igIt) ) ) { toBeIgnored =
true; } }
281 if ( toBeIgnored ) {
continue; }
284 wigNorm = 2.0 * M_PI * sqrt ( 2.0 / (2.0 * bandIter + 1.0 ) ) ;
287 for (
unsigned int ordIter = 0; ordIter < static_cast<unsigned int> ( ( bandIter * 2 ) + 1 ); ordIter++ )
290 if ( ( ordIter-bandIter + bandIter ) % 2 ) { signValue = -1.0 ; }
291 else { signValue = 1.0 ; }
294 for (
unsigned int ordIter2 = 0; ordIter2 < static_cast<unsigned int> ( ( bandIter * 2 ) + 1 ); ordIter2++ )
297 indexO = so3CoefLoc ( ordIter-bandIter, ordIter2-bandIter, bandIter, localComparisonBandLim );
300 _so3Coeffs.at(strIt)[indexO][0] = this->_EMatrices.at(fileUsageIndex).at(bandUsageIndex).at(ordIter).at(ordIter2)[0] * wigNorm * signValue;
301 _so3Coeffs.at(strIt)[indexO][1] = this->_EMatrices.at(fileUsageIndex).at(bandUsageIndex).at(ordIter).at(ordIter2)[1] * wigNorm * signValue;
313 for (
unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt++ )
316 if ( this->_enLevelsDoNotFollow.at(strIt) == 1 ) {
continue; }
317 if ( this->_trSigmaDoNotFollow.at(strIt) == 1 ) {
continue; }
320 localComparisonBandLim = std::min( this->one->_bandwidthLimit, this->all->at(strIt)->_bandwidthLimit );
323 fftw_complex *_so3Workspace1 =
reinterpret_cast<fftw_complex*
> ( fftw_malloc (
sizeof ( fftw_complex ) * 8 * pow( static_cast<double> ( localComparisonBandLim ), 3.0 ) ) );
324 fftw_complex *_so3Workspace2 =
reinterpret_cast<fftw_complex*
> ( fftw_malloc (
sizeof ( fftw_complex ) * 14 * pow( static_cast<double> ( localComparisonBandLim ), 2.0 ) + 48 * localComparisonBandLim ) );
325 double *_so3Workspace3 =
new double [
static_cast<unsigned int> ( 24 * localComparisonBandLim + 2 * pow( static_cast<double> ( localComparisonBandLim ), 2.0 ) )];
328 this->_so3InverseCoeffs.at(strIt) =
reinterpret_cast<fftw_complex*
> ( fftw_malloc (
sizeof ( fftw_complex ) * 8 * pow( static_cast<double> ( localComparisonBandLim ), 3.0 ) ) );
331 fftw_plan inverseSO3;
333 int howmany = 4 * localComparisonBandLim * localComparisonBandLim;
334 int idist = 2 * localComparisonBandLim;
335 int odist = 2 * localComparisonBandLim;
338 int inembed[2], onembed[2];
339 inembed[0] = 2 * localComparisonBandLim;
340 inembed[1] = 4 * localComparisonBandLim * localComparisonBandLim;
341 onembed[0] = 2 * localComparisonBandLim;
342 onembed[1] = 4 * localComparisonBandLim * localComparisonBandLim;
349 na[1] = 2 * localComparisonBandLim;
351 inverseSO3 = fftw_plan_many_dft ( rank,
358 this->_so3InverseCoeffs.at(strIt),
366 Inverse_SO3_Naive_fftw ( localComparisonBandLim,
367 _so3Coeffs.at(strIt),
368 this->_so3InverseCoeffs.at(strIt),
376 fftw_destroy_plan ( inverseSO3 );
377 fftw_free ( _so3Workspace1 );
378 fftw_free ( _so3Workspace2 );
379 fftw_free ( _so3Coeffs.at(strIt) );
380 delete[] _so3Workspace3;
384 this->_so3InvMapComputed =
true;
404 this->_eulerAngles[0] = alpha;
405 this->_eulerAngles[1] = beta;
406 this->_eulerAngles[2] = gamma;
409 this->_eulerAnglesFound =
true;
432 if ( this->_eulerAnglesFound )
434 return ( this->_eulerAngles );
438 if ( !this->_so3InvMapComputed )
440 std::cerr <<
"!!! ProSHADE ERROR !!! Error in file compatison !!! Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). Call the getSO3InverseMap function BEFORE calling the getEulerAngles function." << std::endl;
444 std::stringstream hlpSS;
445 hlpSS <<
"<font color=\"red\">" <<
"Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). This looks like an internal bug, please report this case." <<
"</font>";
446 rvapi_set_text ( hlpSS.str().c_str(),
460 this->_eulerAngles = std::array< double, 3 > { { 0.0, 0.0, 0.0 } };
463 unsigned int maxLoc = 0;
467 for (
unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
469 tmpVal = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
470 if ( maxVal < tmpVal )
478 *correlation = maxVal;
482 ii = std::floor ( maxLoc / ( 4.0 * this->_bandwidthLimit * this->_bandwidthLimit ) );
483 tmp = maxLoc - ( ii * ( 4.0 * this->_bandwidthLimit * this->_bandwidthLimit ) );
484 jj = std::floor ( tmp / ( 2.0 * this->_bandwidthLimit ) );
485 kk = maxLoc - ( ii * ( 4.0 * this->_bandwidthLimit * this->_bandwidthLimit ) ) -
486 jj * ( 2.0 * this->_bandwidthLimit );
488 this->_eulerAngles[2] = ( M_PI * jj / (
static_cast<double> ( this->_bandwidthLimit ) ) );
489 this->_eulerAngles[1] = ( M_PI * ( 2.0 * ii + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
490 this->_eulerAngles[0] = ( M_PI * kk / (
static_cast<double> ( this->_bandwidthLimit ) ) );
493 fftw_free ( this->_so3InvCoeffs );
494 this->_so3InvCoeffs =
nullptr;
497 this->_eulerAnglesFound =
true;
500 return ( this->_eulerAngles );
519 if ( this->_eulerAnglesFound )
521 return ( this->_eulerAngles );
525 if ( !this->_so3InvMapComputed )
527 std::cerr <<
"!!! ProSHADE ERROR !!! Error in file compatison !!! Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). Call the getSO3InverseMap function BEFORE calling the getEulerAngles function." << std::endl;
531 std::stringstream hlpSS;
532 hlpSS <<
"<font color=\"red\">" <<
"Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). This looks like an internal bug, please report this case." <<
"</font>";
533 rvapi_set_text ( hlpSS.str().c_str(),
547 this->_eulerAngles = std::vector< std::array<double,3> > ( this->all->size() );
550 unsigned int maxLoc = 0;
553 unsigned int localComparisonBandLim = 0;
555 for (
unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt++ )
558 if ( this->_enLevelsDoNotFollow.at(strIt) == 1 ) {
continue; }
559 if ( this->_trSigmaDoNotFollow.at(strIt) == 1 ) {
continue; }
567 localComparisonBandLim = std::min ( this->one->_bandwidthLimit, this->all->at(strIt)->_bandwidthLimit );
569 for (
unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( localComparisonBandLim ), 3 ) ); iter++ )
571 tmpVal = pow( this->_so3InverseCoeffs.at(strIt)[iter][0], 2.0 ) + pow( this->_so3InverseCoeffs.at(strIt)[iter][1], 2.0 );
572 if ( maxVal < tmpVal )
581 ii = std::floor ( maxLoc / ( 4.0 * localComparisonBandLim * localComparisonBandLim ) );
582 tmp = maxLoc - ( ii * ( 4.0 * localComparisonBandLim * localComparisonBandLim ) );
583 jj = std::floor ( tmp / ( 2.0 * localComparisonBandLim ) );
584 kk = maxLoc - ( ii * ( 4.0 * localComparisonBandLim * localComparisonBandLim ) ) -
585 jj * ( 2.0 * localComparisonBandLim );
587 this->_eulerAngles.at(strIt)[2] = ( M_PI * jj / (
static_cast<double> ( localComparisonBandLim ) ) );
588 this->_eulerAngles.at(strIt)[1] = ( M_PI * ( 2.0 * ii + 1.0 ) / ( 4.0 * localComparisonBandLim ) ) ;
589 this->_eulerAngles.at(strIt)[0] = ( M_PI * kk / (
static_cast<double> ( localComparisonBandLim ) ) );
592 fftw_free ( this->_so3InverseCoeffs.at(strIt) );
593 this->_so3InverseCoeffs.at(strIt) =
nullptr;
597 this->_eulerAnglesFound =
true;
600 return ( this->_eulerAngles );
629 if ( !this->_so3InvMapComputed )
631 std::cerr <<
"!!! ProSHADE ERROR !!! Error in file compatison !!! Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). Call the getSO3InverseMap function BEFORE calling the getSO3Peaks function." << std::endl;
635 std::stringstream hlpSS;
636 hlpSS <<
"<font color=\"red\">" <<
"Attempted to get Euler angles without pre-computing the required data (SO3 inverse transform). This looks like an internal bug, please report this case." <<
"</font>";
637 rvapi_set_text ( hlpSS.str().c_str(),
650 if ( this->_so3InvCoeffs ==
nullptr )
652 std::cerr <<
"!!! ProSHADE ERROR !!! Error in function getSO3Peaks. The pointer to the SO3 data is empty. Common cause is when you call this function more than once, or when you call the getEulerAngles function before this one. If either is the case, please add \"false\" as an additional parameter to all previous calls of both functions." << std::endl;
656 std::stringstream hlpSS;
657 hlpSS <<
"<font color=\"red\">" <<
"The pointer to the SO3 data is empty. This looks like an internal bug, please report this case." <<
"</font>";
658 rvapi_set_text ( hlpSS.str().c_str(),
676 int xDim = 4 * pow ( this->_bandwidthLimit, 2 );
677 int yDim = 2 * this->_bandwidthLimit;
679 std::array<double,8> hlpArr;
680 std::vector< std::array<double,8> > peakList;
681 std::array<double,4> hlpInterp;
682 std::vector< std::array<double,4> > peakInterpVals;
683 std::vector< std::vector< std::array<double,4> > > peakInterpValsList;
684 std::vector< std::vector<double> > tMat;
685 std::vector< std::vector<double> > avgMat;
686 std::vector< std::vector< std::vector<double> > > uAndVMat;
687 std::vector< std::array<double,4> > rotAxes;
688 std::array<double,4> curAxis;
689 std::vector<double> hlpVec;
690 std::vector< std::array<double,8> > ret;
691 std::vector<double> peakHeights;
694 double angAlpha = 0.0;
695 double angBeta = 0.0;
696 double angGamma = 0.0;
697 double singularityErrTolerance = 0.05;
698 double avgMatWeight = 0.0;
701 double interQuartileRange = 0.0;
707 unsigned int medianPos = 0;
709 bool breakPeak =
false;
712 for (
unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
715 tmpVal = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
718 x = std::floor ( iter / xDim );
719 y = std::floor ( ( iter - x * xDim ) / yDim );
720 z = iter - x * xDim - y * yDim;
724 for (
int xCh = -peakSize; xCh <= +peakSize; xCh++ )
726 if ( breakPeak ) {
break; }
727 for (
int yCh = -peakSize; yCh <= +peakSize; yCh++ )
729 if ( breakPeak ) {
break; }
730 for (
int zCh = -peakSize; zCh <= +peakSize; zCh++ )
732 if ( breakPeak ) {
break; }
733 if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) {
continue; }
735 peakX = x + xCh;
if ( peakX >= yDim ) { peakX -= yDim; };
if ( peakX < 0 ) { peakX += yDim; }
736 peakY = y + yCh;
if ( peakY >= yDim ) { peakY -= yDim; };
if ( peakY < 0 ) { peakY += yDim; }
737 peakZ = z + zCh;
if ( peakZ >= yDim ) { peakZ -= yDim; };
if ( peakZ < 0 ) { peakZ += yDim; }
738 newIter = peakX * xDim + peakY * yDim + peakZ;
740 if ( (pow( this->_so3InvCoeffs[newIter][0], 2.0 ) + pow( this->_so3InvCoeffs[newIter][1], 2.0 )) > tmpVal ) { breakPeak =
true;
break; }
744 if ( breakPeak ) {
continue; }
747 angGamma = ( M_PI * y / (
static_cast<double> ( this->_bandwidthLimit ) ) );
748 angBeta = ( M_PI * ( 2.0 * x + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
749 angAlpha = ( M_PI * z / (
static_cast<double> ( this->_bandwidthLimit ) ) );
752 if ( angAlpha > M_PI ) { angAlpha = angAlpha - M_PI * 2.0; }
753 if ( angBeta > M_PI ) { angBeta = angBeta - M_PI * 2.0; }
754 if ( angGamma > M_PI ) { angGamma = angGamma - M_PI * 2.0; }
757 ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &median, &median, &median, &rotVal,
true );
758 if ( ( rotVal < ( 0.0 + singularityErrTolerance ) ) && ( rotVal > ( 0.0 - singularityErrTolerance ) ) ) {
continue; }
761 peakHeights.emplace_back ( tmpVal );
765 medianPos =
static_cast<unsigned int> ( peakHeights.size() / 2 );
766 std::sort ( peakHeights.begin(), peakHeights.end() );
767 median = peakHeights.at(medianPos);
769 interQuartileRange = peakHeights.at(medianPos+(medianPos/2)) - peakHeights.at(medianPos-(medianPos/2));
770 this->_peakHeightThr = peakHeights.at(medianPos+(medianPos/2)) + ( interQuartileRange * noIQRs );
773 for (
unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
776 tmpVal = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
779 if ( tmpVal < this->_peakHeightThr ) {
continue; }
782 x = std::floor ( iter / xDim );
783 y = std::floor ( ( iter - x * xDim ) / yDim );
784 z = iter - x * xDim - y * yDim;
787 peakInterpVals.clear ( );
791 for (
int xCh = -peakSize; xCh <= +peakSize; xCh++ )
793 if ( breakPeak ) {
break; }
794 for (
int yCh = -peakSize; yCh <= +peakSize; yCh++ )
796 if ( breakPeak ) {
break; }
797 for (
int zCh = -peakSize; zCh <= +peakSize; zCh++ )
799 if ( breakPeak ) {
break; }
800 if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) {
continue; }
802 hlpArr[0] = x + xCh;
if ( hlpArr[0] >= yDim ) { hlpArr[0] -= yDim; };
if ( hlpArr[0] < 0 ) { hlpArr[0] += yDim; }
803 hlpArr[1] = y + yCh;
if ( hlpArr[1] >= yDim ) { hlpArr[1] -= yDim; };
if ( hlpArr[1] < 0 ) { hlpArr[1] += yDim; }
804 hlpArr[2] = z + zCh;
if ( hlpArr[2] >= yDim ) { hlpArr[2] -= yDim; };
if ( hlpArr[2] < 0 ) { hlpArr[2] += yDim; }
805 newIter = hlpArr[0] * xDim + hlpArr[1] * yDim + hlpArr[2];
806 hlpArr[3] = pow( this->_so3InvCoeffs[newIter][0], 2.0 ) + pow( this->_so3InvCoeffs[newIter][1], 2.0 );
807 if ( hlpArr[3] > tmpVal ) { breakPeak =
true;
break; }
808 hlpInterp[0] = hlpArr[0]; hlpInterp[1] = hlpArr[1]; hlpInterp[2] = hlpArr[2]; hlpInterp[3] = hlpArr[3];
809 peakInterpVals.emplace_back ( hlpInterp );
813 if ( breakPeak ) {
continue; }
816 angGamma = ( M_PI * y / (
static_cast<double> ( this->_bandwidthLimit ) ) );
817 angBeta = ( M_PI * ( 2.0 * x + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
818 angAlpha = ( M_PI * z / (
static_cast<double> ( this->_bandwidthLimit ) ) );
820 ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &rotVal, &rotVal, &rotVal, &rotVal,
true );
821 if ( ( rotVal < ( 0.0 + singularityErrTolerance ) ) && ( rotVal > ( 0.0 - singularityErrTolerance ) ) ) {
continue; }
826 hlpArr[0] = x + 0;
if ( hlpArr[0] >= yDim ) { hlpArr[0] -= yDim; };
if ( hlpArr[0] < 0 ) { hlpArr[0] += yDim; }
827 hlpArr[1] = y + 0;
if ( hlpArr[1] >= yDim ) { hlpArr[1] -= yDim; };
if ( hlpArr[1] < 0 ) { hlpArr[1] += yDim; }
828 hlpArr[2] = z + 0;
if ( hlpArr[2] >= yDim ) { hlpArr[2] -= yDim; };
if ( hlpArr[2] < 0 ) { hlpArr[2] += yDim; }
829 newIter = hlpArr[0] * xDim + hlpArr[1] * yDim + hlpArr[2];
830 hlpArr[3] = pow( this->_so3InvCoeffs[newIter][0], 2.0 ) + pow( this->_so3InvCoeffs[newIter][1], 2.0 );
831 peakList.emplace_back ( hlpArr );
832 hlpInterp[0] = hlpArr[0];
833 hlpInterp[1] = hlpArr[1];
834 hlpInterp[2] = hlpArr[2];
835 hlpInterp[3] = hlpArr[3];
836 peakInterpVals.emplace_back ( hlpInterp );
837 peakInterpValsList.emplace_back ( peakInterpVals );
842 std::cout <<
">>>>> " << peakList.size() <<
" peaks detected." << std::endl;
848 fftw_free ( this->_so3InvCoeffs );
849 this->_so3InvCoeffs =
nullptr;
853 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( peakList.size() ); iter++ )
858 hlpVec.emplace_back ( 0.0 );
859 hlpVec.emplace_back ( 0.0 );
860 hlpVec.emplace_back ( 0.0 );
861 avgMat.emplace_back ( hlpVec );
862 avgMat.emplace_back ( hlpVec );
863 avgMat.emplace_back ( hlpVec );
868 for (
unsigned int valIter = 0; valIter < static_cast<unsigned int> ( peakInterpValsList.at(iter).size() ); valIter++ )
871 angGamma = ( M_PI * peakInterpValsList.at(iter).at(valIter)[1] / (
static_cast<double> ( this->_bandwidthLimit ) ) );
872 angBeta = ( M_PI * ( 2.0 * peakInterpValsList.at(iter).at(valIter)[0] + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
873 angAlpha = ( M_PI * peakInterpValsList.at(iter).at(valIter)[2] / (
static_cast<double> ( this->_bandwidthLimit ) ) );
876 if ( angAlpha > M_PI ) { angAlpha = angAlpha - M_PI * 2.0; }
877 if ( angBeta > M_PI ) { angBeta = angBeta - M_PI * 2.0; }
878 if ( angGamma > M_PI ) { angGamma = angGamma - M_PI * 2.0; }
885 hlpVec.emplace_back ( cos ( angAlpha ) * cos ( angBeta ) * cos ( angGamma ) - sin ( angAlpha ) * sin ( angGamma ) );
886 hlpVec.emplace_back ( sin ( angAlpha ) * cos ( angBeta ) * cos ( angGamma ) + cos ( angAlpha ) * sin ( angGamma ) );
887 hlpVec.emplace_back ( -sin ( angBeta ) * cos ( angGamma ) );
888 tMat.emplace_back ( hlpVec );
891 hlpVec.emplace_back ( -cos ( angAlpha ) * cos ( angBeta ) * sin ( angGamma ) - sin ( angAlpha ) * cos ( angGamma ) );
892 hlpVec.emplace_back ( -sin ( angAlpha ) * cos ( angBeta ) * sin ( angGamma ) + cos ( angAlpha ) * cos ( angGamma ) );
893 hlpVec.emplace_back ( sin ( angBeta ) * sin ( angGamma ) );
894 tMat.emplace_back ( hlpVec );
897 hlpVec.emplace_back ( cos ( angAlpha ) * sin ( angBeta ) );
898 hlpVec.emplace_back ( sin ( angAlpha ) * sin ( angBeta ) );
899 hlpVec.emplace_back ( cos ( angBeta ) );
900 tMat.emplace_back ( hlpVec );
903 avgMat.at(0).at(0) += tMat.at(0).at(0) * peakInterpValsList.at(iter).at(valIter)[3];
904 avgMat.at(0).at(1) += tMat.at(0).at(1) * peakInterpValsList.at(iter).at(valIter)[3];
905 avgMat.at(0).at(2) += tMat.at(0).at(2) * peakInterpValsList.at(iter).at(valIter)[3];
907 avgMat.at(1).at(0) += tMat.at(1).at(0) * peakInterpValsList.at(iter).at(valIter)[3];
908 avgMat.at(1).at(1) += tMat.at(1).at(1) * peakInterpValsList.at(iter).at(valIter)[3];
909 avgMat.at(1).at(2) += tMat.at(1).at(2) * peakInterpValsList.at(iter).at(valIter)[3];
911 avgMat.at(2).at(0) += tMat.at(2).at(0) * peakInterpValsList.at(iter).at(valIter)[3];
912 avgMat.at(2).at(1) += tMat.at(2).at(1) * peakInterpValsList.at(iter).at(valIter)[3];
913 avgMat.at(2).at(2) += tMat.at(2).at(2) * peakInterpValsList.at(iter).at(valIter)[3];
916 avgMatWeight += peakInterpValsList.at(iter).at(valIter)[3];
920 avgMat.at(0).at(0) /= avgMatWeight;
921 avgMat.at(0).at(1) /= avgMatWeight;
922 avgMat.at(0).at(2) /= avgMatWeight;
924 avgMat.at(1).at(0) /= avgMatWeight;
925 avgMat.at(1).at(1) /= avgMatWeight;
926 avgMat.at(1).at(2) /= avgMatWeight;
928 avgMat.at(2).at(0) /= avgMatWeight;
929 avgMat.at(2).at(1) /= avgMatWeight;
930 avgMat.at(2).at(2) /= avgMatWeight;
933 uAndVMat = this->getSingularValuesUandVMatrix ( avgMat, 3, settings );
938 hlpVec.emplace_back ( 0.0 );
939 hlpVec.emplace_back ( 0.0 );
940 hlpVec.emplace_back ( 0.0 );
941 tMat.emplace_back ( hlpVec ); tMat.emplace_back ( hlpVec ); tMat.emplace_back ( hlpVec );
944 for (
unsigned int row = 0; row < 3; row++ )
946 for (
unsigned int col = 0; col < 3; col++ )
948 for (
unsigned int inner = 0; inner < 3; inner++ )
950 tMat.at(row).at(col) += uAndVMat.at(0).at(inner).at(row) * uAndVMat.at(1).at(col).at(inner);
957 if ( ( angBeta >= 0.0 ) ) { angBeta = acos ( tMat.at(2).at(2) ); }
958 else { angBeta = -acos ( tMat.at(2).at(2) ); }
962 if ( angAlpha > ( M_PI / 2.0 ) ) { angAlpha = M_PI - asin ( tMat.at(2).at(1) / sin ( angBeta ) ); }
963 else if ( angAlpha < ( -M_PI / 2.0 ) ) { angAlpha = -M_PI - asin ( tMat.at(2).at(1) / sin ( angBeta ) ); }
964 else { angAlpha = asin ( tMat.at(2).at(1) / sin ( angBeta ) ); }
967 if ( angGamma > ( M_PI / 2.0 ) ) { angGamma = M_PI - asin ( tMat.at(1).at(2) / sin ( angBeta ) ); }
968 else if ( angGamma < ( -M_PI / 2.0 ) ) { angGamma = -M_PI - asin ( tMat.at(1).at(2) / sin ( angBeta ) ); }
969 else { angGamma = asin ( tMat.at(1).at(2) / sin ( angBeta ) ); }
972 peakList.at(iter)[0] = angAlpha;
973 peakList.at(iter)[1] = angBeta;
974 peakList.at(iter)[2] = angGamma;
978 std::cout <<
">>>>> Peaks angles interpolated." << std::endl;
982 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( peakList.size() ); iter++ )
985 angGamma = peakList.at(iter)[2];
986 angBeta = peakList.at(iter)[1];
987 angAlpha = peakList.at(iter)[0];
990 ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &curAxis[0], &curAxis[1], &curAxis[2], &curAxis[3] );
993 rotAxes.emplace_back ( curAxis );
997 for (
unsigned int rot = 0; rot < static_cast<unsigned int> ( peakList.size() ); rot++ )
999 hlpArr[0] = peakList.at(rot)[0];
1000 hlpArr[1] = peakList.at(rot)[1];
1001 hlpArr[2] = peakList.at(rot)[2];
1002 hlpArr[3] = peakList.at(rot)[3];
1003 hlpArr[4] = rotAxes.at(rot)[0];
1004 hlpArr[5] = rotAxes.at(rot)[1];
1005 hlpArr[6] = rotAxes.at(rot)[2];
1006 hlpArr[7] = rotAxes.at(rot)[3];
1007 ret.emplace_back ( hlpArr );
1011 std::cout <<
">>>>> Peaks converted to angle-axis." << std::endl;
1015 std::sort ( ret.begin(), ret.end(), [](
const std::array<double,8>& a,
const std::array<double,8>& b) {
return a[3] > b[3]; });
1037 double ProSHADE_internal::ProSHADE_comparePairwise::maxPeakHeightFromAngleAxis (
double X,
1044 if ( this->_so3InvCoeffs ==
nullptr )
1046 std::cerr <<
"!!! ProSHADE ERROR !!! Error in function checkPeakFromAngleAxis. The pointer to the SO3 data is empty. Common cause is when you call this function more than once, or when you call the getEulerAngles function before this one. If either is the case, please add \"false\" as an additional parameter to all previous calls of both functions." << std::endl;
1050 std::stringstream hlpSS;
1051 hlpSS <<
"<font color=\"red\">" <<
"The pointer to the SO3 data is empty. This looks like an internal bug, please report this case." <<
"</font>";
1052 rvapi_set_text ( hlpSS.str().c_str(),
1066 double peakHeight = 0.0;
1070 double angAlpha = 0.0;
1071 double angBeta = 0.0;
1072 double angGamma = 0.0;
1078 int xDim = 4 * pow ( this->_bandwidthLimit, 2 );
1079 int yDim = 2 * this->_bandwidthLimit;
1084 for (
unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
1087 peakHeight = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
1090 if ( peakHeight <= ret ) {
continue; }
1093 x = std::floor ( iter / xDim );
1094 y = std::floor ( ( iter - x * xDim ) / yDim );
1095 z = iter - x * xDim - y * yDim;
1098 angGamma = ( M_PI * y / (
static_cast<double> ( this->_bandwidthLimit ) ) );
1099 angBeta = ( M_PI * ( 2.0 * x + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
1100 angAlpha = ( M_PI * z / (
static_cast<double> ( this->_bandwidthLimit ) ) );
1103 if ( angAlpha > M_PI ) { angAlpha = angAlpha - M_PI * 2.0; }
1104 if ( angBeta > M_PI ) { angBeta = angBeta - M_PI * 2.0; }
1105 if ( angGamma > M_PI ) { angGamma = angGamma - M_PI * 2.0; }
1108 ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &axX, &axY, &axZ, &angle );
1109 if ( !( ( angle < ( Angle + 0.1 ) ) && ( ( angle > ( Angle - 0.1 ) ) ) ) )
1114 if ( !( ( axX < ( X + 0.1 ) ) && ( ( axX > ( X - 0.1 ) ) ) ) ||
1115 !( ( axY < ( Y + 0.1 ) ) && ( ( axY > ( Y - 0.1 ) ) ) ) ||
1116 !( ( axZ < ( Z + 0.1 ) ) && ( ( axZ > ( Z - 0.1 ) ) ) ) )
1122 if ( ret < peakHeight ) { ret = peakHeight; }
1150 if ( this->_so3InvCoeffs ==
nullptr )
1152 std::cerr <<
"!!! ProSHADE ERROR !!! Error in function maxAvgPeakForSymmetry. The pointer to the SO3 data is empty. Common cause is when you call this function more than once, or when you call the getEulerAngles function before this one. If either is the case, please add \"false\" as an additional parameter to all previous calls of both functions." << std::endl;
1156 std::stringstream hlpSS;
1157 hlpSS <<
"<font color=\"red\">" <<
"The pointer to the SO3 data is empty. This looks like an internal bug, please report this case." <<
"</font>";
1158 rvapi_set_text ( hlpSS.str().c_str(),
1175 double angAlpha = 0.0;
1176 double angBeta = 0.0;
1177 double angGamma = 0.0;
1182 int xDim = 4 * pow ( this->_bandwidthLimit, 2 );
1183 int yDim = 2 * this->_bandwidthLimit;
1185 std::vector< std::array<double,2> > peaks;
1186 std::array<double,2> hlpArr;
1187 bool xPos =
true;
if ( X < 0 ) { xPos =
false; }
1188 bool yPos =
true;
if ( Y < 0 ) { yPos =
false; }
1189 bool zPos =
true;
if ( Z < 0 ) { zPos =
false; }
1192 for (
unsigned int iter = 0; iter < ( 8 * pow( static_cast<unsigned int> ( this->_bandwidthLimit ), 3 ) ); iter++ )
1195 x = std::floor ( iter / xDim );
1196 y = std::floor ( ( iter - x * xDim ) / yDim );
1197 z = iter - x * xDim - y * yDim;
1200 angGamma = ( M_PI * y / (
static_cast<double> ( this->_bandwidthLimit ) ) );
1201 angBeta = ( M_PI * ( 2.0 * x + 1.0 ) / ( 4.0 * this->_bandwidthLimit ) ) ;
1202 angAlpha = ( M_PI * z / (
static_cast<double> ( this->_bandwidthLimit ) ) );
1205 ProSHADE_internal_misc::getAxisAngleFromEuler ( angAlpha, angBeta, angGamma, &axX, &axY, &axZ, &angle );
1206 if ( !( ( std::abs( axX ) < ( std::abs( X ) + 0.1 ) ) && ( ( std::abs( axX ) > ( std::abs( X ) - 0.1 ) ) ) ) ||
1207 !( ( std::abs( axY ) < ( std::abs( Y ) + 0.1 ) ) && ( ( std::abs( axY ) > ( std::abs( Y ) - 0.1 ) ) ) ) ||
1208 !( ( std::abs( axZ ) < ( std::abs( Z ) + 0.1 ) ) && ( ( std::abs( axZ ) > ( std::abs( Z ) - 0.1 ) ) ) ) )
1214 if ( ( ( xPos && ( axX > 0.0 ) ) || ( !xPos && ( axX < 0.0 ) ) ) &&
1215 ( ( yPos && ( axY > 0.0 ) ) || ( !yPos && ( axY < 0.0 ) ) ) &&
1216 ( ( zPos && ( axZ > 0.0 ) ) || ( !zPos && ( axZ < 0.0 ) ) ) )
1223 if ( ( ( xPos && ( axX < 0.0 ) ) || ( !xPos && ( axX > 0.0 ) ) ) &&
1224 ( ( yPos && ( axY < 0.0 ) ) || ( !yPos && ( axY > 0.0 ) ) ) &&
1225 ( ( zPos && ( axZ < 0.0 ) ) || ( !zPos && ( axZ > 0.0 ) ) ) )
1238 hlpArr[0] = angle + M_PI;
1239 hlpArr[1] = pow( this->_so3InvCoeffs[iter][0], 2.0 ) + pow( this->_so3InvCoeffs[iter][1], 2.0 );
1240 peaks.emplace_back ( hlpArr );
1244 std::sort ( peaks.begin(), peaks.end(), [](
const std::array<double,2>& a,
const std::array<double,2>& b) {
return a[0] < b[0]; } );
1247 double angStep = 0.1;
1248 std::vector<double> sums ( std::floor ( (2.0*M_PI/angStep) / Angle ) );
1249 double curSum = 0.0;
1250 double maxVal = 0.0;
1251 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( std::floor ( (2.0*M_PI/angStep) / Angle ) ); iter++ )
1254 for (
unsigned int angCmb = 0; angCmb < static_cast<unsigned int> ( Angle ); angCmb++ )
1257 for (
unsigned int angIt = 0; angIt < static_cast<unsigned int> ( peaks.size() ); angIt++ )
1259 if ( peaks.at(angIt)[0] < ( ( iter*angStep ) + ( ( 2.0 * M_PI / Angle ) * angCmb ) ) ) {
continue; }
1260 if ( peaks.at(angIt)[0] > ( ( (iter+1)*angStep ) + ( ( 2.0 * M_PI / Angle ) * angCmb ) ) ) {
break; }
1262 if ( peaks.at(angIt)[1] > maxVal ) { maxVal = peaks.at(angIt)[1]; }
1267 if ( ret < curSum ) { ret = curSum; }
1294 double axisErrorTolerance,
1296 double percAllowedToMiss,
1300 if ( peaks.size() < 1 )
1302 std::cerr <<
"!!! ProSHADE ERROR !!! No peaks in the data - cannot search for symmetry. Either loosen up the criteria for the getSO3Peaks function, or conclude no symmetry." << std::endl;
1306 std::stringstream hlpSS;
1307 hlpSS <<
"<font color=\"red\">" <<
"Found no peaks in the rotation function map. This could signify no symmetry being present, or very low sampling of the map. Please try increasing resolution and bandwidth if you believe symmetry should be present." <<
"</font>";
1308 rvapi_set_text ( hlpSS.str().c_str(),
1322 std::vector< std::array<double,5> > ret;
1323 std::vector< std::array<double,7> > grpPeaks;
1324 std::vector< std::array<double,7> > additionalPeaks;
1325 std::vector<double> angVals;
1326 std::array<double,7> hlpArr;
1327 std::array<double,5> retArr;
1328 bool sameAxisFound =
false;
1329 bool groupContinuous =
false;
1330 bool oneDifferent =
false;
1331 bool matchedThisPeak =
false;
1332 bool additionalPeakFound =
false;
1333 double errTolerance = 0.0;
1334 double angTolerance = 0.0;
1335 double angDivisionRemainder = 0.0;
1336 double angDivisionBasis = 0.0;
1337 double groupAngle = 0.0;
1341 double totVals = 0.0;
1342 double meanPeak = 0.0;
1343 double nextPeakError = ( M_PI * 2.0 ) / ( static_cast<double> ( this->_bandwidthLimit ) * 2.0 );
1344 double nextSymmetryError = 0.0;
1345 double peakVal = 0.0;
1346 double minAngDist = 0.0;
1347 double currentAverage = 0.0;
1348 int peaksNeededMinimum = 0;
1349 int peaksNeededRemaining = 0;
1350 int peaksNeededFound = 0;
1351 std::vector<double> groupAngleHits;
1352 std::vector< std::array<double,2> > groupMembers;
1353 std::array<double,2> hlpArr2;
1354 std::vector<double> missingPeaks;
1355 std::vector<double> angsToTry;
1356 std::vector<unsigned int> thisAxisSyms;
1357 std::vector<unsigned int> delVecElems;
1359 std::vector< std::vector<unsigned int> > uniqueAxisGroups;
1361 if ( axisErrorTolerance == 0.0 ) { errTolerance = 0.1; }
1362 else { errTolerance = axisErrorTolerance; }
1365 for (
unsigned int peakIter = 0; peakIter < static_cast<unsigned int> ( peaks.size() ); peakIter++ )
1368 if ( ( peaks.at(peakIter)[4] == 0.0 ) && ( peaks.at(peakIter)[5] == 0.0 ) && ( peaks.at(peakIter)[6] == 0.0 ) ) {
continue; }
1371 sameAxisFound =
false;
1374 if ( ( std::abs(peaks.at(peakIter)[4]) > std::abs(peaks.at(peakIter)[5]) ) && ( std::abs(peaks.at(peakIter)[4]) > std::abs(peaks.at(peakIter)[6]) ) )
1376 if ( peaks.at(peakIter)[4] < 0 )
1378 peaks.at(peakIter)[4] *= -1.0;
1379 peaks.at(peakIter)[5] *= -1.0;
1380 peaks.at(peakIter)[6] *= -1.0;
1381 peaks.at(peakIter)[7] *= -1.0;
1384 else if ( ( std::abs(peaks.at(peakIter)[5]) > std::abs(peaks.at(peakIter)[4]) ) && ( std::abs(peaks.at(peakIter)[5]) > std::abs(peaks.at(peakIter)[6]) ) )
1386 if ( peaks.at(peakIter)[5] < 0 )
1388 peaks.at(peakIter)[4] *= -1.0;
1389 peaks.at(peakIter)[5] *= -1.0;
1390 peaks.at(peakIter)[6] *= -1.0;
1391 peaks.at(peakIter)[7] *= -1.0;
1394 else if ( ( std::abs(peaks.at(peakIter)[6]) > std::abs(peaks.at(peakIter)[4]) ) && ( std::abs(peaks.at(peakIter)[6]) > std::abs(peaks.at(peakIter)[5]) ) )
1396 if ( peaks.at(peakIter)[6] < 0 )
1398 peaks.at(peakIter)[4] *= -1.0;
1399 peaks.at(peakIter)[5] *= -1.0;
1400 peaks.at(peakIter)[6] *= -1.0;
1401 peaks.at(peakIter)[7] *= -1.0;
1406 if ( ( peaks.at(peakIter)[7] < ( 0.0 + errTolerance ) ) && ( peaks.at(peakIter)[7] > ( 0.0 - errTolerance ) ) )
1412 for (
unsigned int sameAxisGrp = 0; sameAxisGrp < static_cast<unsigned int> ( uniqueAxisGroups.size() ); sameAxisGrp++ )
1415 for (
unsigned int sameAxis = 0; sameAxis < static_cast<unsigned int> ( uniqueAxisGroups.at(sameAxisGrp).size() ); sameAxis++ )
1418 if ( ( ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[4]-errTolerance) <= peaks.at(peakIter)[4] ) &&
1419 ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[4]+errTolerance) > peaks.at(peakIter)[4] ) ) &&
1420 ( ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[5]-errTolerance) <= peaks.at(peakIter)[5] ) &&
1421 ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[5]+errTolerance) > peaks.at(peakIter)[5] ) ) &&
1422 ( ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[6]-errTolerance) <= peaks.at(peakIter)[6] ) &&
1423 ( (peaks.at(uniqueAxisGroups.at(sameAxisGrp).at(sameAxis))[6]+errTolerance) > peaks.at(peakIter)[6] ) ) )
1425 sameAxisFound =
true;
1426 uniqueAxisGroups.at(sameAxisGrp).emplace_back ( peakIter );
1433 if ( sameAxisFound ) {
continue; }
1436 std::vector<unsigned int> hlpVec;
1437 hlpVec.emplace_back ( peakIter );
1438 uniqueAxisGroups.emplace_back ( hlpVec );
1442 for (
unsigned int peakIter = 0; peakIter < static_cast<unsigned int> ( peaks.size() ); peakIter++ )
1445 if ( ( peaks.at(peakIter)[7] < ( 0.0 + errTolerance ) ) && ( peaks.at(peakIter)[7] > ( 0.0 - errTolerance ) ) )
1451 if ( ( peaks.at(peakIter)[4] == 0.0 ) && ( peaks.at(peakIter)[5] == 0.0 ) && ( peaks.at(peakIter)[6] == 0.0 ) )
1453 for (
unsigned int sameAxisGrp = 0; sameAxisGrp < static_cast<unsigned int> ( uniqueAxisGroups.size() ); sameAxisGrp++ )
1455 uniqueAxisGroups.at(sameAxisGrp).emplace_back ( peakIter );
1461 for (
unsigned int axisGroup = 0; axisGroup < static_cast<unsigned int> ( uniqueAxisGroups.size() ); axisGroup++ )
1463 for (
unsigned int ex1 = 0; ex1 < static_cast<unsigned int> ( uniqueAxisGroups.at(axisGroup).size() ); ex1++ )
1465 for (
unsigned int ex2 = 1; ex2 < static_cast<unsigned int> ( uniqueAxisGroups.at(axisGroup).size() ); ex2++ )
1468 if ( ex1 >= ex2 ) {
continue; }
1471 if ( ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[7] < ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[7] + errTolerance ) ) &&
1472 ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[7] > ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[7] - errTolerance ) ) )
1475 oneDifferent =
false;
1476 if ( ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[0] < ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[0] + errTolerance ) ) &&
1477 ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[0] > ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[0] - errTolerance ) ) )
1483 oneDifferent =
true;
1486 if ( ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[1] < ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[1] + errTolerance ) ) &&
1487 ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[1] > ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[1] - errTolerance ) ) )
1493 oneDifferent =
true;
1496 if ( ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[2] < ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[2] + errTolerance ) ) &&
1497 ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex1))[2] > ( peaks.at(uniqueAxisGroups.at(axisGroup).at(ex2))[2] - errTolerance ) ) )
1503 oneDifferent =
true;
1506 if ( !oneDifferent )
1516 std::cout <<
">>>>> Peaks clustered by symmetry axis." << std::endl;
1520 for (
unsigned int axisGroup = 0; axisGroup < static_cast<unsigned int> ( uniqueAxisGroups.size() ); axisGroup++ )
1523 thisAxisSyms.clear ( );
1524 additionalPeaks.clear ( );
1528 additionalPeakFound =
false;
1530 for (
unsigned int grpIter = 0; grpIter < static_cast<unsigned int> ( uniqueAxisGroups.at(axisGroup).size() ); grpIter++ )
1532 hlpArr[0] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[4];
1533 hlpArr[1] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[5];
1534 hlpArr[2] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[6];
1535 hlpArr[3] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[7];
1536 hlpArr[4] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpIter))[3];
1537 hlpArr[5] =
static_cast<double> ( grpIter );
1540 grpPeaks.emplace_back ( hlpArr );
1543 if ( additionalPeaks.size () > 0 )
1545 for (
unsigned int adIt = 0; adIt < static_cast<unsigned int> ( additionalPeaks.size() ); adIt++ )
1547 grpPeaks.emplace_back ( additionalPeaks.at(adIt) );
1553 printf (
">>>>> Symmetry axis group:\n" );
1554 printf (
" Group\tPeak\t x\t y\t z\t Angle\t Peak\n" );
1555 printf (
" Index\tIndex\t elem.\t elem.\t elem.\t (rad)\theight\n" );
1556 for (
int i = 0; i < static_cast<int> ( grpPeaks.size() ); i++ )
1558 std::sort ( grpPeaks.begin(), grpPeaks.end(), [](
const std::array<double,7>& a,
const std::array<double,7>& b) {
return a[4] > b[4]; } );
1560 printf (
" %d\t %d\t%+.3f\t%+.3f\t%+.3f\t%+.3f\t%+.3f\n", axisGroup, i, grpPeaks.at(i)[0], grpPeaks.at(i)[1], grpPeaks.at(i)[2], grpPeaks.at(i)[3], grpPeaks.at(i)[4] );
1567 angsToTry.clear ( );
1570 std::sort ( grpPeaks.begin(), grpPeaks.end(), [](
const std::array<double,7>& a,
const std::array<double,7>& b) {
return a[4] > b[4]; });
1573 if ( grpPeaks.at(0)[4] < ( 0.0 + this->_peakHeightThr ) )
1580 for (
unsigned int pAD = 1; pAD < static_cast<unsigned int> ( grpPeaks.size() ); pAD++ )
1582 if ( ( std::abs( std::abs ( grpPeaks.at(0)[3] ) - std::abs ( grpPeaks.at(pAD)[3] ) ) < minAngDist ) && std::abs( std::abs ( grpPeaks.at(0)[3] ) - std::abs ( grpPeaks.at(pAD)[3] ) ) > 0.1 )
1584 minAngDist = std::abs( std::abs ( grpPeaks.at(0)[3] ) - std::abs ( grpPeaks.at(pAD)[3] ) );
1589 angDivisionRemainder = std::modf ( static_cast<double> ( (2.0*M_PI) / std::abs ( minAngDist ) ), &angDivisionBasis );
1590 if ( angDivisionRemainder > 0.8 )
1592 angDivisionRemainder -= 1.0;
1593 angDivisionBasis += 1.0;
1595 angDivisionRemainder /= angDivisionBasis;
1596 nextSymmetryError = ( M_PI * 2.0 / angDivisionBasis ) - ( M_PI * 2.0 / ( angDivisionBasis + 1.0 ) );
1597 angTolerance = ( nextPeakError / nextSymmetryError );
1599 if ( ( angDivisionRemainder < ( 0.0 + angTolerance ) ) && ( angDivisionRemainder > ( 0.0 - angTolerance ) ) )
1601 if ( angTolerance < 0.5 )
1603 angsToTry.emplace_back ( angDivisionBasis );
1607 angsToTry.emplace_back ( angDivisionBasis - 1.0 );
1608 angsToTry.emplace_back ( angDivisionBasis );
1609 angsToTry.emplace_back ( angDivisionBasis + 1.0 );
1614 angDivisionRemainder = std::modf ( static_cast<double> ( (2.0*M_PI) / std::abs ( grpPeaks.at(0)[3] ) ), &angDivisionBasis );
1616 if ( angDivisionRemainder > 0.8 )
1618 angDivisionRemainder -= 1.0;
1619 angDivisionBasis += 1.0;
1621 angDivisionRemainder /= angDivisionBasis;
1624 nextSymmetryError = ( M_PI * 2.0 / angDivisionBasis ) - ( M_PI * 2.0 / ( angDivisionBasis + 1.0 ) );
1625 angTolerance = ( nextPeakError / nextSymmetryError );
1627 if ( ( angDivisionRemainder < ( 0.0 + angTolerance ) ) && ( angDivisionRemainder > ( 0.0 - angTolerance ) ) )
1629 if ( angTolerance < 0.5 )
1631 angsToTry.emplace_back ( angDivisionBasis );
1635 angsToTry.emplace_back ( angDivisionBasis - 1.0 );
1636 angsToTry.emplace_back ( angDivisionBasis );
1637 angsToTry.emplace_back ( angDivisionBasis + 1.0 );
1642 std::sort ( angsToTry.begin(), angsToTry.end(), std::greater<double>() );
1643 angsToTry.erase ( std::unique( angsToTry.begin(), angsToTry.end() ), angsToTry.end() );
1646 delVecElems.clear();
1647 for (
unsigned int tAS = 0; tAS < static_cast<unsigned int> ( thisAxisSyms.size() ); tAS++ )
1649 for (
unsigned int nTS = 0; nTS < static_cast<unsigned int> ( angsToTry.size() ); nTS++ )
1651 if ( thisAxisSyms.at(tAS) ==
static_cast<unsigned int> ( angsToTry.at(nTS) ) )
1653 delVecElems.emplace_back ( nTS );
1658 std::sort ( delVecElems.begin(), delVecElems.end(), std::greater<double>() );
1659 for (
unsigned int remEl = 0; remEl < static_cast<unsigned int> ( delVecElems.size() ); remEl++ )
1661 angsToTry.erase ( angsToTry.begin() + delVecElems.at(remEl) );
1666 printf (
">>>>>>>>> Testing the following symmetries: " );
1667 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( angsToTry.size() ); iter++ )
1669 printf (
"C%.0f\t", angsToTry.at(iter) );
1674 if ( angsToTry.size () == 0 )
1676 grpPeaks.at(0)[4] = 0.0;
1680 for (
unsigned int angToTry = 0; angToTry < static_cast<unsigned int> ( angsToTry.size() ); angToTry++ )
1683 groupAngle = ( 2.0 * M_PI ) / angsToTry.at(angToTry);
1686 for (
int iter = static_cast<int> ( -(angsToTry.at(angToTry)/2.0 + 1.0) ); iter <= static_cast<int> ( angsToTry.at(angToTry)/2.0 + 1.0 ); iter++ )
1688 angVals.emplace_back ( static_cast<double> ( iter * groupAngle ) );
1692 groupAngleHits.clear ( );
1693 groupMembers.clear ( );
1694 missingPeaks.clear ( );
1695 for (
unsigned int grpAngIt = 0; grpAngIt < static_cast<unsigned int> ( angVals.size() ); grpAngIt++ )
1697 matchedThisPeak =
false;
1698 for (
unsigned int peakIt = 0; peakIt < static_cast<unsigned int> ( grpPeaks.size() ); peakIt++ )
1700 if ( angVals.at(grpAngIt) == 0.0 )
1702 groupAngleHits.emplace_back ( angVals.at(grpAngIt) );
1703 hlpArr2[0] =
static_cast<double> ( peakIt );
1704 hlpArr2[1] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpPeaks.at(peakIt)[5]))[3];
1705 groupMembers.emplace_back ( hlpArr2 );
1706 matchedThisPeak =
true;
1710 if ( ( angVals.at(grpAngIt) < ( grpPeaks.at(peakIt)[3] + errTolerance ) ) &&
1711 ( angVals.at(grpAngIt) > ( grpPeaks.at(peakIt)[3] - errTolerance ) ) )
1713 groupAngleHits.emplace_back ( angVals.at(grpAngIt) );
1714 hlpArr2[0] =
static_cast<double> ( peakIt );
1715 hlpArr2[1] = peaks.at(uniqueAxisGroups.at(axisGroup).at(grpPeaks.at(peakIt)[5]))[3];
1716 groupMembers.emplace_back ( hlpArr2 );
1717 matchedThisPeak =
true;
1722 if ( !matchedThisPeak )
1724 missingPeaks.emplace_back ( angVals.at(grpAngIt) );
1729 groupContinuous =
false;
1730 if ( groupAngleHits.size () > 1 )
1732 groupContinuous =
true;
1733 for (
unsigned int iter = 1; iter < static_cast<unsigned int> ( groupAngleHits.size () ); iter++ )
1735 if ( ( ( groupAngleHits.at(iter-1) + groupAngle ) < ( groupAngleHits.at(iter) + errTolerance ) ) &&
1736 ( ( groupAngleHits.at(iter-1) + groupAngle ) > ( groupAngleHits.at(iter) - errTolerance ) ) )
1742 groupContinuous =
false;
1749 if ( !groupContinuous || ( static_cast<int> ( groupAngleHits.size() ) < static_cast<int> ( angsToTry.at(angToTry) ) ) )
1753 printf (
">>>>>>>>>>>>>> Checking for missing peaks...\n" );
1757 if ( ( ( static_cast<double> ( angsToTry.at(angToTry) ) - static_cast<double> ( groupAngleHits.size() ) ) /
static_cast<double> ( angsToTry.at(angToTry) ) ) < percAllowedToMiss )
1760 peaksNeededMinimum =
static_cast<int> ( angsToTry.at(angToTry) ) - static_cast<int> ( groupAngleHits.size() );
1761 peaksNeededFound = 0;
1762 peaksNeededRemaining =
static_cast<int> ( missingPeaks.size() );
1763 currentAverage = 0.0;
1764 for (
unsigned int grpHlpIt = 0; grpHlpIt < static_cast<unsigned int> ( grpPeaks.size() ); grpHlpIt++ )
1766 currentAverage += grpPeaks.at(grpHlpIt)[4];
1768 currentAverage /=
static_cast<double> ( grpPeaks.size() );
1769 currentAverage = std::max ( currentAverage - this->_peakHeightThr, this->_peakHeightThr );
1770 if ( angsToTry.at(0) < 9 )
1772 currentAverage = std::max ( this->_peakHeightThr, currentAverage / 4.0 );
1776 for (
unsigned int misPeak = 0; misPeak < static_cast<unsigned int> ( missingPeaks.size() ); misPeak++ )
1778 if ( peaksNeededMinimum > ( peaksNeededRemaining + peaksNeededFound ) ) {
break; }
1780 if ( missingPeaks.at(misPeak) > M_PI ) { peaksNeededRemaining -= 1;
continue; }
1781 if ( missingPeaks.at(misPeak) < -M_PI ) { peaksNeededRemaining -= 1;
continue; }
1782 if ( missingPeaks.at(misPeak) == 0.0 ) { peaksNeededRemaining -= 1;
continue; }
1787 bool pExists =
false;
1789 for (
unsigned int grpHlpIt = 0; grpHlpIt < static_cast<unsigned int> ( grpPeaks.size() ); grpHlpIt++ )
1791 newX += grpPeaks.at(grpHlpIt)[0];
1792 newY += grpPeaks.at(grpHlpIt)[1];
1793 newZ += grpPeaks.at(grpHlpIt)[2];
1795 newX /=
static_cast<double> ( grpPeaks.size() );
1796 newY /=
static_cast<double> ( grpPeaks.size() );
1797 newZ /=
static_cast<double> ( grpPeaks.size() );
1799 if ( missingPeaks.at(misPeak) > 0.0 )
1801 peakVal = this->maxPeakHeightFromAngleAxis ( newX, newY, newZ, missingPeaks.at(misPeak), settings );
1802 if ( peakVal >= currentAverage )
1804 peaksNeededFound += 1;
1809 hlpArr[3] = missingPeaks.at(misPeak);
1814 additionalPeaks.emplace_back ( hlpArr );
1815 additionalPeakFound =
true;
1821 peakVal = this->maxPeakHeightFromAngleAxis ( -newX, -newY, -newZ, -missingPeaks.at(misPeak), settings );
1822 if ( peakVal >= currentAverage )
1824 peaksNeededFound += 1;
1829 hlpArr[3] = missingPeaks.at(misPeak);
1834 additionalPeaks.emplace_back ( hlpArr );
1835 additionalPeakFound =
true;
1842 peakVal = this->maxPeakHeightFromAngleAxis ( -newX, -newY, -newZ, -missingPeaks.at(misPeak), settings );
1843 if ( peakVal >= currentAverage )
1845 peaksNeededFound += 1;
1850 hlpArr[3] = missingPeaks.at(misPeak);
1855 additionalPeaks.emplace_back ( hlpArr );
1856 additionalPeakFound =
true;
1862 peakVal = this->maxPeakHeightFromAngleAxis ( newX, newY, newZ, missingPeaks.at(misPeak), settings );
1863 if ( peakVal >= currentAverage )
1865 peaksNeededFound += 1;
1870 hlpArr[3] = missingPeaks.at(misPeak);
1875 additionalPeaks.emplace_back ( hlpArr );
1876 additionalPeakFound =
true;
1882 peaksNeededRemaining -= 1;
1888 if ( additionalPeakFound )
1892 std::cout <<
">>>>>>>>>>>>>>>>> Missing peaks discovered. Redoing the analysis." << std::endl;
1898 if ( !groupContinuous )
1900 groupAngleHits.emplace_back ( 0.0 );
1901 std::sort ( groupAngleHits.begin(), groupAngleHits.end() );
1903 if ( groupAngleHits.size () > 1 )
1905 groupContinuous =
true;
1906 for (
unsigned int iter = 1; iter < static_cast<unsigned int> ( groupAngleHits.size () ); iter++ )
1908 if ( ( ( groupAngleHits.at(iter-1) + groupAngle ) < ( groupAngleHits.at(iter) + errTolerance ) ) &&
1909 ( ( groupAngleHits.at(iter-1) + groupAngle ) > ( groupAngleHits.at(iter) - errTolerance ) ) )
1915 groupContinuous =
false;
1923 thisAxisSyms.emplace_back ( angsToTry.at(angToTry) );
1926 if ( groupContinuous )
1929 if ( static_cast<int> ( groupAngleHits.size() ) >= static_cast<int> ( angsToTry.at(angToTry) ) )
1932 std::sort ( groupMembers.begin(), groupMembers.end(), [](
const std::array<double,2>& a,
const std::array<double,2>& b) {
return a[1] > b[1]; } );
1933 groupMembers.erase ( std::unique( groupMembers.begin(), groupMembers.end() ), groupMembers.end() );
1937 printf (
">>>>>>>>>>> Found symmetry C%+.0f on lines ", angsToTry.at(angToTry) );
1938 for (
unsigned int i = 0; i < static_cast<unsigned int> ( groupMembers.size() ); i++ )
1940 printf (
"%+.0f ", groupMembers.at(i)[0] );
1946 retArr[0] = angsToTry.at(angToTry);
1954 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( groupMembers.size() ); iter++ )
1956 if ( ( grpPeaks.at(groupMembers.at(iter)[0])[0] == 0.0 ) && ( grpPeaks.at(groupMembers.at(iter)[0])[1] == 0.0 ) &&
1957 ( grpPeaks.at(groupMembers.at(iter)[0])[2] == 0.0 ) ) {
continue; }
1958 if ( groupMembers.at(iter)[1] == 0.0 ) {
continue; }
1960 meanX += grpPeaks.at(groupMembers.at(iter)[0])[0];
1961 meanY += grpPeaks.at(groupMembers.at(iter)[0])[1];
1962 meanZ += grpPeaks.at(groupMembers.at(iter)[0])[2];
1963 if ( grpPeaks.at(groupMembers.at(iter)[0])[6] == 0.0 )
1965 meanPeak += groupMembers.at(iter)[1];
1970 if ( ( meanX == 0.0 ) && ( meanY == 0.0 ) && ( meanZ == 0.0 ) )
1975 retArr[4] = groupMembers.at(0)[1];
1982 meanPeak /= totVals;
1987 retArr[4] = meanPeak;
1990 ret.emplace_back ( retArr );
1992 grpPeaks.at(0)[4] = 0.0;
1997 grpPeaks.at(0)[4] = 0.0;
2002 if ( !groupContinuous )
2004 grpPeaks.at(0)[4] = 0.0;
2011 std::cout <<
">>>>> C symmetries detection complete." << std::endl;
2015 if ( ret.size() > 1 )
2017 std::vector<unsigned int> removeThese;
2018 for (
unsigned int symIt = 0; symIt < static_cast<unsigned int> ( ret.size() ); symIt++ )
2020 for (
unsigned int existSym = 0; existSym < static_cast<unsigned int> ( ret.size() ); existSym++ )
2022 if ( symIt >= existSym ) {
continue; }
2025 if ( ( ( ret.at(existSym)[1] < ( ret.at(symIt)[1] + errTolerance ) ) && ( ret.at(existSym)[1] > ( ret.at(symIt)[1] - errTolerance ) ) ) &&
2026 ( ( ret.at(existSym)[2] < ( ret.at(symIt)[2] + errTolerance ) ) && ( ret.at(existSym)[2] > ( ret.at(symIt)[2] - errTolerance ) ) ) &&
2027 ( ( ret.at(existSym)[3] < ( ret.at(symIt)[3] + errTolerance ) ) && ( ret.at(existSym)[3] > ( ret.at(symIt)[3] - errTolerance ) ) ) )
2030 if ( ret.at(existSym)[0] == ret.at(symIt)[0] )
2033 if ( ret.at(existSym)[4] > ret.at(symIt)[4] )
2035 removeThese.emplace_back ( symIt );
2039 removeThese.emplace_back ( existSym );
2047 std::sort ( removeThese.begin(), removeThese.end(), std::greater<unsigned int>() );
2048 removeThese.erase ( std::unique( removeThese.begin(), removeThese.end() ), removeThese.end() );
2050 for (
unsigned int remIt = 0; remIt < static_cast<unsigned int> ( removeThese.size() ); remIt++ )
2052 ret.erase ( ret.begin() + removeThese.at(remIt) );
2057 std::cout <<
">>>>> Duplicate C symmetries removal complete." << std::endl;
2063 fftw_free ( this->_so3InvCoeffs );
2064 this->_so3InvCoeffs =
nullptr;
2068 std::sort ( ret.begin(), ret.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[4] > b[4]; });
2071 this->_CSymmsFound =
true;
2087 if ( this->_so3InvCoeffs !=
nullptr )
2089 fftw_free ( this->_so3InvCoeffs );
2090 this->_so3InvCoeffs =
nullptr;
2110 return ( this->_peakHeightThr );
2134 if ( !this->_CSymmsFound )
2136 std::cerr <<
"!!! ProSHADE ERROR !!! Tried to simplify the C symmetry peaks before detecting them. Terminating..." << std::endl;
2140 std::stringstream hlpSS;
2141 hlpSS <<
"<font color=\"red\">" <<
"Tried to simplify the C symmetry peaks before detecting them. This looks like internal bug, please report this case." <<
"</font>";
2142 rvapi_set_text ( hlpSS.str().c_str(),
2154 if ( CnSymm.size() < 1 )
2156 std::cout <<
"!!! ProSHADE WARNING !!! Tried to simplify the C-symmetries list, but the list has less than 1 entry." << std::endl;
2157 if ( pf !=
nullptr )
2166 std::vector< std::array<double,5> > ret;
2167 std::vector< std::array<double,2> > pHeights;
2168 std::array<double,2> hlpArr;
2171 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( CnSymm.size() ); iter++ )
2173 hlpArr[0] =
static_cast<double> ( iter );
2174 hlpArr[1] = CnSymm.at(iter)[4];
2175 pHeights.emplace_back ( hlpArr );
2179 std::sort ( pHeights.begin(), pHeights.end(), [](
const std::array<double,2>& a,
const std::array<double,2>& b) {
return a[1] > b[1]; });
2180 double maxPeakThres = pHeights.at(0)[1] - ( pHeights.at(0)[1] * maxGap );
2182 ret.emplace_back ( CnSymm.at(pHeights.at(0)[0]) );
2183 for (
unsigned int iter = 1; iter < static_cast<unsigned int> ( pHeights.size() ); iter++ )
2185 if ( pHeights.at(iter)[1] >= maxPeakThres )
2187 ret.emplace_back ( CnSymm.at(pHeights.at(iter)[0]) );
2192 std::sort ( ret.begin(), ret.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[0] > b[0]; });
2212 double axisErrorTolerance )
2215 double dotProduct = 0.0;
2216 bool changeFound =
true;
2217 bool matchedAll =
true;
2218 bool identical =
true;
2219 std::vector< std::vector<std::array<double,6> > > ret;
2220 std::vector< std::vector<std::array<double,6> > > retHlp;
2221 std::vector< std::array<double,6> > hlpVec;
2222 std::array<double,6> hlpArr;
2225 for (
unsigned int cSym1 = 0; cSym1 < static_cast<unsigned int> ( CnSymm.size() ); cSym1++ )
2227 for (
unsigned int cSym2 = 1; cSym2 < static_cast<unsigned int> ( CnSymm.size() ); cSym2++ )
2230 if ( cSym1 >= cSym2 ) {
continue; }
2236 dotProduct = ( CnSymm.at(cSym1)[1] * CnSymm.at(cSym2)[1] +
2237 CnSymm.at(cSym1)[2] * CnSymm.at(cSym2)[2] +
2238 CnSymm.at(cSym1)[3] * CnSymm.at(cSym2)[3] );
2241 if ( std::abs ( dotProduct ) < axisErrorTolerance )
2243 hlpArr[0] = CnSymm.at(cSym1)[0];
2244 hlpArr[1] = CnSymm.at(cSym1)[1];
2245 hlpArr[2] = CnSymm.at(cSym1)[2];
2246 hlpArr[3] = CnSymm.at(cSym1)[3];
2247 hlpArr[4] = CnSymm.at(cSym1)[4];
2248 hlpVec.emplace_back ( hlpArr );
2250 hlpArr[0] = CnSymm.at(cSym2)[0];
2251 hlpArr[1] = CnSymm.at(cSym2)[1];
2252 hlpArr[2] = CnSymm.at(cSym2)[2];
2253 hlpArr[3] = CnSymm.at(cSym2)[3];
2254 hlpArr[4] = CnSymm.at(cSym2)[4];
2255 hlpVec.emplace_back ( hlpArr );
2256 retHlp.emplace_back ( hlpVec );
2262 if ( retHlp.size() == 0 ) { this->_DSymmsFound =
true;
return ( ret ); }
2266 while ( changeFound )
2269 changeFound =
false;
2272 for (
unsigned int cSym1 = 0; cSym1 < static_cast<unsigned int> ( CnSymm.size() ); cSym1++ )
2275 for (
unsigned int eSym = 0; eSym < static_cast<unsigned int> ( retHlp.size() ); eSym++ )
2282 for (
unsigned int mem = 0; mem < static_cast<unsigned int> ( retHlp.at(eSym).size() ); mem++ )
2285 if ( ( CnSymm.at(cSym1)[0] == retHlp.at(eSym).at(mem)[0] ) &&
2286 ( CnSymm.at(cSym1)[1] == retHlp.at(eSym).at(mem)[1] ) &&
2287 ( CnSymm.at(cSym1)[2] == retHlp.at(eSym).at(mem)[2] ) &&
2288 ( CnSymm.at(cSym1)[3] == retHlp.at(eSym).at(mem)[3] ) )
2294 dotProduct = ( CnSymm.at(cSym1)[1] * retHlp.at(eSym).at(mem)[1] +
2295 CnSymm.at(cSym1)[2] * retHlp.at(eSym).at(mem)[2] +
2296 CnSymm.at(cSym1)[3] * retHlp.at(eSym).at(mem)[3] );
2298 if ( std::abs ( dotProduct ) < axisErrorTolerance ) { ; }
2299 else { matchedAll =
false; }
2303 if ( identical ) {
continue; }
2308 hlpArr[0] = CnSymm.at(cSym1)[0];
2309 hlpArr[1] = CnSymm.at(cSym1)[1];
2310 hlpArr[2] = CnSymm.at(cSym1)[2];
2311 hlpArr[3] = CnSymm.at(cSym1)[3];
2312 hlpArr[4] = CnSymm.at(cSym1)[4];
2313 retHlp.at(eSym).emplace_back ( hlpArr );
2322 for (
unsigned int eSym1 = 0; eSym1 < static_cast<unsigned int> ( retHlp.size() ); eSym1++ )
2325 if ( retHlp.at(eSym1).at(0)[0] == -1.0 ) {
continue; }
2327 for (
unsigned int eSym2 = 1; eSym2 < static_cast<unsigned int> ( retHlp.size() ); eSym2++ )
2330 if ( retHlp.at(eSym2).at(0)[0] == -1.0 ) {
continue; }
2333 if ( eSym1 >= eSym2 ) {
continue; }
2337 for (
unsigned int mem1 = 0; mem1 < static_cast<unsigned int> ( retHlp.at(eSym1).size() ); mem1++ )
2340 for (
unsigned int mem2 = 0; mem2 < static_cast<unsigned int> ( retHlp.at(eSym2).size() ); mem2++ )
2342 if ( ( retHlp.at(eSym1).at(mem1)[0] == retHlp.at(eSym2).at(mem2)[0] ) &&
2343 ( retHlp.at(eSym1).at(mem1)[1] == retHlp.at(eSym2).at(mem2)[1] ) &&
2344 ( retHlp.at(eSym1).at(mem1)[2] == retHlp.at(eSym2).at(mem2)[2] ) &&
2345 ( retHlp.at(eSym1).at(mem1)[3] == retHlp.at(eSym2).at(mem2)[3] ) )
2351 if ( !matchedAll ) { identical =
false;
break; }
2356 retHlp.at(eSym2).at(0)[0] = -1.0;
2362 double avgPeak = 0.0;
2363 double sumSym = 0.0;
2365 for (
unsigned int eSym1 = 0; eSym1 < static_cast<unsigned int> ( retHlp.size() ); eSym1++ )
2367 if ( retHlp.at(eSym1).at(0)[0] == -1.0 ) {
continue; }
2371 avgPeak = retHlp.at(eSym1).at(0)[4] * retHlp.at(eSym1).at(0)[0];
2372 sumSym = retHlp.at(eSym1).at(0)[0];
2373 for (
unsigned int mem = 1; mem < static_cast<unsigned int> ( retHlp.at(eSym1).size() ); mem++ )
2375 avgPeak += retHlp.at(eSym1).at(mem)[4] * retHlp.at(eSym1).at(mem)[0];
2376 sumSym += retHlp.at(eSym1).at(mem)[0];
2381 for (
unsigned int mem = 0; mem < static_cast<unsigned int> ( retHlp.at(eSym1).size() ); mem++ )
2383 hlpArr[0] = retHlp.at(eSym1).at(mem)[0];
2384 hlpArr[1] = retHlp.at(eSym1).at(mem)[1];
2385 hlpArr[2] = retHlp.at(eSym1).at(mem)[2];
2386 hlpArr[3] = retHlp.at(eSym1).at(mem)[3];
2387 hlpArr[4] = retHlp.at(eSym1).at(mem)[4];
2388 hlpArr[5] = avgPeak;
2389 hlpVec.emplace_back ( hlpArr );
2391 ret.emplace_back ( hlpVec );
2396 std::sort ( ret.begin(), ret.end(), [](
const std::vector< std::array<double,6> >& a,
const std::vector< std::array<double,6> >& b) {
return a.at(0)[5] > b.at(0)[5]; });
2399 this->_DSymmsFound =
true;
2424 if ( !this->_DSymmsFound )
2426 std::cerr <<
"!!! ProSHADE ERROR !!! Tried to simplify the D symmetry peaks before detecting them. Terminating..." << std::endl;
2430 std::stringstream hlpSS;
2431 hlpSS <<
"<font color=\"red\">" <<
"Tried to simplify the D symmetry peaks before detecting them. This looks like internal bug, please report this case." <<
"</font>";
2432 rvapi_set_text ( hlpSS.str().c_str(),
2444 if ( DnSymm.size() < 1 )
2446 std::cout <<
"!!! ProSHADE WARNING !!! Tried to simplify the D-symmetries list, but the list has less than 1 entry." << std::endl;
2447 if ( pf !=
nullptr )
2451 std::vector< std::vector< std::array<double,6> > > ret;
2456 std::vector< std::vector<std::array<double,6> > > ret;
2457 std::vector< std::array<double,6> > hVec;
2458 std::vector< std::vector<std::array<double,4> > > pHeights;
2459 std::vector< std::array<double,4> > hlpVec;
2460 std::array<double,4> hlpArr;
2463 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( DnSymm.size() ); iter++ )
2466 hlpArr[0] =
static_cast<double> ( iter );
2469 for (
unsigned int it = 0; it < static_cast<unsigned int> ( DnSymm.at(iter).size() ); it++ )
2471 hlpArr[1] =
static_cast<double> ( it );
2472 hlpArr[2] = DnSymm.at(iter).at(it)[0];
2473 hlpArr[3] = DnSymm.at(iter).at(it)[4];
2474 hlpVec.emplace_back ( hlpArr );
2477 pHeights.emplace_back ( hlpVec );
2483 double firstIt = 0.0;
2484 double secondIt = 0.0;
2485 double avgPeak = 0.0;
2486 std::vector< std::array<double,2> > combsEx;
2487 std::array<double,2> hA;
2488 bool uniqCombo =
true;
2489 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( pHeights.size() ); iter++ )
2492 std::sort ( pHeights.at(iter).begin(), pHeights.at(iter).end(), [](
const std::array<double,4>& a,
const std::array<double,4>& b) {
return a[2] > b[2]; });
2495 lSym = pHeights.at(iter).at(0)[2];
2497 for (
unsigned int it = 0; it < static_cast<unsigned int> ( pHeights.at(iter).size() ); it++ )
2499 if ( pHeights.at(iter).at(it)[2] == lSym )
2501 if ( mPeak < pHeights.at(iter).at(it)[3] )
2503 mPeak = pHeights.at(iter).at(it)[3];
2504 firstIt =
static_cast<double> ( it );
2509 secondIt = std::max ( secondIt, pHeights.at(iter).at(it)[2] );
2512 if ( secondIt == 0.0 ) { ; }
2513 else { lSym = secondIt; }
2515 for (
unsigned int it = 0; it < static_cast<unsigned int> ( pHeights.at(iter).size() ); it++ )
2517 if ( pHeights.at(iter).at(it)[2] == lSym )
2519 if ( mPeak < pHeights.at(iter).at(it)[3] )
2521 if ( static_cast<double> ( it ) != firstIt )
2523 mPeak = pHeights.at(iter).at(it)[3];
2524 secondIt =
static_cast<double> ( it );
2533 hVec.emplace_back ( DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1]) );
2534 hVec.emplace_back ( DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1]) );
2535 avgPeak = ( DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1])[4] + DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1])[4] ) / 2.0;
2536 hVec.at(0)[5] = avgPeak;
2537 hVec.at(1)[5] = avgPeak;
2538 ret.emplace_back ( hVec );
2539 hA[0] = DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1])[0];
2540 hA[1] = DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1])[0];
2541 combsEx.emplace_back ( hA );
2546 for (
unsigned int it = 0; it < static_cast<unsigned int> ( combsEx.size() ); it++ )
2548 if ( ( combsEx.at(it)[0] == DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1])[0] ) &&
2549 ( combsEx.at(it)[1] == DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1])[0] ) )
2558 hVec.emplace_back ( DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1]) );
2559 hVec.emplace_back ( DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1]) );
2560 ret.emplace_back ( hVec );
2561 hA[0] = DnSymm.at(iter).at(pHeights.at(iter).at(firstIt)[1])[0];
2562 hA[1] = DnSymm.at(iter).at(pHeights.at(iter).at(secondIt)[1])[0];
2563 combsEx.emplace_back ( hA );
2569 if ( ret.size () < 6 )
2571 std::sort ( ret.begin(), ret.end(), [](
const std::vector< std::array<double,6> >& a,
const std::vector< std::array<double,6> >& b) {
return a.at(0)[0] > b.at(0)[0]; });
2598 double passThreshold )
2612 if ( dist >= passThreshold ) { ret =
true; }
2613 else { ret =
false; }
2634 double *tetrSymmPeakAvg,
2635 double axisErrorTolerance )
2638 std::vector< std::vector<std::array<double,5> > > ret;
2639 std::vector< std::array<double,5> > hlpVec;
2640 std::array<double,5> hlpArr;
2641 std::vector<unsigned int> C3s;
2642 double dotProduct = 0.0;
2643 *tetrSymmPeakAvg = 0.0;
2644 double totPairs = 0.0;
2647 for (
unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2649 if ( CnSymm.at(cSym)[0] == 3 )
2651 C3s.emplace_back ( cSym );
2656 for (
unsigned int c3Sym = 0; c3Sym < static_cast<unsigned int> ( C3s.size() ); c3Sym++ )
2658 for (
unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2660 if ( CnSymm.at(cSym)[0] == 3 )
2662 if ( cSym == C3s.at(c3Sym) ) {
continue; }
2664 dotProduct = ( CnSymm.at(C3s.at(c3Sym))[1] * CnSymm.at(cSym)[1] +
2665 CnSymm.at(C3s.at(c3Sym))[2] * CnSymm.at(cSym)[2] +
2666 CnSymm.at(C3s.at(c3Sym))[3] * CnSymm.at(cSym)[3] );
2668 if ( ( ( 1.0 / 3.0 ) > ( dotProduct - axisErrorTolerance ) ) &&
2669 ( ( 1.0 / 3.0 ) < ( dotProduct + axisErrorTolerance ) ) )
2674 hlpArr[0] = CnSymm.at(C3s.at(c3Sym))[0];
2675 hlpArr[1] = CnSymm.at(C3s.at(c3Sym))[1];
2676 hlpArr[2] = CnSymm.at(C3s.at(c3Sym))[2];
2677 hlpArr[3] = CnSymm.at(C3s.at(c3Sym))[3];
2678 hlpArr[4] = ( CnSymm.at(C3s.at(c3Sym))[4] + CnSymm.at(cSym)[4] ) / 2.0;
2679 hlpVec.emplace_back ( hlpArr );
2681 hlpArr[0] = CnSymm.at(cSym)[0];
2682 hlpArr[1] = CnSymm.at(cSym)[1];
2683 hlpArr[2] = CnSymm.at(cSym)[2];
2684 hlpArr[3] = CnSymm.at(cSym)[3];
2685 hlpVec.emplace_back ( hlpArr );
2686 ret.emplace_back ( hlpVec );
2687 *tetrSymmPeakAvg += hlpArr[4];
2695 *tetrSymmPeakAvg /= totPairs;
2716 double *octaSymmPeakAvg,
2717 double axisErrorTolerance )
2720 std::vector< std::vector< std::array<double,5 > > > ret;
2721 std::vector< std::array<double,5> > hlpVec;
2722 std::array<double,5> hlpArr;
2723 std::vector<unsigned int> C4s;
2724 double dotProduct = 0.0;
2725 *octaSymmPeakAvg = 0.0;
2726 double totPairs = 0.0;
2729 for (
unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2731 if ( CnSymm.at(cSym)[0] == 4 )
2733 C4s.emplace_back ( cSym );
2738 for (
unsigned int c4Sym = 0; c4Sym < static_cast<unsigned int> ( C4s.size() ); c4Sym++ )
2740 for (
unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2742 if ( CnSymm.at(cSym)[0] == 3 )
2744 dotProduct = ( CnSymm.at(C4s.at(c4Sym))[1] * CnSymm.at(cSym)[1] +
2745 CnSymm.at(C4s.at(c4Sym))[2] * CnSymm.at(cSym)[2] +
2746 CnSymm.at(C4s.at(c4Sym))[3] * CnSymm.at(cSym)[3] );
2748 if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( dotProduct - axisErrorTolerance ) ) &&
2749 ( ( 1.0 / sqrt ( 3.0 ) ) < ( dotProduct + axisErrorTolerance ) ) )
2754 hlpArr[0] = CnSymm.at(C4s.at(c4Sym))[0];
2755 hlpArr[1] = CnSymm.at(C4s.at(c4Sym))[1];
2756 hlpArr[2] = CnSymm.at(C4s.at(c4Sym))[2];
2757 hlpArr[3] = CnSymm.at(C4s.at(c4Sym))[3];
2758 hlpArr[4] = ( CnSymm.at(C4s.at(c4Sym))[4] + CnSymm.at(cSym)[4] ) / 2.0;
2759 hlpVec.emplace_back ( hlpArr );
2761 hlpArr[0] = CnSymm.at(cSym)[0];
2762 hlpArr[1] = CnSymm.at(cSym)[1];
2763 hlpArr[2] = CnSymm.at(cSym)[2];
2764 hlpArr[3] = CnSymm.at(cSym)[3];
2765 hlpVec.emplace_back ( hlpArr );
2766 ret.emplace_back ( hlpVec );
2767 *octaSymmPeakAvg += hlpArr[4];
2775 *octaSymmPeakAvg /= totPairs;
2796 double *icosSymmPeakAvg,
2797 double axisErrorTolerance )
2800 std::vector< std::vector< std::array<double,5> > > ret;
2801 std::vector< std::array<double,5> > hlpVec;
2802 std::array<double,5> hlpArr;
2803 std::vector<unsigned int> C5s;
2804 double dotProduct = 0.0;
2805 *icosSymmPeakAvg = 0.0;
2806 double totPairs = 0.0;
2809 for (
unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2811 if ( CnSymm.at(cSym)[0] == 5 )
2813 C5s.emplace_back ( cSym );
2818 for (
unsigned int c5Sym = 0; c5Sym < static_cast<unsigned int> ( C5s.size() ); c5Sym++ )
2820 for (
unsigned int cSym = 0; cSym < static_cast<unsigned int> ( CnSymm.size() ); cSym++ )
2822 if ( CnSymm.at(cSym)[0] == 3 )
2824 dotProduct = ( CnSymm.at(C5s.at(c5Sym))[1] * CnSymm.at(cSym)[1] +
2825 CnSymm.at(C5s.at(c5Sym))[2] * CnSymm.at(cSym)[2] +
2826 CnSymm.at(C5s.at(c5Sym))[3] * CnSymm.at(cSym)[3] );
2828 if ( ( ( sqrt ( 5.0 ) / 3.0 ) > ( dotProduct - axisErrorTolerance ) ) &&
2829 ( ( sqrt ( 5.0 ) / 3.0 ) < ( dotProduct + axisErrorTolerance ) ) )
2834 hlpArr[0] = CnSymm.at(C5s.at(c5Sym))[0];
2835 hlpArr[1] = CnSymm.at(C5s.at(c5Sym))[1];
2836 hlpArr[2] = CnSymm.at(C5s.at(c5Sym))[2];
2837 hlpArr[3] = CnSymm.at(C5s.at(c5Sym))[3];
2838 hlpArr[4] = ( CnSymm.at(C5s.at(c5Sym))[4] + CnSymm.at(cSym)[4] ) / 2.0;
2839 hlpVec.emplace_back ( hlpArr );
2841 hlpArr[0] = CnSymm.at(cSym)[0];
2842 hlpArr[1] = CnSymm.at(cSym)[1];
2843 hlpArr[2] = CnSymm.at(cSym)[2];
2844 hlpArr[3] = CnSymm.at(cSym)[3];
2845 hlpVec.emplace_back ( hlpArr );
2846 ret.emplace_back ( hlpVec );
2847 *icosSymmPeakAvg += hlpArr[4];
2855 *icosSymmPeakAvg /= totPairs;
2881 std::vector< std::array<double,5> > tetrSymms,
2882 std::vector< std::array<double,5> > allCs,
2883 double axisErrorTolerance,
2887 std::vector< std::array<double,5> > ret;
2888 std::array<double,5> hlpArr;
2889 bool allAnglesMet =
false;
2890 double dotProduct = 0.0;
2893 ret.emplace_back ( tetrSymms.at(0) );
2894 ret.emplace_back ( tetrSymms.at(1) );
2898 std::cout <<
">> Searching for tetrahederal symmetry axes." << std::endl;
2902 std::cout <<
">>>>> Searching for the four C3 symmetry axes." << std::endl;
2906 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
2909 if ( ret.size() > 2 )
2915 if ( allCs.at(iter)[0] != 3.0 )
2921 if ( ret.size() == 0 )
2923 ret.emplace_back ( allCs.at(iter) );
2928 allAnglesMet =
true;
2929 for (
unsigned int i = 0; i < static_cast<unsigned int> ( ret.size() ); i++ )
2931 dotProduct = ( ret.at(i)[1] * allCs.at(iter)[1] +
2932 ret.at(i)[2] * allCs.at(iter)[2] +
2933 ret.at(i)[3] * allCs.at(iter)[3] );
2935 if ( ( ( 1.0 / 3.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
2936 ( ( 1.0 / 3.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
2942 allAnglesMet =
false;
2948 ret.emplace_back ( allCs.at(iter) );
2953 if ( ret.size() != 3 )
2957 std::cout <<
">>>>>>>> Found only " << ret.size() <<
" C3 symmetries. Searching for the missing ones." << std::endl;
2961 if ( ret.size() == 2 )
2964 std::vector< std::array<double,3> > axesVec;
2965 std::array<double,3> hArr;
2966 double axNorm, axX, axY, axZ;
2967 for (
double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
2969 for (
double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
2971 for (
double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
2973 if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) {
continue; }
2975 axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
2980 allAnglesMet =
true;
2981 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
2983 dotProduct = ( ret.at(it)[1] * axX +
2984 ret.at(it)[2] * axY +
2985 ret.at(it)[3] * axZ );
2987 if ( ( ( 1.0 / 3.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
2988 ( ( 1.0 / 3.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
2994 allAnglesMet =
false;
3005 allAnglesMet =
true;
3006 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3008 if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3009 ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3010 ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3012 allAnglesMet =
false;
3018 axesVec.emplace_back ( hArr );
3025 std::vector< std::array<double,5> > hRes;
3026 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3029 double res = cmpObj->
maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 3.0, settings );
3032 hlpArr[1] = axesVec.at(axIt)[0];
3033 hlpArr[2] = axesVec.at(axIt)[1];
3034 hlpArr[3] = axesVec.at(axIt)[2];
3036 hRes.emplace_back ( hlpArr );
3039 std::sort ( hRes.begin(), hRes.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[4] > b[4]; });
3043 ret.emplace_back ( hRes.at(0) );
3046 if ( ret.size() != 4 )
3060 std::cout <<
">>>>>>>> All C3 symmetries located." << std::endl;
3065 std::cout <<
">>>>> Searching for the three perpendicular C2 symmetry axes." << std::endl;
3069 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
3072 if ( allCs.at(iter)[0] != 2.0 )
3078 dotProduct = ( ret.at(0)[1] * allCs.at(iter)[1] +
3079 ret.at(0)[2] * allCs.at(iter)[2] +
3080 ret.at(0)[3] * allCs.at(iter)[3] );
3082 if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3083 ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3086 if ( ret.size() == 4 )
3088 ret.emplace_back ( allCs.at(iter) );
3093 allAnglesMet =
true;
3094 for (
unsigned int i = 4; i < static_cast<unsigned int> ( ret.size() ); i++ )
3096 dotProduct = ( ret.at(i)[1] * allCs.at(iter)[1] +
3097 ret.at(i)[2] * allCs.at(iter)[2] +
3098 ret.at(i)[3] * allCs.at(iter)[3] );
3100 if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3101 ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3107 allAnglesMet =
false;
3113 ret.emplace_back ( allCs.at(iter) );
3119 if ( ret.size() != 6 )
3123 std::cout <<
">>>>>>>> Found only " << ret.size() - 4 <<
" C3 symmetries. Searching for the missing ones." << std::endl;
3127 if ( ret.size() == 6 )
3130 std::vector< std::array<double,3> > axesVec;
3131 std::array<double,3> hArr;
3132 double axNorm, axX, axY, axZ;
3133 for (
double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
3135 for (
double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
3137 for (
double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
3139 if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) {
continue; }
3141 axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
3146 allAnglesMet =
true;
3147 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3149 dotProduct = ( ret.at(it)[1] * axX +
3150 ret.at(it)[2] * axY +
3151 ret.at(it)[3] * axZ );
3153 if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3154 ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3160 allAnglesMet =
false;
3171 allAnglesMet =
true;
3172 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3174 if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3175 ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3176 ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3178 allAnglesMet =
false;
3184 axesVec.emplace_back ( hArr );
3191 std::vector< std::array<double,5> > hRes;
3192 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3195 double res = cmpObj->
maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 2.0, settings );
3198 hlpArr[1] = axesVec.at(axIt)[0];
3199 hlpArr[2] = axesVec.at(axIt)[1];
3200 hlpArr[3] = axesVec.at(axIt)[2];
3202 hRes.emplace_back ( hlpArr );
3205 std::sort ( hRes.begin(), hRes.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[4] > b[4]; });
3209 ret.emplace_back ( hRes.at(0) );
3212 if ( ret.size() != 7 )
3226 std::cout <<
">>>>>>>> All C2 symmetries located." << std::endl;
3248 std::vector< std::array<double,5> > ret;
3249 std::array<double,5> hlpArr;
3252 std::cout <<
">> Searching for tetrahederal symmetry elements." << std::endl;
3261 ret.emplace_back ( hlpArr );
3265 std::cout <<
">>>>>>>> Identity element." << std::endl;
3269 for (
unsigned int iter = 0; iter < 3; iter++ )
3272 hlpArr[1] = symmAxes.at(iter)[1];
3273 hlpArr[2] = symmAxes.at(iter)[2];
3274 hlpArr[3] = symmAxes.at(iter)[3];
3276 ret.emplace_back ( hlpArr );
3280 std::cout <<
">>>>>>>> 4 clockwise C3 rotations." << std::endl;
3284 for (
unsigned int iter = 0; iter < 3; iter++ )
3287 hlpArr[1] = symmAxes.at(iter)[1];
3288 hlpArr[2] = symmAxes.at(iter)[2];
3289 hlpArr[3] = symmAxes.at(iter)[3];
3291 ret.emplace_back ( hlpArr );
3295 std::cout <<
">>>>>>>> 4 anti-clockwise C3 rotations." << std::endl;
3299 for (
unsigned int iter = 3; iter < 6; iter++ )
3302 hlpArr[1] = symmAxes.at(iter)[1];
3303 hlpArr[2] = symmAxes.at(iter)[2];
3304 hlpArr[3] = symmAxes.at(iter)[3];
3306 ret.emplace_back ( hlpArr );
3310 std::cout <<
">>>>>>>> 3 clockwise D222 rotations." << std::endl;
3314 std::cout <<
">>>>> Tetrahedral symmetry elements generated." << std::endl;
3341 std::vector< std::array<double,5> > octaSymms,
3342 std::vector< std::array<double,5> > allCs,
3343 double axisErrorTolerance,
3347 std::vector< std::array<double,5> > ret;
3348 std::array<double,5> hlpArr;
3349 bool allAnglesMet =
false;
3350 double dotProduct = 0.0;
3353 if ( octaSymms.at(0)[0] == 4.0 )
3355 ret.emplace_back ( octaSymms.at(0) );
3359 ret.emplace_back ( octaSymms.at(1) );
3364 std::cout <<
">> Searching for (cub)octahedral symmetry elements." << std::endl;
3368 std::cout <<
">>>>> Searching for the three C4 symmetry elements." << std::endl;
3372 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
3375 if ( allCs.at(iter)[0] != 4.0 )
3381 if ( ret.size() == 1 )
3383 dotProduct = ( ret.at(0)[1] * allCs.at(iter)[1] +
3384 ret.at(0)[2] * allCs.at(iter)[2] +
3385 ret.at(0)[3] * allCs.at(iter)[3] );
3387 if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3388 ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3390 ret.emplace_back ( allCs.at(iter) );
3396 if ( ret.size() == 2 )
3398 allAnglesMet =
true;
3399 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3401 dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
3402 ret.at(it)[2] * allCs.at(iter)[2] +
3403 ret.at(it)[3] * allCs.at(iter)[3] );
3405 if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3406 ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3412 allAnglesMet =
false;
3418 ret.emplace_back ( allCs.at(iter) );
3424 if ( ret.size() != 3 )
3428 std::cout <<
">>>>>>>> Found only " << ret.size() <<
" C4 symmetries. Searching for the missing ones." << std::endl;
3432 if ( ret.size() == 2 )
3435 std::vector< std::array<double,3> > axesVec;
3436 std::array<double,3> hArr;
3437 double axNorm, axX, axY, axZ;
3438 for (
double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
3440 for (
double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
3442 for (
double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
3444 if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) {
continue; }
3446 axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
3451 allAnglesMet =
true;
3452 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3454 dotProduct = ( ret.at(it)[1] * axX +
3455 ret.at(it)[2] * axY +
3456 ret.at(it)[3] * axZ );
3458 if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3459 ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3465 allAnglesMet =
false;
3476 allAnglesMet =
true;
3477 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3479 if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3480 ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3481 ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3483 allAnglesMet =
false;
3489 axesVec.emplace_back ( hArr );
3496 std::vector< std::array<double,5> > hRes;
3497 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3500 double res = cmpObj->
maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 4.0, settings );
3503 hlpArr[1] = axesVec.at(axIt)[0];
3504 hlpArr[2] = axesVec.at(axIt)[1];
3505 hlpArr[3] = axesVec.at(axIt)[2];
3507 hRes.emplace_back ( hlpArr );
3510 std::sort ( hRes.begin(), hRes.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[4] > b[4]; });
3514 ret.emplace_back ( hRes.at(0) );
3517 if ( ret.size() != 3 )
3531 std::cout <<
">>>>>>>> All C4 symmetries located." << std::endl;
3536 std::cout <<
">>>>> Searching for the four C3 symmetry elements." << std::endl;
3540 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
3543 if ( allCs.at(iter)[0] != 3.0 )
3549 allAnglesMet =
true;
3550 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3553 if ( ret.at(it)[0] == 3.0 )
3558 dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
3559 ret.at(it)[2] * allCs.at(iter)[2] +
3560 ret.at(it)[3] * allCs.at(iter)[3] );
3562 if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3563 ( ( 1.0 / sqrt ( 3.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3569 allAnglesMet =
false;
3575 ret.emplace_back ( allCs.at(iter) );
3580 if ( ret.size() != 7 )
3584 std::cout <<
">>>>>>>> Found only " << ret.size() - 3 <<
" C3 symmetries. Searching for the missing ones." << std::endl;
3588 if ( ret.size() == 6 )
3591 std::vector< std::array<double,3> > axesVec;
3592 std::array<double,3> hArr;
3593 double axNorm, axX, axY, axZ;
3594 for (
double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
3596 for (
double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
3598 for (
double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
3600 if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) {
continue; }
3602 axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
3607 allAnglesMet =
true;
3608 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3611 if ( ret.at(it)[0] == 3.0 )
3616 dotProduct = ( ret.at(it)[1] * axX +
3617 ret.at(it)[2] * axY +
3618 ret.at(it)[3] * axZ );
3620 if ( ( ( 1.0 / sqrt ( 3.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3621 ( ( 1.0 / sqrt ( 3.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3627 allAnglesMet =
false;
3632 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3635 if ( ret.at(it)[0] == 4.0 )
3640 dotProduct = ( ret.at(it)[1] * axX +
3641 ret.at(it)[2] * axY +
3642 ret.at(it)[3] * axZ );
3644 if ( ( 1.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3645 ( 1.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
3647 allAnglesMet =
false;
3662 allAnglesMet =
true;
3663 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3665 if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3666 ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3667 ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3669 allAnglesMet =
false;
3675 axesVec.emplace_back ( hArr );
3682 std::vector< std::array<double,5> > hRes;
3683 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3686 double res = cmpObj->
maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 3.0, settings );
3689 hlpArr[1] = axesVec.at(axIt)[0];
3690 hlpArr[2] = axesVec.at(axIt)[1];
3691 hlpArr[3] = axesVec.at(axIt)[2];
3693 hRes.emplace_back ( hlpArr );
3696 std::sort ( hRes.begin(), hRes.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[4] > b[4]; });
3700 ret.emplace_back ( hRes.at(0) );
3703 if ( ret.size() != 7 )
3717 std::cout <<
">>>>>>>> All C3 symmetries located." << std::endl;
3722 std::cout <<
">>>>> Searching for the six C2 symmetry elements." << std::endl;
3726 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
3729 if ( allCs.at(iter)[0] != 2.0 )
3735 allAnglesMet =
true;
3736 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3739 if ( ret.at(it)[0] != 4.0 )
3744 dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
3745 ret.at(it)[2] * allCs.at(iter)[2] +
3746 ret.at(it)[3] * allCs.at(iter)[3] );
3748 if ( ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3749 ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) ) ||
3750 ( ( ( 1.0 / sqrt ( 2.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3751 ( ( 1.0 / sqrt ( 2.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) ) )
3757 allAnglesMet =
false;
3763 ret.emplace_back ( allCs.at(iter) );
3768 if ( ret.size() != 13 )
3772 std::cout <<
">>>>>>>> Found only " << ret.size() - 13 <<
" C2 symmetries. Searching for the missing ones." << std::endl;
3776 if ( ret.size() > 10 )
3779 std::vector< std::array<double,3> > axesVec;
3780 std::array<double,3> hArr;
3781 double axNorm, axX, axY, axZ;
3782 for (
double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
3784 for (
double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
3786 for (
double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
3788 if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) {
continue; }
3790 axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
3795 allAnglesMet =
true;
3796 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
3799 if ( ( ret.at(it)[0] == 3.0 ) || ( ret.at(it)[0] == 4.0 ) )
3804 dotProduct = ( ret.at(it)[1] * axX +
3805 ret.at(it)[2] * axY +
3806 ret.at(it)[3] * axZ );
3808 if ( ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3809 ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) ) ||
3810 ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
3811 ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) ) )
3817 allAnglesMet =
false;
3828 allAnglesMet =
true;
3829 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3831 if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
3832 ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
3833 ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
3835 allAnglesMet =
false;
3839 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( ret.size() ); axIt++ )
3841 if ( ( ret.at(axIt)[0] == 3.0 ) || ( ret.at(axIt)[0] == 4.0 ) )
3846 if ( ( ( ret.at(axIt)[1] > ( axX - axisErrorTolerance ) ) && ( ret.at(axIt)[1] < ( axX + axisErrorTolerance ) ) ) &&
3847 ( ( ret.at(axIt)[2] > ( axY - axisErrorTolerance ) ) && ( ret.at(axIt)[2] < ( axY + axisErrorTolerance ) ) ) &&
3848 ( ( ret.at(axIt)[3] > ( axZ - axisErrorTolerance ) ) && ( ret.at(axIt)[3] < ( axZ + axisErrorTolerance ) ) ) )
3850 allAnglesMet =
false;
3856 axesVec.emplace_back ( hArr );
3863 std::vector< std::array<double,5> > hRes;
3864 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
3867 double res = cmpObj->
maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 2.0, settings );
3870 hlpArr[1] = axesVec.at(axIt)[0];
3871 hlpArr[2] = axesVec.at(axIt)[1];
3872 hlpArr[3] = axesVec.at(axIt)[2];
3874 hRes.emplace_back ( hlpArr );
3877 std::sort ( hRes.begin(), hRes.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[4] > b[4]; });
3879 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( hRes.size() ); axIt++ )
3881 if ( ret.size () != 13 )
3885 ret.emplace_back ( hRes.at(axIt) );
3890 if ( ret.size() != 13 )
3904 std::cout <<
">>>>>>>> All C2 symmetries located." << std::endl;
3926 std::vector< std::array<double,5> > ret;
3927 std::array<double,5> hlpArr;
3930 std::cout <<
">> Searching for octahedral symmetry elements." << std::endl;
3939 ret.emplace_back ( hlpArr );
3942 std::cout <<
">>>>>>>> Identity element." << std::endl;
3946 for (
unsigned int iter = 0; iter < 3; iter++ )
3949 hlpArr[1] = symmAxes.at(iter)[1];
3950 hlpArr[2] = symmAxes.at(iter)[2];
3951 hlpArr[3] = symmAxes.at(iter)[3];
3953 ret.emplace_back ( hlpArr );
3957 std::cout <<
">>>>>>>> 3 clockwise 180 C4 rotations." << std::endl;
3961 for (
unsigned int iter = 0; iter < 3; iter++ )
3964 hlpArr[1] = symmAxes.at(iter)[1];
3965 hlpArr[2] = symmAxes.at(iter)[2];
3966 hlpArr[3] = symmAxes.at(iter)[3];
3968 ret.emplace_back ( hlpArr );
3971 hlpArr[1] = symmAxes.at(iter)[1];
3972 hlpArr[2] = symmAxes.at(iter)[2];
3973 hlpArr[3] = symmAxes.at(iter)[3];
3975 ret.emplace_back ( hlpArr );
3979 std::cout <<
">>>>>>>> 6 (anti)clockwise +-90 C4 rotations." << std::endl;
3983 for (
unsigned int iter = 3; iter < 7; iter++ )
3986 hlpArr[1] = symmAxes.at(iter)[1];
3987 hlpArr[2] = symmAxes.at(iter)[2];
3988 hlpArr[3] = symmAxes.at(iter)[3];
3990 ret.emplace_back ( hlpArr );
3993 hlpArr[1] = symmAxes.at(iter)[1];
3994 hlpArr[2] = symmAxes.at(iter)[2];
3995 hlpArr[3] = symmAxes.at(iter)[3];
3997 ret.emplace_back ( hlpArr );
4001 std::cout <<
">>>>>>>> 8 (anti)clockwise +-120 C3 rotations." << std::endl;
4005 for (
unsigned int iter = 7; iter < 13; iter++ )
4008 hlpArr[1] = symmAxes.at(iter)[1];
4009 hlpArr[2] = symmAxes.at(iter)[2];
4010 hlpArr[3] = symmAxes.at(iter)[3];
4012 ret.emplace_back ( hlpArr );
4016 std::cout <<
">>>>>>>> 6 clockwise 120 C2 rotations." << std::endl;
4021 std::cout <<
">>>>> Octahederal symmetry elements generated." << std::endl;
4049 std::vector< std::array<double,5> > icosSymms,
4050 std::vector< std::array<double,5> > allCs,
4051 double axisErrorTolerance,
4055 std::vector< std::array<double,5> > ret;
4056 std::array<double,5> hlpArr;
4057 bool allAnglesMet =
false;
4058 double dotProduct = 0.0;
4061 if ( icosSymms.at(0)[0] == 5.0 )
4063 ret.emplace_back ( icosSymms.at(0) );
4067 ret.emplace_back ( icosSymms.at(1) );
4072 std::cout <<
">> Searching for icosahedral symmetry elements." << std::endl;
4076 std::cout <<
">>>>> Searching for the six C5 symmetry elements." << std::endl;
4079 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
4082 if ( allCs.at(iter)[0] != 5.0 )
4088 allAnglesMet =
true;
4089 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4091 dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
4092 ret.at(it)[2] * allCs.at(iter)[2] +
4093 ret.at(it)[3] * allCs.at(iter)[3] );
4095 if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4096 ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4102 allAnglesMet =
false;
4109 bool alreadyExists =
false;
4110 for (
unsigned int i = 0; i < static_cast<unsigned int> ( ret.size() ); i++ )
4112 if ( ( ( ( ret.at(i)[1] + axisErrorTolerance ) < ( allCs.at(iter)[1] ) ) && ( ( ret.at(i)[1] - axisErrorTolerance ) > ( allCs.at(iter)[1] ) ) ) &&
4113 ( ( ( ret.at(i)[2] + axisErrorTolerance ) < ( allCs.at(iter)[2] ) ) && ( ( ret.at(i)[2] - axisErrorTolerance ) > ( allCs.at(iter)[2] ) ) ) &&
4114 ( ( ( ret.at(i)[3] + axisErrorTolerance ) < ( allCs.at(iter)[3] ) ) && ( ( ret.at(i)[3] - axisErrorTolerance ) > ( allCs.at(iter)[3] ) ) ) )
4116 alreadyExists =
true;
4120 if ( !alreadyExists )
4122 ret.emplace_back ( allCs.at(iter) );
4129 if ( ret.size() != 6 )
4133 std::cout <<
">>>>>>>> Found only " << ret.size() <<
" C5 symmetries. Searching for the missing ones." << std::endl;
4137 if ( ret.size() > 2 )
4140 std::vector< std::array<double,3> > axesVec;
4141 std::array<double,3> hArr;
4142 double axNorm, axX, axY, axZ;
4143 for (
double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
4145 for (
double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
4147 for (
double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
4149 if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) {
continue; }
4151 axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
4156 allAnglesMet =
true;
4157 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4159 dotProduct = ( ret.at(it)[1] * axX +
4160 ret.at(it)[2] * axY +
4161 ret.at(it)[3] * axZ );
4163 if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4164 ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4170 allAnglesMet =
false;
4181 allAnglesMet =
true;
4182 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4184 if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
4185 ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
4186 ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
4188 allAnglesMet =
false;
4194 axesVec.emplace_back ( hArr );
4201 std::vector< std::array<double,5> > hRes;
4202 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4205 double res = cmpObj->
maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 5.0, settings );
4208 hlpArr[1] = axesVec.at(axIt)[0];
4209 hlpArr[2] = axesVec.at(axIt)[1];
4210 hlpArr[3] = axesVec.at(axIt)[2];
4212 hRes.emplace_back ( hlpArr );
4215 std::sort ( hRes.begin(), hRes.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[4] > b[4]; });
4217 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( hRes.size() ); axIt++ )
4219 if ( ret.size() != 6 )
4223 ret.emplace_back ( hRes.at(axIt) );
4228 if ( ret.size() != 6 )
4242 std::cout <<
">>>>>>>> All C5 symmetries located." << std::endl;
4247 std::cout <<
">>>>> Searching for the ten C3 symmetry elements." << std::endl;
4251 double closeC5s = 0.0;
4252 double awayC5s = 0.0;
4253 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
4256 if ( allCs.at(iter)[0] != 3.0 )
4262 allAnglesMet =
true;
4265 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4268 if ( ret.at(it)[0] != 5.0 )
4273 dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
4274 ret.at(it)[2] * allCs.at(iter)[2] +
4275 ret.at(it)[3] * allCs.at(iter)[3] );
4277 if ( ( ( sqrt ( 5.0 ) / 3.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4278 ( ( sqrt ( 5.0 ) / 3.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4283 if ( ( ( 1.0 - ( sqrt ( 5.0 ) / 3.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4284 ( ( 1.0 - ( sqrt ( 5.0 ) / 3.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4291 if ( ( closeC5s == 3.0 ) && ( awayC5s == 3.0 ) )
4293 ret.emplace_back ( allCs.at(iter) );
4298 if ( ret.size() != 16 )
4302 std::cout <<
">>>>>>>> Found only " << ret.size() - 6 <<
" C3 symmetries. Searching for the missing ones." << std::endl;
4306 if ( ret.size() > 5 )
4309 std::vector< std::array<double,3 >> axesVec;
4310 std::array<double,3> hArr;
4311 double axNorm, axX, axY, axZ;
4312 for (
double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
4314 for (
double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
4316 for (
double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
4318 if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) {
continue; }
4320 axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
4327 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4330 if ( ret.at(it)[0] != 5.0 )
4335 dotProduct = ( ret.at(it)[1] * axX +
4336 ret.at(it)[2] * axY +
4337 ret.at(it)[3] * axZ );
4339 if ( ( ( sqrt ( 5.0 ) / 3.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4340 ( ( sqrt ( 5.0 ) / 3.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4345 if ( ( ( 1.0 - ( sqrt ( 5.0 ) / 3.0 ) ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4346 ( ( 1.0 - ( sqrt ( 5.0 ) / 3.0 ) ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4353 if ( ( closeC5s == 3.0 ) && ( awayC5s == 3.0 ) )
4359 allAnglesMet =
true;
4360 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4362 if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
4363 ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
4364 ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
4366 allAnglesMet =
false;
4372 axesVec.emplace_back ( hArr );
4379 std::vector< std::array<double,5> > hRes;
4380 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4383 double res = cmpObj->
maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 3.0, settings );
4386 hlpArr[1] = axesVec.at(axIt)[0];
4387 hlpArr[2] = axesVec.at(axIt)[1];
4388 hlpArr[3] = axesVec.at(axIt)[2];
4390 hRes.emplace_back ( hlpArr );
4393 std::sort ( hRes.begin(), hRes.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[4] > b[4]; });
4395 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( hRes.size() ); axIt++ )
4397 if ( ret.size() != 16 )
4401 ret.emplace_back ( hRes.at(axIt) );
4406 if ( ret.size() != 16 )
4420 std::cout <<
">>>>>>>> All C3 symmetries located." << std::endl;
4425 std::cout <<
">>>>> Searching for the fifteen C2 symmetry elements." << std::endl;
4429 double perpC2s = 0.0;
4430 double sixthC2s = 0.0;
4431 double halfSixtbC2s = 0.0;
4432 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( allCs.size() ); iter++ )
4435 if ( allCs.at(iter)[0] != 2.0 )
4441 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4444 if ( ret.at(it)[0] != 5.0 )
4449 dotProduct = ( ret.at(it)[1] * allCs.at(iter)[1] +
4450 ret.at(it)[2] * allCs.at(iter)[2] +
4451 ret.at(it)[3] * allCs.at(iter)[3] );
4453 if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4454 ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4459 if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4460 ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4465 if ( ( ( sqrt ( 3.0 ) / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4466 ( ( sqrt ( 3.0 ) / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4468 halfSixtbC2s += 1.0;
4473 if ( ( perpC2s == 2.0 ) && ( sixthC2s == 2.0 ) && ( halfSixtbC2s == 2.0 ) )
4475 ret.emplace_back ( allCs.at(iter) );
4480 if ( ret.size() != 31 )
4484 std::cout <<
">>>>>>>> Found only " << ret.size() - 16 <<
" C2 symmetries. Searching for the missing ones." << std::endl;
4488 if ( ret.size() > 15 )
4491 std::vector< std::array<double,3> > axesVec;
4492 std::array<double,3> hArr;
4493 double axNorm, axX, axY, axZ;
4494 for (
double xIt = -1.0; xIt < 1.001; xIt += 0.2 )
4496 for (
double yIt = -1.0; yIt < 1.001; yIt += 0.2 )
4498 for (
double zIt = -1.0; zIt < 1.001; zIt += 0.2 )
4500 if ( xIt == 0.0 && yIt == 0.0 && zIt == 0.0 ) {
continue; }
4502 axNorm = sqrt ( pow ( xIt, 2.0 ) + pow ( yIt, 2.0 ) + pow ( zIt, 2.0 ) );
4510 for (
unsigned int it = 0; it < static_cast<unsigned int> ( ret.size() ); it++ )
4513 if ( ret.at(it)[0] != 5.0 )
4518 dotProduct = ( ret.at(it)[1] * axX +
4519 ret.at(it)[2] * axY +
4520 ret.at(it)[3] * axZ );
4522 if ( ( 0.0 > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4523 ( 0.0 < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4528 if ( ( ( 1.0 / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4529 ( ( 1.0 / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4534 if ( ( ( sqrt ( 3.0 ) / 2.0 ) > ( std::abs ( dotProduct ) - axisErrorTolerance ) ) &&
4535 ( ( sqrt ( 3.0 ) / 2.0 ) < ( std::abs ( dotProduct ) + axisErrorTolerance ) ) )
4537 halfSixtbC2s += 1.0;
4542 if ( ( perpC2s == 2.0 ) && ( sixthC2s == 2.0 ) && ( halfSixtbC2s == 2.0 ) )
4548 allAnglesMet =
true;
4549 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4551 if ( ( ( axesVec.at(axIt)[0] > ( axX - axisErrorTolerance ) ) && ( axesVec.at(axIt)[0] < ( axX + axisErrorTolerance ) ) ) &&
4552 ( ( axesVec.at(axIt)[1] > ( axY - axisErrorTolerance ) ) && ( axesVec.at(axIt)[1] < ( axY + axisErrorTolerance ) ) ) &&
4553 ( ( axesVec.at(axIt)[2] > ( axZ - axisErrorTolerance ) ) && ( axesVec.at(axIt)[2] < ( axZ + axisErrorTolerance ) ) ) )
4555 allAnglesMet =
false;
4561 axesVec.emplace_back ( hArr );
4568 std::vector< std::array<double,5> > hRes;
4569 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( axesVec.size() ); axIt++ )
4572 double res = cmpObj->
maxAvgPeakForSymmetry ( axesVec.at(axIt)[0], axesVec.at(axIt)[1], axesVec.at(axIt)[2], 2.0, settings );
4575 hlpArr[1] = axesVec.at(axIt)[0];
4576 hlpArr[2] = axesVec.at(axIt)[1];
4577 hlpArr[3] = axesVec.at(axIt)[2];
4579 hRes.emplace_back ( hlpArr );
4582 std::sort ( hRes.begin(), hRes.end(), [](
const std::array<double,5>& a,
const std::array<double,5>& b) {
return a[4] > b[4]; });
4584 for (
unsigned int axIt = 0; axIt < static_cast<unsigned int> ( hRes.size() ); axIt++ )
4586 if ( ret.size() != 31 )
4590 ret.emplace_back ( hRes.at(axIt) );
4595 if ( ret.size() != 31 )
4609 std::cout <<
">>>>>>>> All C2 symmetries located." << std::endl;
4631 std::vector< std::array<double,5> > ret;
4632 std::array<double,5> hlpArr;
4635 std::cout <<
">> Searching for icosahedral symmetry elements." << std::endl;
4644 ret.emplace_back ( hlpArr );
4647 std::cout <<
">>>>>>>> Identity element." << std::endl;
4651 for (
unsigned int iter = 0; iter < 6; iter++ )
4654 hlpArr[1] = symmAxes.at(iter)[1];
4655 hlpArr[2] = symmAxes.at(iter)[2];
4656 hlpArr[3] = symmAxes.at(iter)[3];
4658 ret.emplace_back ( hlpArr );
4661 hlpArr[1] = symmAxes.at(iter)[1];
4662 hlpArr[2] = symmAxes.at(iter)[2];
4663 hlpArr[3] = symmAxes.at(iter)[3];
4665 ret.emplace_back ( hlpArr );
4669 std::cout <<
">>>>>>>> 12 (anti)clockwise 144 C5 rotations." << std::endl;
4673 for (
unsigned int iter = 0; iter < 6; iter++ )
4676 hlpArr[1] = symmAxes.at(iter)[1];
4677 hlpArr[2] = symmAxes.at(iter)[2];
4678 hlpArr[3] = symmAxes.at(iter)[3];
4680 ret.emplace_back ( hlpArr );
4683 hlpArr[1] = symmAxes.at(iter)[1];
4684 hlpArr[2] = symmAxes.at(iter)[2];
4685 hlpArr[3] = symmAxes.at(iter)[3];
4687 ret.emplace_back ( hlpArr );
4691 std::cout <<
">>>>>>>> 12 (anti)clockwise 144 C5 rotations." << std::endl;
4695 for (
unsigned int iter = 6; iter < 16; iter++ )
4698 hlpArr[1] = symmAxes.at(iter)[1];
4699 hlpArr[2] = symmAxes.at(iter)[2];
4700 hlpArr[3] = symmAxes.at(iter)[3];
4702 ret.emplace_back ( hlpArr );
4705 hlpArr[1] = symmAxes.at(iter)[1];
4706 hlpArr[2] = symmAxes.at(iter)[2];
4707 hlpArr[3] = symmAxes.at(iter)[3];
4709 ret.emplace_back ( hlpArr );
4713 std::cout <<
">>>>>>>> 20 (anti)clockwise 120 C3 rotations." << std::endl;
4717 for (
unsigned int iter = 16; iter < 31; iter++ )
4720 hlpArr[1] = symmAxes.at(iter)[1];
4721 hlpArr[2] = symmAxes.at(iter)[2];
4722 hlpArr[3] = symmAxes.at(iter)[3];
4724 ret.emplace_back ( hlpArr );
4728 std::cout <<
">>>>>>>> 15 clockwise 180 C2 rotations." << std::endl;
4733 std::cout <<
">>>>> Icosahedral symmetry elements generated." << std::endl;
4761 std::array<double,4> ret = std::array<double,4> { { 0.0, 0.0, 0.0, 0.0 } };
4765 double newXSampling = std::min ( obj1->_xSamplingRate, obj2->_xSamplingRate );
4766 double newYSampling = std::min ( obj1->_ySamplingRate, obj2->_ySamplingRate );
4767 double newZSampling = std::min ( obj1->_zSamplingRate, obj2->_zSamplingRate );
4776 newU = std::floor ( obj1->_xRange / newXSampling );
4777 newV = std::floor ( obj1->_yRange / newYSampling );
4778 newW = std::floor ( obj1->_zRange / newZSampling );
4780 double newXRange =
static_cast<double> ( ( newU ) * newXSampling );
4781 double newYRange =
static_cast<double> ( ( newV ) * newYSampling );
4782 double newZRange =
static_cast<double> ( ( newW ) * newZSampling );
4784 int uDiff = newU - obj1->_maxMapU;
4785 int vDiff = newV - obj1->_maxMapV;
4786 int wDiff = newW - obj1->_maxMapW;
4795 if ( uDiff % 2 == 0 )
4797 newXTo = obj1->_xTo + (uDiff/2);
4798 newXFrom = obj1->_xFrom - (uDiff/2);
4802 newXTo = obj1->_xTo + ((uDiff+1)/2);
4803 newXFrom = obj1->_xFrom - ((uDiff+1)/2);
4806 if ( vDiff % 2 == 0 )
4808 newYTo = obj1->_yTo + (vDiff/2);
4809 newYFrom = obj1->_yFrom - (vDiff/2);
4813 newYTo = obj1->_yTo + ((vDiff+1)/2);
4814 newYFrom = obj1->_yFrom - ((vDiff+1)/2);
4817 if ( wDiff % 2 == 0 )
4819 newZTo = obj1->_zTo + (wDiff/2);
4820 newZFrom = obj1->_zFrom - (wDiff/2);
4824 newZTo = obj1->_zTo + ((wDiff+1)/2);
4825 newZFrom = obj1->_zFrom - ((wDiff+1)/2);
4829 if ( ( ( obj1->_xFrom >= 0 ) && ( newXFrom < 0 ) ) &&
4830 ( ( obj1->_xTo >= 0 ) && ( newXTo > 0 ) ) )
4835 if ( ( ( obj1->_yFrom >= 0 ) && ( newYFrom < 0 ) ) &&
4836 ( ( obj1->_yTo >= 0 ) && ( newYTo > 0 ) ) )
4841 if ( ( ( obj1->_zFrom >= 0 ) && ( newZFrom < 0 ) ) &&
4842 ( ( obj1->_zTo >= 0 ) && ( newZTo > 0 ) ) )
4848 if ( newU % 2 != 0 )
4854 if ( newV % 2 != 0 )
4860 if ( newW % 2 != 0 )
4867 double *hlpMap =
new double [(newU+1) * (newV+1) * (newW+1)];
4877 bool xFound =
false;
4878 bool yFound =
false;
4879 bool zFound =
false;
4881 std::array<double,4> p000, p001, p010, p011, p100, p101, p110, p111;
4882 std::array<double,3> p00, p01, p10, p11;
4883 std::array<double,2> p0, p1;
4884 for (
int uIt = 0; uIt < (newU+1); uIt++ )
4886 for (
int vIt = 0; vIt < (newV+1); vIt++ )
4888 for (
int wIt = 0; wIt < (newW+1); wIt++ )
4891 arrPos = wIt + (newW+1) * ( vIt + (newV+1) * uIt );
4898 newX =
static_cast<double> ( uIt ) * newXSampling;
4899 newY =
static_cast<double> ( vIt ) * newYSampling;
4900 newZ =
static_cast<double> ( wIt ) * newZSampling;
4903 for (
int iter = 0; iter < static_cast<int> (obj1->_maxMapU+1); iter++ )
4905 if ( ( newX >= ( static_cast<double> ( iter ) * obj1->_xSamplingRate ) ) && ( newX < ( static_cast<double> ( iter + 1 ) * obj1->_xSamplingRate ) ) )
4913 if ( ( oldX + 1 ) >=
static_cast<int> (obj1->_maxMapU+1) ) { xOver =
true; }
4914 if ( !xFound ) { hlpMap[arrPos] = 0.0;
continue; }
4916 for (
int iter = 0; iter < static_cast<int> (obj1->_maxMapV+1); iter++ )
4918 if ( ( newY >= ( static_cast<double> ( iter ) * obj1->_ySamplingRate ) ) && ( newY < ( static_cast<double> ( iter + 1 ) * obj1->_ySamplingRate ) ) )
4926 if ( ( oldY + 1 ) >=
static_cast<int> (obj1->_maxMapV+1) ) { yOver =
true; }
4927 if ( !yFound ) { hlpMap[arrPos] = 0.0;
continue; }
4929 for (
int iter = 0; iter < static_cast<int> (obj1->_maxMapW+1); iter++ )
4931 if ( ( newZ >= ( static_cast<double> ( iter ) * obj1->_zSamplingRate ) ) && ( newZ < ( static_cast<double> ( iter + 1 ) * obj1->_zSamplingRate ) ) )
4939 if ( ( oldZ + 1 ) >=
static_cast<int> (obj1->_maxMapW+1) ) { zOver =
true; }
4940 if ( !zFound ) { hlpMap[arrPos] = 0.0;
continue; }
4946 hlpPos = p000[2] + (obj1->_maxMapW+1) * ( p000[1] + (obj1->_maxMapV+1) * p000[0] );
4947 p000[3] = obj1->_densityMapCor[hlpPos];
4952 if ( zOver ) { p001[2] = 0.0; }
4953 hlpPos = p001[2] + (obj1->_maxMapW+1) * ( p001[1] + (obj1->_maxMapV+1) * p001[0] );
4954 p001[3] = obj1->_densityMapCor[hlpPos];
4959 if ( yOver ) { p010[1] = 0.0; }
4960 hlpPos = p010[2] + (obj1->_maxMapW+1) * ( p010[1] + (obj1->_maxMapV+1) * p010[0] );
4961 p010[3] = obj1->_densityMapCor[hlpPos];
4966 if ( yOver ) { p011[1] = 0.0; }
4967 if ( zOver ) { p011[2] = 0.0; }
4968 hlpPos = p011[2] + (obj1->_maxMapW+1) * ( p011[1] + (obj1->_maxMapV+1) * p011[0] );
4969 p011[3] = obj1->_densityMapCor[hlpPos];
4974 if ( xOver ) { p100[0] = 0.0; }
4975 hlpPos = p100[2] + (obj1->_maxMapW+1) * ( p100[1] + (obj1->_maxMapV+1) * p100[0] );
4976 p100[3] = obj1->_densityMapCor[hlpPos];
4981 if ( xOver ) { p101[0] = 0.0; }
4982 if ( zOver ) { p101[2] = 0.0; }
4983 hlpPos = p101[2] + (obj1->_maxMapW+1) * ( p101[1] + (obj1->_maxMapV+1) * p101[0] );
4984 p101[3] = obj1->_densityMapCor[hlpPos];
4989 if ( xOver ) { p110[0] = 0.0; }
4990 if ( yOver ) { p110[1] = 0.0; }
4991 hlpPos = p110[2] + (obj1->_maxMapW+1) * ( p110[1] + (obj1->_maxMapV+1) * p110[0] );
4992 p110[3] = obj1->_densityMapCor[hlpPos];
4997 if ( xOver ) { p111[0] = 0.0; }
4998 if ( yOver ) { p111[1] = 0.0; }
4999 if ( zOver ) { p111[2] = 0.0; }
5000 hlpPos = p111[2] + (obj1->_maxMapW+1) * ( p111[1] + (obj1->_maxMapV+1) * p111[0] );
5001 p111[3] = obj1->_densityMapCor[hlpPos];
5004 xd = 1.0 - ( ( newX - (
static_cast<double> ( oldX ) * obj1->_xSamplingRate ) ) / obj1->_xSamplingRate );
5005 p00[0] = p000[1]; p00[1] = p000[2]; p00[2] = ( xd * p000[3] ) + ( (1.0 - xd) * p100[3] );
5006 p01[0] = p001[1]; p01[1] = p001[2]; p01[2] = ( xd * p001[3] ) + ( (1.0 - xd) * p101[3] );
5007 p10[0] = p010[1]; p10[1] = p010[2]; p10[2] = ( xd * p010[3] ) + ( (1.0 - xd) * p110[3] );
5008 p11[0] = p011[1]; p11[1] = p011[2]; p11[2] = ( xd * p011[3] ) + ( (1.0 - xd) * p111[3] );
5010 yd = 1.0 - ( ( newY - (
static_cast<double> ( oldY ) * obj1->_ySamplingRate ) ) / obj1->_ySamplingRate );
5011 p0[0] = p00[1]; p0[1] = ( yd * p00[2] ) + ( (1.0 - yd) * p10[2] );
5012 p1[0] = p01[1]; p1[1] = ( yd * p01[2] ) + ( (1.0 - yd) * p11[2] );
5014 zd = 1.0 - ( ( newZ - (
static_cast<double> ( oldZ ) * obj1->_zSamplingRate ) ) / obj1->_zSamplingRate );
5015 hlpMap[arrPos] = ( zd * p0[1] ) + ( (1.0 - zd) * p1[1] );
5021 delete[] obj1->_densityMapCor;
5022 obj1->_densityMapCor =
new double [(newU+1) * (newV+1) * (newW+1)];
5023 for (
int iter = 0; iter < static_cast<int> ( (newU+1) * (newV+1) * (newW+1) ); iter++ )
5025 obj1->_densityMapCor[iter] = hlpMap[iter];
5032 int xMov1 = std::round ( ( ( obj1->_xFrom * obj1->_xSamplingRate ) - ( newXFrom * newXSampling ) ) / newXSampling );
5033 int yMov1 = std::round ( ( ( obj1->_yFrom * obj1->_ySamplingRate ) - ( newYFrom * newYSampling ) ) / newYSampling );
5034 int zMov1 = std::round ( ( ( obj1->_zFrom * obj1->_zSamplingRate ) - ( newZFrom * newZSampling ) ) / newZSampling );
5036 obj1->_maxMapU = newU;
5037 obj1->_maxMapV = newV;
5038 obj1->_maxMapW = newW;
5040 obj1->_xFrom = newXFrom + xMov1;
5041 obj1->_yFrom = newYFrom + yMov1;
5042 obj1->_zFrom = newZFrom + zMov1;
5044 obj1->_xTo = newXTo + xMov1;
5045 obj1->_yTo = newYTo + yMov1;
5046 obj1->_zTo = newZTo + zMov1;
5048 obj1->_xRange = newXRange;
5049 obj1->_yRange = newYRange;
5050 obj1->_zRange = newZRange;
5052 obj1->_xSamplingRate = newXSampling;
5053 obj1->_ySamplingRate = newYSampling;
5054 obj1->_zSamplingRate = newZSampling;
5060 newU = std::floor ( obj2->_xRange / newXSampling );
5061 newV = std::floor ( obj2->_yRange / newYSampling );
5062 newW = std::floor ( obj2->_zRange / newZSampling );
5064 double newXRange =
static_cast<double> ( ( newU + 1 ) * newXSampling );
5065 double newYRange =
static_cast<double> ( ( newV + 1 ) * newYSampling );
5066 double newZRange =
static_cast<double> ( ( newW + 1 ) * newZSampling );
5068 int uDiff = newU - obj2->_maxMapU;
5069 int vDiff = newV - obj2->_maxMapV;
5070 int wDiff = newW - obj2->_maxMapW;
5079 if ( uDiff % 2 == 0 )
5081 newXTo = obj2->_xTo + (uDiff/2);
5082 newXFrom = obj2->_xFrom - (uDiff/2);
5086 newXTo = obj2->_xTo + ((uDiff+1)/2);
5087 newXFrom = obj2->_xFrom - ((uDiff+1)/2);
5090 if ( vDiff % 2 == 0 )
5092 newYTo = obj2->_yTo + (vDiff/2);
5093 newYFrom = obj2->_yFrom - (vDiff/2);
5097 newYTo = obj2->_yTo + ((vDiff+1)/2);
5098 newYFrom = obj2->_yFrom - ((vDiff+1)/2);
5101 if ( wDiff % 2 == 0 )
5103 newZTo = obj2->_zTo + (wDiff/2);
5104 newZFrom = obj2->_zFrom - (wDiff/2);
5108 newZTo = obj2->_zTo + ((wDiff+1)/2);
5109 newZFrom = obj2->_zFrom - ((wDiff+1)/2);
5113 if ( ( ( obj2->_xFrom >= 0 ) && ( newXFrom < 0 ) ) &&
5114 ( ( obj2->_xTo >= 0 ) && ( newXTo > 0 ) ) )
5119 if ( ( ( obj2->_yFrom >= 0 ) && ( newYFrom < 0 ) ) &&
5120 ( ( obj2->_yTo >= 0 ) && ( newYTo > 0 ) ) )
5125 if ( ( ( obj2->_zFrom >= 0 ) && ( newZFrom < 0 ) ) &&
5126 ( ( obj2->_zTo >= 0 ) && ( newZTo > 0 ) ) )
5132 if ( newU % 2 == 0 )
5138 if ( newV % 2 == 0 )
5144 if ( newW % 2 == 0 )
5151 double *hlpMap =
new double [(newU+1) * (newV+1) * (newW+1)];
5161 bool xFound =
false;
5162 bool yFound =
false;
5163 bool zFound =
false;
5165 std::array<double,4> p000, p001, p010, p011, p100, p101, p110, p111;
5166 std::array<double,3> p00, p01, p10, p11;
5167 std::array<double,2> p0, p1;
5169 for (
int uIt = 0; uIt < (newU+1); uIt++ )
5171 for (
int vIt = 0; vIt < (newV+1); vIt++ )
5173 for (
int wIt = 0; wIt < (newW+1); wIt++ )
5176 arrPos = wIt + (newW+1) * ( vIt + (newV+1) * uIt );
5183 newX =
static_cast<double> ( uIt ) * newXSampling;
5184 newY =
static_cast<double> ( vIt ) * newYSampling;
5185 newZ =
static_cast<double> ( wIt ) * newZSampling;
5188 for (
int iter = 0; iter < static_cast<int> (obj2->_maxMapU+1); iter++ )
5190 if ( ( newX >= ( static_cast<double> ( iter ) * obj2->_xSamplingRate ) ) && ( newX < ( static_cast<double> ( iter + 1 ) * obj2->_xSamplingRate ) ) )
5198 if ( ( oldX + 1 ) >=
static_cast<int> (obj2->_maxMapU+1) ) { xOver =
true; }
5199 if ( !xFound ) { hlpMap[arrPos] = 0.0;
continue; }
5201 for (
int iter = 0; iter < static_cast<int> (obj2->_maxMapV+1); iter++ )
5203 if ( ( newY >= ( static_cast<double> ( iter ) * obj2->_ySamplingRate ) ) && ( newY < ( static_cast<double> ( iter + 1 ) * obj2->_ySamplingRate ) ) )
5211 if ( ( oldY + 1 ) >=
static_cast<int> (obj2->_maxMapV+1) ) { yOver =
true; }
5212 if ( !yFound ) { hlpMap[arrPos] = 0.0;
continue; }
5214 for (
int iter = 0; iter < static_cast<int> (obj2->_maxMapW+1); iter++ )
5216 if ( ( newZ >= ( static_cast<double> ( iter ) * obj2->_zSamplingRate ) ) && ( newZ < ( static_cast<double> ( iter + 1 ) * obj2->_zSamplingRate ) ) )
5224 if ( ( oldZ + 1 ) >=
static_cast<int> (obj2->_maxMapW+1) ) { zOver =
true; }
5225 if ( !zFound ) { hlpMap[arrPos] = 0.0;
continue; }
5231 hlpPos = p000[2] + (obj2->_maxMapW+1) * ( p000[1] + (obj2->_maxMapV+1) * p000[0] );
5232 p000[3] = obj2->_densityMapCor[hlpPos];
5237 if ( zOver ) { p001[2] = 0.0; }
5238 hlpPos = p001[2] + (obj2->_maxMapW+1) * ( p001[1] + (obj2->_maxMapV+1) * p001[0] );
5239 p001[3] = obj2->_densityMapCor[hlpPos];
5244 if ( yOver ) { p010[1] = 0.0; }
5245 hlpPos = p010[2] + (obj2->_maxMapW+1) * ( p010[1] + (obj2->_maxMapV+1) * p010[0] );
5246 p010[3] = obj2->_densityMapCor[hlpPos];
5251 if ( yOver ) { p011[1] = 0.0; }
5252 if ( zOver ) { p011[2] = 0.0; }
5253 hlpPos = p011[2] + (obj2->_maxMapW+1) * ( p011[1] + (obj2->_maxMapV+1) * p011[0] );
5254 p011[3] = obj2->_densityMapCor[hlpPos];
5259 if ( xOver ) { p100[0] = 0.0; }
5260 hlpPos = p100[2] + (obj2->_maxMapW+1) * ( p100[1] + (obj2->_maxMapV+1) * p100[0] );
5261 p100[3] = obj2->_densityMapCor[hlpPos];
5266 if ( xOver ) { p101[0] = 0.0; }
5267 if ( zOver ) { p101[2] = 0.0; }
5268 hlpPos = p101[2] + (obj2->_maxMapW+1) * ( p101[1] + (obj2->_maxMapV+1) * p101[0] );
5269 p101[3] = obj2->_densityMapCor[hlpPos];
5274 if ( xOver ) { p110[0] = 0.0; }
5275 if ( yOver ) { p110[1] = 0.0; }
5276 hlpPos = p110[2] + (obj2->_maxMapW+1) * ( p110[1] + (obj2->_maxMapV+1) * p110[0] );
5277 p110[3] = obj2->_densityMapCor[hlpPos];
5282 if ( xOver ) { p111[0] = 0.0; }
5283 if ( yOver ) { p111[1] = 0.0; }
5284 if ( zOver ) { p111[2] = 0.0; }
5285 hlpPos = p111[2] + (obj2->_maxMapW+1) * ( p111[1] + (obj2->_maxMapV+1) * p111[0] );
5286 p111[3] = obj2->_densityMapCor[hlpPos];
5289 xd = 1.0 - ( ( newX - (
static_cast<double> ( oldX ) * obj2->_xSamplingRate ) ) / obj2->_xSamplingRate );
5290 p00[0] = p000[1]; p00[1] = p000[2]; p00[2] = ( xd * p000[3] ) + ( (1.0 - xd) * p100[3] );
5291 p01[0] = p001[1]; p01[1] = p001[2]; p01[2] = ( xd * p001[3] ) + ( (1.0 - xd) * p101[3] );
5292 p10[0] = p010[1]; p10[1] = p010[2]; p10[2] = ( xd * p010[3] ) + ( (1.0 - xd) * p110[3] );
5293 p11[0] = p011[1]; p11[1] = p011[2]; p11[2] = ( xd * p011[3] ) + ( (1.0 - xd) * p111[3] );
5295 yd = 1.0 - ( ( newY - (
static_cast<double> ( oldY ) * obj2->_ySamplingRate ) ) / obj2->_ySamplingRate );
5296 p0[0] = p00[1]; p0[1] = ( yd * p00[2] ) + ( (1.0 - yd) * p10[2] );
5297 p1[0] = p01[1]; p1[1] = ( yd * p01[2] ) + ( (1.0 - yd) * p11[2] );
5299 zd = 1.0 - ( ( newZ - (
static_cast<double> ( oldZ ) * obj2->_zSamplingRate ) ) / obj2->_zSamplingRate );
5300 hlpMap[arrPos] = ( zd * p0[1] ) + ( (1.0 - zd) * p1[1] );
5306 delete[] obj2->_densityMapCor;
5307 obj2->_densityMapCor =
new double [(newU+1) * (newV+1) * (newW+1)];
5308 for (
int iter = 0; iter < static_cast<int> ( (newU+1) * (newV+1) * (newW+1) ); iter++ )
5310 obj2->_densityMapCor[iter] = hlpMap[iter];
5317 *ob2XMov = std::round ( ( ( obj1->_xFrom * obj1->_xSamplingRate ) - ( newXFrom * newXSampling ) ) / newXSampling );
5318 *ob2YMov = std::round ( ( ( obj1->_yFrom * obj1->_ySamplingRate ) - ( newYFrom * newYSampling ) ) / newYSampling );
5319 *ob2ZMov = std::round ( ( ( obj1->_zFrom * obj1->_zSamplingRate ) - ( newZFrom * newZSampling ) ) / newZSampling );
5322 obj2->_maxMapU = newU;
5323 obj2->_maxMapV = newV;
5324 obj2->_maxMapW = newW;
5326 obj2->_xFrom = newXFrom + (*ob2XMov);
5327 obj2->_yFrom = newYFrom + (*ob2YMov);
5328 obj2->_zFrom = newZFrom + (*ob2ZMov);
5330 obj2->_xTo = newXTo + (*ob2XMov);
5331 obj2->_yTo = newYTo + (*ob2YMov);
5332 obj2->_zTo = newZTo + (*ob2ZMov);
5334 obj2->_xRange = newXRange;
5335 obj2->_yRange = newYRange;
5336 obj2->_zRange = newZRange;
5338 obj2->_xSamplingRate = newXSampling;
5339 obj2->_ySamplingRate = newYSampling;
5340 obj2->_zSamplingRate = newZSampling;
5342 *ob2XMov *= obj2->_xSamplingRate;
5343 *ob2YMov *= obj2->_ySamplingRate;
5344 *ob2ZMov *= obj2->_zSamplingRate;
5348 newU = std::max ( obj1->_maxMapU, obj2->_maxMapU );
5349 newV = std::max ( obj1->_maxMapV, obj2->_maxMapV );
5350 newW = std::max ( obj1->_maxMapW, obj2->_maxMapW );
5352 if ( ( static_cast<int> ( obj1->_maxMapU ) != newU ) || ( static_cast<int> ( obj1->_maxMapV ) != newV ) || ( static_cast<int> ( obj1->_maxMapW ) != newW ) )
5355 double *hlpMap =
new double [(newU+1) * (newV+1) * (newW+1)];
5358 for (
int uIt = 0; uIt < ( newU + 1 ); uIt++ )
5360 for (
int vIt = 0; vIt < ( newV + 1 ); vIt++ )
5362 for (
int wIt = 0; wIt < ( newW + 1 ); wIt++ )
5364 arrPos = wIt + (newW+1) * ( vIt + (newV+1) * uIt );
5365 hlpPos = wIt + (obj1->_maxMapW+1) * ( vIt + (obj1->_maxMapV+1) * uIt );
5367 if ( ( uIt < ( static_cast<int> ( obj1->_maxMapU ) + 1 ) ) &&
5368 ( vIt < ( static_cast<int> ( obj1->_maxMapV ) + 1 ) ) &&
5369 ( wIt < ( static_cast<int> ( obj1->_maxMapW ) + 1 ) ) )
5371 hlpMap[arrPos] = obj1->_densityMapCor[hlpPos];
5375 hlpMap[arrPos] = 0.0;
5382 delete[] obj1->_densityMapCor;
5383 obj1->_densityMapCor =
new double [(newU+1) * (newV+1) * (newW+1)];
5384 for (
int iter = 0; iter < static_cast<int> ( (newU+1) * (newV+1) * (newW+1) ); iter++ )
5386 obj1->_densityMapCor[iter] = hlpMap[iter];
5393 obj1->_xTo += newU - obj1->_maxMapU;
5394 obj1->_yTo += newV - obj1->_maxMapV;
5395 obj1->_zTo += newW - obj1->_maxMapW;
5397 obj1->_maxMapU = newU;
5398 obj1->_maxMapV = newV;
5399 obj1->_maxMapW = newW;
5401 obj1->_xRange = ( newU + 1 ) * obj1->_xSamplingRate;
5402 obj1->_yRange = ( newV + 1 ) * obj1->_ySamplingRate;
5403 obj1->_zRange = ( newW + 1 ) * obj1->_zSamplingRate;
5406 if ( ( static_cast<int> ( obj2->_maxMapU ) != newU ) || ( static_cast<int> ( obj2->_maxMapV ) != newV ) || ( static_cast<int> ( obj2->_maxMapW ) != newW ) )
5409 double *hlpMap =
new double [(newU+1) * (newV+1) * (newW+1)];
5412 for (
int uIt = 0; uIt < ( newU + 1 ); uIt++ )
5414 for (
int vIt = 0; vIt < ( newV + 1 ); vIt++ )
5416 for (
int wIt = 0; wIt < ( newW + 1 ); wIt++ )
5418 arrPos = wIt + (newW+1) * ( vIt + (newV+1) * uIt );
5419 hlpPos = wIt + (obj2->_maxMapW+1) * ( vIt + (obj2->_maxMapV+1) * uIt );
5421 if ( ( uIt < ( static_cast<int> ( obj2->_maxMapU ) + 1 ) ) &&
5422 ( vIt < ( static_cast<int> ( obj2->_maxMapV ) + 1 ) ) &&
5423 ( wIt < ( static_cast<int> ( obj2->_maxMapW ) + 1 ) ) )
5425 hlpMap[arrPos] = obj2->_densityMapCor[hlpPos];
5429 hlpMap[arrPos] = 0.0;
5436 delete[] obj2->_densityMapCor;
5437 obj2->_densityMapCor =
new double [(newU+1) * (newV+1) * (newW+1)];
5438 for (
int iter = 0; iter < static_cast<int> ( (newU+1) * (newV+1) * (newW+1) ); iter++ )
5440 obj2->_densityMapCor[iter] = hlpMap[iter];
5447 obj2->_xTo += newU - obj2->_maxMapU;
5448 obj2->_yTo += newV - obj2->_maxMapV;
5449 obj2->_zTo += newW - obj2->_maxMapW;
5451 obj2->_maxMapU = newU;
5452 obj2->_maxMapV = newV;
5453 obj2->_maxMapW = newW;
5455 obj2->_xRange = ( newU + 1 ) * obj2->_xSamplingRate;
5456 obj2->_yRange = ( newV + 1 ) * obj2->_ySamplingRate;
5457 obj2->_zRange = ( newW + 1 ) * obj2->_zSamplingRate;
5461 int minU = std::min ( (obj1->_maxMapU+1), (obj2->_maxMapU+1) );
5462 int minV = std::min ( (obj1->_maxMapV+1), (obj2->_maxMapV+1) );
5463 int minW = std::min ( (obj1->_maxMapW+1), (obj2->_maxMapW+1) );
5465 fftw_complex *tmpIn1 =
new fftw_complex[(obj1->_maxMapU+1) * (obj1->_maxMapV+1) * (obj1->_maxMapW+1)];
5466 fftw_complex *tmpOut1 =
new fftw_complex[(obj1->_maxMapU+1) * (obj1->_maxMapV+1) * (obj1->_maxMapW+1)];
5467 fftw_complex *tmpIn2 =
new fftw_complex[(obj2->_maxMapU+1) * (obj2->_maxMapV+1) * (obj2->_maxMapW+1)];
5468 fftw_complex *tmpOut2 =
new fftw_complex[(obj2->_maxMapU+1) * (obj2->_maxMapV+1) * (obj2->_maxMapW+1)];
5469 fftw_complex *resIn =
new fftw_complex[minU * minV * minW];
5470 fftw_complex *resOut =
new fftw_complex[minU * minV * minW];
5473 fftw_plan forwardFourierObj1;
5474 fftw_plan forwardFourierObj2;
5475 fftw_plan inverseFourierCombo;
5477 forwardFourierObj1 = fftw_plan_dft_3d ( (obj1->_maxMapU+1), (obj1->_maxMapV+1), (obj1->_maxMapW+1), tmpIn1, tmpOut1, FFTW_FORWARD, FFTW_ESTIMATE );
5478 forwardFourierObj2 = fftw_plan_dft_3d ( (obj2->_maxMapU+1), (obj2->_maxMapV+1), (obj2->_maxMapW+1), tmpIn2, tmpOut2, FFTW_FORWARD, FFTW_ESTIMATE );
5479 inverseFourierCombo = fftw_plan_dft_3d ( minU, minV, minW, resOut, resIn, FFTW_BACKWARD, FFTW_ESTIMATE );
5482 for (
int iter = 0; iter < static_cast<int> ( (obj1->_maxMapU+1) * (obj1->_maxMapV+1) * (obj1->_maxMapW+1) ); iter++ )
5484 tmpIn1[iter][0] = obj1->_densityMapCor[iter];
5485 tmpIn1[iter][1] = 0.0;
5488 for (
int iter = 0; iter < static_cast<int> ( (obj2->_maxMapU+1) * (obj2->_maxMapV+1) * (obj2->_maxMapW+1) ); iter++ )
5490 tmpIn2[iter][0] = obj2->_densityMapCor[iter];
5491 tmpIn2[iter][1] = 0.0;
5495 fftw_execute ( forwardFourierObj1 );
5496 fftw_execute ( forwardFourierObj2 );
5499 double normFactor =
static_cast<double> ( minU * minV * minW );
5500 std::array<double,2> hlpArr = std::array<double,2> { { 2.0, 2.0 } };
5502 int uo2 = (obj1->_maxMapU+1) / 2;
5503 int vo2 = (obj1->_maxMapV+1) / 2;
5504 int wo2 = (obj1->_maxMapW+1) / 2;
5505 int muo2 = (minU+1) / 2;
5506 int mvo2 = (minV+1) / 2;
5507 int mwo2 = (minW+1) / 2;
5509 for (
int uIt = 0; uIt < static_cast<int> ( obj1->_maxMapU+1 ); uIt++ )
5511 for (
int vIt = 0; vIt < static_cast<int> ( obj1->_maxMapV+1 ); vIt++ )
5513 for (
int wIt = 0; wIt < static_cast<int> ( obj1->_maxMapW+1 ); wIt++ )
5515 if ( uIt > uo2 ) { h1 = uIt - (obj1->_maxMapU+1); }
else { h1 = uIt; }
5516 if ( vIt > vo2 ) { k1 = vIt - (obj1->_maxMapV+1); }
else { k1 = vIt; }
5517 if ( wIt > wo2 ) { l1 = wIt - (obj1->_maxMapW+1); }
else { l1 = wIt; }
5519 if ( ( std::abs ( h1 ) <= muo2 ) && ( std::abs ( k1 ) <= mvo2 ) && ( std::abs ( l1 ) <= mwo2 ) )
5521 if ( h1 < 0 ) { h1 += minU; }
5522 if ( k1 < 0 ) { k1 += minV; }
5523 if ( l1 < 0 ) { l1 += minW; }
5525 hlpPos = l1 + minW * ( k1 + minV * h1 );
5526 arrPos = wIt + (obj1->_maxMapW+1) * ( vIt + (obj1->_maxMapV+1) * uIt );
5527 resIn[hlpPos][0] = tmpOut1[arrPos][0];
5528 resIn[hlpPos][1] = tmpOut1[arrPos][1];
5534 uo2 = (obj2->_maxMapU+1) / 2;
5535 vo2 = (obj2->_maxMapV+1) / 2;
5536 wo2 = (obj2->_maxMapW+1) / 2;
5537 for (
int uIt = 0; uIt < static_cast<int> ( obj2->_maxMapU+1 ); uIt++ )
5539 for (
int vIt = 0; vIt < static_cast<int> ( obj2->_maxMapV+1 ); vIt++ )
5541 for (
int wIt = 0; wIt < static_cast<int> ( obj2->_maxMapW+1 ); wIt++ )
5543 if ( uIt > uo2 ) { h1 = uIt - (obj2->_maxMapU+1); }
else { h1 = uIt; }
5544 if ( vIt > vo2 ) { k1 = vIt - (obj2->_maxMapV+1); }
else { k1 = vIt; }
5545 if ( wIt > wo2 ) { l1 = wIt - (obj2->_maxMapW+1); }
else { l1 = wIt; }
5547 if ( ( std::abs ( h1 ) <= muo2 ) && ( std::abs ( k1 ) <= mvo2 ) && ( std::abs ( l1 ) <= mwo2 ) )
5549 if ( h1 < 0 ) { h1 += minU; }
5550 if ( k1 < 0 ) { k1 += minV; }
5551 if ( l1 < 0 ) { l1 += minW; }
5553 hlpPos = l1 + minW * ( k1 + minV * h1 );
5554 arrPos = wIt + (obj2->_maxMapW+1) * ( vIt + (obj2->_maxMapV+1) * uIt );
5555 hlpArr = ProSHADE_internal_misc::complexMultiplicationConjug ( &resIn[hlpPos][0],
5557 &tmpOut2[arrPos][0],
5558 &tmpOut2[arrPos][1] );
5559 resOut[hlpPos][0] = hlpArr[0] / normFactor;
5560 resOut[hlpPos][1] = hlpArr[1] / normFactor;
5567 fftw_execute ( inverseFourierCombo );
5570 double mapPeak = 0.0;
5571 for (
int uIt = 0; uIt < minU; uIt++ )
5573 for (
int vIt = 0; vIt < minV; vIt++ )
5575 for (
int wIt = 0; wIt < minW; wIt++ )
5577 arrPos = wIt + minW * ( vIt + minV * uIt );
5578 if ( resIn[arrPos][0] > mapPeak )
5580 mapPeak = resIn[arrPos][0];
5588 ret[3] = mapPeak / ( ( obj2->_maxMapU+1 ) * ( obj2->_maxMapV+1 ) * ( obj2->_maxMapW+1 ) );
5591 if ( ret[0] > muo2 ) { ret[0] = ret[0] - minU; }
5592 if ( ret[1] > mvo2 ) { ret[1] = ret[1] - minV; }
5593 if ( ret[2] > mwo2 ) { ret[2] = ret[2] - minW; }
5596 fftw_destroy_plan ( forwardFourierObj1 );
5597 fftw_destroy_plan ( forwardFourierObj2 );
5598 fftw_destroy_plan ( inverseFourierCombo );
5607 ret[0] *= obj2->_xSamplingRate;
5608 ret[1] *= obj2->_ySamplingRate;
5609 ret[2] *= obj2->_zSamplingRate;
5612 if ( ret[0] > ( obj2->_xRange / 2.0 ) ) { ret[0] -= obj2->_xRange; }
5613 if ( ret[1] > ( obj2->_yRange / 2.0 ) ) { ret[1] -= obj2->_yRange; }
5614 if ( ret[2] > ( obj2->_zRange / 2.0 ) ) { ret[2] -= obj2->_zRange; }
5616 if ( ret[0] < ( -obj2->_xRange / 2.0 ) ) { ret[0] += obj2->_xRange; }
5617 if ( ret[1] < ( -obj2->_yRange / 2.0 ) ) { ret[1] += obj2->_yRange; }
5618 if ( ret[2] < ( -obj2->_zRange / 2.0 ) ) { ret[2] += obj2->_zRange; }
5643 std::string saveName,
5649 if ( ( this->_eulerAngles[0] == 0.0 ) && ( this->_eulerAngles[1] == 0.0 ) && ( this->_eulerAngles[2] == 0.0 ) && ( !internalUse ) )
5651 cmpObj1->writeMap ( saveName, cmpObj1->_densityMapCor, axOrd );
5655 if ( ( this->_eulerAngles[0] == 0.0 ) && ( this->_eulerAngles[1] == 0.0 ) && ( this->_eulerAngles[2] == 0.0 ) && ( internalUse ) )
5661 std::array<double,2> hlpArr;
5664 double **realSHCoeffsRes =
new double* [cmpObj1->_noShellsWithData];
5665 double **imagSHCoeffsRes =
new double* [cmpObj1->_noShellsWithData];
5666 for (
unsigned int i = 0; i < cmpObj1->_noShellsWithData; i++ )
5668 realSHCoeffsRes[i] =
new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5669 imagSHCoeffsRes[i] =
new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5671 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( cmpObj1->_oneDimmension * cmpObj1->_oneDimmension ); iter++ )
5673 realSHCoeffsRes[i][iter] = 0.0;
5674 imagSHCoeffsRes[i][iter] = 0.0;
5683 std::cout <<
">>>>>>>> Wigner matrices for given Euler angles obtained." << std::endl;
5688 for (
int shell = 0; shell < static_cast<int> ( cmpObj1->_noShellsWithData ); shell++ )
5691 for (
int bandIter = 0; bandIter < static_cast<int> ( this->_bandwidthLimit ); bandIter++ )
5694 for (
int ord1 = 0; ord1 < (bandIter * 2 + 1); ord1++ )
5696 arrPos2 = seanindex ( ord1-bandIter, bandIter, this->_bandwidthLimit );
5699 for (
int ord2 = 0; ord2 < (bandIter * 2 + 1); ord2++ )
5701 arrPos = seanindex ( ord2-bandIter, bandIter, this->_bandwidthLimit );
5702 hlpArr = ProSHADE_internal_misc::complexMultiplication ( &cmpObj1->_realSHCoeffs[shell][arrPos],
5703 &cmpObj1->_imagSHCoeffs[shell][arrPos],
5704 &this->_wignerMatrices.at(bandIter).at(ord1).at(ord2)[0],
5705 &this->_wignerMatrices.at(bandIter).at(ord1).at(ord2)[1] );
5707 realSHCoeffsRes[shell][arrPos2] += hlpArr[0];
5708 imagSHCoeffsRes[shell][arrPos2] += hlpArr[1];
5715 std::cout <<
"Rotated all coefficients." << std::endl;
5720 fftw_plan idctPlan, ifftPlan;
5721 int rank, howmany_rank ;
5722 fftw_iodim dims[1], howmany_dims[1];
5724 double *sigR =
new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5725 double *sigI =
new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5726 double *rcoeffs =
new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5727 double *icoeffs =
new double [cmpObj1->_oneDimmension * cmpObj1->_oneDimmension];
5728 double *weights =
new double [4 * this->_bandwidthLimit];
5729 double *workspace =
new double [( 10 * this->_bandwidthLimit * this->_bandwidthLimit ) + ( 24 * this->_bandwidthLimit )];
5732 double **newShellMappedData =
new double*[cmpObj1->_noShellsWithData];
5733 for (
unsigned int sh = 0; sh < cmpObj1->_noShellsWithData; sh++ )
5735 newShellMappedData[sh] =
new double[
static_cast<unsigned int> ( 4 * this->_bandwidthLimit * this->_bandwidthLimit )];
5739 std::cout <<
">>>>> Memory allocation complete." << std::endl;
5743 idctPlan = fftw_plan_r2r_1d ( 2 * this->_bandwidthLimit,
5750 dims[0].n = 2 * this->_bandwidthLimit;
5751 dims[0].is = 2 * this->_bandwidthLimit;
5754 howmany_dims[0].n = 2 * this->_bandwidthLimit;
5755 howmany_dims[0].is = 1;
5756 howmany_dims[0].os = 2 * this->_bandwidthLimit;
5759 ifftPlan = fftw_plan_guru_split_dft ( rank,
5770 for (
int shell = 0; shell < static_cast<int> ( cmpObj1->_noShellsWithData ); shell++ )
5773 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( cmpObj1->_oneDimmension * cmpObj1->_oneDimmension ); iter++ )
5775 rcoeffs[iter] = realSHCoeffsRes[shell][iter];
5776 icoeffs[iter] = imagSHCoeffsRes[shell][iter];
5780 InvFST_semi_fly ( rcoeffs,
5784 this->_bandwidthLimit,
5787 this->_bandwidthLimit,
5792 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( cmpObj1->_oneDimmension * cmpObj1->_oneDimmension ); iter++ )
5794 newShellMappedData[shell][iter] = sigR[iter];
5799 std::cout <<
">> Inverted all rotated spherical harmonics coefficients." << std::endl;
5811 for (
unsigned int i = 0; i < cmpObj1->_noShellsWithData; i++ )
5813 delete[] realSHCoeffsRes[i];
5814 delete[] imagSHCoeffsRes[i];
5816 delete[] realSHCoeffsRes;
5817 delete[] imagSHCoeffsRes;
5821 double *densityMapRotated =
new double [(cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1)];
5823 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ((cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1)); iter++ )
5825 densityMapRotated[iter] = 0.0;
5829 std::vector<double> lonCO ( cmpObj1->_thetaAngle + 1 );
5830 for (
unsigned int iter = 0; iter <= cmpObj1->_thetaAngle; iter++ ) { lonCO.at(iter) =
static_cast<double> ( iter ) * ( ( static_cast<double> ( M_PI ) * 2.0 ) /
static_cast<double> ( this->_thetaAngle ) ) - (
static_cast<double> ( M_PI ) ); }
5831 std::vector<double> latCO ( cmpObj1->_phiAngle + 1 );
5832 for (
unsigned int iter = 0; iter <= cmpObj1->_phiAngle; iter++ ) { latCO.at(iter) = (
static_cast<double> ( iter ) * ( static_cast<double> ( M_PI ) /
static_cast<double> ( this->_phiAngle ) ) - (
static_cast<double> ( M_PI ) / 2.0 ) ); }
5841 unsigned int lowerLon = 0;
5842 unsigned int upperLon = 0;
5843 unsigned int lowerLat = 0;
5844 unsigned int upperLat = 0;
5845 unsigned int lowerShell = 0;
5846 unsigned int upperShell = 0;
5851 double dist00 = 0.0;
5852 double dist01 = 0.0;
5853 double dist10 = 0.0;
5854 double dist11 = 0.0;
5855 double wCoeff = 0.0;
5859 for (
int uIt = 0; uIt < static_cast<int> (cmpObj1->_maxMapU+1); uIt++ )
5861 for (
int vIt = 0; vIt < static_cast<int> (cmpObj1->_maxMapV+1); vIt++ )
5863 for (
int wIt = 0; wIt < static_cast<int> (cmpObj1->_maxMapW+1); wIt++ )
5867 newU =
static_cast<double> ( uIt - (
static_cast<int> (cmpObj1->_maxMapU+1) / 2 ) );
5868 newV =
static_cast<double> ( vIt - (
static_cast<int> (cmpObj1->_maxMapV+1) / 2 ) );
5869 newW =
static_cast<double> ( wIt - (
static_cast<int> (cmpObj1->_maxMapW+1) / 2 ) );
5872 if ( ( newU == 0.0 ) && ( newV == 0.0 ) && ( newW == 0.0 ) )
5874 arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
5875 densityMapRotated[arrPos] = cmpObj1->_densityMapCor[arrPos];
5880 rad = sqrt ( pow( (newU*cmpObj1->_xSamplingRate), 2.0 ) +
5881 pow( (newV*cmpObj1->_ySamplingRate), 2.0 ) +
5882 pow( (newW*cmpObj1->_zSamplingRate), 2.0 ) );
5883 lon = atan2 ( (newV*cmpObj1->_ySamplingRate) , (newU*cmpObj1->_xSamplingRate) );
5884 lat = asin ( (newW*cmpObj1->_zSamplingRate) / rad );
5887 if ( rad != rad ) { rad = 0.0; }
5888 if ( lon != lon ) { lon = 0.0; }
5889 if ( lat != lat ) { lat = 0.0; }
5892 for (
unsigned int iter = 0; iter < cmpObj1->_thetaAngle; iter++ )
5894 if ( ( std::trunc(10000. * lonCO.at(iter)) <= std::trunc(10000. * lon) ) && ( std::trunc(10000. * lonCO.at(iter+1)) > std::trunc(10000. * lon) ) )
5901 if ( upperLon == cmpObj1->_thetaAngle ) { upperLon = 0; }
5903 for (
unsigned int iter = 0; iter < cmpObj1->_phiAngle; iter++ )
5905 if ( ( std::trunc(10000. * latCO.at(iter)) <= std::trunc(10000. * lat) ) && ( std::trunc(10000. * latCO.at(iter+1)) > std::trunc(10000. * lat) ) )
5912 if ( upperLat == cmpObj1->_phiAngle ) { upperLat = 0; }
5917 for (
unsigned int iter = 0; iter < (cmpObj1->_noShellsWithData-1); iter++ )
5919 if ( ( cmpObj1->_shellPlacement.at(iter) <= rad ) && ( cmpObj1->_shellPlacement.at(iter+1) > rad ) )
5922 upperShell = iter+1;
5926 if ( upperShell == 0 )
5928 arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
5929 densityMapRotated[arrPos] = 0.0;
5934 x00 = newShellMappedData[lowerShell][
static_cast<int> ( lowerLat * cmpObj1->_thetaAngle + lowerLon )];
5935 dist00 = sqrt( pow( (lon - lonCO.at(lowerLon)), 2.0 ) +
5936 pow( (lat - latCO.at(lowerLat)), 2.0 ) );
5937 x01 = newShellMappedData[lowerShell][
static_cast<int> ( lowerLat * cmpObj1->_thetaAngle + upperLon )];
5938 dist01 = sqrt( pow( (lon - lonCO.at(upperLon)), 2.0 ) +
5939 pow( (lat - latCO.at(lowerLat)), 2.0 ) );
5940 x10 = newShellMappedData[lowerShell][
static_cast<int> ( upperLat * cmpObj1->_thetaAngle + lowerLon )];
5941 dist10 = sqrt( pow( (lon - lonCO.at(lowerLon)), 2.0 ) +
5942 pow( (lat - latCO.at(upperLat)), 2.0 ) );
5943 x11 = newShellMappedData[lowerShell][
static_cast<int> ( upperLat * cmpObj1->_thetaAngle + upperLon )];
5944 dist11 = sqrt( pow( (lon - lonCO.at(upperLon)), 2.0 ) +
5945 pow( (lat - latCO.at(upperLat)), 2.0 ) );
5947 if ( !( dist00 == 0.0 ) || !( dist01 == 0.0 ) || !( dist10 == 0.0 ) || !( dist11 == 0.0 ) )
5949 if ( ( dist00 == 0.0 ) && ( dist01 == 0.0 ) && ( dist10 == 0.0 ) && ( dist11 == 0.0 ) )
5951 arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
5952 densityMapRotated[arrPos] = 0.0;
5957 if ( dist00 == 0.0 ) { dist00 = std::max ( dist01, std::max ( dist10, dist11 ) ) * 0.0000001; }
5958 if ( dist01 == 0.0 ) { dist01 = std::max ( dist00, std::max ( dist10, dist11 ) ) * 0.0000001; }
5959 if ( dist10 == 0.0 ) { dist10 = std::max ( dist01, std::max ( dist00, dist11 ) ) * 0.0000001; }
5960 if ( dist11 == 0.0 ) { dist11 = std::max ( dist01, std::max ( dist10, dist00 ) ) * 0.0000001; }
5964 wCoeff = 1.0/( 1.0/dist00 + 1.0/dist01 + 1.0/dist10 + 1.0/dist11 );
5965 lsVal = (x00 * (( 1.0/dist00 )*wCoeff)) + (x01 * ((1.0/dist01)*wCoeff)) + (x10 * ((1.0/dist10)*wCoeff)) + (x11 * ((1.0/dist11)*wCoeff));
5968 x00 = newShellMappedData[upperShell][
static_cast<int> ( lowerLat * cmpObj1->_thetaAngle + lowerLon )];
5969 x01 = newShellMappedData[upperShell][
static_cast<int> ( lowerLat * cmpObj1->_thetaAngle + upperLon )];
5970 x10 = newShellMappedData[upperShell][
static_cast<int> ( upperLat * cmpObj1->_thetaAngle + lowerLon )];
5971 x11 = newShellMappedData[upperShell][
static_cast<int> ( upperLat * cmpObj1->_thetaAngle + upperLon )];
5973 upVal = (x00 * (( 1.0/dist00 )*wCoeff)) + (x01 * ((1.0/dist01)*wCoeff)) + (x10 * ((1.0/dist10)*wCoeff)) + (x11 * ((1.0/dist11)*wCoeff));
5976 dist00 = cmpObj1->_shellPlacement.at(upperShell) - rad;
5977 dist11 = rad - cmpObj1->_shellPlacement.at(lowerShell);
5979 arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
5980 densityMapRotated[arrPos] = upVal*dist11 + lsVal*dist00;
5988 std::vector<double> rMapVals ( (cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1) );
5989 std::vector<double> oMapVals ( (cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1) );
5995 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( (cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1) ); iter++ )
6000 if ( densityMapRotated[iter] == densityMapRotated[iter] ) { rMapVals.at(iter) = densityMapRotated[iter]; }
6001 else { rMapVals.at(iter) = 0.0; }
6004 oMapVals.at(iter) = cmpObj1->_densityMapCor[iter];
6007 double rMapMean = std::accumulate ( rMapVals.begin(), rMapVals.end(), 0.0 ) / static_cast<double> ( rMapVals.size() );
6008 double rMapSdev = std::sqrt ( static_cast<double> ( std::inner_product ( rMapVals.begin(), rMapVals.end(), rMapVals.begin(), 0.0 ) ) /
static_cast<double> ( rMapVals.size() ) - rMapMean * rMapMean );
6010 double oMapMean = std::accumulate ( oMapVals.begin(), oMapVals.end(), 0.0 ) / static_cast<double> ( oMapVals.size() );
6011 double oMapSdev = std::sqrt ( static_cast<double> ( std::inner_product ( oMapVals.begin(), oMapVals.end(), oMapVals.begin(), 0.0 ) ) /
static_cast<double> ( oMapVals.size() ) - oMapMean * oMapMean );
6014 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( (cmpObj1->_maxMapU+1)*(cmpObj1->_maxMapV+1)*(cmpObj1->_maxMapW+1) ); iter++ )
6016 densityMapRotated[iter] = ( rMapVals.at(iter) - rMapMean ) / rMapSdev;
6017 cmpObj1->_densityMapCor[iter] = ( cmpObj1->_densityMapCor[iter] - oMapMean ) / oMapSdev;
6022 double dist = std::max ( (cmpObj1->_maxMapU+1)/2.0, std::max ( (cmpObj1->_maxMapV+1)/2.0, (cmpObj1->_maxMapW+1)/2.0 ) );
6024 for (
int uIt = 0; uIt < static_cast<int> (cmpObj1->_maxMapU+1); uIt++ )
6026 for (
int vIt = 0; vIt < static_cast<int> (cmpObj1->_maxMapV+1); vIt++ )
6028 for (
int wIt = 0; wIt < static_cast<int> (cmpObj1->_maxMapW+1); wIt++ )
6030 arrPos = wIt + (cmpObj1->_maxMapW + 1) * ( vIt + (cmpObj1->_maxMapV + 1) * uIt );
6032 if ( sqrt ( pow( uIt-((cmpObj1->_maxMapU+1)/2.0), 2.0 ) +
6033 pow( vIt-((cmpObj1->_maxMapV+1)/2.0), 2.0 ) +
6034 pow( wIt-((cmpObj1->_maxMapW+1)/2.0), 2.0 ) ) > dist )
6036 densityMapRotated[arrPos] = 0.0;
6037 cmpObj1->_densityMapCor[arrPos] = 0.0;
6044 for (
unsigned int sh = 0; sh < cmpObj1->_noShellsWithData; sh++ )
6046 delete[] newShellMappedData[sh];
6048 delete[] newShellMappedData;
6053 cmpObj1->writeMap ( saveName, densityMapRotated, axOrd );
6055 delete[] cmpObj1->_densityMapCor;
6056 cmpObj1->_densityMapCor =
nullptr;
6061 for (
unsigned int iter = 0; iter < static_cast<unsigned int> ( (cmpObj1->_maxMapU+1) * (cmpObj1->_maxMapV+1) * (cmpObj1->_maxMapW+1) ); iter++ )
6063 cmpObj1->_densityMapCor[iter] = densityMapRotated[iter];
6068 delete[] densityMapRotated;
6088 if ( ( this->all->size() < 1 ) || ( ( this->all->size() % 2 ) != 0 ) )
6090 std::cerr <<
"!!! ProSHADE ERROR !!! Error in file comparison !!! The compare against part has incorrect number of entries." << std::endl;
6094 int checkValue =
static_cast<double> ( this->one->_keepOrRemove );
6095 for (
unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt += 2 )
6097 if ( checkValue != static_cast<double> ( this->all->at(strIt)->_keepOrRemove ) )
6099 std::cerr <<
"!!! ProSHADE ERROR !!! Error in file comparison between " << this->one->_inputFileName <<
" AND " << this->all->at(strIt)->_inputFileName <<
" !!! The phase treatment is different, use the same phase treatment (i.e. remove all, or keep all) to make the strucutres comparable." << std::endl;
6104 checkValue =
static_cast<double> ( this->two->_keepOrRemove );
6105 for (
unsigned int strIt = 1; strIt < static_cast<unsigned int> ( this->all->size() ); strIt += 2 )
6107 if ( checkValue != static_cast<double> ( this->all->at(strIt)->_keepOrRemove ) )
6109 std::cerr <<
"!!! ProSHADE ERROR !!! Error in file comparison between " << this->one->_inputFileName <<
" AND " << this->all->at(strIt)->_inputFileName <<
" !!! The phase treatment is different, use the same phase treatment (i.e. remove all, or keep all) to make the strucutres comparable." << std::endl;
6115 std::vector<ProSHADE_internal::ProSHADE_data*> newAll;
6116 std::vector<double> pattCorrelationVec;
6117 std::vector<std::array<double,4>> translationVec;
6119 bool reverseOrder =
false;
6122 for (
unsigned int strIt = 0; strIt < static_cast<unsigned int> ( this->all->size() ); strIt += 2 )
6124 if ( this->one->_keepOrRemove ==
false )
6127 if ( ( strIt == 0 ) && ( this->all->at(strIt)->_keepOrRemove != false ) )
6130 reverseOrder =
true;
6135 this->all->at(strIt),
6142 pattCorrelationVec.emplace_back ( 0.0 );
6148 std::array<double,3> euAngs = cmpObj->
getEulerAngles ( settings, &pattCorrelationVec.at(totStrIt) );
6157 this->all->at(strIt-1),
6166 this->all->at(strIt+1),
6175 double mat22 = cos ( euAngs[1] );
6176 double mat02 = -sin ( euAngs[1] ) * cos ( euAngs[2] );
6177 double mat12 = sin ( euAngs[1] ) * sin ( euAngs[2] );
6178 double mat20 = cos ( euAngs[0] ) * sin ( euAngs[1] );
6179 double mat21 = sin ( euAngs[0] ) * sin ( euAngs[1] );
6182 double rMat02 = mat20;
6183 double rMat20 = mat02;
6184 double rMat21 = mat12;
6185 double rMat12 = mat21;
6188 euAngs[0] = atan2 ( rMat21, rMat20 );
6189 euAngs[1] = acos ( mat22 );
6190 euAngs[2] = atan2 ( rMat12, -rMat02 );
6192 if ( euAngs[0] < 0.0 ) { euAngs[0]= 2.0 * M_PI + euAngs[0]; }
6193 if ( euAngs[1] < 0.0 ) { euAngs[1]= M_PI + euAngs[1]; }
6194 if ( euAngs[2] < 0.0 ) { euAngs[2]= 2.0 * M_PI + euAngs[2]; }
6209 translationVec.emplace_back ( cmpObj->
getTranslationFunctionMap ( &str1Copy, this->all->at(strIt-1), &xMov1, &yMov1, &zMov1 ) );
6211 this->all->at(strIt-1)->translateMap ( -translationVec.at(totStrIt)[0],
6212 -translationVec.at(totStrIt)[1],
6213 -translationVec.at(totStrIt)[2] );
6225 translationVec.emplace_back ( cmpObj->
getTranslationFunctionMap ( &str1Copy, this->all->at(strIt+1), &xMov1, &yMov1, &zMov1 ) );
6227 this->all->at(strIt+1)->translateMap ( -translationVec.at(totStrIt)[0],
6228 -translationVec.at(totStrIt)[1],
6229 -translationVec.at(totStrIt)[2] );
6235 if ( ( strIt == 0 ) && ( this->all->at(strIt)->_keepOrRemove != false ) )
6238 reverseOrder =
true;
6243 this->all->at(strIt),
6250 pattCorrelationVec.emplace_back ( 0.0 );
6256 std::array<double,3> euAngs = cmpObj->
getEulerAngles ( settings, &pattCorrelationVec.at(totStrIt) );
6265 this->all->at(strIt-1),
6274 this->all->at(strIt+1),
6283 double mat22 = cos ( euAngs[1] );
6284 double mat02 = -sin ( euAngs[1] ) * cos ( euAngs[2] );
6285 double mat12 = sin ( euAngs[1] ) * sin ( euAngs[2] );
6286 double mat20 = cos ( euAngs[0] ) * sin ( euAngs[1] );
6287 double mat21 = sin ( euAngs[0] ) * sin ( euAngs[1] );
6290 double rMat02 = mat20;
6291 double rMat20 = mat02;
6292 double rMat21 = mat12;
6293 double rMat12 = mat21;
6296 euAngs[0] = atan2 ( rMat21, rMat20 );
6297 euAngs[1] = acos ( mat22 );
6298 euAngs[2] = atan2 ( rMat12, -rMat02 );
6300 if ( euAngs[0] < 0.0 ) { euAngs[0]= 2.0 * M_PI + euAngs[0]; }
6301 if ( euAngs[1] < 0.0 ) { euAngs[1]= M_PI + euAngs[1]; }
6302 if ( euAngs[2] < 0.0 ) { euAngs[2]= 2.0 * M_PI + euAngs[2]; }
6317 translationVec.emplace_back ( cmpObj->
getTranslationFunctionMap ( &str1Copy, this->all->at(strIt-1), &xMov1, &yMov1, &zMov1 ) );
6319 this->all->at(strIt-1)->translateMap ( -translationVec.at(totStrIt)[0],
6320 -translationVec.at(totStrIt)[1],
6321 -translationVec.at(totStrIt)[2] );
6333 translationVec.emplace_back ( cmpObj->
getTranslationFunctionMap ( &str1Copy, this->all->at(strIt+1), &xMov1, &yMov1, &zMov1 ) );
6335 this->all->at(strIt+1)->translateMap ( -translationVec.at(totStrIt)[0],
6336 -translationVec.at(totStrIt)[1],
6337 -translationVec.at(totStrIt)[2] );
6345 for (
signed int strIt = static_cast<signed int> ( this->all->size() - 1 ); strIt >= 0 ; strIt-- )
6347 if ( this->all->at(strIt)->_keepOrRemove == false )
6349 delete this->all->at(strIt);
6350 this->all->at(strIt) =
nullptr;
6351 this->all->erase ( this->all->begin ( ) + strIt );
This class deals with reading in the data and computing structure specific information including the ...
void generateWignerMatrices(ProSHADE::ProSHADE_settings *settings)
This function is responsible for computing the Wigner D matrices for symmetry detection purposes...
std::vector< std::array< double, 8 > > getSO3Peaks(ProSHADE::ProSHADE_settings *settings, double noIQRs=1.5, bool freeMem=true, int peakSize=1, double realDist=0.4, int verbose=1)
This function detects 'true' peaks from the background of the SO3 inverse transform map...
std::vector< std::vector< std::array< double, 5 > > > findTetrSymmetry(std::vector< std::array< double, 5 > > CnSymm, double *tetrSymmPeakAvg, double axisErrorTolerance=0.1)
This function detects the tetrahedral symmetries in the structure.
void rotateStructure(ProSHADE_data *cmpObj1, ProSHADE::ProSHADE_settings *settings, std::string saveName, int verbose=0, std::string axOrd="xyz", bool internalUse=false)
Rotates the density map be given Angle-Axis rotation using the Wigner matrices and inverse spharical ...
std::string clearMapFile
If map features are to be extracted, should the clear map be saved (then give file name here)...
std::vector< std::vector< std::array< double, 5 > > > findIcosSymmetry(std::vector< std::array< double, 5 > > CnSymm, double *icosSymmPeakAvg, double axisErrorTolerance=0.1)
This function detects the icosahedral symmetries in the structure.
void alignDensities(ProSHADE::ProSHADE_settings *settings)
Takes the internal objects with and without phases and aligns them to all the other objects...
bool htmlReport
Should HTML report for the run be created?
std::vector< std::array< double, 5 > > findCnSymmetry(std::vector< std::array< double, 8 > > peaks, ProSHADE::ProSHADE_settings *settings, double axisErrorTolerance=0.0, bool freeMem=true, double percAllowedToMiss=0.33, int verbose=1)
This function attempts to detect the C-symmetries present in the list of overlay peaks.
void setEulerAngles(double alpha, double beta, double gamma)
This function is used to set the Euler angles for further processing.
int verbose
Should the software report on the progress, or just be quiet? Value between 0 (quiet) and 4 (loud) ...
std::vector< std::array< double, 5 > > generateTetrElements(std::vector< std::array< double, 5 > > symmAxes, ProSHADE::ProSHADE_settings *settings, int verbose=0)
This function generates the 12 unique tetrahedral symmetry group elements from the already detected a...
The main header file containing all declarations the user of the library needs.
double getRotCoeffDistance(ProSHADE::ProSHADE_settings *settings)
This function computes the full rotation function descriptor distances.
void precomputeTrSigmaDescriptor()
This function computes the E matrices required for the trace sigma descriptor, the rotation function ...
std::vector< std::array< double, 5 > > generateTetrAxes(ProSHADE_comparePairwise *cmpObj, ProSHADE::ProSHADE_settings *settings, std::vector< std::array< double, 5 > > tetrSymm, std::vector< std::array< double, 5 > > allCs, double axisErrorTolerance=0.1, int verbose=0)
This function generates all the tetrahedral symmetry axes from a pair detected by findTetrSymmetry fu...
std::vector< std::vector< std::array< double, 5 > > > findOctaSymmetry(std::vector< std::array< double, 5 > > CnSymm, double *octaSymmPeakAvg, double axisErrorTolerance=0.1)
This function detects the octahedral symmetries in the structure.
void getSO3InverseMap(ProSHADE::ProSHADE_settings *settings)
This function is responsible for computing the SO3 inverse transform.
double maxAvgPeakForSymmetry(double X, double Y, double Z, double Angle, ProSHADE::ProSHADE_settings *settings)
This function takes angle and axis and checks the SO(3) map for this specific symmetry only...
double getPeakThreshold()
This is an accessor function to the peak height threshold.
std::vector< int > ignoreLs
This vector lists all the bandwidth values which should be ignored and not part of the computations...
std::vector< std::array< double, 5 > > generateOctaElements(std::vector< std::array< double, 5 > > symmAxes, ProSHADE::ProSHADE_settings *settings, int verbose=0)
This function generates the 24 octahedral symmetry group elements from the axes generated by generate...
std::vector< std::vector< std::array< double, 6 > > > findDnSymmetry(std::vector< std::array< double, 5 > > CnSymm, double axisErrorTolerance=0.1)
This function detects dihedral (D) symmetries from the list of already detected C symmetries...
std::vector< std::array< double, 5 > > generateOctaAxes(ProSHADE_comparePairwise *cmpObj, ProSHADE::ProSHADE_settings *settings, std::vector< std::array< double, 5 > > octaSymm, std::vector< std::array< double, 5 > > allCs, double axisErrorTolerance=0.1, int verbose=0)
This function generates all the ctahedral symmetry elements from a pair detected by findOctaSymmetry ...
This is the executive class responsible for comparing strictly two structures.
std::vector< std::array< double, 5 > > findCnSymmetryClear(std::vector< std::array< double, 5 > > CnSymm, ProSHADE::ProSHADE_settings *settings, double maxGap=0.2, bool *pf=nullptr)
This function takes the detected C-symmetries list, removes low probability symmetries and sorts as t...
std::vector< std::array< double, 5 > > generateIcosAxes(ProSHADE_comparePairwise *cmpObj, ProSHADE::ProSHADE_settings *settings, std::vector< std::array< double, 5 > > icosSymm, std::vector< std::array< double, 5 > > allCs, double axisErrorTolerance=0.1, int verbose=0)
This function generates all the icosahedral symmetry elements from a pair detected by findIcosSymmetr...
void getSO3InverseMap(ProSHADE::ProSHADE_settings *settings)
This function is responsible for computing the SO3 inverse transform.
std::vector< std::array< double, 5 > > generateIcosElements(std::vector< std::array< double, 5 > > symmAxes, ProSHADE::ProSHADE_settings *settings, int verbose=0)
This function generates the 60 icosahedral symmetry group elements from the axes generated by generat...
std::string axisOrder
A string specifying the order of the axis. Must have three characters and any permutation of 'x'...
The main header file containing all declarations for the innter workings of the library.
bool checkPeakExistence(double alpha, double beta, double gamma, ProSHADE::ProSHADE_settings *settings, double passThreshold=0.3)
This function checks if a peak Euler angles actually produce something the descriptors would agree is...
int htmlReportLineProgress
Iterator for current HTML line in the progress bar.
ProSHADE_comparePairwise(ProSHADE_data *cmpObj1, ProSHADE_data *cmpObj2, double mPower, std::vector< int > ignoreL, unsigned int order, ProSHADE::ProSHADE_settings *settings)
Contructor for the ProSHADE_comparePairwise class.
std::vector< std::vector< std::array< double, 6 > > > findDnSymmetryClear(std::vector< std::vector< std::array< double, 6 > > > DnSymm, ProSHADE::ProSHADE_settings *settings, double maxGap=0.2, bool *pf=nullptr)
This function sorts the D symmetries and removed low probability ones.
std::vector< std::array< double, 3 > > getEulerAngles(ProSHADE::ProSHADE_settings *settings)
This function finds the highest peak in the SO3 inverse transform map and sets it as the optimal over...
This class stores all the settings and is passed to the executive classes instead of multitude of par...
void freeInvMap()
This function frees the SOFT inverse map.
double mPower
This parameter determines the scaling for trace sigma descriptor.
std::array< double, 4 > getTranslationFunctionMap(ProSHADE_data *obj1, ProSHADE_data *obj2, double *ob2XMov=nullptr, double *ob2YMov=nullptr, double *ob2ZMov=nullptr)
Computes the translation function for the two objects and returns the position of the highest peak...
This file contains the ProSHADE_internal_misc namespace and its miscellaneous functions.
std::array< double, 3 > getEulerAngles(ProSHADE::ProSHADE_settings *settings, double *correlation=nullptr)
This function finds the highest peak in the SO3 inverse transform map and sets it as the optimal over...
unsigned int glIntegOrder
This parameter controls the Gauss-Legendre integration order and so the radial resolution.