#include <reconstruct_fourier_starpu.h>

Inheritance diagram for ProgRecFourierStarPU:

Collaboration diagram for ProgRecFourierStarPU:

Classes
struct	BatchProvider

struct	ComputeConstants

struct	ComputeStarPUResult

struct	Params

struct	SimpleBatchProvider

Public Member Functions
	ProgRecFourierStarPU ()

void	defineParams () override

void	readParams () override

void	show () const override

void	run () override

void	setIO (const FileName &fn_in, const FileName &fn_out) override

Public Member Functions inherited from ProgReconsBase
virtual	~ProgReconsBase ()

Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)

const char *	getParam (const char param, const char subparam, int arg=0)

int	getIntParam (const char *param, int arg=0)

int	getIntParam (const char param, const char subparam, int arg=0)

double	getDoubleParam (const char *param, int arg=0)

double	getDoubleParam (const char param, const char subparam, int arg=0)

float	getFloatParam (const char *param, int arg=0)

float	getFloatParam (const char param, const char subparam, int arg=0)

void	getListParam (const char *param, StringVector &list)

int	getCountParam (const char *param)

bool	checkParam (const char *param)

bool	existsParam (const char *param)

void	addParamsLine (const String &line)

void	addParamsLine (const char *line)

ParamDef *	getParamDef (const char *param) const

virtual void	quit (int exit_code=0) const

virtual int	tryRun ()

void	initProgress (size_t total, size_t stepBin=60)

void	setProgress (size_t value=0)

void	endProgress ()

void	processDefaultComment (const char param, const char left)

void	setDefaultComment (const char param, const char comment)

virtual void	initComments ()

void	setProgramName (const char *name)

void	addUsageLine (const char *line, bool verbatim=false)

void	clearUsage ()

void	addExampleLine (const char *example, bool verbatim=true)

void	addSeeAlsoLine (const char *seeAlso)

void	addKeywords (const char *keywords)

const char *	name () const

virtual void	usage (int verb=0) const

virtual void	usage (const String &param, int verb=2)

int	version () const

virtual void	read (int argc, const char **argv, bool reportErrors=true)

virtual void	read (int argc, char **argv, bool reportErrors=true)

void	read (const String &argumentsLine)

	XmippProgram ()

	XmippProgram (int argc, const char **argv)

virtual	~XmippProgram ()

Static Protected Member Functions
static void	prepareMetaData (const FileName &fn_in, MetaDataVec &SF)

static uint32_t	computeBatchCount (const Params &params, const MetaData &SF)

static void	prepareConstants (const Params &params, const MetaData &SF, const FileName &fn_sym, ComputeConstants &constants)

static void	initStarPU ()

static ComputeStarPUResult	computeStarPU (const Params &params, const MetaData &SF, const ComputeConstants &computeConstants, BatchProvider &batches, bool verbose)

static void	shutdownStarPU ()

static void	postProcessAndSave (const Params &params, const ComputeConstants &computeConstants, const FileName &fn_out, std::complex< float > *tempVolume, float *tempWeights)

Protected Attributes
FileName	fn_in

FileName	fn_out

FileName	fn_sym

struct ProgRecFourierStarPU::Params	params

MetaDataVec	SF

struct ProgRecFourierStarPU::ComputeConstants	computeConstants

Protected Attributes inherited from XmippProgram
int	errorCode

ProgramDef *	progDef
	Program definition and arguments parser. More...

std::map< String, CommentList >	defaultComments

int	argc
	Original command line arguments. More...

const char **	argv

Additional Inherited Members
Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Detailed Description

Definition at line 44 of file reconstruct_fourier_starpu.h.

Constructor & Destructor Documentation

◆ ProgRecFourierStarPU()

ProgRecFourierStarPU::ProgRecFourierStarPU ( )

inline

Definition at line 118 of file reconstruct_fourier_starpu.h.

                            {
         // Check that there is no logical problem in the defines used by program. If found, error is thrown.
         if (!((TILE == 1) || (BLOCK_DIM % TILE == 0))) {
             REPORT_ERROR(ERR_PARAM_INCORRECT,"TILE has to be set to 1(one) or BLOCK_DIM has to be a multiple of TILE");
         }
         if ((SHARED_IMG == 1) && (TILE != 1)) {
             REPORT_ERROR(ERR_PARAM_INCORRECT,"TILE cannot be used while SHARED_IMG is active");
         }
         if (TILE >= BLOCK_DIM) {
             REPORT_ERROR(ERR_PARAM_INCORRECT,"TILE must be smaller than BLOCK_DIM");
         }
         if ((SHARED_BLOB_TABLE == 1) && (PRECOMPUTE_BLOB_VAL == 0)) {
             REPORT_ERROR(ERR_PARAM_INCORRECT,"PRECOMPUTE_BLOB_VAL must be set to 1(one) when SHARED_BLOB_TABLE is active");
         }
     }

Member Function Documentation

◆ computeBatchCount()

uint32_t ProgRecFourierStarPU::computeBatchCount	(	const Params &	params,
		const MetaData &	SF
	)

staticprotected

Definition at line 118 of file reconstruct_fourier_starpu.cpp.

                                                                                                              {
     return (static_cast<uint32_t>(SF.size()) + params.batchSize - 1) / params.batchSize;
 }

◆ computeStarPU()

ProgRecFourierStarPU::ComputeStarPUResult ProgRecFourierStarPU::computeStarPU	(	const Params &	params,
		const MetaData &	SF,
		const ComputeConstants &	computeConstants,
		BatchProvider &	batches,
		bool	verbose
	)

staticprotected

Definition at line 317 of file reconstruct_fourier_starpu.cpp.

                       {
 
     uint32_t maxVolumeIndex = computeConstants.maxVolumeIndex;
     uint32_t paddedImgSize = computeConstants.paddedImgSize;
 
     // Initialize GPU
     padded_image_to_fft_cuda_initialize(paddedImgSize);
     reconstruct_cuda_initialize_constants(maxVolumeIndex, maxVolumeIndex,
                                           static_cast<float>(params.blob.radius),
                                           static_cast<float>(params.blob.alpha),
                                           computeConstants.iDeltaSqrt, static_cast<float>(computeConstants.iw0),
                                           1.f / getBessiOrderAlpha(params.blob));
 
     std::vector<size_t> selFileObjectIds;
     SF.findObjects(selFileObjectIds);
     const uint32_t fftSizeX = maxVolumeIndex / 2;
     const uint32_t fftSizeY = maxVolumeIndex;
 
     PaddedImageToFftArgs imageToFftArg = {
             params.maxResolution * params.maxResolution,
             paddedImgSize,
             fftSizeX, fftSizeY
     };
 
     ReconstructFftArgs reconstructFftArg = {
             static_cast<float>(params.blob.radius),
             maxVolumeIndex,
             params.fastLateBlobbing,
             params.blob.order,
             static_cast<float>(params.blob.alpha),
             static_cast<uint32_t>(computeConstants.R_symmetries.size()),
             fftSizeX, fftSizeY
     };
 
     starpu_data_handle_t resultVolumeHandle = {0};
     starpu_data_handle_t resultWeightsHandle = {0};
 
     ComputeStarPUResult result;
 
     {
         /*
          * Allocate 3D array (continuous) of given size^3.
          * Allocated array is cleared (to zero)
          */
         uint32_t dim = static_cast<uint32_t>(maxVolumeIndex + 1);
         uint32_t dim3 = align(dim * dim * dim, 4); // Make the array size always a multiple of 4, for faster redux operations
         size_t volumeDataSize = dim3 * sizeof(std::complex<float>);
         size_t weightDataSize = dim3 * sizeof(float);
         CHECK_STARPU(starpu_malloc((void **) &result.volumeData, volumeDataSize));
         CHECK_STARPU(starpu_malloc((void **) &result.weightsData, weightDataSize));
         // As a final step, StarPU will redux the data created on GPU into this
         memset(result.volumeData, 0, volumeDataSize);
         memset(result.weightsData, 0, weightDataSize);
 
         starpu_vector_data_register(&resultVolumeHandle, STARPU_MAIN_RAM, (uintptr_t) result.volumeData, dim3, sizeof(std::complex<float>));
         starpu_vector_data_register(&resultWeightsHandle, STARPU_MAIN_RAM, (uintptr_t) result.weightsData, dim3, sizeof(float));
         starpu_data_set_reduction_methods(resultVolumeHandle, &codelets.redux_sum_volume, &codelets.redux_init_volume);
         starpu_data_set_reduction_methods(resultWeightsHandle, &codelets.redux_sum_weights, &codelets.redux_init_weights);
         starpu_data_set_name(resultVolumeHandle, "Result Volume Data");
         starpu_data_set_name(resultWeightsHandle, "Result Weights Data");
     }
 
     starpu_data_handle_t blobTableSquaredHandle = {0};
     if (params.fastLateBlobbing) {
         // Blob Table is not used on GPU, but we need to register empty regardless
         starpu_variable_data_register(&blobTableSquaredHandle, -1, 0, 1);
     } else {
         starpu_variable_data_register(&blobTableSquaredHandle, STARPU_MAIN_RAM,
                                       reinterpret_cast<uintptr_t>(&computeConstants.blobTableSqrt), sizeof(computeConstants.blobTableSqrt));
     }
     starpu_data_set_name(blobTableSquaredHandle, "Blob Table Squared");
 
     const uint32_t totalImages = static_cast<uint32_t>(SF.size());
     const uint32_t maxBatches = batches.maxBatches();
 
     std::vector<LoadProjectionArgs> loadProjectionArgs;
     // WARNING: Backing arrays of this vector must never reallocate, as we use pointers into it for codelet arguments
     loadProjectionArgs.reserve(maxBatches);
 
     // Used in a STARPU_SCRATCH mode, so the handle can be shared between all relevant tasks,
     // as the all will get a different memory to work with.
     starpu_data_handle_t fftScratchMemoryHandle = {0};
     // As documented in http://fftw.org/fftw3_doc/Multi_002dDimensional-DFTs-of-Real-Data.html#Multi_002dDimensional-DFTs-of-Real-Data
     uint32_t fftResultX = paddedImgSize/2+1;
     starpu_matrix_data_register(&fftScratchMemoryHandle, -1, 0, fftResultX /* no padding */, fftResultX, paddedImgSize, sizeof(float) * 2);
     starpu_data_set_name(fftScratchMemoryHandle, "Scratch for raw FFT");
 
     int32_t batch;
     while ((batch = batches.nextBatch()) != -1) {
         const uint32_t batchStart = batch * params.batchSize;
         const uint32_t batchEnd = XMIPP_MIN(batchStart + params.batchSize, totalImages);
         const uint32_t currentBatchSize = static_cast<uint32_t>(batchEnd - batchStart);
 
         starpu_data_handle_t amountLoadedHandle = {0};
         starpu_variable_data_register(&amountLoadedHandle, -1, 0, sizeof(LoadedImagesBuffer));
         starpu_data_set_name(amountLoadedHandle, "Batch Meta Data");
 
         starpu_data_handle_t paddedImagesDataHandle = {0};
         starpu_vector_data_register(&paddedImagesDataHandle, -1, 0, currentBatchSize, align(paddedImgSize * paddedImgSize * sizeof(float), ALIGNMENT));
         starpu_data_set_name(paddedImagesDataHandle, "Batch Padded Image Data");
 
         starpu_data_handle_t traverseSpacesHandle = {0};
         starpu_matrix_data_register(&traverseSpacesHandle, -1, 0, computeConstants.R_symmetries.size(), computeConstants.R_symmetries.size(), currentBatchSize, sizeof(RecFourierProjectionTraverseSpace));
         starpu_data_set_name(traverseSpacesHandle, "Batch Traverse Spaces");
 
         // Submit the task to load the projections
         {
             const size_t argIndex = loadProjectionArgs.size();
 
             // Create new LoadProjectionArgs for this batch
             loadProjectionArgs.push_back(LoadProjectionArgs {
                     batchStart, batchEnd,
                     SF,
                     selFileObjectIds,
                     params.do_weights,
                     computeConstants.R_symmetries,
                     maxVolumeIndex, maxVolumeIndex,
                     static_cast<float>(params.blob.radius),
                     params.fastLateBlobbing,
                     paddedImgSize,
                     fftSizeX, fftSizeY
             });
 
             const LoadProjectionArgs& loadProjectionArg = loadProjectionArgs[argIndex];
 
             starpu_task *loadProjectionsTask = starpu_task_create();
             loadProjectionsTask->name = "LoadProjections";
             loadProjectionsTask->cl = &codelets.load_projections;
             loadProjectionsTask->cl_arg = (void *) &loadProjectionArg;
             loadProjectionsTask->cl_arg_size = sizeof(loadProjectionArg);
             loadProjectionsTask->cl_arg_free = 0; // Do not free! (probably default)
             loadProjectionsTask->handles[0] = amountLoadedHandle;
             loadProjectionsTask->handles[1] = paddedImagesDataHandle;
             loadProjectionsTask->handles[2] = traverseSpacesHandle;
             loadProjectionsTask->synchronous = DEBUG_SYNCHRONOUS_TASKS;
 
             CHECK_STARPU(starpu_task_submit(loadProjectionsTask));
         }
 
         // Compute FFT of padded image data and adjust it (crop, shift, normalize)
         starpu_data_handle_t fftHandle = {0};
         starpu_vector_data_register(&fftHandle, -1, 0, currentBatchSize, fftSizeX * fftSizeY * 2 * sizeof(float));
         starpu_data_set_name(fftHandle, "Batch FFT Data");
 
         {
             starpu_task* paddedImageToFftTask = starpu_task_create();
             paddedImageToFftTask->name = "PaddedImageToFFT";
             paddedImageToFftTask->cl = &codelets.padded_image_to_fft;
             paddedImageToFftTask->handles[0] = paddedImagesDataHandle;
             paddedImageToFftTask->handles[1] = fftHandle;
             paddedImageToFftTask->handles[2] = fftScratchMemoryHandle;
             paddedImageToFftTask->handles[3] = amountLoadedHandle;
             paddedImageToFftTask->cl_arg = &imageToFftArg;
             paddedImageToFftTask->cl_arg_size = sizeof(PaddedImageToFftArgs);
             paddedImageToFftTask->cl_arg_free = 0;
             paddedImageToFftTask->synchronous = DEBUG_SYNCHRONOUS_TASKS;
             CHECK_STARPU(starpu_task_submit(paddedImageToFftTask));
         }
 
         starpu_data_unregister_submit(paddedImagesDataHandle);
 
         {// Submit the actual reconstruction
             starpu_task* reconstructFftTask = starpu_task_create();
             reconstructFftTask->name = "ReconstructFFT";
             reconstructFftTask->cl = &codelets.reconstruct_fft;
             reconstructFftTask->handles[0] = fftHandle;
             reconstructFftTask->handles[1] = traverseSpacesHandle;
             reconstructFftTask->handles[2] = blobTableSquaredHandle;
             reconstructFftTask->handles[3] = resultVolumeHandle;
             reconstructFftTask->handles[4] = resultWeightsHandle;
             reconstructFftTask->handles[5] = amountLoadedHandle;
             reconstructFftTask->cl_arg = &reconstructFftArg;
             reconstructFftTask->cl_arg_size = sizeof(ReconstructFftArgs);
             reconstructFftTask->cl_arg_free = 0;
             reconstructFftTask->callback_func = BatchProvider::batchCompleted;
             reconstructFftTask->callback_arg = &batches;
             reconstructFftTask->callback_arg_free = 0;
             reconstructFftTask->synchronous = DEBUG_SYNCHRONOUS_TASKS;
 
             CHECK_STARPU(starpu_task_submit(reconstructFftTask));
         }
 
         // Mark buffers shared between tasks of a batch for removal
         starpu_data_unregister_submit(fftHandle);
         starpu_data_unregister_submit(traverseSpacesHandle);
         starpu_data_unregister_submit(amountLoadedHandle);
     }
 
     // Mark helper buffers that are shared between all tasks for removal
     starpu_data_unregister_submit(blobTableSquaredHandle);
     starpu_data_unregister_submit(fftScratchMemoryHandle);
 
     // Complete all processing
     /*while (starpu_task_nsubmitted() > 0) {
         starpu_do_schedule();
         usleep(100000);//100ms
     }*/
     CHECK_STARPU(starpu_task_wait_for_all());
 
     // Release result buffers synchronously, implicitly causes them to be copied to original memory location
     starpu_data_unregister(resultVolumeHandle);
     starpu_data_unregister(resultWeightsHandle);
 
     padded_image_to_fft_cuda_deinitialize();
 
     return result;
 }

◆ defineParams()

void ProgRecFourierStarPU::defineParams ( )

overridevirtual

Specify supported command line arguments

Reimplemented from XmippProgram.

Reimplemented in ProgRecFourierMpiStarPU.

Definition at line 48 of file reconstruct_fourier_starpu.cpp.

                                         {
     //usage
     addUsageLine("Generate 3D reconstructions from projections using direct Fourier interpolation with arbitrary geometry.");
     addUsageLine("Kaisser-windows are used for interpolation in Fourier space.");
     //params
     addParamsLine("   -i <md_file>                  : Metadata file with input projections");
     addParamsLine("  [-o <volume_file=\"rec_fourier.vol\">] : Filename for output volume");
     addParamsLine("  [--sym <symfile=c1>]           : Enforce symmetry in projections");
     addParamsLine("  [--padding <proj=2.0> <vol=2.0>] : Padding used for projections and volume");
     addParamsLine("  [--max_resolution <p=0.5>]     : Max resolution (Nyquist=0.5)");
     addParamsLine("  [--weight]                     : Use weights stored in the image metadata");
     addParamsLine("  [--blob <radius=1.9> <order=0> <alpha=15>] : Blob parameters");
     addParamsLine("                                 : radius in pixels, order of Bessel function in blob and parameter alpha");
     addParamsLine("  [--fast]                       : Do the blobbing at the end of the computation.");
     addParamsLine("                                 : Gives slightly different results, but is faster.");
     addParamsLine("  [--batchSize <size=25>]        : Amount of images in single batch. Too many won't fit in memory and will not be able to run as fast, too few will have large overhead.");
     addParamsLine("                                 : Each additional image in batch will require at least imageSide^2 * 32 bytes, e.g.");
     addParamsLine("                                 : cca 52MB per batch of 25 256x256 images or 210MB for 512x512 images.");
     addParamsLine("  [--fourierThreads <num=0>]     : Number of threads used for final fourier transformation. Zero means all available.");
     addExampleLine("For reconstruct enforcing i3 symmetry and using stored weights:", false);
     addExampleLine("   xmipp_starpu_reconstruct_fourier  -i reconstruction.sel --sym i3 --weight");
 
     /* Buffer sizes:
      *  paddedImages:
      *      batchSize * paddedImgSize * paddedImgSize * sizeof(float)
      *      batchSize * symmetryCount * sizeof(TraverseSpace)       <- Typically in order of few kB per batch
      *  fft:
      *      batchSize * fftSizeX * fftSizeY * 2 * sizeof(float)
      *  fftWorkSpace:
      *      paddedImgSize * paddedImgSize/2 * sizeof(float) * 2     <- Does not scale with batchSize
      */
 }

◆ initStarPU()

void ProgRecFourierStarPU::initStarPU ( )

staticprotected

Initialize StarPU. Separate method to allow MPI override

Definition at line 290 of file reconstruct_fourier_starpu.cpp.

                                       {
     // Request more workers per CUDA-capable GPU
     setenv("STARPU_NWORKER_PER_CUDA",
            "2" /* seems to work best, 1 leaves GPU idle for fractions of a second (but shouldn't) */,
            false /* don't overwrite user specified value */);
     // Be more reluctant to discard calibration data
     setenv("STARPU_HISTORY_MAX_ERROR", "200" /*%*/, false);
     // Do not let users to specify scheduling policy, because we are using our own and StarPU REALLY doesn't like this conflict
     setenv("STARPU_SCHED", "", true);
 
     starpu_conf starpu_configuration;
     starpu_conf_init(&starpu_configuration);
     starpu_configuration.sched_policy = &schedulers.reconstruct_fourier;
     CHECK_STARPU(starpu_init(&starpu_configuration));
 
     int cpuWorkers[STARPU_NMAXWORKERS];
     unsigned cpuWorkerCount = starpu_worker_get_ids_by_type(starpu_worker_archtype::STARPU_CPU_WORKER, cpuWorkers, STARPU_NMAXWORKERS);
     if (cpuWorkerCount > 1) {
         int combinedWorker = starpu_combined_worker_assign_workerid(cpuWorkerCount, cpuWorkers);
         starpu_sched_ctx_add_workers(&combinedWorker, 1, 0);
 
         schedulers.reconstruct_fourier.add_workers(0, nullptr, 0);// Just force a refresh, not a very clean solution
     }
 
     starpu_malloc_set_align(ALIGNMENT);
 }

◆ postProcessAndSave()

void ProgRecFourierStarPU::postProcessAndSave	(	const Params &	params,
		const ComputeConstants &	computeConstants,
		const FileName &	fn_out,
		std::complex< float > ***	tempVolume,
		float ***	tempWeights
	)

staticprotected

Method will take data stored in tempVolume and tempWeights, crops it in X axis, calculates IFFT and stores the result to final destination.

Definition at line 550 of file reconstruct_fourier_starpu.cpp.

                                                                {
 
     const uint32_t maxVolumeIndex = computeConstants.maxVolumeIndex;
 
     // remove complex conjugate of the intermediate result
     const uint32_t maxVolumeIndexX = maxVolumeIndex/2;// just half of the space is necessary, the rest is complex conjugate
     const uint32_t maxVolumeIndexYZ = maxVolumeIndex;
     mirrorAndCropTempSpaces(tempVolume, tempWeights, maxVolumeIndexX, maxVolumeIndexYZ);
 
     if (params.fastLateBlobbing) {
         tempVolume = applyBlob(tempVolume, (float) params.blob.radius, (float*)computeConstants.blobTableSqrt, computeConstants.iDeltaSqrt, (int)maxVolumeIndexX, (int)maxVolumeIndexYZ);
         tempWeights = applyBlob(tempWeights, (float) params.blob.radius, (float*)computeConstants.blobTableSqrt, computeConstants.iDeltaSqrt, (int)maxVolumeIndexX, (int)maxVolumeIndexYZ);
     }
 
     const uint32_t paddedImgSize = computeConstants.paddedImgSize;
     const uint32_t imgSize = computeConstants.imgSize;
 
     forceHermitianSymmetry(tempVolume, tempWeights, maxVolumeIndexYZ);
     processWeights(tempVolume, tempWeights, maxVolumeIndexX, maxVolumeIndexYZ, params.padding_factor_proj, params.padding_factor_vol, imgSize);
     release(tempWeights, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1);
     MultidimArray< std::complex<double> > VoutFourier(paddedImgSize, paddedImgSize, paddedImgSize/2 + 1); // allocate space for output Fourier transformation
 
     convertToExpectedSpace(tempVolume, maxVolumeIndexYZ, VoutFourier);
     release(tempVolume, maxVolumeIndexYZ+1, maxVolumeIndexYZ+1);
 
     // Output volume
     Image<double> Vout(paddedImgSize, paddedImgSize, paddedImgSize);
     MultidimArray<double> &mVout = Vout.data;
 
     {
         FourierTransformer transformerVol;
         transformerVol.setThreadsNumber(params.fourierTransformThreads);
         transformerVol.fReal = &mVout;
         transformerVol.setFourierAlias(VoutFourier);
         transformerVol.recomputePlanR2C();
         transformerVol.inverseFourierTransform();
     }
 
     CenterFFT(mVout, false);
 
     // Correct by the Fourier transform of the blob
     mVout.setXmippOrigin();
     mVout.selfWindow(FIRST_XMIPP_INDEX(imgSize),FIRST_XMIPP_INDEX(imgSize),
                       FIRST_XMIPP_INDEX(imgSize),LAST_XMIPP_INDEX(imgSize),
                       LAST_XMIPP_INDEX(imgSize),LAST_XMIPP_INDEX(imgSize));
     double pad_relation = params.padding_factor_proj / params.padding_factor_vol;
     pad_relation = (pad_relation * pad_relation * pad_relation);
 
     double ipad_relation = 1.0 / pad_relation;
     double meanFactor2 = 0;
     FOR_ALL_ELEMENTS_IN_ARRAY3D(mVout) {
         double radius = sqrt((double)(k*k + i*i + j*j));
         double aux = radius * computeConstants.iDeltaFourier;
         double factor = computeConstants.fourierBlobTable[ROUND(aux)];
         double factor2 = pow(Sinc(radius / (2*imgSize)), 2);
         A3D_ELEM(mVout,k,i,j) /= (ipad_relation * factor2 * factor);
         meanFactor2 += factor2;
     }
     meanFactor2 /= MULTIDIM_SIZE(mVout);
     FOR_ALL_ELEMENTS_IN_ARRAY3D(mVout) {
         A3D_ELEM(mVout, k, i, j) *= meanFactor2;
     }
     Vout.write(fn_out);
 }

◆ prepareConstants()

void ProgRecFourierStarPU::prepareConstants	(	const Params &	params,
		const MetaData &	SF,
		const FileName &	fn_sym,
		ComputeConstants &	constants
	)

staticprotected

Load or compute constants needed for the computation.

Definition at line 122 of file reconstruct_fourier_starpu.cpp.

                                                                                                                                          {
     // Ask for memory for the output volume and its Fourier transform
     size_t imageSize;
     {
         size_t objId = SF.firstRowId();
         FileName fnImg;
         SF.getValue(MDL_IMAGE, fnImg, objId);
         Image<double> I;
         I.read(fnImg, HEADER);
 
         imageSize = I().xdim;
         if (imageSize != I().ydim)
             REPORT_ERROR(ERR_MULTIDIM_SIZE, "This algorithm only works for squared images");
     }
 
     constants.imgSize = static_cast<int>(imageSize);
     constants.paddedImgSize = static_cast<uint32_t>(imageSize * params.padding_factor_vol);
     {
         uint32_t conserveRows = (uint32_t) ceil(2.0f * constants.paddedImgSize * params.maxResolution);
         // Round up to nearest even number (i.e. divisible by 2)
         constants.maxVolumeIndex = conserveRows + (conserveRows & 1); // second term is 1 iff conserveRows is odd, 0 otherwise
     }
 
     // Build a table of blob values
     blobtype blobFourier = params.blob;
     blobtype blobNormalized = params.blob;
     blobFourier.radius /= params.padding_factor_vol * imageSize;
     blobNormalized.radius /= params.padding_factor_proj / params.padding_factor_vol;
 
     double deltaSqrt     = (params.blob.radius * params.blob.radius) / (BLOB_TABLE_SIZE_SQRT - 1);
     double deltaFourier  = (sqrt(3.0) * imageSize / 2.0) / (BLOB_TABLE_SIZE_SQRT-1);
     constants.iDeltaSqrt    = static_cast<float>(1 / deltaSqrt);
     constants.iDeltaFourier = static_cast<float>(1 / deltaFourier);
 
     // The interpolation kernel must integrate to 1
     constants.iw0 = 1.0 / blob_Fourier_val(0.0, blobNormalized);
     double padXdim3 = params.padding_factor_vol * imageSize;
     padXdim3 = padXdim3 * padXdim3 * padXdim3;
     double blobTableSize = params.blob.radius * sqrt(1.0 / (BLOB_TABLE_SIZE_SQRT-1));
     for (int i = 0; i < BLOB_TABLE_SIZE_SQRT; i++) {
         //use a r*r sample instead of r
         constants.blobTableSqrt[i] = static_cast<float>(blob_val(blobTableSize*sqrt((double)i), params.blob) * constants.iw0);
         constants.fourierBlobTable[i] = blob_Fourier_val(deltaFourier * i, blobFourier) * padXdim3 * constants.iw0;
     }
 
     {// Get symmetries
         Matrix2D<double> Identity(3, 3);
         Identity.initIdentity();
         constants.R_symmetries.push_back(Identity);
         if (!fn_sym.isEmpty()) {
             SymList SL;
             SL.readSymmetryFile(fn_sym);
             constants.R_symmetries.reserve(constants.R_symmetries.size() + SL.symsNo());
             for (int isym = 0; isym < SL.symsNo(); isym++) {
                 Matrix2D<double> L(4, 4), R(4, 4);
                 SL.getMatrices(isym, L, R);
                 R.resize(3, 3);
                 constants.R_symmetries.push_back(R);
             }
         }
     }
 }

◆ prepareMetaData()

void ProgRecFourierStarPU::prepareMetaData	(	const FileName &	fn_in,
		MetaDataVec &	SF
	)

staticprotected

Load meta data, selection file with a list of images to process.

Returns: total batch count

Definition at line 112 of file reconstruct_fourier_starpu.cpp.

                                                                                  {
     // Read the input images
     SF.read(fn_in);
     SF.removeDisabled();
 }

◆ readParams()

void ProgRecFourierStarPU::readParams ( )

overridevirtual

Read arguments from command line

Reimplemented from XmippProgram.

Reimplemented in ProgRecFourierMpiStarPU.

Definition at line 82 of file reconstruct_fourier_starpu.cpp.

                                       {
     fn_in = getParam("-i");
     fn_out = getParam("-o");
     fn_sym = getParam("--sym");
 
     params.padding_factor_proj = getDoubleParam("--padding", 0);
     params.padding_factor_vol = getDoubleParam("--padding", 1);
     params.maxResolution = (float)getDoubleParam("--max_resolution");
     params.do_weights = checkParam("--weight");
     params.blob.radius = getDoubleParam("--blob", 0);
     params.blob.order  = getIntParam("--blob", 1);
     params.blob.alpha  = getDoubleParam("--blob", 2);
     params.fastLateBlobbing = checkParam("--fast");
 
     int batchSize = getIntParam("--batchSize");
     params.batchSize = batchSize <= 0 ? 25 : static_cast<uint32_t>(batchSize);
 
     {// Number of threads of final fourier transformation
         int fourierThreads = getIntParam("--fourierThreads");
         if (fourierThreads <= 0) {
             fourierThreads = (int)sysconf(_SC_NPROCESSORS_ONLN); // Default to all
         }
         if (fourierThreads <= 0) {
             fourierThreads = 1;
             std::cerr << "--fourierThreads: Cannot obtain number of cores. Defaulting to one." << std::endl;
         }
         params.fourierTransformThreads = static_cast<uint32_t>(fourierThreads);
     }
 }

◆ run()

void ProgRecFourierStarPU::run ( )

overridevirtual

Run the image processing. Method will load data, process them and store result to final destination.

Reimplemented from XmippProgram.

Reimplemented in ProgRecFourierMpiStarPU.

Definition at line 264 of file reconstruct_fourier_starpu.cpp.

                                {
     if (verbose) {
         show();
     }
 
     prepareMetaData(fn_in, SF);
     prepareConstants(params, SF, fn_sym, computeConstants);
     initStarPU();
 
     SimpleBatchProvider batchProvider(computeBatchCount(params, SF), (bool) verbose);
 
     ComputeStarPUResult result = computeStarPU(params, SF, computeConstants, batchProvider, (bool) verbose);
 
     // We won't need StarPU anymore
     shutdownStarPU();
 
     // TODO(jp): This step seems unnecessary (later steps would have to be rewritten to work with flat arrays) (also in MPI version)
     // Convert flat volume and weight arrays into multidimensional arrays and destroy originals
     std::complex<float>*** tempVolume = result.createXmippStyleVolume(computeConstants.maxVolumeIndex);
     float*** tempWeights = result.createXmippStyleWeights(computeConstants.maxVolumeIndex);
     result.destroy();
 
     // Adjust and save the resulting volume
     postProcessAndSave(params, computeConstants, fn_out, tempVolume, tempWeights);
 }

◆ setIO()

void ProgRecFourierStarPU::setIO	(	const FileName &	fn_in,
		const FileName &	fn_out
	)

overridevirtual

Functions of common reconstruction interface

Implements ProgReconsBase.

Definition at line 185 of file reconstruct_fourier_starpu.cpp.

                                                                               {
     this->fn_in = fn_in;
     this->fn_out = fn_out;
 }

◆ show()

void ProgRecFourierStarPU::show ( ) const

overridevirtual

Show basic info to standard output

Reimplemented from XmippProgram.

Definition at line 191 of file reconstruct_fourier_starpu.cpp.

                                       {
     std::cout << " =====================================================================" << std::endl;
     std::cout << " Direct 3D reconstruction method using Kaiser windows as interpolators" << std::endl;
     std::cout << " =====================================================================" << std::endl;
     std::cout << " Input selfile             : "  << fn_in << std::endl;
     std::cout << " padding_factor_proj       : "  << params.padding_factor_proj << std::endl;
     std::cout << " padding_factor_vol        : "  << params.padding_factor_vol << std::endl;
     std::cout << " Output volume             : "  << fn_out << std::endl;
     if (!fn_sym.isEmpty())
         std::cout << " Symmetry file for projections : "  << fn_sym << std::endl;
     if (params.do_weights)
         std::cout << " Use weights stored in the image headers or doc file" << std::endl;
     else
         std::cout << " Do NOT use weights" << std::endl;
     std::cout << "\n Interpolation Function"
               << "\n   blrad                 : "  << params.blob.radius
               << "\n   blord                 : "  << params.blob.order
               << "\n   blalpha               : "  << params.blob.alpha
               << "\n max_resolution          : "  << params.maxResolution
               << "\n -----------------------------------------------------------------" << std::endl;
 }

◆ shutdownStarPU()

void ProgRecFourierStarPU::shutdownStarPU ( )

staticprotected

Shutdown StarPU. Separate method to allow MPI override

Definition at line 546 of file reconstruct_fourier_starpu.cpp.

                                           {
     starpu_shutdown();
 }

Member Data Documentation

◆ computeConstants

struct ProgRecFourierStarPU::ComputeConstants ProgRecFourierStarPU::computeConstants

protected

◆ fn_in

FileName ProgRecFourierStarPU::fn_in

protected

Input file name

Definition at line 48 of file reconstruct_fourier_starpu.h.

◆ fn_out

FileName ProgRecFourierStarPU::fn_out

protected

Output file name

Definition at line 51 of file reconstruct_fourier_starpu.h.

◆ fn_sym

FileName ProgRecFourierStarPU::fn_sym

protected

File with symmetries

Definition at line 54 of file reconstruct_fourier_starpu.h.

◆ params

struct ProgRecFourierStarPU::Params ProgRecFourierStarPU::params

protected

◆ SF

MetaDataVec ProgRecFourierStarPU::SF

protected

SelFile containing all projections

Definition at line 85 of file reconstruct_fourier_starpu.h.

The documentation for this class was generated from the following files:

xmipp/legacy/libraries/reconstruction_starpu/reconstruct_fourier_starpu.h
xmipp/legacy/libraries/reconstruction_starpu/reconstruct_fourier_starpu.cpp

Classes

Public Member Functions

Static Protected Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ ProgRecFourierStarPU()

Member Function Documentation

◆ computeBatchCount()

◆ computeStarPU()

◆ defineParams()

◆ initStarPU()

◆ postProcessAndSave()

◆ prepareConstants()

◆ prepareMetaData()

◆ readParams()

◆ run()

◆ setIO()

◆ show()

◆ shutdownStarPU()

Member Data Documentation

◆ computeConstants

◆ fn_in

◆ fn_out

◆ fn_sym

◆ params

◆ SF