37 gemmi::Structure structure = gemmi::read_structure ( gemmi::MaybeGzipped ( fName ) );
39 catch ( std::runtime_error& e )
62 gemmi::Ccp4<float> map;
65 map.read_ccp4 ( gemmi::MaybeGzipped (fName.c_str() ) );
67 catch ( std::runtime_error& e )
109 void ProSHADE_internal_io::readInMapHeader ( gemmi::Ccp4<float> *map, proshade_unsign *xDimInds, proshade_unsign *yDimInds, proshade_unsign *zDimInds, proshade_single *xDim, proshade_single *yDim, proshade_single *zDim, proshade_single *aAng, proshade_single *bAng, proshade_single *cAng, proshade_signed *xFrom, proshade_signed *yFrom, proshade_signed *zFrom, proshade_signed *xAxOrigin, proshade_signed *yAxOrigin, proshade_signed *zAxOrigin, proshade_unsign *xAxOrder, proshade_unsign *yAxOrder, proshade_unsign *zAxOrder, proshade_unsign *xGridInds, proshade_unsign *yGridInds, proshade_unsign *zGridInds )
112 *xDimInds =
static_cast<proshade_unsign
> ( map->header_i32 ( 1 ) );
113 *yDimInds =
static_cast<proshade_unsign
> ( map->header_i32 ( 2 ) );
114 *zDimInds =
static_cast<proshade_unsign
> ( map->header_i32 ( 3 ) );
116 *xFrom =
static_cast<proshade_signed
> ( map->header_i32 ( 5 ) );
117 *yFrom =
static_cast<proshade_signed
> ( map->header_i32 ( 6 ) );
118 *zFrom =
static_cast<proshade_signed
> ( map->header_i32 ( 7 ) );
120 *xDim =
static_cast<proshade_single
> ( map->header_float ( 11 ) );
121 *yDim =
static_cast<proshade_single
> ( map->header_float ( 12 ) );
122 *zDim =
static_cast<proshade_single
> ( map->header_float ( 13 ) );
124 *aAng =
static_cast<proshade_single
> ( map->header_float ( 14 ) );
125 *bAng =
static_cast<proshade_single
> ( map->header_float ( 15 ) );
126 *cAng =
static_cast<proshade_single
> ( map->header_float ( 16 ) );
128 *xAxOrigin =
static_cast<proshade_signed
> ( map->header_i32 ( 50 ) ) + (*xFrom);
129 *yAxOrigin =
static_cast<proshade_signed
> ( map->header_i32 ( 51 ) ) + (*yFrom);
130 *zAxOrigin =
static_cast<proshade_signed
> ( map->header_i32 ( 52 ) ) + (*zFrom);
132 *xAxOrder =
static_cast<proshade_unsign
> ( map->header_i32 ( 17 ) );
133 *yAxOrder =
static_cast<proshade_unsign
> ( map->header_i32 ( 18 ) );
134 *zAxOrder =
static_cast<proshade_unsign
> ( map->header_i32 ( 19 ) );
136 *xGridInds =
static_cast<proshade_unsign
> ( map->header_i32 ( 8 ) );
137 *yGridInds =
static_cast<proshade_unsign
> ( map->header_i32 ( 9 ) );
138 *zGridInds =
static_cast<proshade_unsign
> ( map->header_i32 ( 10 ) );
141 if ( *xGridInds != *xDimInds )
143 *xDim = *xDim * (
static_cast<proshade_single
> ( *xDimInds ) /
static_cast<proshade_single
> ( *xGridInds ) );
144 *xGridInds = *xDimInds;
147 if ( *yGridInds != *yDimInds )
149 *yDim = *yDim * (
static_cast<proshade_single
> ( *yDimInds ) /
static_cast<proshade_single
> ( *yGridInds ) );
150 *yGridInds = *yDimInds;
153 if ( *zGridInds != *zDimInds )
155 *zDim = *zDim * (
static_cast<proshade_single
> ( *zDimInds ) /
static_cast<proshade_single
> ( *zGridInds ) );
156 *zGridInds = *zDimInds;
178 void ProSHADE_internal_io::readInMapData ( gemmi::Ccp4<float> *gemmiMap, proshade_double*& map, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_unsign xAxOrder, proshade_unsign yAxOrder, proshade_unsign zAxOrder )
181 proshade_unsign *axOrdArr =
new proshade_unsign[3];
182 proshade_unsign *axDimArr =
new proshade_unsign[3];
183 proshade_unsign arrPos = 0;
188 axDimArr[0] = xDimInds;
189 axDimArr[1] = yDimInds;
190 axDimArr[2] = zDimInds;
193 map =
new proshade_double [xDimInds * yDimInds * zDimInds];
197 for ( axOrdArr[0] = 0; axOrdArr[0] < axDimArr[xAxOrder-1]; axOrdArr[0]++ )
199 for ( axOrdArr[1] = 0; axOrdArr[1] < axDimArr[yAxOrder-1]; axOrdArr[1]++ )
201 for ( axOrdArr[2] = 0; axOrdArr[2] < axDimArr[zAxOrder-1]; axOrdArr[2]++ )
203 arrPos = axOrdArr[2] + axDimArr[zAxOrder-1] * ( axOrdArr[1] + axDimArr[yAxOrder-1] * axOrdArr[0] );
204 map[arrPos] =
static_cast< proshade_double
> ( gemmiMap->grid.get_value_q(
static_cast< int > ( axOrdArr[xAxOrder-1] ),
205 static_cast< int > ( axOrdArr[yAxOrder-1] ),
206 static_cast< int > ( axOrdArr[zAxOrder-1] ) ) );
234 void ProSHADE_internal_io::readInMapData ( gemmi::Ccp4<int8_t> *gemmiMap, proshade_double*& map, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_unsign xAxOrder, proshade_unsign yAxOrder, proshade_unsign zAxOrder )
237 proshade_unsign *axOrdArr =
new proshade_unsign[3];
238 proshade_unsign *axDimArr =
new proshade_unsign[3];
239 proshade_unsign arrPos = 0;
244 axDimArr[0] = xDimInds;
245 axDimArr[1] = yDimInds;
246 axDimArr[2] = zDimInds;
249 map =
new proshade_double [xDimInds * yDimInds * zDimInds];
253 for ( axOrdArr[0] = 0; axOrdArr[0] < axDimArr[xAxOrder-1]; axOrdArr[0]++ )
255 for ( axOrdArr[1] = 0; axOrdArr[1] < axDimArr[yAxOrder-1]; axOrdArr[1]++ )
257 for ( axOrdArr[2] = 0; axOrdArr[2] < axDimArr[zAxOrder-1]; axOrdArr[2]++ )
259 arrPos = axOrdArr[2] + axDimArr[zAxOrder-1] * ( axOrdArr[1] + axDimArr[yAxOrder-1] * axOrdArr[0] );
260 map[arrPos] =
static_cast< proshade_double
> ( gemmiMap->grid.get_value_q(
static_cast< int > ( axOrdArr[xAxOrder-1] ),
261 static_cast< int > ( axOrdArr[yAxOrder-1] ),
262 static_cast< int > ( axOrdArr[zAxOrder-1] ) ) );
294 void ProSHADE_internal_io::applyMask ( proshade_double*& map, std::string maskFile, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_signed verbose, proshade_double* maskArray, proshade_unsign maXInds, proshade_unsign maYInds, proshade_unsign maZInds )
297 std::stringstream hlpSS;
298 hlpSS <<
"Reading mask file " << maskFile;
302 if ( ( maskArray !=
nullptr ) && ( maXInds != 0 ) && ( maYInds != 0 ) && ( maZInds != 0 ) )
310 if ( maskFile ==
"" ) {
return; }
313 gemmi::Ccp4<float> mask;
314 mask.read_ccp4 ( gemmi::MaybeGzipped ( maskFile.c_str() ) );
317 mask.setup ( gemmi::GridSetup::ReorderOnly, 0 );
320 proshade_unsign xDI, yDI, zDI, xAOR, yAOR, zAOR, xGI, yGI, zGI;
321 proshade_single xDS, yDS, zDS, aA, bA, cA;
322 proshade_signed xF, yF, zF, xAO, yAO, zAO;
333 proshade_double* internalMask =
nullptr;
340 delete[] internalMask;
365 void ProSHADE_internal_io::applyMaskFromArray ( proshade_double*& map, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_double*& mask, proshade_unsign xDimIndsMsk, proshade_unsign yDimIndsMsk, proshade_unsign zDimIndsMsk, proshade_signed verbose )
368 size_t origVolume = xDimInds * yDimInds * zDimInds;
369 size_t newVolume = xDimIndsMsk * yDimIndsMsk * zDimIndsMsk;
370 proshade_double* maskFinal;
373 if ( ( xDimIndsMsk != xDimInds ) || ( yDimIndsMsk != yDimInds ) || ( zDimIndsMsk != zDimInds ) )
376 fftw_complex* origCoeffs =
new fftw_complex [newVolume ];
377 fftw_complex* origCoeffsHKL =
new fftw_complex [newVolume ];
378 fftw_complex* modifCoeffs =
new fftw_complex [origVolume];
379 fftw_complex* modifCoeffsHKL =
new fftw_complex [origVolume];
380 fftw_complex* inMap =
new fftw_complex [newVolume ];
381 fftw_complex* outMap =
new fftw_complex [origVolume];
392 for (
size_t iter = 0; iter < origVolume; iter++ ) { modifCoeffsHKL[iter][0] = 0.0; modifCoeffsHKL[iter][1] = 0.0; }
395 for (
size_t iter = 0; iter < newVolume; iter++ ) { inMap[iter][0] = mask[iter]; inMap[iter][1] = 0.0; }
398 fftw_plan planForwardFourier = fftw_plan_dft_3d (
static_cast< int > ( xDimIndsMsk ),
static_cast< int > ( yDimIndsMsk ),
static_cast< int > ( zDimIndsMsk ), inMap, origCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
399 fftw_plan inverseFoourier = fftw_plan_dft_3d (
static_cast< int > ( xDimInds ),
static_cast< int > ( yDimInds ),
static_cast< int > ( zDimInds ), modifCoeffs, outMap, FFTW_BACKWARD, FFTW_ESTIMATE );
402 proshade_signed xPre, yPre, zPre;
403 xPre = std::abs ( (
static_cast< proshade_signed
> ( xDimInds ) -
static_cast< proshade_signed
> ( xDimIndsMsk ) ) / 2 );
404 yPre = std::abs ( (
static_cast< proshade_signed
> ( yDimInds ) -
static_cast< proshade_signed
> ( yDimIndsMsk ) ) / 2 );
405 zPre = std::abs ( (
static_cast< proshade_signed
> ( zDimInds ) -
static_cast< proshade_signed
> ( zDimIndsMsk ) ) / 2 );
407 if ( ( (
static_cast< proshade_signed
> ( xDimInds ) -
static_cast< proshade_signed
> ( xDimIndsMsk ) ) % 2 ) == 1 ) { xPre -= 1; }
408 if ( ( (
static_cast< proshade_signed
> ( yDimInds ) -
static_cast< proshade_signed
> ( yDimIndsMsk ) ) % 2 ) == 1 ) { yPre -= 1; }
409 if ( ( (
static_cast< proshade_signed
> ( zDimInds ) -
static_cast< proshade_signed
> ( zDimIndsMsk ) ) % 2 ) == 1 ) { zPre -= 1; }
412 fftw_execute ( planForwardFourier );
415 proshade_signed maskMapIndex = 0;
416 proshade_signed densMapIndex = 0;
417 proshade_signed xMaskPos, yMaskPos, zMaskPos, xDensPos, yDensPos, zDensPos;
418 proshade_signed maskH, maskK, maskL;
421 for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xDimIndsMsk ); xIt++ )
423 for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yDimIndsMsk ); yIt++ )
425 for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zDimIndsMsk ); zIt++ )
428 maskH = xIt +
static_cast< proshade_signed
> ( (xDimIndsMsk+1) / 2 );
if ( maskH >=
static_cast< proshade_signed
> ( xDimIndsMsk ) ) { maskH -= xDimIndsMsk; }
429 maskK = yIt +
static_cast< proshade_signed
> ( (yDimIndsMsk+1) / 2 );
if ( maskK >=
static_cast< proshade_signed
> ( yDimIndsMsk ) ) { maskK -= yDimIndsMsk; }
430 maskL = zIt +
static_cast< proshade_signed
> ( (zDimIndsMsk+1) / 2 );
if ( maskL >=
static_cast< proshade_signed
> ( zDimIndsMsk ) ) { maskL -= zDimIndsMsk; }
433 maskMapIndex = zIt +
static_cast< proshade_signed
> ( zDimIndsMsk ) * ( yIt +
static_cast< proshade_signed
> ( yDimIndsMsk ) * xIt );
434 densMapIndex = maskL +
static_cast< proshade_signed
> ( zDimIndsMsk ) * ( maskK +
static_cast< proshade_signed
> ( yDimIndsMsk ) * maskH );
437 origCoeffsHKL[densMapIndex][0] = origCoeffs[maskMapIndex][0];
438 origCoeffsHKL[densMapIndex][1] = origCoeffs[maskMapIndex][1];
444 for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xDimInds ); xIt++ )
446 for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yDimInds ); yIt++ )
448 for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zDimInds ); zIt++ )
451 if ( xDimIndsMsk >= xDimInds ) { xMaskPos = xIt + xPre; }
452 else { xMaskPos = xIt - xPre; }
456 if ( yDimIndsMsk >= yDimInds ) { yMaskPos = yIt + yPre; }
457 else { yMaskPos = yIt - yPre; }
461 if ( zDimIndsMsk >= zDimInds ) { zMaskPos = zIt + zPre; }
462 else { zMaskPos = zIt - zPre; }
466 if ( ( xMaskPos < 0 ) || ( xMaskPos >=
static_cast< proshade_signed
> ( xDimIndsMsk ) ) ) {
continue; }
467 if ( ( yMaskPos < 0 ) || ( yMaskPos >=
static_cast< proshade_signed
> ( yDimIndsMsk ) ) ) {
continue; }
468 if ( ( zMaskPos < 0 ) || ( zMaskPos >=
static_cast< proshade_signed
> ( zDimIndsMsk ) ) ) {
continue; }
471 maskMapIndex = zMaskPos +
static_cast< proshade_signed
> ( zDimIndsMsk ) * ( yMaskPos +
static_cast< proshade_signed
> ( yDimIndsMsk ) * xMaskPos );
472 densMapIndex = zDensPos +
static_cast< proshade_signed
> ( zDimInds ) * ( yDensPos +
static_cast< proshade_signed
> ( yDimInds ) * xDensPos );
475 modifCoeffsHKL[densMapIndex][0] = origCoeffsHKL[maskMapIndex][0];
476 modifCoeffsHKL[densMapIndex][1] = origCoeffsHKL[maskMapIndex][1];
482 for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xDimInds ); xIt++ )
484 for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yDimInds ); yIt++ )
486 for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zDimInds ); zIt++ )
489 maskH = xIt +
static_cast< proshade_signed
> ( xDimInds / 2 );
if ( maskH >=
static_cast< proshade_signed
> ( xDimInds ) ) { maskH -= xDimInds; }
490 maskK = yIt +
static_cast< proshade_signed
> ( yDimInds / 2 );
if ( maskK >=
static_cast< proshade_signed
> ( yDimInds ) ) { maskK -= yDimInds; }
491 maskL = zIt +
static_cast< proshade_signed
> ( zDimInds / 2 );
if ( maskL >=
static_cast< proshade_signed
> ( zDimInds ) ) { maskL -= zDimInds; }
494 maskMapIndex = zIt +
static_cast< proshade_signed
> ( zDimInds ) * ( yIt +
static_cast< proshade_signed
> ( yDimInds ) * xIt );
495 densMapIndex = maskL +
static_cast< proshade_signed
> ( zDimInds ) * ( maskK +
static_cast< proshade_signed
> ( yDimInds ) * maskH );
498 modifCoeffs[densMapIndex][0] = modifCoeffsHKL[maskMapIndex][0];
499 modifCoeffs[densMapIndex][1] = modifCoeffsHKL[maskMapIndex][1];
505 fftw_execute ( inverseFoourier );
508 maskFinal =
new proshade_double [origVolume];
512 for (
size_t iter = 0; iter < origVolume; iter++ ) { maskFinal[iter] = outMap[iter][0]; }
515 fftw_destroy_plan ( planForwardFourier );
516 fftw_destroy_plan ( inverseFoourier );
518 delete[] modifCoeffs;
519 delete[] origCoeffsHKL;
520 delete[] modifCoeffsHKL;
526 maskFinal =
new proshade_double [origVolume];
528 for (
size_t iter = 0; iter < origVolume; iter++ ) { maskFinal[iter] = mask[iter]; }
532 for (
size_t iter = 0; iter < static_cast< size_t > ( xDimInds * yDimInds * zDimInds ); iter++ ) { map[iter] *= maskFinal[iter]; }
563 void ProSHADE_internal_io::applyWeights ( proshade_double*& map, std::string weightsFile, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_signed verbose, proshade_double* weightsArray, proshade_unsign waXInds, proshade_unsign waYInds, proshade_unsign waZInds )
566 std::stringstream hlpSS;
570 if ( ( weightsArray !=
nullptr ) && ( waXInds != 0 ) && ( waYInds != 0 ) && ( waZInds != 0 ) )
578 if ( weightsFile ==
"" ) {
return; }
581 gemmi::Ccp4<float> weights;
582 weights.read_ccp4 ( gemmi::MaybeGzipped ( weightsFile.c_str() ) );
585 weights.setup ( gemmi::GridSetup::ReorderOnly, 0 );
588 proshade_unsign xDI, yDI, zDI, xAOR, yAOR, zAOR, xGI, yGI, zGI;
589 proshade_single xDS, yDS, zDS, aA, bA, cA;
590 proshade_signed xF, yF, zF, xAO, yAO, zAO;
601 proshade_double* internalWeights =
nullptr;
608 delete[] internalWeights;
633 void ProSHADE_internal_io::applyWeightsFromArray ( proshade_double*& map, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_double*& weights, proshade_unsign xDimIndsWgh, proshade_unsign yDimIndsWgh, proshade_unsign zDimIndsWgh, proshade_signed verbose )
636 proshade_double* weightsFinal;
637 size_t origVolume = xDimInds * yDimInds * zDimInds;
638 size_t newVolume = xDimIndsWgh * yDimIndsWgh * zDimIndsWgh;
641 if ( ( xDimIndsWgh != xDimInds ) || ( yDimIndsWgh != yDimInds ) || ( zDimIndsWgh != zDimInds ) )
644 fftw_complex* origCoeffs =
new fftw_complex [newVolume ];
645 fftw_complex* origCoeffsHKL =
new fftw_complex [newVolume ];
646 fftw_complex* modifCoeffs =
new fftw_complex [origVolume];
647 fftw_complex* modifCoeffsHKL =
new fftw_complex [origVolume];
648 fftw_complex* inMap =
new fftw_complex [newVolume ];
649 fftw_complex* outMap =
new fftw_complex [origVolume];
660 for (
size_t iter = 0; iter < origVolume; iter++ ) { modifCoeffsHKL[iter][0] = 0.0; modifCoeffsHKL[iter][1] = 0.0; }
663 for (
size_t iter = 0; iter < newVolume; iter++ ) { inMap[iter][0] = weights[iter]; inMap[iter][1] = 0.0; }
666 fftw_plan planForwardFourier = fftw_plan_dft_3d (
static_cast< int > ( xDimIndsWgh ),
static_cast< int > ( yDimIndsWgh ),
static_cast< int > ( zDimIndsWgh ), inMap, origCoeffs, FFTW_FORWARD, FFTW_ESTIMATE );
667 fftw_plan inverseFoourier = fftw_plan_dft_3d (
static_cast< int > ( xDimInds ),
static_cast< int > ( yDimInds ),
static_cast< int > ( zDimInds ), modifCoeffs, outMap, FFTW_BACKWARD, FFTW_ESTIMATE );
670 proshade_signed xPre, yPre, zPre;
671 xPre = std::abs ( (
static_cast< proshade_signed
> ( xDimInds ) -
static_cast< proshade_signed
> ( xDimIndsWgh ) ) / 2 );
672 yPre = std::abs ( (
static_cast< proshade_signed
> ( yDimInds ) -
static_cast< proshade_signed
> ( yDimIndsWgh ) ) / 2 );
673 zPre = std::abs ( (
static_cast< proshade_signed
> ( zDimInds ) -
static_cast< proshade_signed
> ( zDimIndsWgh ) ) / 2 );
675 if ( ( (
static_cast< proshade_signed
> ( xDimInds ) -
static_cast< proshade_signed
> ( xDimIndsWgh ) ) % 2 ) == 1 ) { xPre -= 1; }
676 if ( ( (
static_cast< proshade_signed
> ( yDimInds ) -
static_cast< proshade_signed
> ( yDimIndsWgh ) ) % 2 ) == 1 ) { yPre -= 1; }
677 if ( ( (
static_cast< proshade_signed
> ( zDimInds ) -
static_cast< proshade_signed
> ( zDimIndsWgh ) ) % 2 ) == 1 ) { zPre -= 1; }
680 fftw_execute ( planForwardFourier );
683 proshade_signed maskMapIndex = 0;
684 proshade_signed densMapIndex = 0;
685 proshade_signed xMaskPos, yMaskPos, zMaskPos, xDensPos, yDensPos, zDensPos;
686 proshade_signed maskH, maskK, maskL;
689 for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xDimIndsWgh ); xIt++ )
691 for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yDimIndsWgh ); yIt++ )
693 for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zDimIndsWgh ); zIt++ )
696 maskH = xIt +
static_cast< proshade_signed
> ( (xDimIndsWgh+1) / 2 );
if ( maskH >=
static_cast< proshade_signed
> ( xDimIndsWgh ) ) { maskH -= xDimIndsWgh; }
697 maskK = yIt +
static_cast< proshade_signed
> ( (yDimIndsWgh+1) / 2 );
if ( maskK >=
static_cast< proshade_signed
> ( yDimIndsWgh ) ) { maskK -= yDimIndsWgh; }
698 maskL = zIt +
static_cast< proshade_signed
> ( (zDimIndsWgh+1) / 2 );
if ( maskL >=
static_cast< proshade_signed
> ( zDimIndsWgh ) ) { maskL -= zDimIndsWgh; }
701 maskMapIndex = zIt +
static_cast< proshade_signed
> ( zDimIndsWgh ) * ( yIt +
static_cast< proshade_signed
> ( yDimIndsWgh ) * xIt );
702 densMapIndex = maskL +
static_cast< proshade_signed
> ( zDimIndsWgh ) * ( maskK +
static_cast< proshade_signed
> ( yDimIndsWgh ) * maskH );
705 origCoeffsHKL[densMapIndex][0] = origCoeffs[maskMapIndex][0];
706 origCoeffsHKL[densMapIndex][1] = origCoeffs[maskMapIndex][1];
712 for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xDimInds ); xIt++ )
714 for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yDimInds ); yIt++ )
716 for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zDimInds ); zIt++ )
719 if ( xDimIndsWgh >= xDimInds ) { xMaskPos = xIt + xPre; }
720 else { xMaskPos = xIt - xPre; }
724 if ( yDimIndsWgh >= yDimInds ) { yMaskPos = yIt + yPre; }
725 else { yMaskPos = yIt - yPre; }
729 if ( zDimIndsWgh >= zDimInds ) { zMaskPos = zIt + zPre; }
730 else { zMaskPos = zIt - zPre; }
734 if ( ( xMaskPos < 0 ) || ( xMaskPos >=
static_cast< proshade_signed
> ( xDimIndsWgh ) ) ) {
continue; }
735 if ( ( yMaskPos < 0 ) || ( yMaskPos >=
static_cast< proshade_signed
> ( yDimIndsWgh ) ) ) {
continue; }
736 if ( ( zMaskPos < 0 ) || ( zMaskPos >=
static_cast< proshade_signed
> ( zDimIndsWgh ) ) ) {
continue; }
739 maskMapIndex = zMaskPos +
static_cast< proshade_signed
> ( zDimIndsWgh ) * ( yMaskPos +
static_cast< proshade_signed
> ( yDimIndsWgh ) * xMaskPos );
740 densMapIndex = zDensPos +
static_cast< proshade_signed
> ( zDimInds ) * ( yDensPos +
static_cast< proshade_signed
> ( yDimInds ) * xDensPos );
743 modifCoeffsHKL[densMapIndex][0] = origCoeffsHKL[maskMapIndex][0];
744 modifCoeffsHKL[densMapIndex][1] = origCoeffsHKL[maskMapIndex][1];
750 for ( proshade_signed xIt = 0; xIt < static_cast< proshade_signed > ( xDimInds ); xIt++ )
752 for ( proshade_signed yIt = 0; yIt < static_cast< proshade_signed > ( yDimInds ); yIt++ )
754 for ( proshade_signed zIt = 0; zIt < static_cast< proshade_signed > ( zDimInds ); zIt++ )
757 maskH = xIt +
static_cast< proshade_signed
> ( xDimInds / 2 );
if ( maskH >=
static_cast< proshade_signed
> ( xDimInds ) ) { maskH -= xDimInds; }
758 maskK = yIt +
static_cast< proshade_signed
> ( yDimInds / 2 );
if ( maskK >=
static_cast< proshade_signed
> ( yDimInds ) ) { maskK -= yDimInds; }
759 maskL = zIt +
static_cast< proshade_signed
> ( zDimInds / 2 );
if ( maskL >=
static_cast< proshade_signed
> ( zDimInds ) ) { maskL -= zDimInds; }
762 maskMapIndex = zIt +
static_cast< proshade_signed
> ( zDimInds ) * ( yIt +
static_cast< proshade_signed
> ( yDimInds ) * xIt );
763 densMapIndex = maskL +
static_cast< proshade_signed
> ( zDimInds ) * ( maskK +
static_cast< proshade_signed
> ( yDimInds ) * maskH );
766 modifCoeffs[densMapIndex][0] = modifCoeffsHKL[maskMapIndex][0];
767 modifCoeffs[densMapIndex][1] = modifCoeffsHKL[maskMapIndex][1];
773 fftw_execute ( inverseFoourier );
776 weightsFinal =
new proshade_double [origVolume];
780 for (
size_t iter = 0; iter < origVolume; iter++ ) { weightsFinal[iter] = outMap[iter][0]; }
783 fftw_destroy_plan ( planForwardFourier );
784 fftw_destroy_plan ( inverseFoourier );
786 delete[] modifCoeffs;
787 delete[] origCoeffsHKL;
788 delete[] modifCoeffsHKL;
794 weightsFinal =
new proshade_double [origVolume];
796 for (
size_t iter = 0; iter < origVolume; iter++ ) { weightsFinal[iter] = weights[iter]; }
800 fftw_complex* inMap =
new fftw_complex [origVolume];
801 fftw_complex* outMap =
new fftw_complex [origVolume];
804 fftw_plan planForwardFourier = fftw_plan_dft_3d (
static_cast< int > ( xDimInds ),
static_cast< int > ( yDimInds ),
static_cast< int > ( zDimInds ), inMap, outMap, FFTW_FORWARD, FFTW_ESTIMATE );
805 fftw_plan inverseFoourier = fftw_plan_dft_3d (
static_cast< int > ( xDimInds ),
static_cast< int > ( yDimInds ),
static_cast< int > ( zDimInds ), outMap, inMap, FFTW_BACKWARD, FFTW_ESTIMATE );
808 for (
size_t iter = 0; iter < static_cast< size_t > ( origVolume ); iter++ ) { inMap[iter][0] = map[iter]; inMap[iter][1] = 0.0; }
811 fftw_execute ( planForwardFourier );
814 proshade_double normFactor =
static_cast<proshade_double
> ( origVolume );
815 for (
size_t iter = 0; iter < static_cast< size_t > ( origVolume ); iter++ ) { outMap[iter][0] *= weightsFinal[iter] / normFactor; outMap[iter][1] *= weightsFinal[iter] / normFactor; }
818 fftw_execute ( inverseFoourier );
821 for (
size_t iter = 0; iter < static_cast< size_t > ( xDimInds * yDimInds * zDimInds ); iter++ ) { map[iter] = inMap[iter][0]; }
824 delete[] weightsFinal;
827 fftw_destroy_plan ( planForwardFourier );
828 fftw_destroy_plan ( inverseFoourier );
868 void ProSHADE_internal_io::writeOutMapHeader ( gemmi::Ccp4<float> *map, proshade_unsign xDimInds, proshade_unsign yDimInds, proshade_unsign zDimInds, proshade_single xDim, proshade_single yDim, proshade_single zDim, proshade_single aAng, proshade_single bAng, proshade_single cAng, proshade_signed xFrom, proshade_signed yFrom, proshade_signed zFrom, proshade_signed xAxOrigin, proshade_signed yAxOrigin, proshade_signed zAxOrigin, proshade_unsign xAxOrder, proshade_unsign yAxOrder, proshade_unsign zAxOrder, proshade_unsign xGridInds, proshade_unsign yGridInds, proshade_unsign zGridInds, std::string title,
int mode )
871 map->set_header_i32 ( 1 ,
static_cast<int32_t
> ( xDimInds ) );
872 map->set_header_i32 ( 2 ,
static_cast<int32_t
> ( yDimInds ) );
873 map->set_header_i32 ( 3 ,
static_cast<int32_t
> ( zDimInds ) );
874 map->set_header_i32 ( 4 ,
static_cast<int32_t
> ( mode ) );
875 map->set_header_i32 ( 5 ,
static_cast<int32_t
> ( xFrom ) );
876 map->set_header_i32 ( 6 ,
static_cast<int32_t
> ( yFrom ) );
877 map->set_header_i32 ( 7 ,
static_cast<int32_t
> ( zFrom ) );
878 map->set_header_i32 ( 8 ,
static_cast<int32_t
> ( xGridInds ) );
879 map->set_header_i32 ( 9 ,
static_cast<int32_t
> ( yGridInds ) );
880 map->set_header_i32 ( 10,
static_cast<int32_t
> ( zGridInds ) );
881 map->set_header_float ( 11,
static_cast<float> ( xDim ) );
882 map->set_header_float ( 12,
static_cast<float> ( yDim ) );
883 map->set_header_float ( 13,
static_cast<float> ( zDim ) );
884 map->set_header_float ( 14,
static_cast<float> ( aAng ) );
885 map->set_header_float ( 15,
static_cast<float> ( bAng ) );
886 map->set_header_float ( 16,
static_cast<float> ( cAng ) );
887 map->set_header_i32 ( 17,
static_cast<int32_t
> ( xAxOrder ) );
888 map->set_header_i32 ( 18,
static_cast<int32_t
> ( yAxOrder ) );
889 map->set_header_i32 ( 19,
static_cast<int32_t
> ( zAxOrder ) );
890 if ( map->grid.spacegroup ) { map->set_header_i32 ( 23,
static_cast<int32_t
> ( map->grid.spacegroup->ccp4 ) ); }
891 else { map->set_header_i32 ( 23,
static_cast<int32_t
> ( 1 ) ); }
892 map->set_header_i32 ( 24,
static_cast<int32_t
> ( map->grid.spacegroup->operations().order() * 80 ) );
893 map->set_header_str ( 27,
"CCP4" );
894 map->set_header_i32 ( 28,
static_cast<int32_t
> ( 20140 ) );
895 map->set_header_i32 ( 50,
static_cast<int32_t
> ( xAxOrigin ) );
896 map->set_header_i32 ( 51,
static_cast<int32_t
> ( yAxOrigin ) );
897 map->set_header_i32 ( 52,
static_cast<int32_t
> ( zAxOrigin ) );
898 map->set_header_str ( 53,
"MAP" );
899 if ( gemmi::is_little_endian() ) { map->set_header_i32 ( 54,
static_cast<int32_t
> ( 0x00004144 ) ); }
900 else { map->set_header_i32 ( 54,
static_cast<int32_t
> ( 0x11110000 ) ); }
901 map->set_header_i32 ( 56,
static_cast<int32_t
> ( 1 ) );
902 std::memset (
reinterpret_cast<void*
> ( &(map->ccp4_header.at( 56 )) ),
' ',
static_cast< size_t > ( 800 + map->grid.spacegroup->operations().order() * 80 ) );
903 map->set_header_str ( 57, title );
960 void ProSHADE_internal_io::writeRotationTranslationJSON ( proshade_double trsX1, proshade_double trsY1, proshade_double trsZ1, proshade_double eulA, proshade_double eulB, proshade_double eulG, proshade_double trsX2, proshade_double trsY2, proshade_double trsZ2, std::string fileName )
963 std::ofstream jsonFile;
964 jsonFile.open ( fileName );
967 if ( !jsonFile.is_open( ) )
969 throw ProSHADE_exception (
"Failed to open JSON output file.",
"E000056", __FILE__, __LINE__, __func__,
"Failed to open json file to which the rotation and\n : translation would be written into. Most likely cause is\n : lack of rights to write in the current folder." );
973 proshade_double* rotMat =
new proshade_double[9];
979 jsonFile <<
" \"translationToOrigin\" : [ " << trsX1 <<
", " << trsY1 <<
", " << trsZ1 <<
" ], \n";
981 jsonFile <<
" \"rotationMatrix:\" : [ " << rotMat[0] <<
", " << rotMat[1] <<
", " << rotMat[2] <<
", \n";
982 jsonFile <<
" " << rotMat[3] <<
", " << rotMat[4] <<
", " << rotMat[5] <<
", \n";
983 jsonFile <<
" " << rotMat[6] <<
", " << rotMat[7] <<
", " << rotMat[8] <<
"], \n";
985 jsonFile <<
" \"translationFromRotCenToOverlay\" : [ " << trsX2 <<
", " << trsY2 <<
", " << trsZ2 <<
" ] \n";