46 ProSHADE_internal_spheres::ProSHADE_sphere::ProSHADE_sphere ( proshade_unsign xDimMax, proshade_unsign yDimMax, proshade_unsign zDimMax, proshade_single xSize, proshade_single ySize, proshade_single zSize, proshade_unsign shOrder, std::vector<proshade_single>* spherePos,
bool progressiveMapping, proshade_unsign band, proshade_double* map, proshade_unsign* maxShellBand )
49 this->shellOrder = shOrder;
50 this->sphereWidth = ( spherePos->at(0) + spherePos->at(1) ) / 2.0;
51 this->sphereRadius = spherePos->at(shOrder);
54 proshade_double maxDist = 0.0;
55 if ( shOrder ==
static_cast<proshade_unsign
> ( spherePos->size() - 1 ) ) { maxDist =
static_cast<proshade_double
> ( spherePos->at(spherePos->size()-1) + ( spherePos->at(1) - spherePos->at(0) ) ); }
56 else { maxDist =
static_cast<proshade_double
> ( ( spherePos->at(shOrder) + spherePos->at(shOrder+1) ) / 2.0 ); }
59 this->maxSphereRange =
static_cast<proshade_single
> ( 2.0 * maxDist );
62 this->xDimSampling = xSize /
static_cast<proshade_single
> (xDimMax);
63 this->yDimSampling = ySize /
static_cast<proshade_single
> (yDimMax);
64 this->zDimSampling = zSize /
static_cast<proshade_single
> (zDimMax);
67 proshade_unsign maxCircumference = this->
getMaxCircumference ( xDimMax, yDimMax, zDimMax, this->maxSphereRange, xSize, ySize, zSize );
70 if ( progressiveMapping )
73 this->localAngRes = this->localBandwidth * 2;
77 this->localBandwidth = band;
78 this->localAngRes = this->localBandwidth * 2;
82 if ( *maxShellBand < this->localBandwidth ) { *maxShellBand = this->localBandwidth; }
85 this->mappedData =
new proshade_double[this->localAngRes * this->localAngRes];
89 this->mappedDataRot = NULL;
92 this->
mapData ( map, xDimMax, yDimMax, zDimMax );
104 delete[] this->mappedData;
106 if ( this->mappedDataRot != NULL )
108 delete[] this->mappedDataRot;
127 proshade_signed xFromInd =
static_cast<proshade_signed
> ( xDimMax / 2 ) -
128 static_cast<proshade_signed
> ( (maxRange/2) / this->xDimSampling );
129 proshade_signed yFromInd =
static_cast<proshade_signed
> ( yDimMax / 2 ) -
130 static_cast<proshade_signed
> ( (maxRange/2) / this->yDimSampling );
131 proshade_signed zFromInd =
static_cast<proshade_signed
> ( zDimMax / 2 ) -
132 static_cast<proshade_signed
> ( (maxRange/2) / this->zDimSampling );
134 proshade_signed xToInd =
static_cast<proshade_signed
> ( xDimMax / 2 ) +
135 static_cast<proshade_signed
> ( (maxRange/2) / this->xDimSampling );
136 proshade_signed yToInd =
static_cast<proshade_signed
> ( yDimMax / 2 ) +
137 static_cast<proshade_signed
> ( (maxRange/2) / this->yDimSampling );
138 proshade_signed zToInd =
static_cast<proshade_signed
> ( zDimMax / 2 ) +
139 static_cast<proshade_signed
> ( (maxRange/2) / this->zDimSampling );
142 if ( xFromInd < 0 ) { xFromInd = 0; }
143 if ( yFromInd < 0 ) { yFromInd = 0; }
144 if ( zFromInd < 0 ) { zFromInd = 0; }
145 if ( xToInd >
static_cast<proshade_signed
> ( xDimMax ) ) { xToInd = xDimMax; }
146 if ( yToInd >
static_cast<proshade_signed
> ( yDimMax ) ) { yToInd = yDimMax; }
147 if ( zToInd >
static_cast<proshade_signed
> ( zDimMax ) ) { zToInd = zDimMax; }
150 proshade_unsign xDimSZ = xToInd - xFromInd;
151 proshade_unsign yDimSZ = yToInd - yFromInd;
152 proshade_unsign zDimSZ = zToInd - zFromInd;
155 if ( ( xDimSZ == 0 ) && ( yDimSZ == 0 ) && ( zDimSZ == 0 ) ) { xDimSZ = 1; yDimSZ = 1; zDimSZ = 1; }
158 std::vector<proshade_unsign> dimSizes;
162 std::sort ( dimSizes.begin(), dimSizes.end() );
163 proshade_unsign maxDim = dimSizes.at(2);
164 proshade_unsign midDim = dimSizes.at(1);
167 return ( maxDim + midDim );
184 proshade_double x, y, z, xRelative, yRelative, zRelative;
185 proshade_signed xBottom, yBottom, zBottom, xTop, yTop, zTop;
186 std::vector<proshade_double> lonCO ( this->localAngRes + 1 );
187 std::vector<proshade_double> latCO ( this->localAngRes + 1 );
188 std::vector<proshade_double> c000 = std::vector<proshade_double> ( 4 );
189 std::vector<proshade_double> c001 = std::vector<proshade_double> ( 4 );
190 std::vector<proshade_double> c010 = std::vector<proshade_double> ( 4 );
191 std::vector<proshade_double> c011 = std::vector<proshade_double> ( 4 );
192 std::vector<proshade_double> c100 = std::vector<proshade_double> ( 4 );
193 std::vector<proshade_double> c101 = std::vector<proshade_double> ( 4 );
194 std::vector<proshade_double> c110 = std::vector<proshade_double> ( 4 );
195 std::vector<proshade_double> c111 = std::vector<proshade_double> ( 4 );
196 std::vector<proshade_double> c00 = std::vector<proshade_double> ( 4 );
197 std::vector<proshade_double> c01 = std::vector<proshade_double> ( 4 );
198 std::vector<proshade_double> c10 = std::vector<proshade_double> ( 4 );
199 std::vector<proshade_double> c11 = std::vector<proshade_double> ( 4 );
200 std::vector<proshade_double> c0 = std::vector<proshade_double> ( 4 );
201 std::vector<proshade_double> c1 = std::vector<proshade_double> ( 4 );
204 this->getLongitudeCutoffs ( &lonCO );
205 this->getLattitudeCutoffs ( &latCO );
208 for (
unsigned int thIt = 0; thIt < this->localAngRes; thIt++ )
210 for (
unsigned int phIt = 0; phIt < this->localAngRes; phIt++ )
213 this->getInterpolationXYZ ( &x, &y, &z, thIt, &lonCO, phIt, &latCO );
216 this->getXYZTopBottoms ( xDimMax, yDimMax, zDimMax, x, y, z, &xBottom, &yBottom, &zBottom, &xTop, &yTop, &zTop );
219 if ( !getMapPoint ( map, xDimMax, yDimMax, zDimMax, xBottom, yBottom, zBottom, &c000 ) ) { this->mappedData[phIt * this->localAngRes + thIt] = 0.0;
continue; }
220 if ( !getMapPoint ( map, xDimMax, yDimMax, zDimMax, xBottom, yBottom, zTop , &c001 ) ) { this->mappedData[phIt * this->localAngRes + thIt] = 0.0;
continue; }
221 if ( !getMapPoint ( map, xDimMax, yDimMax, zDimMax, xBottom, yTop , zBottom, &c010 ) ) { this->mappedData[phIt * this->localAngRes + thIt] = 0.0;
continue; }
222 if ( !getMapPoint ( map, xDimMax, yDimMax, zDimMax, xBottom, yTop , zTop , &c011 ) ) { this->mappedData[phIt * this->localAngRes + thIt] = 0.0;
continue; }
223 if ( !getMapPoint ( map, xDimMax, yDimMax, zDimMax, xTop , yBottom, zBottom, &c100 ) ) { this->mappedData[phIt * this->localAngRes + thIt] = 0.0;
continue; }
224 if ( !getMapPoint ( map, xDimMax, yDimMax, zDimMax, xTop , yBottom, zTop , &c101 ) ) { this->mappedData[phIt * this->localAngRes + thIt] = 0.0;
continue; }
225 if ( !getMapPoint ( map, xDimMax, yDimMax, zDimMax, xTop , yTop , zBottom, &c110 ) ) { this->mappedData[phIt * this->localAngRes + thIt] = 0.0;
continue; }
226 if ( !getMapPoint ( map, xDimMax, yDimMax, zDimMax, xTop , yTop , zTop , &c111 ) ) { this->mappedData[phIt * this->localAngRes + thIt] = 0.0;
continue; }
229 xRelative = ( x - ( ( xBottom -
static_cast<proshade_signed
> ( ( ( xDimMax ) / 2 ) ) ) * this->xDimSampling ) ) / this->xDimSampling;
230 this->interpolateAlongFirst ( c000, c001, c010, c011, c100, c101, c110, c111, &c00, &c01, &c10, &c11, xRelative );
233 yRelative = ( y - ( ( yBottom -
static_cast<proshade_signed
> ( ( ( yDimMax ) / 2 ) ) ) * this->yDimSampling ) ) / this->yDimSampling;
234 this->interpolateAlongSecond ( c00, c01, c10, c11, &c0, &c1, yRelative );
237 zRelative = ( z - ( ( zBottom -
static_cast<proshade_signed
> ( ( ( zDimMax ) / 2 ) ) ) * this->zDimSampling ) ) / this->zDimSampling;
238 this->mappedData[phIt * this->localAngRes + thIt] = ( c0.at(3) * ( 1.0 - zRelative ) ) + ( c1.at(3) * zRelative );
266 proshade_signed posIter = zPos + zDimMax * ( yPos + yDimMax * xPos );
269 if ( ( xPos < 0 ) || ( xPos >=
static_cast<proshade_signed
> ( xDimMax ) ) ) {
return (
false ); }
270 if ( ( yPos < 0 ) || ( yPos >=
static_cast<proshade_signed
> ( yDimMax ) ) ) {
return (
false ); }
271 if ( ( zPos < 0 ) || ( zPos >=
static_cast<proshade_signed
> ( zDimMax ) ) ) {
return (
false ); }
274 interpVec->at(0) =
static_cast<proshade_double
> ( xPos * this->xDimSampling );
275 interpVec->at(1) =
static_cast<proshade_double
> ( yPos * this->yDimSampling );
276 interpVec->at(2) =
static_cast<proshade_double
> ( zPos * this->zDimSampling );
277 interpVec->at(3) = map[posIter];
293 for ( proshade_unsign iter = 0; iter <= this->localAngRes; iter++ )
295 lonCO->at(iter) =
static_cast<proshade_double
> ( iter ) *
296 ( (
static_cast<proshade_double
> ( M_PI ) * 2.0 ) /
static_cast<proshade_double
> ( this->localAngRes ) ) -
297 (
static_cast<proshade_double
> ( M_PI ) );
314 for ( proshade_unsign iter = 0; iter <= this->localAngRes; iter++ )
316 latCO->at(iter) = (
static_cast<proshade_double
> ( iter ) *
317 (
static_cast<proshade_double
> ( M_PI ) /
static_cast<proshade_double
> ( this->localAngRes ) ) -
318 (
static_cast<proshade_double
> ( M_PI ) / 2.0 ) );
342 *x = this->sphereRadius * cos( ( lonCO->at(thetaIt) + lonCO->at(thetaIt+1) ) / 2.0 ) *
343 cos( ( latCO->at(phiIt) + latCO->at(phiIt+1) ) / 2.0 );
344 *y = this->sphereRadius * sin( ( lonCO->at(thetaIt) + lonCO->at(thetaIt+1) ) / 2.0 ) *
345 cos( ( latCO->at(phiIt) + latCO->at(phiIt+1) ) / 2.0 );
346 *z = this->sphereRadius * sin( ( latCO->at(phiIt) + latCO->at(phiIt+1) ) / 2.0 );
370 void ProSHADE_internal_spheres::ProSHADE_sphere::getXYZTopBottoms ( proshade_unsign xDimMax, proshade_unsign yDimMax, proshade_unsign zDimMax, proshade_double x, proshade_double y, proshade_double z, proshade_signed* xBottom, proshade_signed* yBottom, proshade_signed* zBottom, proshade_signed* xTop, proshade_signed* yTop, proshade_signed* zTop )
373 *xBottom = std::floor ( (x / this->xDimSampling) ) + (xDimMax/2);
374 *yBottom = std::floor ( (y / this->yDimSampling) ) + (yDimMax/2);
375 *zBottom = std::floor ( (z / this->zDimSampling) ) + (zDimMax/2);
377 *xTop = *xBottom + 1;
378 *yTop = *yBottom + 1;
379 *zTop = *zBottom + 1;
396 return ( this->localBandwidth );
409 return ( this->localAngRes );
422 return ( this->mappedData );
445 void ProSHADE_internal_spheres::ProSHADE_sphere::interpolateAlongFirst ( std::vector<proshade_double> c000, std::vector<proshade_double> c001, std::vector<proshade_double> c010, std::vector<proshade_double> c011, std::vector<proshade_double> c100, std::vector<proshade_double> c101, std::vector<proshade_double> c110, std::vector<proshade_double> c111, std::vector<proshade_double>* c00, std::vector<proshade_double>* c01, std::vector<proshade_double>* c10, std::vector<proshade_double>* c11, proshade_double xd )
448 c00->at(0) = ( this->xDimSampling * xd ) + c000.at(0);
449 c00->at(1) = c000.at(1);
450 c00->at(2) = c000.at(2);
451 c00->at(3) = ( c000.at(3) * ( 1.0 - xd ) ) + ( c100.at(3) * xd );
454 c01->at(0) = ( this->xDimSampling * xd ) + c001.at(0);
455 c01->at(1) = c001.at(1);
456 c01->at(2) = c001.at(2);
457 c01->at(3) = ( c001.at(3) * ( 1.0 - xd ) ) + ( c101.at(3) * xd );
460 c10->at(0) = ( this->xDimSampling * xd ) + c010.at(0);
461 c10->at(1) = c010.at(1);
462 c10->at(2) = c010.at(2);
463 c10->at(3) = ( c010.at(3) * ( 1.0 - xd ) ) + ( c110.at(3) * xd );
466 c11->at(0) = ( this->xDimSampling * xd ) + c011.at(0);
467 c11->at(1) = c011.at(1);
468 c11->at(2) = c011.at(2);
469 c11->at(3) = ( c011.at(3) * ( 1.0 - xd ) ) + ( c111.at(3) * xd );
494 c0->at(0) = c00.at(0);
495 c0->at(1) = ( this->yDimSampling * yd ) + c00.at(1);
496 c0->at(2) = c00.at(2);
497 c0->at(3) = ( c00.at(3) * ( 1.0 - yd ) ) + ( c10.at(3) * yd );
500 c1->at(0) = c01.at(0);
501 c1->at(1) = ( this->yDimSampling * yd ) + c01.at(1);
502 c1->at(2) = c01.at(2);
503 c1->at(3) = ( c01.at(3) * ( 1.0 - yd ) ) + ( c11.at(3) * yd );
521 if (
static_cast<proshade_unsign
> ( std::ceil ( circumference / 2 ) ) % 2 == 0 )
523 return (
static_cast<proshade_unsign
> ( std::ceil ( circumference / 2 ) ) );
527 return (
static_cast<proshade_unsign
> ( std::ceil ( circumference / 2 ) ) + 1 );
543 proshade_single ret =
static_cast<proshade_single
> ( resolution / 2.0 );
546 while ( std::floor ( maxMapRange / ret ) < 10 )
567 proshade_double sphereDistanceAsFractionOfTotal =
static_cast<proshade_double
> ( sphereDist ) / ( maxMapRange / 2.0 );
568 proshade_unsign ret = 0;
571 for ( proshade_unsign iter = 2; iter < static_cast<proshade_unsign> ( 10000 ); iter++ )
573 if ( ProSHADE_internal_precomputedVals::glIntMaxDists[iter] >= sphereDistanceAsFractionOfTotal )
593 return ( this->sphereRadius );
602 this->mappedDataRot = NULL;
603 this->mappedDataRot =
new proshade_double[this->localAngRes * this->localAngRes];
619 this->mappedDataRot[pos] = value;
633 return ( this->mappedDataRot[pos] );
652 this->angularDim = dim;
653 this->radiusMin = this->radius - ( radRange / 2.0 );
654 this->radiusMax = this->radius + ( radRange / 2.0 );
655 this->representedAngle = repAng;
656 this->sphereNumber = sphNo;
659 this->axesValues =
new proshade_double[dim*dim];
663 for ( proshade_unsign iter = 0; iter < ( dim * dim ); iter++ ) { this->axesValues[iter] = 0.0; }
674 if ( this->axesValues !=
nullptr )
676 delete[] this->axesValues;
687 return ( this->radius );
697 return ( this->radiusMax );
708 return ( this->angularDim );
719 return ( this->radiusMin );
730 return ( this->representedAngle );
741 return ( this->sphereNumber );
752 return ( this->peaks );
768 proshade_double lonSampling = ( M_PI ) /
static_cast<proshade_double
> ( this->angularDim );
769 proshade_double latSampling = ( M_PI * 2.0 ) /
static_cast<proshade_double
> ( this->angularDim );
771 proshade_double lat, lon, cX, cY, cZ, c000, c001, c010, c011, c100, c101, c110, c111, c00, c01, c10, c11, c0, c1, xRelative, yRelative, zRelative, eulerAlpha, eulerBeta, eulerGamma, mapX, mapY, mapZ;
772 proshade_signed xBottom, xTop, yBottom, yTop, zBottom, zTop, mapIndex;
775 for ( proshade_signed lonIt = 0; lonIt < static_cast<proshade_signed> ( this->angularDim ); lonIt++ )
777 for ( proshade_signed latIt = 0; latIt < static_cast<proshade_signed> ( this->angularDim ); latIt++ )
780 lon =
static_cast<proshade_double
> ( lonIt ) * lonSampling;
781 lat =
static_cast<proshade_double
> ( latIt ) * latSampling;
782 cX = 1.0 * std::sin ( lon ) * std::cos ( lat );
783 cY = 1.0 * std::sin ( lon ) * std::sin ( lat );
784 cZ = 1.0 * std::cos ( lon );
793 xBottom = std::floor ( mapX );
if ( xBottom < 0.0 ) { xBottom += this->angularDim; }
if ( xBottom >=
static_cast<proshade_signed
> ( this->angularDim ) ) { xBottom -=
static_cast<proshade_signed
> ( this->angularDim ); }
794 yBottom = std::floor ( mapY );
if ( yBottom < 0.0 ) { yBottom += this->angularDim; }
if ( yBottom >=
static_cast<proshade_signed
> ( this->angularDim ) ) { yBottom -=
static_cast<proshade_signed
> ( this->angularDim ); }
795 zBottom = std::floor ( mapZ );
if ( zBottom < 0.0 ) { zBottom += this->angularDim; }
if ( zBottom >=
static_cast<proshade_signed
> ( this->angularDim ) ) { zBottom -=
static_cast<proshade_signed
> ( this->angularDim ); }
796 xTop = std::ceil ( mapX );
if ( xTop < 0.0 ) { xTop += this->angularDim; }
if ( xTop >=
static_cast<proshade_signed
> ( this->angularDim ) ) { xTop -=
static_cast<proshade_signed
> ( this->angularDim ); }
797 yTop = std::ceil ( mapY );
if ( yTop < 0.0 ) { yTop += this->angularDim; }
if ( yTop >=
static_cast<proshade_signed
> ( this->angularDim ) ) { yTop -=
static_cast<proshade_signed
> ( this->angularDim ); }
798 zTop = std::ceil ( mapZ );
if ( zTop < 0.0 ) { zTop += this->angularDim; }
if ( zTop >=
static_cast<proshade_signed
> ( this->angularDim ) ) { zTop -=
static_cast<proshade_signed
> ( this->angularDim ); }
801 mapIndex = zBottom + this->angularDim * ( yBottom + this->angularDim * xBottom );
802 c000 = pow( rotFun[mapIndex][0], 2.0 ) + pow( rotFun[mapIndex][1], 2.0 );
805 mapIndex = zTop + this->angularDim * ( yBottom + this->angularDim * xBottom );
806 c001 = pow( rotFun[mapIndex][0], 2.0 ) + pow( rotFun[mapIndex][1], 2.0 );
809 mapIndex = zBottom + this->angularDim * ( yTop + this->angularDim * xBottom );
810 c010 = pow( rotFun[mapIndex][0], 2.0 ) + pow( rotFun[mapIndex][1], 2.0 );
813 mapIndex = zTop + this->angularDim * ( yTop + this->angularDim * xBottom );
814 c011 = pow( rotFun[mapIndex][0], 2.0 ) + pow( rotFun[mapIndex][1], 2.0 );
817 mapIndex = zBottom + this->angularDim * ( yBottom + this->angularDim * xTop );
818 c100 = pow( rotFun[mapIndex][0], 2.0 ) + pow( rotFun[mapIndex][1], 2.0 );
821 mapIndex = zTop + this->angularDim * ( yBottom + this->angularDim * xTop );
822 c101 = pow( rotFun[mapIndex][0], 2.0 ) + pow( rotFun[mapIndex][1], 2.0 );
825 mapIndex = zBottom + this->angularDim * ( yTop + this->angularDim * xTop );
826 c110 = pow( rotFun[mapIndex][0], 2.0 ) + pow( rotFun[mapIndex][1], 2.0 );
829 mapIndex = zTop + this->angularDim * ( yTop + this->angularDim * xTop );
830 c111 = pow( rotFun[mapIndex][0], 2.0 ) + pow( rotFun[mapIndex][1], 2.0 );
833 xRelative = mapX - std::floor( mapX );
834 c00 = ( c000 * ( 1.0 - xRelative ) ) + ( c100 * xRelative );
835 c01 = ( c001 * ( 1.0 - xRelative ) ) + ( c101 * xRelative );
836 c10 = ( c010 * ( 1.0 - xRelative ) ) + ( c110 * xRelative );
837 c11 = ( c011 * ( 1.0 - xRelative ) ) + ( c111 * xRelative );
840 yRelative = mapY - std::floor( mapY );
841 c0 = ( c00 * ( 1.0 - yRelative ) ) + ( c10 * yRelative );
842 c1 = ( c01 * ( 1.0 - yRelative ) ) + ( c11 * yRelative );
845 zRelative = mapZ - std::floor( mapZ );
848 mapIndex = lonIt + ( latIt *
static_cast<proshade_unsign
> ( this->angularDim ) );
849 this->axesValues[mapIndex] = ( c0 * ( 1.0 - zRelative ) ) + ( c1 * zRelative );
867 return ( this->axesValues[longitude + ( lattitude *
static_cast<proshade_unsign
> ( this->angularDim ) )] );
881 proshade_double c00, c01, c10, c11, c0, c1, latRelative, lonRelative;
882 proshade_signed latTop, latBottom, lonTop, lonBottom, gridIndex;
885 latBottom = std::floor ( lattitude );
if ( latBottom < 0.0 ) { latBottom += this->angularDim; }
if ( latBottom >=
static_cast<proshade_signed
> ( this->angularDim ) ) { latBottom -= this->angularDim; }
886 lonBottom = std::floor ( longitude );
if ( lonBottom < 0.0 ) { lonBottom += this->angularDim; }
if ( lonBottom >=
static_cast<proshade_signed
> ( this->angularDim ) ) { lonBottom -= this->angularDim; }
887 latTop = std::ceil ( lattitude );
if ( latTop < 0.0 ) { latTop += this->angularDim; }
if ( latTop >=
static_cast<proshade_signed
> ( this->angularDim ) ) { latTop -= this->angularDim; }
888 lonTop = std::ceil ( longitude );
if ( lonTop < 0.0 ) { lonTop += this->angularDim; }
if ( lonTop >=
static_cast<proshade_signed
> ( this->angularDim ) ) { lonTop -= this->angularDim; }
891 gridIndex = lonBottom + ( latBottom *
static_cast<proshade_unsign
> ( this->angularDim ) );
892 c00 = this->axesValues[gridIndex];
894 gridIndex = lonBottom + ( latTop *
static_cast<proshade_unsign
> ( this->angularDim ) );
895 c01 = this->axesValues[gridIndex];
897 gridIndex = lonTop + ( latBottom *
static_cast<proshade_unsign
> ( this->angularDim ) );
898 c10 = this->axesValues[gridIndex];
900 gridIndex = lonTop + ( latTop *
static_cast<proshade_unsign
> ( this->angularDim ) );
901 c11 = this->axesValues[gridIndex];
904 lonRelative = longitude - std::floor( longitude );
905 c0 = ( c00 * ( 1.0 - lonRelative ) ) + ( c10 * lonRelative );
906 c1 = ( c01 * ( 1.0 - lonRelative ) ) + ( c11 * lonRelative );
909 latRelative = lattitude - std::floor ( lattitude );
910 proshade_double res = ( c0 * ( 1.0 - latRelative ) ) + ( c1 * latRelative );
927 proshade_double currentHeight;
928 proshade_signed nbLat, nbLon;
932 for ( proshade_signed latIt = 0; latIt < static_cast<proshade_signed> ( this->angularDim ); latIt++ )
934 for ( proshade_signed lonIt = 0; lonIt < static_cast<proshade_signed> ( this->angularDim ); lonIt++ )
937 currentHeight = this->getSphereLatLonPosition ( latIt, lonIt );
941 for ( proshade_signed latRound = -noSmNeighbours; latRound <= noSmNeighbours; latRound++ )
943 for ( proshade_signed lonRound = -noSmNeighbours; lonRound <= noSmNeighbours; lonRound++ )
946 if ( latRound == 0 && lonRound == 0 ) {
continue; }
949 nbLat = latIt + latRound;
950 nbLon = lonIt + lonRound;
951 if ( nbLat < 0 ) { nbLat += this->angularDim; }
if ( nbLat >=
static_cast<proshade_signed
> ( this->angularDim ) ) { nbLat -= this->angularDim; }
952 if ( nbLon < 0 ) { nbLon += this->angularDim; }
if ( nbLon >=
static_cast<proshade_signed
> ( this->angularDim ) ) { nbLon -= this->angularDim; }
955 if ( this->getSphereLatLonPosition ( nbLat, nbLon ) > currentHeight ) { isPeak =
false;
break; }
958 if ( !isPeak ) {
break; }
964 this->peaks.emplace_back ( std::pair<proshade_unsign,proshade_unsign> ( latIt, lonIt ) );
989 proshade_double curHeight;
990 std::vector< proshade_unsign > dels ( 0, this->peaks.size() );
993 for ( proshade_unsign peakIt = 0; peakIt < static_cast<proshade_unsign> ( this->peaks.size() ); peakIt++ )
996 curHeight = this->getSphereLatLonPosition ( this->peaks.at(peakIt).first, this->peaks.at(peakIt).second );
999 if ( curHeight < peakThres )
1006 std::sort ( dels.begin(), dels.end(), std::greater <proshade_unsign>() );
1009 for ( proshade_unsign delIt = 0; delIt < static_cast< proshade_unsign > ( dels.size() ); delIt++ )
1011 this->peaks.erase ( this->peaks.begin() + dels.at(delIt) );
1032 this->dimension = angDim;
1033 this->lonSampling = ( M_PI ) /
static_cast<proshade_double
> ( this->dimension );
1034 this->latSampling = ( M_PI * 2.0 ) /
static_cast<proshade_double
> ( this->dimension );
1037 this->latFrom =
static_cast<proshade_double
> ( lat ) * this->latSampling;
1038 this->latTo =
static_cast<proshade_double
> ( lat ) * this->latSampling;
1039 this->lonFrom =
static_cast<proshade_double
> ( lon ) * this->lonSampling;
1040 this->lonTo =
static_cast<proshade_double
> ( lon ) * this->lonSampling;
1041 this->latFromInds = lat;
1042 this->latToInds = lat;
1043 this->lonFromInds = lon;
1044 this->lonToInds = lon;
1048 this->latMinLonMinXYZ =
new proshade_double[3];
1049 this->latMaxLonMinXYZ =
new proshade_double[3];
1050 this->latMinLonMaxXYZ =
new proshade_double[3];
1051 this->latMaxLonMaxXYZ =
new proshade_double[3];
1058 this->computeCornerPositions ( );
1071 delete[] this->latMinLonMinXYZ;
1072 delete[] this->latMaxLonMinXYZ;
1073 delete[] this->latMinLonMaxXYZ;
1074 delete[] this->latMaxLonMaxXYZ;
1082 this->latMinLonMinXYZ[0] = 1.0 * std::sin ( this->lonFrom ) * std::cos ( this->latFrom );
1083 this->latMinLonMinXYZ[1] = 1.0 * std::sin ( this->lonFrom ) * std::sin ( this->latFrom );
1084 this->latMinLonMinXYZ[2] = 1.0 * std::cos ( this->lonFrom );
1086 this->latMaxLonMinXYZ[0] = 1.0 * std::sin ( this->lonFrom ) * std::cos ( this->latTo );
1087 this->latMaxLonMinXYZ[1] = 1.0 * std::sin ( this->lonFrom ) * std::sin ( this->latTo );
1088 this->latMaxLonMinXYZ[2] = 1.0 * std::cos ( this->lonFrom );
1090 this->latMinLonMaxXYZ[0] = 1.0 * std::sin ( this->lonTo ) * std::cos ( this->latFrom );
1091 this->latMinLonMaxXYZ[1] = 1.0 * std::sin ( this->lonTo ) * std::sin ( this->latFrom );
1092 this->latMinLonMaxXYZ[2] = 1.0 * std::cos ( this->lonTo );
1094 this->latMaxLonMaxXYZ[0] = 1.0 * std::sin ( this->lonTo ) * std::cos ( this->latTo );
1095 this->latMaxLonMaxXYZ[1] = 1.0 * std::sin ( this->lonTo ) * std::sin ( this->latTo );
1096 this->latMaxLonMaxXYZ[2] = 1.0 * std::cos ( this->lonTo );
1112 proshade_signed smallerAngul = newAngul - this->dimension;
1113 proshade_signed largerAngul = newAngul + this->dimension;
1116 proshade_signed ret = std::min ( std::abs ( currentAngul - newAngul ), std::min ( std::abs ( currentAngul - smallerAngul ), std::abs ( currentAngul - largerAngul ) ) );
1140 bool peakAdded =
false;
1141 std::stringstream hlpSS;
1142 std::stringstream hlpSS2;
1145 proshade_double xPos = 1.0 * std::sin ( lon * this->lonSampling ) * std::cos ( lat * this->latSampling );
1146 proshade_double yPos = 1.0 * std::sin ( lon * this->lonSampling ) * std::sin ( lat * this->latSampling );
1147 proshade_double zPos = 1.0 * std::cos ( lon * this->lonSampling );
1148 hlpSS2 <<
"Peak " << xPos <<
" ; " << yPos <<
" ; " << zPos <<
" is close enough to group with corner ";
1151 if (
ProSHADE_internal_maths::vectorOrientationSimilaritySameDirection ( xPos, yPos, zPos, this->latMaxLonMinXYZ[0], this->latMaxLonMinXYZ[1], this->latMaxLonMinXYZ[2], cosTol ) && !peakAdded ) { peakAdded =
true; hlpSS2 << this->latMaxLonMinXYZ[0] <<
" ; " << this->latMaxLonMinXYZ[1] <<
" ; " << this->latMaxLonMinXYZ[2]; }
1152 if (
ProSHADE_internal_maths::vectorOrientationSimilaritySameDirection ( xPos, yPos, zPos, this->latMinLonMaxXYZ[0], this->latMinLonMaxXYZ[1], this->latMinLonMaxXYZ[2], cosTol ) && !peakAdded ) { peakAdded =
true; hlpSS2 << this->latMinLonMaxXYZ[0] <<
" ; " << this->latMinLonMaxXYZ[1] <<
" ; " << this->latMinLonMaxXYZ[2]; }
1153 if (
ProSHADE_internal_maths::vectorOrientationSimilaritySameDirection ( xPos, yPos, zPos, this->latMaxLonMaxXYZ[0], this->latMaxLonMaxXYZ[1], this->latMaxLonMaxXYZ[2], cosTol ) && !peakAdded ) { peakAdded =
true; hlpSS2 << this->latMaxLonMaxXYZ[0] <<
" ; " << this->latMaxLonMaxXYZ[1] <<
" ; " << this->latMaxLonMaxXYZ[2]; }
1159 hlpSS <<
"Peak group dimensions changed from LAT " << this->latFromInds <<
" - " << this->latToInds <<
" and LON " << this->lonFromInds <<
" - " << this->lonToInds <<
" to ";
1163 proshade_unsign largerCorner, smallerCorner;
1164 bool latCornersDone =
false;
1165 bool lonCornersDone =
false;
1168 if ( ( this->latFromInds <= this->latToInds ) && !( ( lat >= this->latFromInds ) && ( lat <= this->latToInds ) ) )
1171 smallerCorner = angularDistanceWithBorders ( lat, this->latFromInds );
1172 largerCorner = angularDistanceWithBorders ( lat, this->latToInds );
1174 if ( smallerCorner < largerCorner )
1176 this->latFromInds = lat;
1177 latCornersDone =
true;
1179 if ( smallerCorner > largerCorner )
1181 this->latToInds = lat;
1182 latCornersDone =
true;
1184 if ( smallerCorner == largerCorner )
1186 if ( lat < this->latFromInds ) { this->latFromInds = lat; latCornersDone =
true; }
1187 else if ( lat > this->latToInds ) { this->latToInds = lat; latCornersDone =
true; }
1191 if ( ( this->latFromInds > this->latToInds ) && !( ( lat >= this->latFromInds ) || ( lat <= this->latToInds ) ) )
1194 smallerCorner = angularDistanceWithBorders ( lat, this->latFromInds );
1195 largerCorner = angularDistanceWithBorders ( lat, this->latToInds );
1197 if ( smallerCorner < largerCorner )
1199 this->latFromInds = lat;
1200 latCornersDone =
true;
1202 if ( smallerCorner > largerCorner )
1204 this->latToInds = lat;
1205 latCornersDone =
true;
1207 if ( smallerCorner == largerCorner )
1209 if ( lat < this->latFromInds ) { this->latFromInds = lat; latCornersDone =
true; }
1210 else if ( lat > this->latToInds ) { this->latToInds = lat; latCornersDone =
true; }
1216 if ( ( this->lonFromInds <= this->lonToInds ) && !( ( lon >= this->lonFromInds ) && ( lon <= this->lonToInds ) ) )
1219 smallerCorner = angularDistanceWithBorders ( lon, this->lonFromInds );
1220 largerCorner = angularDistanceWithBorders ( lon, this->lonToInds );
1222 if ( smallerCorner < largerCorner )
1224 this->lonFromInds = lon;
1225 lonCornersDone =
true;
1227 if ( smallerCorner > largerCorner )
1229 this->lonToInds = lon;
1230 lonCornersDone =
true;
1232 if ( smallerCorner == largerCorner )
1234 if ( lon < this->lonFromInds ) { this->lonFromInds = lon; lonCornersDone =
true; }
1235 else if ( lon > this->lonToInds ) { this->lonToInds = lon; lonCornersDone =
true; }
1239 if ( ( this->lonFromInds > this->lonToInds ) && !( ( lon >= this->lonFromInds ) || ( lon <= this->lonToInds ) ) )
1242 smallerCorner = angularDistanceWithBorders ( lon, this->lonFromInds );
1243 largerCorner = angularDistanceWithBorders ( lon, this->lonToInds );
1245 if ( smallerCorner < largerCorner )
1247 this->lonFromInds = lon;
1248 lonCornersDone =
true;
1250 if ( smallerCorner > largerCorner )
1252 this->lonToInds = lon;
1253 lonCornersDone =
true;
1255 if ( smallerCorner == largerCorner )
1257 if ( lon < this->lonFromInds ) { this->lonFromInds = lon; lonCornersDone =
true; }
1258 else if ( lon > this->lonToInds ) { this->lonToInds = lon; lonCornersDone =
true; }
1263 if ( latCornersDone )
1265 this->latFrom =
static_cast<proshade_double
> ( this->latFromInds ) * this->latSampling;
1266 this->latTo =
static_cast<proshade_double
> ( this->latToInds ) * this->latSampling;
1269 if ( lonCornersDone )
1271 this->lonFrom =
static_cast<proshade_double
> ( this->lonFromInds ) * this->lonSampling;
1272 this->lonTo =
static_cast<proshade_double
> ( this->lonToInds ) * this->lonSampling;
1276 this->computeCornerPositions ( );
1277 hlpSS <<
"LAT " << this->latFromInds <<
" - " << this->latToInds <<
" and LON " << this->lonFromInds <<
" - " << this->lonToInds <<
" ( peak position LAT " << lat <<
" LON " << lon <<
" )";
1281 bool isSphereNew =
true;
1282 for ( proshade_unsign iter = 0; iter < static_cast<proshade_unsign> ( this->spherePositions.size() ); iter++ ) {
if ( this->spherePositions.at(iter) == sphPos ) { isSphereNew =
false; } }
1287 return ( peakAdded );
1298 return ( this->latFromInds );
1309 return ( this->latToInds );
1320 return ( this->lonFromInds );
1331 return ( this->lonToInds );
1342 return ( this->spherePositions );
1369 if ( ( fold - 1 ) != spherePositions.size() ) { return ; }
1372 proshade_double bestPosVal, bestLatInd, bestLonInd;
1373 std::vector< proshade_unsign > spheresFormingFold;
1376 for ( proshade_unsign shIt = 0; shIt < static_cast<proshade_unsign> ( sphereVals.size() ); shIt++ )
1382 this->getBestIndexForFold ( &bestPosVal, &bestLatInd, &bestLonInd, &spheresFormingFold, sphereVals );
1385 if ( bicubicInterp )
1391 proshade_double* detectedSymmetry =
new proshade_double[6];
1394 detectedSymmetry[0] =
static_cast<proshade_double
> ( fold );
1395 detectedSymmetry[1] = 1.0 * std::sin ( bestLonInd * this->lonSampling ) * std::cos ( bestLatInd * this->latSampling );
1396 detectedSymmetry[2] = 1.0 * std::sin ( bestLonInd * this->lonSampling ) * std::sin ( bestLatInd * this->latSampling );
1397 detectedSymmetry[3] = 1.0 * std::cos ( bestLonInd * this->lonSampling );
1398 detectedSymmetry[4] = ( 2.0 * M_PI ) / detectedSymmetry[0];
1399 detectedSymmetry[5] = ( bestPosVal - 1.0 ) / ( detectedSymmetry[0] - 1 );
1402 if ( ( ( std::max ( std::abs ( detectedSymmetry[1] ), std::max ( std::abs ( detectedSymmetry[2] ), std::abs ( detectedSymmetry[3] ) ) ) == std::abs ( detectedSymmetry[1] ) ) && ( detectedSymmetry[1] < 0.0 ) ) ||
1403 ( ( std::max ( std::abs ( detectedSymmetry[1] ), std::max ( std::abs ( detectedSymmetry[2] ), std::abs ( detectedSymmetry[3] ) ) ) == std::abs ( detectedSymmetry[2] ) ) && ( detectedSymmetry[2] < 0.0 ) ) ||
1404 ( ( std::max ( std::abs ( detectedSymmetry[1] ), std::max ( std::abs ( detectedSymmetry[2] ), std::abs ( detectedSymmetry[3] ) ) ) == std::abs ( detectedSymmetry[3] ) ) && ( detectedSymmetry[3] < 0.0 ) ) )
1406 detectedSymmetry[1] *= -1.0;
1407 detectedSymmetry[2] *= -1.0;
1408 detectedSymmetry[3] *= -1.0;
1409 detectedSymmetry[4] *= -1.0;
1416 std::stringstream hlpSS;
1417 hlpSS <<
"Detected group with fold " << detectedSymmetry[0] <<
" along axis " << detectedSymmetry[1] <<
" ; " << detectedSymmetry[2] <<
" ; " << detectedSymmetry[3] <<
" and with peak height " << detectedSymmetry[5];
1433 std::vector< proshade_double > angs;
1436 for ( proshade_unsign shPos = 0; shPos < static_cast<proshade_unsign> ( this->spherePositions.size() ); shPos++ )
1442 for ( proshade_unsign ang1It = 0; ang1It < static_cast<proshade_unsign> ( angs.size() ); ang1It++ )
1444 for ( proshade_unsign ang2It = 1; ang2It < static_cast<proshade_unsign> ( angs.size() ); ang2It++ )
1447 if ( ang1It >= ang2It ) {
continue; }
1455 std::sort ( (*angDiffs).begin(), (*angDiffs).end() );
1456 (*angDiffs).erase ( std::unique ( (*angDiffs).begin(), (*angDiffs).end() ), (*angDiffs).end() );
1476 proshade_double divRem, divBasis, symmErr, angTolerance, angToleranceNext;
1477 proshade_double peakErr = ( M_PI * 2.0 ) / (
static_cast<proshade_double
> ( this->dimension ) );
1480 for ( proshade_unsign diffIt = 0; diffIt < static_cast<proshade_unsign> ( angDiffs->size() ); diffIt++ )
1483 divRem = std::modf (
static_cast<proshade_double
> ( ( 2.0 * M_PI ) / std::abs ( angDiffs->at(diffIt) ) ), &divBasis );
1493 if ( divBasis == 1 ) {
continue; }
1496 if (
static_cast<proshade_unsign
> ( this->spherePositions.size() ) < ( divBasis - 1 ) ) {
continue; }
1499 symmErr = divRem * ( ( 2.0 * M_PI ) /
static_cast<proshade_double
> ( divBasis ) );
1500 angTolerance = std::abs( symmErr / peakErr );
1501 angToleranceNext = ( ( ( 2.0 * M_PI ) /
static_cast<proshade_double
> ( divBasis ) ) - ( ( 2.0 * M_PI ) /
static_cast<proshade_double
> ( divBasis + 1 ) ) ) / peakErr;
1504 if ( angTolerance < std::max ( 3.0, ( 0.1 / peakErr ) ) )
1507 if ( angToleranceNext < std::max ( 1.5, ( 0.1 / peakErr ) ) )
1520 std::sort ( (*foldsToTry).begin(), (*foldsToTry).end(), std::greater <proshade_unsign>() );
1521 foldsToTry->erase ( std::unique ( (*foldsToTry).begin(), (*foldsToTry).end() ), (*foldsToTry).end() );
1538 proshade_double soughtAngle, minSphereVal;
1539 proshade_unsign minSpherePos;
1542 for ( proshade_double fIt = 1.0; fIt < static_cast<proshade_double> ( foldToTry ); fIt += 1.0 )
1545 minSphereVal = 999.9;
1546 soughtAngle = fIt * ( 2.0 * M_PI /
static_cast<proshade_double
> ( foldToTry ) );
1549 for ( proshade_unsign angsIt = 0; angsIt < static_cast<proshade_unsign> ( this->spherePositions.size() ); angsIt++ )
1551 if ( std::abs ( sphereVals.at(this->spherePositions.at(angsIt))->getRepresentedAngle() - soughtAngle ) < sphereAngleTolerance )
1553 if ( minSphereVal > 1.0 - std::cos ( std::abs ( sphereVals.at(this->spherePositions.at(angsIt))->getRepresentedAngle() - soughtAngle ) ) )
1555 minSphereVal = 1.0 - std::cos ( std::abs ( sphereVals.at(this->spherePositions.at(angsIt))->getRepresentedAngle() - soughtAngle ) );
1556 minSpherePos = angsIt;
1562 if ( minSphereVal == 999.9 ) {
break; }
1585 proshade_double curPosVal;
1587 if ( this->latFromInds > this->latToInds ) { this->latToInds += this->dimension; }
1588 if ( this->lonFromInds > this->lonToInds ) { this->lonToInds += this->dimension; }
1591 for ( proshade_unsign latIt = this->latFromInds; latIt <= this->latToInds; latIt++ )
1594 if ( latIt >= this->dimension ) { latIt -= this->dimension; this->latToInds -= this->dimension; }
1596 for ( proshade_unsign lonIt = this->lonFromInds; lonIt <= this->lonToInds; lonIt++ )
1599 if ( lonIt >= this->dimension ) { lonIt -= this->dimension; this->lonToInds -= this->dimension; }
1605 for ( proshade_unsign sphIt = 0; sphIt < static_cast<proshade_unsign> ( spheresFormingFold->size() ); sphIt++ )
1607 curPosVal += sphereVals.at(spheresFormingFold->at(sphIt))->getSphereLatLonPosition ( latIt, lonIt );
1611 if ( curPosVal > *bestPosVal )
1613 *bestPosVal = curPosVal;
1614 *bestLatInd = latIt;
1615 *bestLonInd = lonIt;