#include <adjust_volume_grey_levels.h>

Inheritance diagram for ProgAdjustVolume:

Collaboration diagram for ProgAdjustVolume:

Public Member Functions
void	run ()

double	mismatching (double a, double b)

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	show () 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 ()

Public Attributes
Image< double >	ImIn

MultidimArray< double >	V

MetaDataVec	SF

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Protected Member Functions
void	defineParams ()

void	readParams ()

void	show ()

void	apply (MultidimArray< float > &output_volume)

Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Attributes
FileName	fn_vol
	Filename with the input volume. More...

FileName	fn_sel
	Filename with the input projections. More...

FileName	fn_out

bool	optimize
	Optimize. More...

double	probb_eval
	Probability of being evaluated. More...

bool	tempFile

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

Detailed Description

Definition at line 37 of file adjust_volume_grey_levels.h.

Member Function Documentation

◆ apply()

void ProgAdjustVolume::apply ( MultidimArray< float > & output_volume )

protected

Apply. This is the function that really does the job

Definition at line 169 of file adjust_volume_grey_levels.cpp.

 {
     // Compute the average power and average value of all the projections
     double sum = 0, sum2 = 0, N = 0;
     std::cerr << "Computing first estimate of the linear transformation ...\n";
     int imgno = SF.size();
     init_progress_bar(imgno);
     int i = 0;
 
     int projXdim, projYdim;
     if(SF.isEmpty())
     {
         std::cerr << "Empty inputFile File\n";
         exit(1);
     }
     Image<double> I;
     for (size_t objId : SF.ids())
     {
         // Read image
         I.readApplyGeo(SF, objId);
         projXdim = XSIZE(I());
         projYdim = YSIZE(I());
 
         // Compute the image statistics
         double avg=0., stddev=0., min, max;
         I().computeStats(avg, stddev, min, max);
         double Ni = projXdim * projYdim;
         sum += avg;
         sum2 += stddev * stddev;
         N += Ni;
 
         // End of loop
         i++;
         if (i % 10 == 0)
             progress_bar(i);
     }
     progress_bar(imgno);
     std::cout << std::endl;
 
     // Statistics of the volume
     double avg0=0., stddev0=0., min0, max0;
     V.computeStats(avg0, stddev0, min0, max0);
 
     // First guess of the transformation parameters a*(x-vm)+b
     // r is the average length of a ray
     // The equations to solve are
     // v_i follows currently N(avg0,stddev0)
     // we want it to follow N(avgF,stddevF)
     // such that sum_{i=1}^r{v_i} follows N(avg_pict,stddev_pict)
     // On the other hand we know that
     //     sum_{i=1}^r{v_i} follows N(r avgF, sqrt(r)*stddevF)
     // Thus,
     //     avgF=avg_pict/r
     //     stddevF=stddev_pict/sqrt(r)
     // the equations of the transformation are
     //    y=ax+b
     //    avg_y     =a*avg_x+b=avgF
     //    stddev_y  =a*stddev_y=stddevF
     // Therefore
     //    a=stddevF/stddev0
     //    b=avgF-a*avg0
     double r = pow(MULTIDIM_SIZE(V), 1.0 / 3.0);
     double avg_pict = sum / imgno;
     double stddev_pict = sqrt(sum2 / imgno);
     double avgF = avg_pict / r;
     double stddevF = stddev_pict / sqrt(r);
     double a = stddevF / stddev0;
     double b = avgF - a * avg0;
     std::cout << "First Linear transformation: y=" << a << "*x+"  << b << std::endl;
 
     // Optimize
     if (optimize)
     {
         Matrix1D<double> p(2), steps(2);
         p(0) = a;
         p(1) = b;
         steps.initConstant(1);
         double ftol = 0.01, fret;
         int iter;
         globalAdjustVolumeProg = this;
         powellOptimizer(p, 1, 2, &projectionMismatching, nullptr,
                         ftol, fret, iter, steps, true);
         a = p(0);
         b = p(1);
     }
 
     // Apply the transformation
 //    out = V;
     FOR_ALL_ELEMENTS_IN_ARRAY3D(V) out(k, i, j) = a * V(k, i, j) + b;
 }

◆ defineParams()

void ProgAdjustVolume::defineParams ( )

protectedvirtual

Function in which the param of each Program are defined.

Reimplemented from XmippProgram.

Definition at line 31 of file adjust_volume_grey_levels.cpp.

 {
     // Usage
     addUsageLine("Search for a linear transformation of the gray values of a volume so that the");
     addUsageLine("error between the theoretically projected images and the experimental images");
     addUsageLine("is minimized. This program must be used before computing the Volumetric SSNR");
     addUsageLine("if the reconstruction algorithm scales the output volume differently.");
     // See also
 
     // Examples
     addExampleLine("Adjust a volume to a set of images:", false);
     addExampleLine("xmipp_adjust_volume_grey_levels -i input_volume.vol -m experimental.sel -o output_volume.vol ");
 
     // Parameters
     addParamsLine(" -i <volume_file>      : Volume to adjust its range.");
     addParamsLine(" alias --input;");
     addParamsLine(" -m  <metadata_file>   : Set of projections of the volume.");
     addParamsLine(" alias --metadata;");
     addParamsLine(" [-o <volume_file=\"\">]    : Output adjusted volume. By default, the input one.");
     addParamsLine(" alias --output;");
     addParamsLine(" [--optimize]          : Optimize the linear transformation. By default, keep the initially computed.");
     addParamsLine(" [--probb_eval <p=0.2>] : The goal function is evaluated each time from a random set of projections.");
     addParamsLine("                       : Each image has a probability of p of being evaluated.");
 }

◆ mismatching()

double ProgAdjustVolume::mismatching	(	double	a,
		double	b
	)

Mismatching. This function returns the overall mismatiching between the experimental projections and the theoretical projections of the current volume.

Definition at line 103 of file adjust_volume_grey_levels.cpp.

 {
     if (a <= 0)
         return 1e38;
 
     // Transform the volume
     MultidimArray<double> aux = V;
     FOR_ALL_ELEMENTS_IN_ARRAY3D(aux) aux(k, i, j) = a * aux(k, i, j) + b;
 
     // Compute the mismatching
     double retval = 0;
     //SF.go_first_ACTIVE();
     int N = 0;
     //while (!SF.eof())
     if(SF.isEmpty())
     {
         std::cerr << "Empty inputFile File\n";
         exit(1);
     }
     Image<double> I;
     ApplyGeoParams params;
     params.only_apply_shifts=true;
     for (size_t objId : SF.ids())
     {
         // Skip randomly some images
         double x = rnd_unif(0, 1);
         if (x > probb_eval)
             continue;
         N++;
 
         I.readApplyGeo(SF, objId, params);
         I().setXmippOrigin();
 
         // Project the auxiliary volume in the same direction
         Projection P;
         projectVolume(aux, P, YSIZE(I()), XSIZE(I()),
                       I.rot(), I.tilt(), I.psi());
 
         // Compute the difference
         MultidimArray<double> diff;
         diff = I() - P();
         retval += diff.sum2();
 
 #ifdef DEBUG
 
         I.write("PPPexp.xmp");
         I() = P();
         I.write("PPPtheo.xmp");
         std::cout << "Difference=" << diff.sum2() << std::endl;
         std::cout << "Press any key\n";
         char c;
         std::cin >> c;
 #endif
 
     }
     //TODO Review this
     //    while (SF.nextObject()!= NO_MORE_OBJECTS)
     //        ;
     if (0 == N) {
         REPORT_ERROR(ERR_NUMERICAL, "N is zero (0), which would lead to division by zero");
     }
     return retval / N;
 }

◆ readParams()

void ProgAdjustVolume::readParams ( )

protectedvirtual

Function in which each program will read parameters that it need. If some error occurs the usage will be printed out.

Reimplemented from XmippProgram.

Definition at line 56 of file adjust_volume_grey_levels.cpp.

 {
     fn_vol = getParam("-i");
     fn_sel = getParam("-m");
     fn_out = (checkParam("-o"))? getParam("-o") : fn_vol;
     optimize = checkParam("--optimize");
     probb_eval = getDoubleParam("--probb_eval");
     verbose = (checkParam("-v"))? getIntParam("-v"): false;
     tempFile = (STR_EQUAL(fn_vol.c_str(), fn_out.c_str()));
 }

◆ run()

void ProgAdjustVolume::run ( )

virtual

Run. Calls apply and save the result.

Reimplemented from XmippProgram.

Definition at line 67 of file adjust_volume_grey_levels.cpp.

 {
     if (verbose)
         show();
 
     // Read input volume
     ImIn.read(fn_vol);
     V.alias(ImIn());
     V.setXmippOrigin();
 
     int Xdim, Ydim, Zdim;
     Xdim = XSIZE(V);
     Ydim = YSIZE(V);
     Zdim = ZSIZE(V);
 
     Image<float> ImOut;
     ImOut.mapFile2Write(Xdim, Ydim, Zdim, fn_out, tempFile);
     ImOut().setXmippOrigin();
 
     // Read input metadataFile
     SF.read(fn_sel,nullptr);
 
     apply(ImOut());
 
     ImIn.clear();
     ImOut.write();
 }

◆ show()

void ProgAdjustVolume::show ( )

protected

Show parameters.

Definition at line 261 of file adjust_volume_grey_levels.cpp.

 {
     std::cout << "Input Volume:  " << fn_vol   << std::endl
     << "Input MetaDAtaFile: " << fn_sel   << std::endl
     << "Output Volume: " << fn_out   << std::endl
     << "Optimize:      " << optimize << std::endl;
 }