#include <volume_to_pseudoatoms.h>

Inheritance diagram for ProgVolumeToPseudoatoms:

Collaboration diagram for ProgVolumeToPseudoatoms:

Public Member Functions
void	readParams ()
	Read parameters from command line. More...

void	show () const
	show parameters More...

void	defineParams ()
	define parameters More...

void	produceSideInfo ()
	Prepare side info. More...

void	run ()
	Run. More...

void	placeSeeds (int Nseeds)
	Place seeds. More...

void	removeSeeds (int Nseeds)
	Remove seeds. More...

void	removeTooCloseSeeds ()
	Remove too close seeds. More...

double	computeAverage (int k, int i, int j, MultidimArray< double > &V)
	Compute average of a volume. More...

void	drawGaussian (double k, double i, double j, MultidimArray< double > &V, double intensity)
	Draw a Gaussian on a volume. More...

void	drawApproximation ()
	Draw approximation. More...

void	extractRegion (int idxGaussian, MultidimArray< double > &region, bool extended=false) const
	Extract region around a Gaussian. More...

void	insertRegion (const MultidimArray< double > &region)
	Insert region. More...

double	evaluateRegion (const MultidimArray< double > &region) const
	Evaluate region. More...

void	optimizeCurrentAtoms ()
	Optimize current atoms. More...

void	writeResults ()
	Write results. More...

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 Public Member Functions
static void *	optimizeCurrentAtomsThread (void *threadArgs)
	Optimize current atoms (thread) More...

Public Attributes
FileName	fnVol
	Volume to convert. More...

FileName	fnOut
	Output volume. More...

Mask	mask_prm

bool	useMask

double	sigma
	Sigma. More...

double	targetError
	Maximum error (as a percentage) More...

double	stop
	Stop criterion. More...

int	initialSeeds
	Initial seeds. More...

double	growSeeds
	Grow seeds. More...

bool	allowMovement
	Allow gaussians to move. More...

bool	allowIntensity
	Allow gaussians to vary intensity. More...

double	intensityFraction

std::string	intensityColumn
	Column for the intensity (if any) More...

double	minDistance
	Mindistance. More...

double	penalty
	Penalization. More...

int	numThreads
	Number of threads. More...

double	sampling
	Sampling rate. More...

size_t	Nclosest
	N closest atoms for the distance histogram. More...

bool	dontScale
	Don't scale the atom weights at the end. More...

bool	binarize
	Binarize. More...

double	threshold
	Threshold for the binarization. More...

Image< double >	Vin

Image< double >	Vcurrent

double	energyDiff

double	percentageDiff

double	energyOriginal

std::vector< PseudoAtom >	atoms

double	sigma3

double	percentil1

double	range

double	smallAtom

MultidimArray< double >	gaussianTable

barrier_t	barrier

int	threadOpCode

Prog_Convert_Vol2Pseudo_ThreadParams *	threadArgs

pthread_t *	threadIds

FourierFilter	Filter

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Additional Inherited Members
Protected Member Functions inherited from XmippProgram
void	defineCommons ()

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 72 of file volume_to_pseudoatoms.h.

Member Function Documentation

◆ computeAverage()

double ProgVolumeToPseudoatoms::computeAverage	(	int	k,
		int	i,
		int	j,
		MultidimArray< double > &	V
	)

Compute average of a volume.

Definition at line 586 of file volume_to_pseudoatoms.cpp.

 {
     int k0=std::max(STARTINGZ(V),(int)floor(k-sigma3));
     int i0=std::max(STARTINGY(V),(int)floor(i-sigma3));
     int j0=std::max(STARTINGX(V),(int)floor(j-sigma3));
     int kF=std::min(FINISHINGZ(V),(int)ceil(k+sigma3));
     int iF=std::min(FINISHINGY(V),(int)ceil(i+sigma3));
     int jF=std::min(FINISHINGX(V),(int)ceil(j+sigma3));
     double sum=0;
     for (int kk=k0; kk<=kF; kk++)
         for (int ii=i0; ii<=iF; ii++)
             for (int jj=j0; jj<=jF; jj++)
                 sum+=V(kk,ii,jj);
     return sum/((kF-k0+1)*(iF-i0+1)*(jF-j0+1));
 }

◆ defineParams()

void ProgVolumeToPseudoatoms::defineParams ( )

virtual

define parameters

Reimplemented from XmippProgram.

Definition at line 111 of file volume_to_pseudoatoms.cpp.

 {
     addUsageLine("Creates a set of pseudoatoms representing the density of an EM volume. ");
     addUsageLine("+This is useful for the vector quantization process needed in problems ");
     addUsageLine("+like docking, approximation of structures by Alpha Shapes, Normal Mode ");
     addUsageLine("+Analysis, etc.");
     addUsageLine("+");
     addUsageLine("+The volume is approximated by Gaussians of a desired size. The user can ");
     addUsageLine("+specify whether the Gaussians can have different intensities or not, as well ");
     addUsageLine("+as the level of precision with which she desires to approximate the input ");
     addUsageLine("+EM volume.");
     addParamsLine("   -i <volume>                       : Input volume");
     addParamsLine("  [-o <rootname=\"\">]               : Output rootname. If not given, the rootname of the input volume is taken.");
     addParamsLine("                                     : The output of the program is: [rootname].pdb");
     addParamsLine("                                     : (PDB file with the pseudo atoms).");
     addParamsLine("                                     :+If verbose is set to 2, then also [rootname].vol ");
     addParamsLine("                                     :+(the approximation volume), [rootname].hist ");
     addParamsLine("                                     :+(histogram of the Gaussian intensities), ");
     addParamsLine("                                     :+[rootname]_rawDiff.vol (difference between the ");
     addParamsLine("                                     :+input volume and its approximation), ");
     addParamsLine("                                     :+[rootname]_relativeDiff.vol (the raw difference ");
     addParamsLine("                                     :+divided by the input volume at that location; this ");
     addParamsLine("                                     :+gives an idea of how much the error represents with ");
     addParamsLine("                                     :+respect to the input");
     addParamsLine("  [--sigma <s=1.5>]                  : Sigma of gaussians (in Angstroms)");
     addParamsLine("                                     : It should be comparable to the sampling rate");
     addParamsLine("  [--initialSeeds+ <N=300>]          : Initial number of pseudoatoms");
     addParamsLine("  [--growSeeds+ <percentage=30>]     : Percentage of growth");
     addParamsLine("                                     :+At each iteration the smallest percentage/2 ");
     addParamsLine("                                     :+pseudoatoms will be removed, and percentage new pseudoatoms will be created.");
     addParamsLine("  [--stop+ <p=0.001>]                : Stop criterion (0<p<1) for inner iterations");
     addParamsLine("                                     :+At each iteration the current number of gaussians will be optimized until ");
     addParamsLine("                                     :+the average error does not decrease at least this amount relative to the previous iteration.");
     addParamsLine("  [--targetError+ <e=2>]             : Finish when the average representation");
     addParamsLine("                                     : error is below this threshold (in percentage; by default, 2%)");
     addParamsLine("  [--dontAllowMovement+]             : Don't allow pseudoatoms to move");
     addParamsLine("  [--dontAllowIntensity+ <f=0.01>]   : Don't allow pseudoatoms to change intensity");
     addParamsLine("                                     : f determines the fraction of intensity");
     addParamsLine("                                     : held by each pseudoatom");
     addParamsLine("  [--intensityColumn+ <s=Bfactor>]   : Where to write the intensity in the PDB file");
     addParamsLine("                   where <s>");
     addParamsLine("                         occupancy");
     addParamsLine("                         Bfactor");
     addParamsLine("  [--Nclosest+ <N=3>]                : N closest atoms, it is used only for the");
     addParamsLine("                                     : distance histogram");
     addParamsLine("  [--minDistance+ <d=0.001>]         : Minimum distance between two pseudoatoms (in Angstroms)");
     addParamsLine("                                     : Set it to -1 to disable");
     addParamsLine("  [--penalty+ <p=10>]                : Penalty for overshooting");
     addParamsLine("  [--sampling_rate <Ts=1>]           : Sampling rate Angstroms/pixel");
     addParamsLine("  [--dontScale+]                     : Don't scale atom weights in the PDB");
     addParamsLine("  [--binarize+ <threshold>]          : Binarize the volume for a more uniform distribution");
     addParamsLine("  [--thr <n=1>]                      : Number of threads");
     mask_prm.defineParams(this,INT_MASK,nullptr,"Statistics restricted to the mask area.",true);
     addExampleLine("Convert volume to pseudoatoms, each pseudoatom with a different weight",false);
     addExampleLine("xmipp_volume_to_pseudoatoms -i volume.vol -o pseudoatoms");
     addExampleLine("Convert volume to pseudoatoms, all pseudoatoms with the same weight",false);
     addExampleLine("xmipp_volume_to_pseudoatoms -i volume.vol -o pseudoatoms --dontAllowIntensity");
 }

◆ drawApproximation()

void ProgVolumeToPseudoatoms::drawApproximation ( )

Draw approximation.

Definition at line 554 of file volume_to_pseudoatoms.cpp.

 {
     Vcurrent().initZeros(Vin());
     int nmax=atoms.size();
     for (int n=0; n<nmax; n++)
         drawGaussian(atoms[n].location(0),atoms[n].location(1),
                      atoms[n].location(2),Vcurrent(),atoms[n].intensity);
 
     energyDiff=0;
     double N=0;
     percentageDiff=0;
     const MultidimArray<int> &iMask3D=mask_prm.get_binary_mask();
     const MultidimArray<double> &mVcurrent=Vcurrent();
     const MultidimArray<double> &mVin=Vin();
     FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(mVcurrent)
     {
         if (useMask && DIRECT_MULTIDIM_ELEM(iMask3D,n)==0)
             continue;
         double Vinv=DIRECT_MULTIDIM_ELEM(mVin,n);
         if (Vinv<=0)
             continue;
         double vdiff=Vinv-DIRECT_MULTIDIM_ELEM(mVcurrent,n);
         double vperc=fabs(vdiff);
         energyDiff+=vdiff*vdiff;
         percentageDiff+=vperc;
         N++;
     }
     energyDiff/=N;
     percentageDiff/=(N*range);
 }

◆ drawGaussian()

void ProgVolumeToPseudoatoms::drawGaussian	(	double	k,
		double	i,
		double	j,
		MultidimArray< double > &	V,
		double	intensity
	)

Draw a Gaussian on a volume.

Definition at line 603 of file volume_to_pseudoatoms.cpp.

 {
     int k0=CEIL(XMIPP_MAX(STARTINGZ(V),k-sigma3));
     int i0=CEIL(XMIPP_MAX(STARTINGY(V),i-sigma3));
     int j0=CEIL(XMIPP_MAX(STARTINGX(V),j-sigma3));
     int kF=FLOOR(XMIPP_MIN(FINISHINGZ(V),k+sigma3));
     int iF=FLOOR(XMIPP_MIN(FINISHINGY(V),i+sigma3));
     int jF=FLOOR(XMIPP_MIN(FINISHINGX(V),j+sigma3));
     for (int kk=k0; kk<=kF; kk++)
     {
         double aux=kk-k;
         double diffkk2=aux*aux;
         for (int ii=i0; ii<=iF; ii++)
         {
             aux=ii-i;
             double diffiikk2=aux*aux+diffkk2;
             for (int jj=j0; jj<=jF; jj++)
             {
                 aux=jj-j;
                 double r=sqrt(diffiikk2+aux*aux);
                 aux=r*1000;
                 long iaux=lround(aux);
                 A3D_ELEM(V,kk,ii,jj)+=intensity*DIRECT_A1D_ELEM(gaussianTable,iaux);
             }
         }
     }
 }

◆ evaluateRegion()

double ProgVolumeToPseudoatoms::evaluateRegion ( const MultidimArray< double > & region ) const

Evaluate region.

Definition at line 662 of file volume_to_pseudoatoms.cpp.

 {
     double avgDiff=0;
     double N=0;
     const MultidimArray<int> &iMask3D=mask_prm.get_binary_mask();
     const MultidimArray<double> &mVin=Vin();
 #ifdef DEBUG
     Image<double> save, save2;
     save().initZeros(region);
     save2().initZeros(region);
 #endif
     FOR_ALL_ELEMENTS_IN_ARRAY3D(region)
     {
         double Vinv=A3D_ELEM(mVin,k,i,j);
         if (Vinv<=0 || (useMask && A3D_ELEM(iMask3D,k,i,j)==0))
             continue;
 #ifdef DEBUG
         save(k,i,j)=Vinv;
         save2(k,i,j)=A3D_ELEM(region,k,i,j);
 #endif
         double vdiff=A3D_ELEM(region,k,i,j)-Vinv;
         double vperc=(vdiff<0)?-vdiff:penalty*vdiff;
         avgDiff+=vperc;
 #ifdef DEBUG
         std::cout << "(k,i,j)=(" << k << "," << i << "," << j << ") toSimulate=" << Vinv << " simulated=" << A3D_ELEM(region,k,i,j) << " vdiff=" << vdiff << " vperc=" << vperc << std::endl;
 #endif
         ++N;
     }
 #ifdef DEBUG
     save.write("PPPtoSimulate.vol");
     save2.write("PPPsimulated.vol");
     std::cout << "Error=" << avgDiff/(N*range) << std::endl;
     std::cout << "Press any key\n";
     char c; std::cin >> c;
 #endif
     return avgDiff/(N*range);
 }

◆ extractRegion()

void ProgVolumeToPseudoatoms::extractRegion	(	int	idxGaussian,
		MultidimArray< double > &	region,
		bool	extended = `false`
	)		const

Extract region around a Gaussian.

Definition at line 632 of file volume_to_pseudoatoms.cpp.

 {
     double k=atoms[idxGaussian].location(0);
     double i=atoms[idxGaussian].location(1);
     double j=atoms[idxGaussian].location(2);
 
     double sigma3ToUse=sigma3;
     if (extended)
         sigma3ToUse+=1.5;
 
     int k0=CEIL(XMIPP_MAX(STARTINGZ(Vcurrent()),k-sigma3ToUse));
     int i0=CEIL(XMIPP_MAX(STARTINGY(Vcurrent()),i-sigma3ToUse));
     int j0=CEIL(XMIPP_MAX(STARTINGX(Vcurrent()),j-sigma3ToUse));
     int kF=FLOOR(XMIPP_MIN(FINISHINGZ(Vcurrent()),k+sigma3ToUse));
     int iF=FLOOR(XMIPP_MIN(FINISHINGY(Vcurrent()),i+sigma3ToUse));
     int jF=FLOOR(XMIPP_MIN(FINISHINGX(Vcurrent()),j+sigma3ToUse));
 
     region.resizeNoCopy(kF-k0+1,iF-i0+1,jF-j0+1);
     STARTINGZ(region)=k0;
     STARTINGY(region)=i0;
     STARTINGX(region)=j0;
     const MultidimArray<double> &mVcurrent=Vcurrent();
     for (int k=k0; k<=kF; k++)
         for (int i=i0; i<=iF; i++)
             for (int j=j0; j<=jF; j++)
                 A3D_ELEM(region,k,i,j)=A3D_ELEM(mVcurrent,k,i,j);
 }

◆ insertRegion()

void ProgVolumeToPseudoatoms::insertRegion ( const MultidimArray< double > & region )

Insert region.

Definition at line 702 of file volume_to_pseudoatoms.cpp.

 {
     const MultidimArray<double> &mVcurrent=Vcurrent();
     FOR_ALL_ELEMENTS_IN_ARRAY3D(region)
     A3D_ELEM(mVcurrent,k,i,j)=A3D_ELEM(region,k,i,j);
 }

◆ optimizeCurrentAtoms()

void ProgVolumeToPseudoatoms::optimizeCurrentAtoms ( )

Optimize current atoms.

Definition at line 837 of file volume_to_pseudoatoms.cpp.

 {
     if (!allowIntensity && !allowMovement)
         return;
     bool finished=false;
     int iter=0;
     do
     {
         double oldError=percentageDiff;
 
         threadOpCode=WORKTHREAD;
         // Launch workers
         barrier_wait(&barrier);
         // Wait for workers to finish
         barrier_wait(&barrier);
 
         // Retrieve results
         int Nintensity=0;
         int Nmovement=0;
         for (int i=0; i<numThreads; i++)
         {
             Nintensity+=threadArgs[i].Nintensity;
             Nmovement+=threadArgs[i].Nmovement;
         }
 
         // Remove all the removed atoms
         int nmax=atoms.size();
         for (int n=nmax-1; n>=0; n--)
             if (atoms[n].intensity==0)
                 atoms.erase(atoms.begin()+n);
 
         drawApproximation();
         if (verbose>0)
             std::cout << "Iteration " << iter << " error= " << percentageDiff
             << " Natoms= " << atoms.size()
             << " Intensity= " << Nintensity
             << " Location= " << Nmovement
             << std::endl;
 
         if (iter>0)
             if ((oldError-percentageDiff)/oldError<stop)
                 finished=true;
         iter++;
     }
     while (!finished);
 }

◆ optimizeCurrentAtomsThread()

void * ProgVolumeToPseudoatoms::optimizeCurrentAtomsThread ( void * threadArgs )

static

Optimize current atoms (thread)

Definition at line 713 of file volume_to_pseudoatoms.cpp.

 {
     auto *myArgs=(Prog_Convert_Vol2Pseudo_ThreadParams *) threadArgs;
     ProgVolumeToPseudoatoms *parent=myArgs->parent;
     std::vector< PseudoAtom > &atoms=parent->atoms;
     bool allowIntensity=parent->allowIntensity;
     bool allowMovement=parent->allowMovement;
     MultidimArray<double> region, regionBackup;
 
     barrier_t *barrier=&(parent->barrier);
     do
     {
         barrier_wait( barrier );
         if (parent->threadOpCode==KILLTHREAD)
             return nullptr;
 
         myArgs->Nintensity=0;
         myArgs->Nmovement=0;
         int nmax=atoms.size();
         for (int n=0; n<nmax; n++)
         {
             if ((n+1)%parent->numThreads!=myArgs->myThreadID)
                 continue;
 
             parent->extractRegion(n,region,true);
             double currentRegionEval=parent->evaluateRegion(region);
             parent->drawGaussian(atoms[n].location(0), atoms[n].location(1),
                                  atoms[n].location(2),region,-atoms[n].intensity);
             regionBackup=region;
 
 #ifdef DEBUG
             std::cout << "Atom n=" << n << " current intensity=" << atoms[n].intensity << " -> " << currentRegionEval << std::endl;
 #endif
             // Change intensity
             if (allowIntensity)
             {
                 // Try with a Gaussian that is of different intensity
                 double tryCoeffs[8]={0, 0.1, 0.2, 0.5, 0.9, 0.99, 1.01, 1.1};
                 double bestRed=0;
                 int bestT=-1;
                 for (int t=0; t<8; t++)
                 {
                     region=regionBackup;
                     parent->drawGaussian(atoms[n].location(0),
                                          atoms[n].location(1), atoms[n].location(2),region,
                                          tryCoeffs[t]*atoms[n].intensity);
                     double trialRegionEval=parent->evaluateRegion(region);
                     double reduction=trialRegionEval-currentRegionEval;
                     if (reduction<bestRed)
                     {
                         bestRed=reduction;
                         bestT=t;
 #ifdef DEBUG
                         std::cout << "    better -> " << trialRegionEval << " (factor=" << tryCoeffs[t]  << ")" << std::endl;
 #endif
                     }
                 }
                 if (bestT!=-1)
                 {
                     atoms[n].intensity*=tryCoeffs[bestT];
                     region=regionBackup;
                     parent->drawGaussian(atoms[n].location(0), atoms[n].location(1),
                                          atoms[n].location(2),region,atoms[n].intensity);
                     pthread_mutex_lock(&mutexUpdateVolume);
                     parent->insertRegion(region);
                     pthread_mutex_unlock(&mutexUpdateVolume);
                     currentRegionEval=parent->evaluateRegion(region);
                     parent->drawGaussian(atoms[n].location(0),
                                          atoms[n].location(1), atoms[n].location(2),region,
                                          -atoms[n].intensity);
                     regionBackup=region;
 #ifdef DEBUG
                     std::cout << "    finally -> " << currentRegionEval << " (intensity=" << atoms[n].intensity  << ")" << std::endl;
 #endif
                     myArgs->Nintensity++;
                 }
             }
 
             // Change location
             if (allowMovement && atoms[n].intensity>0)
             {
                 double tryX[6]={-0.45,0.5, 0.0 ,0.0, 0.0 ,0.0};
                 double tryY[6]={ 0.0 ,0.0,-0.45,0.5, 0.0 ,0.0};
                 double tryZ[6]={ 0.0 ,0.0, 0.0 ,0.0,-0.45,0.5};
                 double bestRed=0;
                 int bestT=-1;
                 for (int t=0; t<6; t++)
                 {
                     region=regionBackup;
                     parent->drawGaussian(atoms[n].location(0)+tryZ[t],
                                          atoms[n].location(1)+tryY[t],
                                          atoms[n].location(2)+tryX[t],
                                          region,atoms[n].intensity);
                     double trialRegionEval=parent->evaluateRegion(region);
                     double reduction=trialRegionEval-currentRegionEval;
                     if (reduction<bestRed)
                     {
                         bestRed=reduction;
                         bestT=t;
                     }
                 }
                 if (bestT!=-1)
                 {
                     atoms[n].location(0)+=tryZ[bestT];
                     atoms[n].location(1)+=tryY[bestT];
                     atoms[n].location(2)+=tryX[bestT];
                     region=regionBackup;
                     parent->drawGaussian(atoms[n].location(0),
                                          atoms[n].location(1), atoms[n].location(2),region,
                                          atoms[n].intensity);
                     pthread_mutex_lock(&mutexUpdateVolume);
                     parent->insertRegion(region);
                     pthread_mutex_unlock(&mutexUpdateVolume);
                     myArgs->Nmovement++;
                 }
             }
         }
 
         barrier_wait( barrier );
     }
     while (true);
 }

◆ placeSeeds()

void ProgVolumeToPseudoatoms::placeSeeds ( int Nseeds )

Place seeds.

Definition at line 329 of file volume_to_pseudoatoms.cpp.

 {
     // Convolve the difference with the Gaussian to know
     // where it would be better to put a Gaussian
     MultidimArray<double> Vdiff=Vin();
     Vdiff-=Vcurrent();
     Filter.applyMaskSpace(Vdiff);
 
     // Look for the Nseeds
     const MultidimArray<int> &iMask3D=mask_prm.get_binary_mask();
     size_t idx=0;
     std::list<SeedCandidate> candidateList;
     size_t listSize=0;
     double lastValue=0;
     FOR_ALL_ELEMENTS_IN_ARRAY3D(Vdiff)
     {
         if (!useMask || useMask && DIRECT_MULTIDIM_ELEM(iMask3D,idx))
         {
             double v=A3D_ELEM(Vdiff,k,i,j);
             if (listSize==0 || v>lastValue)
             {
                 // Check if there is any other candidate around
                 bool found=false;
                 for (std::list<SeedCandidate>::iterator iter=candidateList.begin(); iter!=candidateList.end(); iter++)
                     if (std::abs(iter->k-k)<sigma && std::abs(iter->i-i)<sigma && std::abs(iter->j-j)<sigma)
                     {
                         found=true;
                         break;
                     }
                 if (!found)
                 {
                     SeedCandidate aux;
                     aux.v=v;
                     aux.k=k;
                     aux.i=i;
                     aux.j=j;
                     candidateList.push_back(aux);
                     candidateList.sort();
                     listSize++;
                     if (listSize>Nseeds)
                     {
                         candidateList.pop_back();
                         listSize--;
                     }
                     lastValue=candidateList.back().v;
                 }
             }
         }
         idx++;
     }
 
     // Place atoms
     for (std::list<SeedCandidate>::iterator iter=candidateList.begin(); iter!=candidateList.end(); iter++)
     {
         PseudoAtom a;
         VEC_ELEM(a.location,0)=iter->k;
         VEC_ELEM(a.location,1)=iter->i;
         VEC_ELEM(a.location,2)=iter->j;
         if (allowIntensity)
             a.intensity=iter->v;
         else
         {
             if (iter->v<smallAtom)
                 break;
             a.intensity=smallAtom;
         }
         atoms.push_back(a);
 
         // Remove this density from the difference
         drawGaussian(iter->k,iter->i,iter->j,Vdiff,-a.intensity);
     }
 }

◆ produceSideInfo()

void ProgVolumeToPseudoatoms::produceSideInfo ( )

Prepare side info.

Definition at line 170 of file volume_to_pseudoatoms.cpp.

 {
     sigma/=sampling;
     minDistance/=sampling;
 
     if (intensityColumn!="occupancy" && intensityColumn!="Bfactor")
         REPORT_ERROR(ERR_VALUE_INCORRECT,(std::string)"Unknown column: "+intensityColumn);
 
     Vin.read(fnVol);
     Vin().setXmippOrigin();
     if (binarize)
         Vin().binarize(threshold,0);
 
     if (fnOut=="")
         fnOut=fnVol.withoutExtension();
 
     Vcurrent().initZeros(Vin());
     mask_prm.generate_mask(Vin());
 
     sigma3=3*sigma;
     gaussianTable.resize(CEIL(sigma3*sqrt(3.0)*1000));
     FOR_ALL_ELEMENTS_IN_ARRAY1D(gaussianTable)
     gaussianTable(i)=gaussian1D(i/1000.0,sigma);
 
     energyOriginal=0;
     double N=0;
     double minval=1e38, maxval=-1e38;
     const MultidimArray<int> &iMask3D=mask_prm.get_binary_mask();
     FOR_ALL_ELEMENTS_IN_ARRAY3D(Vin())
     {
         if (useMask && iMask3D(k,i,j)==0)
             continue;
         double v=Vin(k,i,j);
         energyOriginal+=v*v;
         minval=XMIPP_MIN(minval,v);
         maxval=XMIPP_MAX(maxval,v);
         N++;
     }
     energyOriginal/=N;
 
     Histogram1D hist;
     if (useMask)
         compute_hist_within_binary_mask(iMask3D, Vin(), hist,
                                         minval, maxval, 200);
     else
         compute_hist(Vin(), hist, minval, maxval, 200);
     percentil1=hist.percentil(1);
     if (percentil1<=0)
         percentil1=maxval/500;
     range = hist.percentil(99)-percentil1;
 
     if (XMIPP_EQUAL_ZERO(range))
         REPORT_ERROR(ERR_VALUE_INCORRECT, "Range cannot be zero.");
 
     smallAtom=range*intensityFraction;
 
     // Create threads
     barrier_init(&barrier,numThreads+1);
     threadIds=(pthread_t *)malloc(numThreads*sizeof(pthread_t));
     threadArgs=(Prog_Convert_Vol2Pseudo_ThreadParams *)
                malloc(numThreads*sizeof(Prog_Convert_Vol2Pseudo_ThreadParams));
     for (int i=0; i<numThreads; i++)
     {
         threadArgs[i].myThreadID=i;
         threadArgs[i].parent=this;
         pthread_create( (threadIds+i), nullptr, optimizeCurrentAtomsThread,
                         (void *) (threadArgs+i));
     }
 
     // Filter for the difference volume
     Filter.FilterBand=LOWPASS;
     Filter.FilterShape=REALGAUSSIAN;
     Filter.w1=sigma;
     Filter.generateMask(Vin());
     Filter.do_generate_3dmask=false;
 }

◆ readParams()

void ProgVolumeToPseudoatoms::readParams ( )

virtual

Read parameters from command line.

Reimplemented from XmippProgram.

Definition at line 51 of file volume_to_pseudoatoms.cpp.

 {
     fnVol = getParam("-i");
     fnOut = getParam("-o");
     mask_prm.allowed_data_types = INT_MASK;
     if ((useMask = checkParam("--mask")))
         mask_prm.readParams(this);
     sigma = getDoubleParam("--sigma");
     targetError = getDoubleParam("--targetError")/100.0;
     stop = getDoubleParam("--stop");
     initialSeeds = getIntParam("--initialSeeds");
     growSeeds = getDoubleParam("--growSeeds");
     allowMovement = !checkParam("--dontAllowMovement");
     allowIntensity = !checkParam("--dontAllowIntensity");
     if (!allowIntensity)
         intensityFraction = getDoubleParam("--dontAllowIntensity",1);
     intensityColumn = getParam("--intensityColumn");
     Nclosest = getIntParam("--Nclosest");
     minDistance = getDoubleParam("--minDistance");
     penalty = getDoubleParam("--penalty");
     numThreads = getIntParam("--thr");
     sampling = getDoubleParam("--sampling_rate");
     dontScale = checkParam("--dontScale");
     binarize = checkParam("--binarize");
     if (binarize)
         threshold=getDoubleParam("--binarize");
     else
         threshold=0;
 }

◆ removeSeeds()

void ProgVolumeToPseudoatoms::removeSeeds ( int Nseeds )

Remove seeds.

Definition at line 403 of file volume_to_pseudoatoms.cpp.

 {
     int fromNegative=ROUND(Nseeds*0.5);
     int fromSmall=Nseeds-fromNegative;
 
     if (allowIntensity)
     {
         // Remove too small atoms
         std::sort(atoms.begin(),atoms.end());
         atoms.erase(atoms.begin(),atoms.begin()+fromSmall);
     }
     else
     {
         fromNegative=Nseeds;
         fromSmall=0;
     }
 
     // Remove atoms from regions in which the error is too negative
     MultidimArray<double> Vdiff=Vin();
     Vdiff-=Vcurrent();
     int alreadyRemoved=0;
     double vmin=Vdiff.computeMin();
     if (vmin<0)
     {
         for (double v=vmin+vmin/20; v<0; v-=vmin/20)
         {
             size_t oldListSize;
             do
             {
                 oldListSize=atoms.size();
 
                 // Search for a point within a negative region
                 bool found=false;
                 int kneg=0, ineg=0, jneg=0;
                 const MultidimArray<int> &iMask3D=mask_prm.get_binary_mask();
                 for (int k=STARTINGZ(Vdiff); k<=FINISHINGZ(Vdiff) && !found; k++)
                     for (int i=STARTINGY(Vdiff); i<=FINISHINGY(Vdiff) && !found; i++)
                         for (int j=STARTINGX(Vdiff); j<=FINISHINGX(Vdiff) && !found; j++)
                         {
                             if (useMask && iMask3D(k,i,j)==0)
                                 continue;
                             if (A3D_ELEM(Vdiff,k,i,j)<v)
                             {
                                 kneg=k;
                                 ineg=i;
                                 jneg=j;
                                 A3D_ELEM(Vdiff,k,i,j)=0;
                                 found=true;
                             }
                         }
 
                 // If found such a point, search for a nearby atom
                 if (found)
                 {
                     // Search for atom
                     int nmax=atoms.size();
                     for (int n=0; n<nmax; n++)
                     {
                         double r=
                             (kneg-atoms[n].location(0))*(kneg-atoms[n].location(0))+
                             (ineg-atoms[n].location(1))*(ineg-atoms[n].location(1))+
                             (jneg-atoms[n].location(2))*(jneg-atoms[n].location(2));
                         r=sqrt(r);
                         if (r<sigma3)
                         {
                             drawGaussian(atoms[n].location(0),
                                          atoms[n].location(1),
                                          atoms[n].location(2),
                                          Vdiff,
                                          atoms[n].intensity);
                             atoms.erase(atoms.begin()+n);
                             alreadyRemoved++;
                             break;
                         }
                     }
                 }
             }
             while (oldListSize>atoms.size() && alreadyRemoved<fromNegative);
             if (alreadyRemoved==fromNegative)
                 break;
         }
     }
 
     removeTooCloseSeeds();
 }

◆ removeTooCloseSeeds()

void ProgVolumeToPseudoatoms::removeTooCloseSeeds ( )

Remove too close seeds.

Definition at line 489 of file volume_to_pseudoatoms.cpp.

 {
     // Remove atoms that are too close to each other
     if (minDistance>0 && allowIntensity)
     {
         std::vector<int> toRemove;
         int nmax=atoms.size();
         double minDistance2=minDistance*minDistance;
         for (int n1=0; n1<nmax; n1++)
         {
             bool found=false;
             int nn=0, nnmax=toRemove.size();
             while (nn<nnmax)
             {
                 if (toRemove[nn]==n1)
                 {
                     found=true;
                     break;
                 }
                 else if (toRemove[nn]>n1)
                     break;
                 nn++;
             }
             if (found)
                 continue;
             for (int n2=n1+1; n2<nmax; n2++)
             {
                 nn=0;
                 found=false;
                 while (nn<nnmax)
                 {
                     if (toRemove[nn]==n2)
                     {
                         found=true;
                         break;
                     }
                     else if (toRemove[nn]>n2)
                         break;
                     nn++;
                 }
                 if (found)
                     continue;
                 double diffZ=atoms[n1].location(0)-atoms[n2].location(0);
                 double diffY=atoms[n1].location(1)-atoms[n2].location(1);
                 double diffX=atoms[n1].location(2)-atoms[n2].location(2);
                 double d2=diffZ*diffZ+diffY*diffY+diffX*diffX;
                 if (d2<minDistance2)
                 {
                     if (atoms[n1].intensity<atoms[n2].intensity)
                     {
                         toRemove.push_back(n1);
                         break;
                     }
                     else
                         toRemove.push_back(n2);
                     std::sort(toRemove.begin(),toRemove.end());
                 }
             }
         }
         for (int n=toRemove.size()-1; n>=0; n--)
             atoms.erase(atoms.begin()+toRemove[n]);
     }
 }

◆ run()

void ProgVolumeToPseudoatoms::run ( )

virtual

Run.

Reimplemented from XmippProgram.

Definition at line 966 of file volume_to_pseudoatoms.cpp.

 {
     produceSideInfo();
     int iter=0;
     double previousNAtoms=0;
     percentageDiff=1;
     double actualGrowSeeds=0.;
     do
     {
         // Place seeds
         if (iter==0)
             placeSeeds(initialSeeds);
         else
         {
             double Natoms=atoms.size();
             actualGrowSeeds=growSeeds*std::min(1.0,0.1+(percentageDiff-targetError)/targetError);
             removeSeeds(FLOOR(Natoms*(actualGrowSeeds/2)/100));
             placeSeeds(FLOOR(Natoms*actualGrowSeeds/100));
         }
         drawApproximation();
 
         if (iter==0 && verbose>0)
             std::cout << "Initial error with " << atoms.size()
             << " pseudo-atoms " << percentageDiff << std::endl;
 
         // Optimize seeds until convergence
         optimizeCurrentAtoms();
         if (verbose>0)
             std::cout << "Error with " << atoms.size() << " pseudo-atoms "
             << percentageDiff << std::endl;
         writeResults();
         iter++;
         
         if (fabs(previousNAtoms-atoms.size())/atoms.size()<0.01*actualGrowSeeds/100)
         {
             std::cout << "The required precision cannot be attained\n"
             << "Suggestion: Reduce sigma and/or minDistance\n"
             << "Writing best approximation with current parameters\n";
 
             break;
         }
         previousNAtoms=atoms.size();
     }
     while (percentageDiff>targetError);
     removeTooCloseSeeds();
     writeResults();
 
     // Kill threads
     threadOpCode=KILLTHREAD;
     barrier_wait(&barrier);
     free(threadIds);
     free(threadArgs);
 }

◆ show()

void ProgVolumeToPseudoatoms::show ( ) const

virtual

show parameters

Reimplemented from XmippProgram.

Definition at line 81 of file volume_to_pseudoatoms.cpp.

 {
     if (verbose==0)
         return;
     std::cout << "Input volume:   " << fnVol             << std::endl
     << "Output volume:  " << fnOut             << std::endl
     << "Sigma:          " << sigma             << std::endl
     << "Initial seeds:  " << initialSeeds      << std::endl
     << "Grow seeds:     " << growSeeds         << std::endl
     << "Target error:   " << targetError       << std::endl
     << "Stop:           " << stop              << std::endl
     << "AllowMovement:  " << allowMovement     << std::endl
     << "AllowIntensity: " << allowIntensity    << std::endl
     << "Intensity Frac: " << intensityFraction << std::endl
     << "Intensity Col:  " << intensityColumn   << std::endl
     << "Nclosest:       " << Nclosest          << std::endl
     << "Min. Distance:  " << minDistance       << std::endl
     << "Penalty:        " << penalty           << std::endl
     << "Threads:        " << numThreads        << std::endl
     << "Sampling Rate:  " << sampling          << std::endl
     << "Don't scale:    " << dontScale         << std::endl
     << "Binarize:       " << binarize          << std::endl
     << "Threshold:      " << threshold         << std::endl
     ;
     if (useMask)
         mask_prm.show();
     else
         std::cout << "No mask\n";
 }

◆ writeResults()

void ProgVolumeToPseudoatoms::writeResults ( )

Write results.

Definition at line 885 of file volume_to_pseudoatoms.cpp.

 {
     // Compute the histogram of intensities
     MultidimArray<double> intensities;
     intensities.initZeros(atoms.size());
     FOR_ALL_ELEMENTS_IN_ARRAY1D(intensities)
     intensities(i)=atoms[i].intensity;
     Histogram1D hist;
     compute_hist(intensities, hist, 0, intensities.computeMax(), 100);
 
     if (verbose>=2)
     {
         Vcurrent.write(fnOut+"_approximation.vol");
         hist.write(fnOut+"_approximation.hist");
 
         // Save the difference
         Image<double> Vdiff;
         Vdiff()=Vin()-Vcurrent();
         const MultidimArray<int> &iMask3D=mask_prm.get_binary_mask();
         if (useMask && XSIZE(iMask3D)!=0)
             FOR_ALL_ELEMENTS_IN_ARRAY3D(Vdiff())
             if (!iMask3D(k,i,j))
                 Vdiff(k,i,j)=0;
         Vdiff.write(fnOut+"_rawDiff.vol");
 
         Vdiff()/=range;
         Vdiff.write(fnOut+"_relativeDiff.vol");
     }
 
     // Write the PDB
     double minIntensity=intensities.computeMin();
     double maxIntensity=intensities.computeMax();
     double a=0.99/(maxIntensity-minIntensity);
     if (dontScale)
         a=1;
 
     FILE *fhOut=nullptr;
     fhOut=fopen((fnOut+".pdb").c_str(),"w");
     if (!fhOut)
         REPORT_ERROR(ERR_IO_NOWRITE,fnOut+".pdb");
     int nmax=atoms.size();
     int col=1;
     if (intensityColumn=="Bfactor")
         col=2;
     fprintf(fhOut,"REMARK xmipp_volume_to_pseudoatoms\n");
     fprintf(fhOut,"REMARK fixedGaussian %f\n",sigma*sampling);
     fprintf(fhOut,"REMARK intensityColumn %s\n",intensityColumn.c_str());
     for (int n=0; n<nmax; n++)
     {
         double intensity=1.0;
         if (allowIntensity)
             intensity=0.01+ROUND(100*a*(atoms[n].intensity-minIntensity))/100.0;
         if (col==1)
             fprintf(fhOut,
                     "ATOM  %5d DENS DENS%5d    %8.3f%8.3f%8.3f%6.2f     1      DENS\n",
                     n+1,n+1,
                     (float)(atoms[n].location(2)*sampling),
                     (float)(atoms[n].location(1)*sampling),
                     (float)(atoms[n].location(0)*sampling),
                     (float)intensity);
         else
             fprintf(fhOut,
                     "ATOM  %5d DENS DENS%5d    %8.3f%8.3f%8.3f     1%6.2f      DENS\n",
                     n+1,n+1,
                     (float)(atoms[n].location(2)*sampling),
                     (float)(atoms[n].location(1)*sampling),
                     (float)(atoms[n].location(0)*sampling),
                     (float)intensity);
     }
     fclose(fhOut);
 
     if (verbose>=2)
     {
         PDBPhantom pdb;
         pdb.read(fnOut+".pdb");
         distanceHistogramPDB(pdb,Nclosest,-1,200,hist);
         hist.write(fnOut+"_distance.hist");
     }
 }

Member Data Documentation

◆ allowIntensity

bool ProgVolumeToPseudoatoms::allowIntensity

Allow gaussians to vary intensity.

Definition at line 106 of file volume_to_pseudoatoms.h.

◆ allowMovement

bool ProgVolumeToPseudoatoms::allowMovement

Allow gaussians to move.

Definition at line 103 of file volume_to_pseudoatoms.h.

◆ atoms

std::vector< PseudoAtom > ProgVolumeToPseudoatoms::atoms

Definition at line 210 of file volume_to_pseudoatoms.h.

◆ barrier

barrier_t ProgVolumeToPseudoatoms::barrier

Definition at line 228 of file volume_to_pseudoatoms.h.

◆ binarize

bool ProgVolumeToPseudoatoms::binarize

Binarize.

Definition at line 136 of file volume_to_pseudoatoms.h.

◆ dontScale

bool ProgVolumeToPseudoatoms::dontScale

Don't scale the atom weights at the end.

Definition at line 133 of file volume_to_pseudoatoms.h.

◆ energyDiff

double ProgVolumeToPseudoatoms::energyDiff

Definition at line 201 of file volume_to_pseudoatoms.h.

◆ energyOriginal

double ProgVolumeToPseudoatoms::energyOriginal

Definition at line 207 of file volume_to_pseudoatoms.h.

◆ Filter

FourierFilter ProgVolumeToPseudoatoms::Filter

Definition at line 242 of file volume_to_pseudoatoms.h.

◆ fnOut

FileName ProgVolumeToPseudoatoms::fnOut

Output volume.

Definition at line 79 of file volume_to_pseudoatoms.h.

◆ fnVol

FileName ProgVolumeToPseudoatoms::fnVol

Volume to convert.

Definition at line 76 of file volume_to_pseudoatoms.h.

◆ gaussianTable

MultidimArray<double> ProgVolumeToPseudoatoms::gaussianTable

Definition at line 225 of file volume_to_pseudoatoms.h.

◆ growSeeds

double ProgVolumeToPseudoatoms::growSeeds

Grow seeds.

Definition at line 100 of file volume_to_pseudoatoms.h.

◆ initialSeeds

int ProgVolumeToPseudoatoms::initialSeeds

Initial seeds.

Definition at line 97 of file volume_to_pseudoatoms.h.

◆ intensityColumn

std::string ProgVolumeToPseudoatoms::intensityColumn

Column for the intensity (if any)

Definition at line 115 of file volume_to_pseudoatoms.h.

◆ intensityFraction

double ProgVolumeToPseudoatoms::intensityFraction

Intensity fraction. In case intensity is not allowed to change, this fraction is multiplied by the intensity range and all atoms will have this intensity value.

Definition at line 112 of file volume_to_pseudoatoms.h.

◆ mask_prm

Mask ProgVolumeToPseudoatoms::mask_prm

Definition at line 82 of file volume_to_pseudoatoms.h.

◆ minDistance

double ProgVolumeToPseudoatoms::minDistance

Mindistance.

Definition at line 118 of file volume_to_pseudoatoms.h.

◆ Nclosest

size_t ProgVolumeToPseudoatoms::Nclosest

N closest atoms for the distance histogram.

Definition at line 130 of file volume_to_pseudoatoms.h.

◆ numThreads

int ProgVolumeToPseudoatoms::numThreads

Number of threads.

Definition at line 124 of file volume_to_pseudoatoms.h.

◆ penalty

double ProgVolumeToPseudoatoms::penalty

Penalization.

Definition at line 121 of file volume_to_pseudoatoms.h.

◆ percentageDiff

double ProgVolumeToPseudoatoms::percentageDiff

Definition at line 204 of file volume_to_pseudoatoms.h.

◆ percentil1

double ProgVolumeToPseudoatoms::percentil1

Definition at line 216 of file volume_to_pseudoatoms.h.

◆ range

double ProgVolumeToPseudoatoms::range

Definition at line 219 of file volume_to_pseudoatoms.h.

◆ sampling

double ProgVolumeToPseudoatoms::sampling

Sampling rate.

Definition at line 127 of file volume_to_pseudoatoms.h.

◆ sigma

double ProgVolumeToPseudoatoms::sigma

Sigma.

Definition at line 88 of file volume_to_pseudoatoms.h.

◆ sigma3

double ProgVolumeToPseudoatoms::sigma3

Definition at line 213 of file volume_to_pseudoatoms.h.

◆ smallAtom

double ProgVolumeToPseudoatoms::smallAtom

Definition at line 222 of file volume_to_pseudoatoms.h.

◆ stop

double ProgVolumeToPseudoatoms::stop

Stop criterion.

Definition at line 94 of file volume_to_pseudoatoms.h.

◆ targetError

double ProgVolumeToPseudoatoms::targetError

Maximum error (as a percentage)

Definition at line 91 of file volume_to_pseudoatoms.h.

◆ threadArgs

Prog_Convert_Vol2Pseudo_ThreadParams* ProgVolumeToPseudoatoms::threadArgs

Definition at line 236 of file volume_to_pseudoatoms.h.

◆ threadIds

pthread_t* ProgVolumeToPseudoatoms::threadIds

Definition at line 239 of file volume_to_pseudoatoms.h.

◆ threadOpCode

int ProgVolumeToPseudoatoms::threadOpCode

Definition at line 233 of file volume_to_pseudoatoms.h.

◆ threshold

double ProgVolumeToPseudoatoms::threshold

Threshold for the binarization.

Definition at line 139 of file volume_to_pseudoatoms.h.

◆ useMask

bool ProgVolumeToPseudoatoms::useMask

Definition at line 85 of file volume_to_pseudoatoms.h.

◆ Vcurrent

Image<double> ProgVolumeToPseudoatoms::Vcurrent

Definition at line 198 of file volume_to_pseudoatoms.h.

◆ Vin

Image<double> ProgVolumeToPseudoatoms::Vin

Definition at line 195 of file volume_to_pseudoatoms.h.

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

xmipp/libraries/reconstruction/volume_to_pseudoatoms.h
xmipp/libraries/reconstruction/volume_to_pseudoatoms.cpp

Public Member Functions

Static Public Member Functions

Public Attributes

Additional Inherited Members

Detailed Description

Member Function Documentation

◆ computeAverage()

◆ defineParams()

◆ drawApproximation()

◆ drawGaussian()

◆ evaluateRegion()

◆ extractRegion()

◆ insertRegion()

◆ optimizeCurrentAtoms()

◆ optimizeCurrentAtomsThread()

◆ placeSeeds()

◆ produceSideInfo()

◆ readParams()

◆ removeSeeds()

◆ removeTooCloseSeeds()

◆ run()

◆ show()

◆ writeResults()

Member Data Documentation

◆ allowIntensity

◆ allowMovement

◆ atoms

◆ barrier

◆ binarize

◆ dontScale

◆ energyDiff

◆ energyOriginal

◆ Filter

◆ fnOut

◆ fnVol

◆ gaussianTable

◆ growSeeds

◆ initialSeeds

◆ intensityColumn

◆ intensityFraction

◆ mask_prm

◆ minDistance

◆ Nclosest

◆ numThreads

◆ penalty

◆ percentageDiff

◆ percentil1

◆ range

◆ sampling

◆ sigma

◆ sigma3

◆ smallAtom

◆ stop

◆ targetError

◆ threadArgs

◆ threadIds

◆ threadOpCode

◆ threshold

◆ useMask

◆ Vcurrent

◆ Vin