41 std::vector< proshade_double* > ret;
42 std::vector< proshade_double > nonPeakVals;
43 proshade_double* retHlp = NULL;
44 proshade_double pointHeight = 0.0;
45 proshade_signed x, y, z, peakX, peakY, peakZ, newIter, ptrIter;
46 proshade_signed xDim = pow ( dim, 2 );
47 proshade_signed yDim = dim;
51 for ( proshade_unsign iter = 0; iter < ( pow( static_cast<proshade_unsign> ( dim ), 3 ) ); iter++ )
54 pointHeight = pow( map[iter][0], 2.0 ) + pow( map[iter][1], 2.0 );
57 x = std::floor ( iter / xDim );
58 y = std::floor ( ( iter - x * xDim ) / yDim );
59 z = iter - x * xDim - y * yDim;
64 retHlp =
new proshade_double[
static_cast<proshade_unsign
>( pow( ((peakSize*2)+1), 3) * 4 )];
71 for ( proshade_signed xCh = -peakSize; xCh <= +peakSize; xCh++ )
73 if ( breakPeak ) {
break; }
74 for ( proshade_signed yCh = -peakSize; yCh <= +peakSize; yCh++ )
76 if ( breakPeak ) {
break; }
77 for ( proshade_signed zCh = -peakSize; zCh <= +peakSize; zCh++ )
79 if ( breakPeak ) {
break; }
80 if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) {
continue; }
83 peakX = x + xCh;
if ( peakX >= yDim ) { peakX = yDim - 1; };
if ( peakX < 0 ) { peakX = 0; }
84 peakY = y + yCh;
if ( peakY >= yDim ) { peakY = yDim - 1; };
if ( peakY < 0 ) { peakY = 0; }
85 peakZ = z + zCh;
if ( peakZ >= yDim ) { peakZ = yDim - 1; };
if ( peakZ < 0 ) { peakZ = 0; }
86 newIter = peakX * xDim + peakY * yDim + peakZ;
89 if ( ( pow ( map[newIter][0], 2.0 ) + pow( map[newIter][1], 2.0 ) ) > pointHeight ) { breakPeak =
true;
break; }
92 retHlp[ptrIter] =
static_cast<proshade_double
> ( peakX );
93 retHlp[ptrIter+1] =
static_cast<proshade_double
> ( peakY );
94 retHlp[ptrIter+2] =
static_cast<proshade_double
> ( peakZ );
95 retHlp[ptrIter+3] = pow ( map[newIter][0], 2.0 ) + pow( map[newIter][1], 2.0 );
106 retHlp[3] = pointHeight;
115 if ( retHlp != NULL ) {
delete[] retHlp; }
135 proshade_double backgroundThreshold = std::min ( std::max ( medianIQR[0] + ( medianIQR[1] * noIQRs ), 0.05 ), 0.3 );
138 for ( proshade_signed iter =
static_cast<proshade_signed
> ( pointVec->size()-1 ); iter >= 0 ; iter-- )
140 if ( pointVec->at(iter)[3] <= backgroundThreshold )
142 delete[] pointVec->at(iter);
143 pointVec->erase ( pointVec->begin() + iter );
164 avgMat =
new proshade_double [9];
165 hlpMat =
new proshade_double [9];
166 eA =
new proshade_double;
167 eB =
new proshade_double;
168 eG =
new proshade_double;
169 uAV =
new proshade_double [18];
224 proshade_double *avgMat, *hlpMat, *eulA, *eulB, *eulG, *uAndV;
225 proshade_unsign noPoints = pow( ( ( peakSize * 2 ) + 1 ), 3 ) * 4;
226 proshade_double matWeight;
232 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( pointVec->size() ); iter++ )
235 for ( proshade_unsign i = 0; i < 9; i++ ) { avgMat[i] = 0.0; }
239 for ( proshade_unsign it = 0; it < noPoints; it += 4 )
248 for ( proshade_unsign i = 0; i < 9; i++ ) { avgMat[i] += hlpMat[i] * pointVec->at(iter)[it+3]; }
251 matWeight += pointVec->at(iter)[it+3];
255 for ( proshade_unsign i = 0; i < 9; i++ ) { avgMat[i] /= matWeight; }
259 if ( uAndV[0] == -777.7 )
263 matWeight = pointVec->at(iter)[3];
264 pointVec->at(iter) =
new proshade_double [4];
266 pointVec->at(iter)[0] = *eulA;
267 pointVec->at(iter)[1] = *eulB;
268 pointVec->at(iter)[2] = *eulG;
269 pointVec->at(iter)[3] = matWeight;
275 for ( proshade_unsign i = 0; i < 9; i++ ) { avgMat[i] = 0.0; }
282 matWeight = pointVec->at(iter)[3];
283 delete[] pointVec->at(iter);
284 pointVec->at(iter) =
new proshade_double [4];
286 pointVec->at(iter)[0] = *eulA;
287 pointVec->at(iter)[1] = *eulB;
288 pointVec->at(iter)[2] = *eulG;
289 pointVec->at(iter)[3] = matWeight;
318 std::vector< proshade_double* > allHigherIndices;
319 proshade_double* nonPeakMedianIQR =
new proshade_double[2];
330 delete[] nonPeakMedianIQR;
333 return ( allHigherIndices );
360 std::stringstream hlpSSP;
361 hlpSSP <<
"Found " << allPeaks.size() <<
" possible peaks.";
365 if ( allPeaks.size() == 0 )
374 proshade_double highestPeak = 0.0;
375 proshade_unsign highestPeakIndex = 0;
376 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign>( allPeaks.size() ); iter++ )
378 if ( allPeaks.at(iter)[3] > highestPeak ) { highestPeak = allPeaks.at(iter)[3]; highestPeakIndex = iter; }
382 *eulA = allPeaks.at(highestPeakIndex)[0];
383 *eulB = allPeaks.at(highestPeakIndex)[1];
384 *eulG = allPeaks.at(highestPeakIndex)[2];
387 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( allPeaks.size() ); iter++ )
389 delete[] allPeaks.at(iter);
393 std::stringstream hlpSS;
394 hlpSS <<
"Optimal Euler angles are " << *eulA <<
" ; " << *eulB <<
" ; " << *eulG <<
" with peak height " << highestPeak;
418 void ProSHADE_internal_peakSearch::allocateSmoothingZScoreMemory ( proshade_unsign dim, proshade_double*& scoreOverVals, proshade_signed*& signals, proshade_double*& filteredY, proshade_double*& avgFilter, proshade_double*& stdFilter, proshade_double*& subVec, proshade_double*& medianIQR, proshade_double*& YZMap, proshade_double*& XZMap, proshade_double*& XYMap, proshade_unsign smLag )
421 signals =
new proshade_signed[dim];
422 scoreOverVals =
new proshade_double[dim];
423 filteredY =
new proshade_double[dim];
424 avgFilter =
new proshade_double[dim];
425 stdFilter =
new proshade_double[dim];
426 subVec =
new proshade_double[smLag];
427 medianIQR =
new proshade_double[2];
428 YZMap =
new proshade_double[dim * dim * dim];
429 XZMap =
new proshade_double[dim * dim * dim];
430 XYMap =
new proshade_double[dim * dim * dim];
462 void ProSHADE_internal_peakSearch::releaseSmoothingZScoreMemory ( proshade_double*& scoreOverVals, proshade_signed*& signals, proshade_double*& filteredY, proshade_double*& avgFilter, proshade_double*& stdFilter, proshade_double*& subVec, proshade_double*& medianIQR, proshade_double*& YZMap, proshade_double*& XZMap, proshade_double*& XYMap )
465 delete[] scoreOverVals;
500 void ProSHADE_internal_peakSearch::getSmoothedZScore1D ( proshade_unsign dim, proshade_unsign smoothingLag, proshade_double ZScoreThreshold, proshade_signed* signals, proshade_double* filteredY, proshade_double* avgFilter, proshade_double* stdFilter, proshade_double* subVec, proshade_double* medianIQR, proshade_double* scoreOverVals )
503 for ( proshade_unsign i = 0; i < dim; i++ ) { signals[i] = 0; avgFilter[i] = 0.0; stdFilter[i] = 0.0; filteredY[i] = 0.0; }
504 for ( proshade_unsign i = 0; i < smoothingLag; i++ ) { subVec[i] = scoreOverVals[i-smoothingLag+dim]; }
506 avgFilter[0] = medianIQR[0];
507 stdFilter[0] = medianIQR[1];
510 for ( proshade_unsign i = 0; i < dim; i++)
513 if ( std::abs ( scoreOverVals[i] - avgFilter[i]) > ZScoreThreshold * stdFilter[i] )
515 if ( scoreOverVals[i] > avgFilter[i] )
525 if ( i != 0 ) { filteredY[i] = 0.5 * scoreOverVals[i] + (1 - 0.5) * filteredY[i - 1]; }
526 else { filteredY[i] = 0.5 * scoreOverVals[i]; }
531 filteredY[i] = scoreOverVals[i];
535 for ( proshade_signed subIt = 0; subIt < static_cast<proshade_signed> ( smoothingLag ); subIt++ )
537 if ( (
static_cast<proshade_signed
>(i) +
static_cast<proshade_signed
>(subIt) -
static_cast<proshade_signed
>(smoothingLag) + 1 ) < 0 )
539 subVec[subIt] = scoreOverVals[(
static_cast<proshade_signed
>(i) +
static_cast<proshade_signed
>(subIt) -
540 static_cast<proshade_signed
>(smoothingLag) + 1 ) + dim];
544 subVec[subIt] = filteredY[(
static_cast<proshade_signed
>(i) +
static_cast<proshade_signed
>(subIt) -
545 static_cast<proshade_signed
>(smoothingLag) + 1 )];
549 avgFilter[i+1] = medianIQR[0];
550 stdFilter[i+1] = medianIQR[1];
577 void ProSHADE_internal_peakSearch::getXAxisArraysSmoothedZScorePeaks ( proshade_unsign dim, proshade_unsign smoothingLag, proshade_double ZScoreThreshold, proshade_signed* signals, proshade_double* filteredY, proshade_double* avgFilter, proshade_double* stdFilter, proshade_double* subVec, proshade_double* medianIQR, proshade_double* scoreOverVals, proshade_complex* map, proshade_double* YZMap )
580 for ( proshade_unsign yIt = 0; yIt < dim; yIt++ )
582 for ( proshade_unsign zIt = 0; zIt < dim; zIt++ )
585 for ( proshade_unsign xIt = 0; xIt < dim; xIt++ )
587 scoreOverVals[xIt] = pow ( map[xIt *
static_cast<proshade_unsign
>(4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt][0], 2.0 ) +
588 pow ( map[xIt *
static_cast<proshade_unsign
>(4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt][1], 2.0 );
592 getSmoothedZScore1D ( dim, smoothingLag, ZScoreThreshold, signals, filteredY, avgFilter, stdFilter, subVec, medianIQR, scoreOverVals );
595 for ( proshade_unsign xIt = 0; xIt < dim; xIt++ )
597 YZMap[xIt *
static_cast<proshade_unsign
> (4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt] = signals[xIt];
626 void ProSHADE_internal_peakSearch::getYAxisArraysSmoothedZScorePeaks ( proshade_unsign dim, proshade_unsign smoothingLag, proshade_double ZScoreThreshold, proshade_signed* signals, proshade_double* filteredY, proshade_double* avgFilter, proshade_double* stdFilter, proshade_double* subVec, proshade_double* medianIQR, proshade_double* scoreOverVals, proshade_complex* map, proshade_double* XZMap )
629 for ( proshade_unsign xIt = 0; xIt < dim; xIt++ )
631 for ( proshade_unsign zIt = 0; zIt < dim; zIt++ )
634 for ( proshade_unsign yIt = 0; yIt < dim; yIt++ )
636 scoreOverVals[yIt] = pow ( map[xIt *
static_cast<proshade_unsign
>(4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt][0], 2.0 ) +
637 pow ( map[xIt *
static_cast<proshade_unsign
>(4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt][1], 2.0 );
641 getSmoothedZScore1D ( dim, smoothingLag, ZScoreThreshold, signals, filteredY, avgFilter, stdFilter, subVec, medianIQR, scoreOverVals );
644 for ( proshade_unsign yIt = 0; yIt < dim; yIt++ )
646 XZMap[xIt *
static_cast<proshade_unsign
> (4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt] = signals[xIt];
675 void ProSHADE_internal_peakSearch::getZAxisArraysSmoothedZScorePeaks ( proshade_unsign dim, proshade_unsign smoothingLag, proshade_double ZScoreThreshold, proshade_signed* signals, proshade_double* filteredY, proshade_double* avgFilter, proshade_double* stdFilter, proshade_double* subVec, proshade_double* medianIQR, proshade_double* scoreOverVals, proshade_complex* map, proshade_double* XYMap )
678 for ( proshade_unsign xIt = 0; xIt < dim; xIt++ )
680 for ( proshade_unsign yIt = 0; yIt < dim; yIt++ )
683 for ( proshade_unsign zIt = 0; zIt < dim; zIt++ )
685 scoreOverVals[yIt] = pow ( map[xIt *
static_cast<proshade_unsign
>(4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt][0], 2.0 ) +
686 pow ( map[xIt *
static_cast<proshade_unsign
>(4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt][1], 2.0 );
690 getSmoothedZScore1D ( dim, smoothingLag, ZScoreThreshold, signals, filteredY, avgFilter, stdFilter, subVec, medianIQR, scoreOverVals );
693 for ( proshade_unsign zIt = 0; zIt < dim; zIt++ )
695 XYMap[xIt *
static_cast<proshade_unsign
> (4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt] = signals[xIt];
720 void ProSHADE_internal_peakSearch::findAllPointNeighbours ( proshade_double* YZMap, proshade_double* XZMap, proshade_double* XYMap, proshade_unsign* visitedMap, proshade_signed dim, proshade_signed x, proshade_signed y, proshade_signed z, std::vector< proshade_unsign >* retVals )
723 proshade_signed newIter = 0;
724 proshade_signed peakX, peakY, peakZ;
725 proshade_signed xDim =
static_cast<proshade_signed
>(4 * pow ( ( dim / 2 ), 2 ));
728 for ( proshade_signed xCh = -1; xCh <= +1; xCh++ )
730 for ( proshade_signed yCh = -1; yCh <= +1; yCh++ )
732 for ( proshade_signed zCh = -1; zCh <= +1; zCh++ )
735 if ( ( xCh == 0 ) && ( yCh == 0 ) && ( zCh == 0 ) ) {
continue; }
738 peakX = x + xCh;
if ( peakX >= dim ) { peakX -= dim; };
if ( peakX < 0 ) { peakX += dim; }
739 peakY = y + yCh;
if ( peakY >= dim ) { peakY -= dim; };
if ( peakY < 0 ) { peakY += dim; }
740 peakZ = z + zCh;
if ( peakZ >= dim ) { peakZ -= dim; };
if ( peakZ < 0 ) { peakZ += dim; }
741 newIter = peakX * xDim + peakY * dim + peakZ;
744 if ( visitedMap[newIter] == 1 ) {
continue; }
745 else { visitedMap[newIter] = 1; }
748 if ( ( YZMap[newIter] + XZMap[newIter] + XYMap[newIter] ) != 3 ) {
continue; }
781 std::vector< proshade_unsign > ret;
782 std::vector< proshade_unsign > thisIsland;
783 proshade_unsign* visitedMap =
new proshade_unsign[dim*dim*dim];
785 for ( proshade_unsign i = 0; i < ( dim * dim * dim ); i++ ) { visitedMap[i] = 0; }
786 proshade_unsign index, maxIndex;
787 proshade_double maxHeight;
790 for ( proshade_unsign xIt = 0; xIt < dim; xIt++ )
792 for ( proshade_unsign yIt = 0; yIt < dim; yIt++ )
794 for ( proshade_unsign zIt = 0; zIt < dim; zIt++ )
797 index = xIt *
static_cast<proshade_unsign
>(4 * pow ( ( dim / 2 ), 2 )) + yIt * dim + zIt;
800 if ( visitedMap[index] == 1 ) {
continue; }
801 else { visitedMap[index] = 1; }
804 if ( ( YZMap[index] + XZMap[index] + XYMap[index] ) != 3 ) {
continue; }
807 thisIsland.clear ( );
811 findAllPointNeighbours ( YZMap, XZMap, XYMap, visitedMap,
static_cast<proshade_signed
>(dim),
static_cast<proshade_signed
>(xIt),
static_cast<proshade_signed
>(yIt),
static_cast<proshade_signed
>(zIt), &thisIsland );
816 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( thisIsland.size() ); iter++ )
818 if ( ( pow( map[thisIsland.at(iter)][0], 2.0 ) + pow( map[thisIsland.at(iter)][1], 2.0) ) > maxHeight )
820 maxHeight = pow( map[thisIsland.at(iter)][0], 2.0 ) + pow( map[thisIsland.at(iter)][1], 2.0 );
821 maxIndex = thisIsland.at(iter);
856 proshade_signed x, y, z, peakX, peakY, peakZ, newIter, ptrIter;
857 proshade_unsign noNeighbours = pow( ((peakSize*2)+1), 3) * 4;
860 std::vector < proshade_unsign > allPeaksRepresentatives;
864 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( allPeaksRepresentatives.size() ); iter++ )
870 z = ( allPeaksRepresentatives.at(iter) % (dim*dim) ) % dim;
871 y = ( allPeaksRepresentatives.at(iter) - z ) % (dim*dim) / dim;
872 x = ( allPeaksRepresentatives.at(iter) - z - ( y * dim ) ) / (dim*dim);
875 proshade_double* neighbours =
new proshade_double[noNeighbours];
880 for ( proshade_signed xCh = -peakSize; xCh <= +peakSize; xCh++ )
882 for ( proshade_signed yCh = -peakSize; yCh <= +peakSize; yCh++ )
884 for ( proshade_signed zCh = -peakSize; zCh <= +peakSize; zCh++ )
887 peakX = x + xCh;
if ( peakX >= dim ) { peakX -= dim; };
if ( peakX < 0 ) { peakX += dim; }
888 peakY = y + yCh;
if ( peakY >= dim ) { peakY -= dim; };
if ( peakY < 0 ) { peakY += dim; }
889 peakZ = z + zCh;
if ( peakZ >= dim ) { peakZ -= dim; };
if ( peakZ < 0 ) { peakZ += dim; }
890 newIter = peakX * (dim*dim) + peakY * dim + peakZ;
893 allPeaksWithNeighbours->at(iter)[ptrIter] =
static_cast<proshade_double
> ( peakX );
894 allPeaksWithNeighbours->at(iter)[ptrIter+1] =
static_cast<proshade_double
> ( peakY );
895 allPeaksWithNeighbours->at(iter)[ptrIter+2] =
static_cast<proshade_double
> ( peakZ );
896 allPeaksWithNeighbours->at(iter)[ptrIter+3] = pow ( map[newIter][0], 2.0 ) + pow( map[newIter][1], 2.0 );
925 std::vector< proshade_double* > allHigherIndices;
926 proshade_unsign smoothingLag = std::floor ( smoothingFraction * dim );
927 proshade_double ZScoreThreshold = noIQRs;
928 proshade_signed *signals;
929 proshade_double *scoreOverVals, *filteredY, *avgFilter, *stdFilter, *subVec, *medianIQR, *YZMap, *XZMap, *XYMap;
932 if ( dim <= smoothingLag + 2)
938 allocateSmoothingZScoreMemory ( dim, scoreOverVals, signals, filteredY, avgFilter, stdFilter, subVec, medianIQR, YZMap, XZMap, XYMap, smoothingLag );
941 getXAxisArraysSmoothedZScorePeaks ( dim, smoothingLag, ZScoreThreshold, signals, filteredY, avgFilter, stdFilter, subVec, medianIQR, scoreOverVals, map, YZMap );
944 getYAxisArraysSmoothedZScorePeaks ( dim, smoothingLag, ZScoreThreshold, signals, filteredY, avgFilter, stdFilter, subVec, medianIQR, scoreOverVals, map, XZMap );
947 getZAxisArraysSmoothedZScorePeaks ( dim, smoothingLag, ZScoreThreshold, signals, filteredY, avgFilter, stdFilter, subVec, medianIQR, scoreOverVals, map, XYMap );
959 return ( allHigherIndices );
978 std::vector< proshade_double* > allPeaks =
getAllPeaksSmoothedZ ( map, dim, smoothingFraction, noIQRs, peakSize );
981 proshade_double highestPeak = 0.0;
982 proshade_unsign highestPeakIndex = 0;
983 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign>( allPeaks.size() ); iter++ )
985 if ( allPeaks.at(iter)[4] > highestPeak ) { highestPeak = allPeaks.at(iter)[4]; highestPeakIndex = iter; }
989 proshade_double* rotMat =
new proshade_double[9];
996 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( allPeaks.size() ); iter++ )
998 delete[] allPeaks.at(iter);