XRayPSF Class Reference

#include <psf_xr.h>

Collaboration diagram for XRayPSF:

Public Types
enum	operMode { GENERATE_PSF, PSF_FROM_FILE }

Public Member Functions
	XRayPSF ()
	XRayPSF (const XRayPSF &)=delete
	XRayPSF (const XRayPSF &&)=delete
XRayPSF &	operator= (const XRayPSF &)=delete
XRayPSF &	operator= (const XRayPSF &&)=delete
	~XRayPSF ()
void	init ()
void	clear ()
	Clear. More...
void	readParams (XmippProgram *program)
void	read (const FileName &fn, bool readVolume=true)
void	write (const FileName &fn)
void	show ()
	Show the microscope parameters. More...
void	setFocalShift (double zShift)
	Add focal shift to previously read psf zshift. More...
void	calculateParams (double _dxo, double _dzo=-1, double threshold=0.)
	Produce Side information. More...
void	reducePSF2Slabs (double threshold)
	Calculate the width of the slabs to reduce computing time and the mean PSF for each. More...
void	applyOTF (MultidimArray< double > &Im, const double sliceOffset) const
	Apply the OTF to the image, by means of the convolution. More...
void	generateOTF (MultidimArray< std::complex< double > > &OTF, double Zpos) const
	Generate the Optical Transfer Function (OTF) for a slice according to Microscope and Im parameters. More...
void	generatePSF ()
	Generate the 3D Point Spread Function (PSF) according to Microscope parameters. More...
void	adjustParam ()
	Calculate if a resize of the X-Y plane is needed to avoid the Nyquist Limit. More...
void	adjustParam (MultidimArray< double > &Vol)

Static Public Member Functions
static void	defineParams (XmippProgram *program)

Public Attributes
operMode	mode
PsfType	type
	Define the selected PSF generation algorithm. More...
ImageGeneric	psfGen
Image< double >	psfVol
	3D PSF read from file More...
MultidimArray< double >	PSF
	Working PSF with nonlinear zdim whose slices are the mean PSF for slabs. More...
double	slabThr
	Threshold to separate The volume into slabs to use the same PSF. More...
std::vector< int >	slabIndex
	Z positions in the original PSF Volume to determine de slabs. More...
Matrix2D< double >	T
MultidimArray< double > *	mask
	Lens shape Mask. More...
double	Rlens
	Lens Aperture Radius. More...
double	Zo
	Object plane on Focus (Reference) More...
double	Zi
	Object plane. More...
double	DoF
	Depth of focus. Only for information purposes. More...
size_t	Nox
	X size of the input image (object plane size) More...
size_t	Noy
	Y size of the input image (object plane size) More...
size_t	Noz
	Z size of the input image (object plane size) More...
double	dxiMax
double	dxl
	Pixel size in X-dim in lens plane. More...
double	dyl
	Pixel size in Y-dim in lens plane. More...
double	deltaZMaxX
	Z limits around Zo in the psf generation due to Nyquist Limit. More...
double	deltaZMaxY
double	deltaZMinX
double	deltaZMinY
PsfxrAdjust	AdjustType
	Parameters to change image size to avoid Nyquist limit. More...
double	pupileSizeMin
	Minimum diameter size of the microscope pupile in the lens plane, measured in pixels. More...
double	lambda
	Lambda of illumination. More...
double	Flens
	Focal length. More...
double	Nzp
	Number of zones in zone plate. More...
double	deltaR
	Outermost zone width. More...
double	Ms
	Magnification. More...
double	DeltaZo
	Z axis global shift. More...
double	dxo
	object space XY-plane sampling rate More...
double	dxi
	Image space XY-plane sampling rate. More...
double	dzo
	object space Z sampling rate More...
size_t	Nix
	Size of the image in image plane, to be rescaled if needed. More...
size_t	Niy
double	dxoPSF
	object space XY-plane sampling rate of the PSF Volume More...
double	dzoPSF
	object space Z sampling rate of the PSF Volume More...
int	verbose
	Switch to control verbose mode. More...
int	nThr

Protected Member Functions
void	generatePSFIdealLens (MultidimArray< double > &PSFi, double Zpos) const
	Generate the PSF for a single plane according to a ideal lens. More...

Friends
std::ostream &	operator<< (std::ostream &out, const XRayPSF &psf)
	Show. More...

Detailed Description

X-ray PSF class. Here goes how to filter an image with the Point Spread Function of a X-ray microscope optics

#include <data/psf_xr.h>


int main(int argc, char **argv)
{
 FileName fnPSF, fn_input,fn_output;
 XRayPSF psf;


if (checkParameter(argc, argv, "-psf"))
{
 fnPSF = getParameter(argc, argv, "-psf");
 psf.read(fnPSF);
}
else
 psf.clear(); 


fn_input = getParameter(argc, argv, "-i");

psf.produceSideInfo();

Image<double>   inputVol, imOut;

inputVol.read(fn_input);
psf.adjustParam(inputVol);
project_xr(psf, inputVol, imOut);

if (checkParameter(argc, argv, "-out"))
 fn_output = getParameter(argc, argv, "-out");
else
 fn_output = file_name.without_extension().add_extension("out").add_extension("spi");

imOut.write(fn_output);
}

Definition at line 103 of file psf_xr.h.

Member Enumeration Documentation

operMode

enum XRayPSF::operMode

Enumerator
GENERATE_PSF
PSF_FROM_FILE

Definition at line 108 of file psf_xr.h.

    {
        GENERATE_PSF,
        PSF_FROM_FILE
    } operMode;

Constructor & Destructor Documentation

XRayPSF() [1/3]

XRayPSF::XRayPSF

(

)

Empty constructor.

Definition at line 33 of file psf_xr.cpp.

{
    init();
}

XRayPSF() [2/3]

XRayPSF::XRayPSF ( const XRayPSF &  )

delete

XRayPSF() [3/3]

XRayPSF::XRayPSF ( const XRayPSF &&  )

delete

~XRayPSF()

XRayPSF::~XRayPSF

(

)

Definition at line 38 of file psf_xr.cpp.

{
    clear();
}

Member Function Documentation

adjustParam() [1/2]

void XRayPSF::adjustParam

(

)

Calculate if a resize of the X-Y plane is needed to avoid the Nyquist Limit.

Lens Radius in pixels higher than image

For indices in standard fashion

Definition at line 673 of file psf_xr.cpp.

{
    Nix = Nox;
    Niy = Noy;
    AdjustType = PSFXR_STD;

    if (mode == GENERATE_PSF)
    {
        switch (type)
        {
        case IDEAL_FRESNEL_LENS:
            {
                // Reset the image plane sampling value
                dxi = dxo * Ms;

                dxl = lambda*Zi / (Nox * dxi); // Pixel X-size en the plane of lens aperture
                dyl = lambda*Zi / (Noy * dxi); // Pixel Y-size en the plane of lens aperture

                deltaZMaxX = Zo*((Rlens*2*dxl)/(Rlens*2*dxl - Zo*lambda) - 1);
                deltaZMaxY = Zo*((Rlens*2*dyl)/(Rlens*2*dyl - Zo*lambda) - 1);
                deltaZMinX = Zo*((Rlens*2*dxl)/(Rlens*2*dxl + Zo*lambda) - 1);
                deltaZMinY = Zo*((Rlens*2*dyl)/(Rlens*2*dyl + Zo*lambda) - 1);


                if (verbose > 1)
                {
                    std::cout << std::endl;
                    std::cout << "----------------------" << std::endl;
                    std::cout << "XrayPSF::Param adjust:" << std::endl;
                    std::cout << "----------------------" << std::endl;
                    std::cout << "(Nox,Noy,Nz) = (" << Nox << "," << Noy << "," << Noz << ")" << std::endl;
                    std::cout << "Larger volume Z plane to be calculated = " << (ABS(DeltaZo)+(Noz-1)/2*dzo)*1e6 << " um" << std::endl;
                    std::cout << "Larger allowed discrete Z plane (x,y) = (" << deltaZMaxX*1e6 << ", " << deltaZMaxY*1e6 << ") um" << std::endl;
                }

                if (dxi>dxiMax) 
                {
                    Nix = (size_t)ceil(Nox * dxi/dxiMax);
                    Niy = (size_t)ceil(Noy * dxi/dxiMax);

                    dxi *= Nox/Nix;

                    AdjustType = PSFXR_INT;

                    if (verbose > 1)
                        std::cout << "XrayPSF: Image plane sampling too small: increasing resolution" << std::endl;

                }
                else
                {
                    if (dxi < pupileSizeMin/Nox * dxiMax)
                    {
                        Nix = (size_t)ceil(pupileSizeMin * dxiMax/dxi);
                        AdjustType = PSFXR_ZPAD;
                    }
                    if (dxi < pupileSizeMin/Noy * dxiMax)
                    {
                        Niy = (size_t)ceil(pupileSizeMin * dxiMax/dxi);
                        AdjustType = PSFXR_ZPAD;
                    }
                    if (DeltaZo + Noz/2*dzo > deltaZMaxX)
                    {
                        Nix = std::max(Nix,(size_t)ceil(Zi*Rlens*2*fabs(DeltaZo+Noz/2*dzo)/(Zo*dxi*(Zo+DeltaZo+Noz/2*dzo))));
                        AdjustType = PSFXR_ZPAD;
                    }
                    if (DeltaZo - Noz/2*dzo < deltaZMinX)
                    {
                        Nix = std::max(Nix,(size_t)ceil(Zi*Rlens*2*fabs(DeltaZo-Noz/2*dzo)/(Zo*dxi*(Zo+DeltaZo-Noz/2*dzo))));
                        AdjustType = PSFXR_ZPAD;
                    }
                    if (DeltaZo + Noz/2*dzo > deltaZMaxY)
                    {
                        Niy = std::max(Niy,(size_t)ceil(Zi*Rlens*2*fabs(DeltaZo+Noz/2*dzo)/(Zo*dxi*(Zo+DeltaZo+Noz/2*dzo))));
                        AdjustType = PSFXR_ZPAD;
                    }
                    if ( DeltaZo - Noz/2*dzo < deltaZMinY)
                    {
                        Niy = std::max(Niy,(size_t)ceil(Zi*Rlens*2*fabs(DeltaZo-Noz/2*dzo)/(Zo*dxi*(Zo+DeltaZo-Noz/2*dzo))));
                        AdjustType = PSFXR_ZPAD;
                    }
                }

                dxl = lambda*Zi / (Nix * dxi); // Pixel X-size en the plane of lens aperture
                dyl = lambda*Zi / (Niy * dxi); // Pixel Y-size en the plane of lens aperture

                if (verbose > 1)
                {
                    if (AdjustType!=PSFXR_STD)
                    {
                        if (AdjustType==PSFXR_ZPAD)
                            std::cout << "XrayPSF: Image plane size too small: increasing size" << std::endl;
                        std::cout << "New slice dimensions:  " <<  "(Nix, Niy) = (" << Nix << ", " << Niy << ")" << std::endl;
                    }
                    std::cout << "(dxl, dyl) = (" << dxl*1e6 << ", " << dyl*1e6 << ") um" << std::endl;
                    std::cout << "Pupile Diameter in pixels (NpX, NpY) = (" << ceil(2*Rlens / dxl)
                    << ", " << ceil(2*Rlens / dyl) << ")"  << std::endl;
                    std::cout << std::endl;
                }

                // Calculate the mask to be applied when generating PSFIdealLens

                mask->initZeros(Niy,Nix);
                mask->setXmippOrigin();

                double Rlens2 = Rlens*Rlens;
                //                double auxY = dyl*(1 - (int)Niy);
                //                double auxX = dxl*(1 - (int)Nix);

                for (int i = STARTINGY(*mask); i <= FINISHINGY(*mask); ++i)
                {
                    double y = dyl * i;
                    double y2 = y * y;
                    for (int  j = STARTINGX(*mask); j <= FINISHINGX(*mask); ++j)// Circular mask
                    {
                        double x = dxl * j;
                        if (x*x + y2 <= Rlens2)
                            A2D_ELEM(*mask,i,j) = 1;
                    }
                }


                break;
            }
        case ANALYTIC_ZP:
            {
                REPORT_ERROR(ERR_NOT_IMPLEMENTED, "Zone Plate algorithm not implemented yet.");
                break;
            }
        }
    }
    else
    {
        //        OTF.clear(); // Only for each projection
    }

    if (verbose > 1)
        --verbose;
}

adjustParam() [2/2]

void XRayPSF::adjustParam

(

MultidimArray< double > &

Vol

)

Definition at line 664 of file psf_xr.cpp.

{
    Nox = XSIZE(vol);
    Noy = YSIZE(vol);
    Noz = ZSIZE(vol);

    adjustParam();
}

applyOTF()

void XRayPSF::applyOTF	(	MultidimArray< double > &	Im,
		const double	sliceOffset
	)		const

Apply the OTF to the image, by means of the convolution.

Definition at line 384 of file psf_xr.cpp.

{
    //    double Z;
    //    MultidimArray< std::complex<double> > OTF;

    if (mode == GENERATE_PSF)
        REPORT_ERROR(ERR_VALUE_NOTSET, "XrayPSF::applyOTF: a PSF from file must be read.");

    double Z = Zo + DeltaZo + sliceOffset;

    MultidimArray<std::complex<double> > OTF;
    generateOTF(OTF, Z);

    MultidimArray<std::complex<double> > ImFT;
    FourierTransformer FTransAppOTF(FFTW_BACKWARD);

    //#define DEBUG
#ifdef DEBUG

    Image<double> _Im;
    _Im() = Im;

    _Im.write(("psfxr-ImIn.spi"));
#endif

    FTransAppOTF.FourierTransform(Im, ImFT, false);

#ifdef DEBUG

    _Im() = Im;

    _Im.write(("psfxr-ImIn_after.spi"));

    //    Image<double> _Im;
    _Im().clear();
    _Im().resizeNoCopy(OTF);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(OTF)
    dAij(_Im(),i,j) = abs(dAij(OTF,i,j));
    _Im.write("psfxr-otf.spi");

    _Im().resizeNoCopy(ImFT);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(ImFT)
    dAij(_Im(),i,j) = abs(dAij(ImFT,i,j));
    _Im.write(("psfxr-imft1.spi"));
#endif

    //    double normSize = MULTIDIM_SIZE(OTF);

    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(ImFT)
    dAij(ImFT,i,j) *= dAij(OTF,i,j);

#ifdef DEBUG

    _Im().resize(ImFT);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(ImFT)
    dAij(_Im(),i,j) = abs(dAij(ImFT,i,j));
    _Im.write(("psfxr-imft2.spi"));
#endif

    FTransAppOTF.inverseFourierTransform();

#ifdef DEBUG

    _Im() = Im;

    //    _Im().resize(Im);
    //    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(Im)
    //    dAij(_Im(),i,j) = abs(dAij(Im,i,j));
    _Im.write(("psfxr-imout.spi"));
#endif
#undef DEBUG
}

calculateParams()

void XRayPSF::calculateParams	(	double	_dxo,
		double	_dzo = -1,
		double	threshold = 0.
	)

Produce Side information.

Calculation of the rest of microscope parameters

TODO: The scale in Z may not be necessary to save memory and should be used the nearest PSF slice for each Z plane

Definition at line 279 of file psf_xr.cpp.

{

    dxo = (_dxo > 0)? _dxo : dxoPSF;
    dzo = (_dzo > 0)? _dzo : dxo;


    Rlens = (4 * Nzp * deltaR)*0.5; // Eq. 9.13 from Soft X-rays and .... David Attwood
    Flens = (2 * Rlens * deltaR)/lambda;

    Zo = (1 + 1 / Ms) * Flens;
    Zi = Zo * Ms;
    dxi = dxo * Ms;
    dxiMax = lambda * Zi / (2 * Rlens);
    DoF = 4*deltaR*deltaR/lambda;

    if(verbose > 0)
    {
        show();
        verbose++; // Show only once the param adjust info
    }

    if (mode == PSF_FROM_FILE)
    {
        T.initIdentity(4);
        double scaleFactor = dxo/dxoPSF;
        //        double scaleFactorZ = dxo/dzoPSF;
        double scaleFactorZ = 1;

        dMij(T, 0, 0) = scaleFactor;
        dMij(T, 1, 1) = scaleFactor;
        dMij(T, 2, 2) = scaleFactorZ;
        MultidimArray<double> &mdaPsfVol = psfVol();

        mdaPsfVol.resize(1, (size_t)(Noz/scaleFactorZ),
                         (size_t)(Noy/scaleFactor),
                         (size_t)(Nox/scaleFactor), false);

        mdaPsfVol.setXmippOrigin();

        applyGeometry(xmipp_transformation::LINEAR, mdaPsfVol, psfGen(), T,
                      xmipp_transformation::IS_INV, xmipp_transformation::DONT_WRAP);

        psfGen.clear(); // Free mem in case image is not mapped

        // Rescale the PSF values to keep the normalization in each X-Y plane
        mdaPsfVol *= scaleFactor*scaleFactor;

        // Reset the sampling values so now psfVol matches input phantom's sampling
        dxoPSF = dxo;
        //        dzoPSF = dzo;
        T.initIdentity(4);

        if (threshold != 0.)
            reducePSF2Slabs(threshold);
    }
    else // Mask is used when creating PSF on demand
    {
        mask = new MultidimArray<double>;
    }
}

clear()

void XRayPSF::clear

(

)

Clear.

Definition at line 267 of file psf_xr.cpp.

{
    psfVol.clear();
    init();
}

defineParams()

void XRayPSF::defineParams ( XmippProgram *  program )

static

Definition at line 45 of file psf_xr.cpp.

{
    //    program->addParamsLine(" [-file <psf_param_file>]       : Read X-ray psf parameters from file.");
    program->addParamsLine(" [--lambda <v=2.43>]             : X-ray wavelength (nm).");
    program->addParamsLine(" [--out_width <deltaR=40>]       : Outermost zone width of the X-ray Fresnel lens (nm).");
    program->addParamsLine(" [--zones <N=560>]               : Number of zones of the X-ray Fresnel lens.");
    program->addParamsLine(" [--mag <Ms=2304>]               : Magnification of the X-ray microscope.");
    program->addParamsLine(" [--sampling <dxy=10> <dz=dxy>]  : Sampling rate in X-Y plane and Z axis (nm).");
    program->addParamsLine(" [--zshift <deltaZ=0>]           : Longitudinal displacement along Z axis in microns.");
    program->addParamsLine(" [--size <x> <y=x> <z=x>]        : Size of the X-ray PSF volume.");
    program->addParamsLine(" [--type  <lens_type=ideal>]     : Lens type to generate the PSF.");
    program->addParamsLine("           where <lens_type>");
    program->addParamsLine("                  ideal          : Ideal phase Fresnel lens.");
    program->addParamsLine("                  zp             : Fresnel Zone Plate lens.");

}

generateOTF()

void XRayPSF::generateOTF	(	MultidimArray< std::complex< double > > &	OTF,
		double	Zpos
	)		const

Generate the Optical Transfer Function (OTF) for a slice according to Microscope and Im parameters.

Definition at line 458 of file psf_xr.cpp.

{

    MultidimArray<double> PSFi;
    FourierTransformer ftGenOTF(FFTW_BACKWARD);
    switch (mode)
    {
    case GENERATE_PSF:
        {
            if (type == IDEAL_FRESNEL_LENS)
                generatePSFIdealLens(PSFi, Zpos);
            break;
        }
    case PSF_FROM_FILE:
        {
            PSFi.resizeNoCopy(Niy, Nix);
            PSFi.setXmippOrigin();

            const MultidimArray<double> &mPsfVol = MULTIDIM_ARRAY(psfVol);

            double zIndexPSF = (Zo - Zpos)/dzoPSF; // Slice index for Zpos

            if (zIndexPSF < STARTINGZ(mPsfVol) ||
                zIndexPSF > FINISHINGZ(mPsfVol))
                PSFi.initConstant(1/YXSIZE(PSFi));
            else
            {   /* Actually T transform matrix is set to identity, as sampling is the same for both psfVol and phantom
                                                                                                                 * It is only missing to set Z shift to select the slice */
                dMij(T, 2, 3) = zIndexPSF; // Distance from the focal plane
                applyGeometry(xmipp_transformation::LINEAR, PSFi, mPsfVol, T,
                              xmipp_transformation::IS_INV, xmipp_transformation::DONT_WRAP, dAkij(mPsfVol,0,0,0));
                CenterFFT(PSFi, true);
            }
            break;
        }
    }

    ftGenOTF.FourierTransform(PSFi, OTF);

    //#define DEBUG
#ifdef DEBUG

    Image<double> _Im;
    //    _Im() = PSFi;
    CenterFFT(PSFi, false);
    _Im().alias(PSFi);
    _Im.write("psfxr-psfi.spi");
    //    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(PSFi)
    //    dAij(_Im(),i,j) = abs(dAij(PSFi,i,j));
    //    _Im.write("psfxr-psfi.spi");
    //    CenterOriginFFT(OTFTemp,1);
    //    _Im().resize(OTFTemp);
    //    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(OTFTemp)
    //    dAij(_Im(),i,j) = abs(dAij(OTFTemp,i,j));
    //    _Im.write("psfxr-otf1.spi");
    //    CenterOriginFFT(OTFTemp,0);
    _Im.clear();
    _Im().resize(OTF);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(OTF)
    dAij(_Im(),i,j) = abs(dAij(OTF,i,j));
    _Im.write("psfxr-otf2.spi");

#endif
#undef DEBUG
}

generatePSF()

void XRayPSF::generatePSF

(

)

Generate the 3D Point Spread Function (PSF) according to Microscope parameters.

Definition at line 525 of file psf_xr.cpp.

{
    calculateParams(dxoPSF, dzoPSF);
    adjustParam();

    //    PSF = new Image<double>(Nox, Noy, Noz);
    //    VOLMATRIX(*PSF).setXmippOrigin();

    MultidimArray<double> &mPsfVol = MULTIDIM_ARRAY(psfVol);
    mPsfVol.resize(1,Noz, Niy, Nix, false);
    mPsfVol.setXmippOrigin();

    MultidimArray<double> PSFi;

    init_progress_bar(ZSIZE(mPsfVol));

    for (int k = STARTINGZ(mPsfVol), k2 = 0; k <= FINISHINGZ(mPsfVol); k++, k2++)
    {
        /* We keep sign of Z, Zo and DeltaZo positives in object space for the sake of simplicity in calculations,
         * it is, they increase opposite to Zdim in Xmipp reference system, so this is the reason to calculate
         * the plane Z by subtracting k*dzoPSF which keeps positive in XMIPP system.
         */

        double Z = Zo + DeltaZo - k*dzoPSF;

        switch (type)
        {
        case IDEAL_FRESNEL_LENS:
            {
                generatePSFIdealLens(PSFi, Z);
                break;
            }
        case ANALYTIC_ZP:
            {
                REPORT_ERROR(ERR_NOT_IMPLEMENTED, "Zone plates algorithm not implemented yet.");
                break;
            }
        }

        CenterFFT(PSFi, false);
        mPsfVol.setSlice(k, PSFi);

        progress_bar(k2+1);
    }

    /* When PSF is generated from IDEAL_FRESNEL_LENS, dimensions are adjusted to avoid subsampling,
     * in that case it is not recommended to crop XY plane to avoid the loss of normalization in PSF plane. */

    //    if (Nix != Nox || Niy != Noy)
    //        mPsfVol.selfWindow(STARTINGZ(mPsfVol),
    //                           FIRST_XMIPP_INDEX(Noy),FIRST_XMIPP_INDEX(Nox),
    //                           FINISHINGZ(mPsfVol),
    //                           FIRST_XMIPP_INDEX(Noy)+Noy-1,FIRST_XMIPP_INDEX(Nox)+Nox-1);
}

generatePSFIdealLens()

void XRayPSF::generatePSFIdealLens ( MultidimArray< double > &  PSFi, double  Zpos  ) const

protected

Generate the PSF for a single plane according to a ideal lens.

Definition at line 581 of file psf_xr.cpp.

{
    double focalEquiv = 1/(1/Zpos - 1/Zo); // inverse of defocus = 1/Z - 1/Zo

    MultidimArray< std::complex<double> > OTFTemp(Niy, Nix);
    MultidimArray< std::complex<double> > PSFiTemp;
    //    Mask_Params mask_prm; TODO do we have to include masks using this class?

    lensPD(OTFTemp, focalEquiv, lambda, dxl, dyl);

    // Apply the Shape mask
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(OTFTemp)
    {
        if (dAi(*mask,n) == 0)
            dAi(OTFTemp,n) = 0;
    }

    //#define DEBUG
#ifdef DEBUG
    Image<double> _Im;
    _Im().resize(OTFTemp);
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(OTFTemp)
    dAij(_Im(),i,j) = arg(dAij(OTFTemp,i,j));
    _Im.write("phase_lens.spi");
    FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY2D(OTFTemp)
    dAij(_Im(),i,j) = abs(dAij(OTFTemp,i,j));
    _Im.write("abs_lens.spi");
#endif
#undef DEBUG

    FourierTransformer transformer(FFTW_BACKWARD);
    transformer.FourierTransform(OTFTemp, PSFiTemp, false);
    double norm=0;
    double iNorm;

    PSFi.resizeNoCopy(PSFiTemp);

#ifdef DEBUG

    _Im.data.resizeNoCopy(PSFiTemp);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(PSFiTemp)
    DIRECT_MULTIDIM_ELEM(_Im.data,n) = abs(DIRECT_MULTIDIM_ELEM(PSFiTemp,n));
    _Im.write("psfitemp.spi");
#endif

    double aux;
    double aux2;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(PSFi)
    {
        aux=abs(DIRECT_MULTIDIM_ELEM(PSFiTemp,n));
        DIRECT_MULTIDIM_ELEM(PSFi,n)=aux2=aux*aux;
        norm += aux2;
    }

    iNorm = 1/norm;
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(PSFi)
    DIRECT_MULTIDIM_ELEM(PSFi,n) *= iNorm;

    //#define DEBUG
#ifdef DEBUG

    _Im() = PSFi;
    CenterFFT(_Im(),0);
    _Im.write("generate-psfi.spi");
#endif
#undef DEBUG


#ifdef DEBUG

    transformer.inverseFourierTransform();

    CenterOriginFFT(OTFTemp,0);
    FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(OTFTemp)
    {
        DIRECT_MULTIDIM_ELEM(_Im.data,n) = abs(DIRECT_MULTIDIM_ELEM(OTFTemp,n));
    }
    _Im.write("otftemp.spi");
#endif

#undef DEBUG
}

init()

void XRayPSF::init

(

)

Definition at line 250 of file psf_xr.cpp.

{
    lambda = 2.43e-9;
    //    Flens = 1.4742e-3;
    deltaR = 40e-9;
    Nzp = 560;
    Ms = 2304;
    dzo = dxo = 1e-9;
    DeltaZo = 0;
    pupileSizeMin = 5;

    mode = GENERATE_PSF;
    type = IDEAL_FRESNEL_LENS;

    //    ftGenOTF.setNormalizationSign(FFTW_BACKWARD);
}

operator=() [1/2]

XRayPSF& XRayPSF::operator= ( const XRayPSF &  )

delete

operator=() [2/2]

XRayPSF& XRayPSF::operator= ( const XRayPSF &&  )

delete

read()

void XRayPSF::read	(	const FileName &	fn,
		bool	readVolume = true
	)

Read from file. An exception is thrown if the file cannot be open.

Definition at line 82 of file psf_xr.cpp.

{

    if (fn == "")
    {
        clear();
        return;
    }

    if (fn.isMetaData())
    {
        MetaDataVec MD;
        MD.read(fn);
        size_t id = MD.firstRowId();

        FileName fnPSF;
        if ( MD.getValue(MDL_IMAGE,fnPSF, id) && readVolume )
        {
            mode = PSF_FROM_FILE;
            psfGen.readMapped(fn.getDir()+fnPSF);
            ArrayDim aDim;
            psfGen.getDimensions(aDim);
            //            MULTIDIM_ARRAY(psfVol).getDimensions(aDim);
            Nox = aDim.xdim;
            Noy = aDim.ydim;
            Noz = aDim.zdim;

            psfGen().setXmippOrigin();

            //            MULTIDIM_ARRAY(psfVol).setXmippOrigin();
        }
        else
        {
            std::vector<double> dimV;
            if (!MD.getValue(MDL_CTF_DIMENSIONS,dimV, id))
                REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_CTF_DIMENSIONS) + " argument not present.");
            Nox = (int)dimV[0];
            Noy = (int)dimV[1];
            Noz = (int)dimV[2];
        }
        String typeS;
        if (!MD.getValue(MDL_CTF_XRAY_LENS_TYPE, typeS, id))
            REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_CTF_XRAY_LENS_TYPE) + " argument not present.");
        if (typeS == "ZP")
            type = ANALYTIC_ZP;
        else
            type = IDEAL_FRESNEL_LENS;

        if (!MD.getValue(MDL_CTF_LAMBDA,lambda, id))
            REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_CTF_LAMBDA) + " argument not present.");
        lambda *= 1e-9;
        if (!MD.getValue(MDL_CTF_XRAY_OUTER_ZONE_WIDTH,deltaR, id))
            REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_CTF_XRAY_OUTER_ZONE_WIDTH) + " argument not present.");
        deltaR *= 1e-9;
        if (!MD.getValue(MDL_CTF_XRAY_ZONES_NUMBER,Nzp, id))
            REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_CTF_XRAY_ZONES_NUMBER) + " argument not present.");
        if (!MD.getValue(MDL_MAGNIFICATION,Ms, id))
            REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_MAGNIFICATION) + " argument not present.");
        if (!MD.getValue(MDL_CTF_SAMPLING_RATE,dxoPSF, id))
            REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_CTF_SAMPLING_RATE) + " argument not present.");
        dxoPSF *= 1e-9;
        if (!MD.getValue(MDL_CTF_SAMPLING_RATE_Z,dzoPSF, id))
            REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_CTF_SAMPLING_RATE_Z) + " argument not present.");
        dzoPSF *= 1e-9;
        if (!MD.getValue(MDL_CTF_LONGITUDINAL_DISPLACEMENT,DeltaZo, id))
            REPORT_ERROR(ERR_ARG_MISSING, MDL::label2Str(MDL_CTF_LONGITUDINAL_DISPLACEMENT) + " argument not present.");
        DeltaZo *= 1e-6;
    }
    else
    {
        FILE *fh_param;
        if ((fh_param = fopen(fn.c_str(), "r")) == nullptr)
            REPORT_ERROR(ERR_IO_NOTOPEN,
                         (std::string)"XmippXROTF::read: There is a problem "
                         "opening the file " + fn);

        try
        {
            lambda = textToFloat(getParameter(fh_param, "lambda", 0, "2.43"))* 1e-9;
            deltaR = textToFloat(getParameter(fh_param, "outer_zone_width", 0, "40")) * 1e-9;
            Nzp = textToFloat(getParameter(fh_param, "zones_number", 0, "560"));
            Ms = textToFloat(getParameter(fh_param, "magnification", 0, "2304"));
            dxoPSF = textToFloat(getParameter(fh_param, "sampling_rate", 0, "1")) *1e-9;
            if (checkParameter(fh_param, "z_sampling_rate"))
                dzoPSF = textToFloat(getParameter(fh_param, "z_sampling_rate", 0)) *1e-9;
            else
                dzoPSF = dxoPSF;
            DeltaZo = textToFloat(getParameter(fh_param, "z_axis_shift", 0, "0")) *1e-6;

            Nox = textToInteger(getParameter(fh_param, "x_dim", 0, "0"));

        }
        catch (XmippError &XE)
        {
            std::cout << XE << std::endl;
            REPORT_ERROR(ERR_IO_NOREAD, (std::string)"There is an error reading " + fn);
        }
        fclose(fh_param);
    }
}

readParams()

void XRayPSF::readParams

(

XmippProgram *

program

)

Definition at line 64 of file psf_xr.cpp.

{
    lambda  = program->getDoubleParam("--lambda")*1e-9;
    deltaR  = program->getDoubleParam("--out_width")*1e-9;
    Nzp     = program->getDoubleParam("--zones");
    Ms      = program->getDoubleParam("--mag");
    dxoPSF  = program->getDoubleParam("--sampling",0)*1e-9;
    dzoPSF  = STR_EQUAL(program->getParam("--sampling", 1), "dxy") ? dxoPSF : program->getDoubleParam("--sampling", 1)*1e-9;
    DeltaZo = program->getDoubleParam("--zshift")*1e-6;
    Nox     = program->getIntParam("--size",0);
    Noy     = STR_EQUAL(program->getParam("--size", 1),"x") ? Nox : program->getIntParam("--size", 1);
    Noz     = STR_EQUAL(program->getParam("--size", 2),"x") ? Nox : program->getIntParam("--size", 2);
    type    = STR_EQUAL(program->getParam("--type"), "zp") ? ANALYTIC_ZP : IDEAL_FRESNEL_LENS;

    verbose = program->getIntParam("-v");
}

reducePSF2Slabs()

void XRayPSF::reducePSF2Slabs

(

double

threshold

)

Calculate the width of the slabs to reduce computing time and the mean PSF for each.

Definition at line 345 of file psf_xr.cpp.

{
    // find max position
    ArrayCoord maxPos;
    MultidimArray<double> &mdaPsfVol = psfVol();

    //MULTIDIM_ARRAY_GENERIC(PSFGen).maxIndex(maxPos);
    mdaPsfVol.maxIndex(maxPos);
    double maxValue = A3D_ELEM(mdaPsfVol,maxPos.z,maxPos.y,maxPos.x);
    threshold *= maxValue;
    slabIndex.clear();

    double maxTemp = maxValue;
    for (int k = 0; k > STARTINGZ(mdaPsfVol); --k)
    {
        double tempV = A3D_ELEM(mdaPsfVol,k,maxPos.y,maxPos.x);
        if ( maxTemp - tempV >  threshold)
        {
            slabIndex.insert(slabIndex.begin(), k);
            maxTemp = tempV;
        }
    }
    slabIndex.insert(slabIndex.begin(), STARTINGZ(mdaPsfVol));

    maxTemp = maxValue;
    for (int k = 0; k < FINISHINGZ(mdaPsfVol); ++k)
    {
        double tempV = A3D_ELEM(mdaPsfVol,k,maxPos.y,maxPos.x);
        if ( maxTemp - tempV >  threshold)
        {
            slabIndex.push_back(k);
            maxTemp = tempV;
        }
    }
    slabIndex.push_back(FINISHINGZ(mdaPsfVol));

} //reducePSF2Slabs

setFocalShift()

void XRayPSF::setFocalShift

(

double

zShift

)

Add focal shift to previously read psf zshift.

Definition at line 273 of file psf_xr.cpp.

{
    DeltaZo -= zShift;
}

show()

void XRayPSF::show

(

)

Show the microscope parameters.

Definition at line 244 of file psf_xr.cpp.

{
    std::cout << *this << std::endl;
}

write()

void XRayPSF::write

(

const FileName &

fn

)

Write to file. An exception is thrown if the file cannot be open.

Definition at line 184 of file psf_xr.cpp.

{
    MetaDataVec MD;
    MD.setColumnFormat(false);
    size_t id = MD.addObject();

    FileName fnPSF = fn.withoutExtension().addExtension("vol");

    MD.setValue(MDL_IMAGE, fnPSF, id);
    String typeS = (type==ANALYTIC_ZP)? "ZP": "ideal";
    MD.setValue(MDL_CTF_XRAY_LENS_TYPE, typeS, id);
    MD.setValue(MDL_CTF_LAMBDA,lambda*1e9, id);
    MD.setValue(MDL_CTF_XRAY_OUTER_ZONE_WIDTH,deltaR*1e9, id);
    MD.setValue(MDL_CTF_XRAY_ZONES_NUMBER,Nzp, id);
    MD.setValue(MDL_MAGNIFICATION,Ms, id);
    MD.setValue(MDL_CTF_SAMPLING_RATE,dxoPSF*1e9, id);
    MD.setValue(MDL_CTF_SAMPLING_RATE_Z,dzoPSF*1e9, id);
    MD.setValue(MDL_CTF_LONGITUDINAL_DISPLACEMENT,DeltaZo*1e6, id);
    std::vector<double> dimV(3);
    dimV[0] = Nox;
    dimV[1] = Noy;
    dimV[2] = Noz;
    MD.setValue(MDL_CTF_DIMENSIONS,dimV, id);

    MD.write(fn.withoutExtension().addExtension("xmd"));
    //    PSF->write(fnPSF);
    psfVol.write(fnPSF);
}

Friends And Related Function Documentation

operator<<

std::ostream& operator<< ( std::ostream &  out, const XRayPSF &  psf  )

friend

Show.

Definition at line 214 of file psf_xr.cpp.

{

    out     << std::endl
    << "--------------------------------------" << std::endl
    << "XrayPSF:: X-Ray Microscope parameters:" << std::endl
    << "--------------------------------------" << std::endl
    << "type =             " << ((psf.type==ANALYTIC_ZP)? "ZP": "ideal") << std::endl
    << "lambda =           " << psf.lambda * 1e9 << " nm" << std::endl
    << "zones_number =     " << psf.Nzp << std::endl
    << "outer_zone_width = " << psf.deltaR * 1e9 << " nm" << std::endl
    << "magnification =    " << psf.Ms << std::endl
    << "sampling_rate =    " << psf.dxo * 1e9 << " nm" << std::endl
    << "z_sampling_rate =  " << psf.dzo * 1e9 << " nm" << std::endl
    << "z_axis_shift =     " << psf.DeltaZo * 1e6 << " um" << std::endl
    << std::endl
    << "focal_length =     " << psf.Flens * 1e3 << " mm" << std::endl
    << "Lens Radius =      " << psf.Rlens * 1e6 << " um" << std::endl
    << "Zo =               " << psf.Zo * 1e3 << " mm" << std::endl
    << "Zi =               " << psf.Zi * 1e3 << " mm" << std::endl
    << "Depth_of_Focus =   " << psf.DoF * 1e6 << " um" << std::endl
    << std::endl
    << "dxi =              " << psf.dxi * 1e6 << " um" << std::endl
    << "dxiMax =           " << psf.dxiMax * 1e6 << " um" << std::endl
    ;

    return out;
}

Member Data Documentation

AdjustType

PsfxrAdjust XRayPSF::AdjustType

Parameters to change image size to avoid Nyquist limit.

Definition at line 171 of file psf_xr.h.

deltaR

double XRayPSF::deltaR

Outermost zone width.

Definition at line 188 of file psf_xr.h.

deltaZMaxX

double XRayPSF::deltaZMaxX

Z limits around Zo in the psf generation due to Nyquist Limit.

Definition at line 165 of file psf_xr.h.

deltaZMaxY

double XRayPSF::deltaZMaxY

Definition at line 166 of file psf_xr.h.

deltaZMinX

double XRayPSF::deltaZMinX

Definition at line 167 of file psf_xr.h.

deltaZMinY

double XRayPSF::deltaZMinY

Definition at line 168 of file psf_xr.h.

DeltaZo

double XRayPSF::DeltaZo

Z axis global shift.

Definition at line 192 of file psf_xr.h.

DoF

double XRayPSF::DoF

Depth of focus. Only for information purposes.

Definition at line 148 of file psf_xr.h.

dxi

double XRayPSF::dxi

Image space XY-plane sampling rate.

Definition at line 199 of file psf_xr.h.

dxiMax

double XRayPSF::dxiMax

Definition at line 159 of file psf_xr.h.

dxl

double XRayPSF::dxl

Pixel size in X-dim in lens plane.

Definition at line 161 of file psf_xr.h.

dxo

double XRayPSF::dxo

object space XY-plane sampling rate

Definition at line 197 of file psf_xr.h.

dxoPSF

double XRayPSF::dxoPSF

object space XY-plane sampling rate of the PSF Volume

Definition at line 207 of file psf_xr.h.

dyl

double XRayPSF::dyl

Pixel size in Y-dim in lens plane.

Definition at line 163 of file psf_xr.h.

dzo

double XRayPSF::dzo

object space Z sampling rate

Definition at line 201 of file psf_xr.h.

dzoPSF

double XRayPSF::dzoPSF

object space Z sampling rate of the PSF Volume

Definition at line 209 of file psf_xr.h.

Flens

double XRayPSF::Flens

Focal length.

Definition at line 184 of file psf_xr.h.

lambda

double XRayPSF::lambda

Lambda of illumination.

Definition at line 180 of file psf_xr.h.

mask

MultidimArray<double>* XRayPSF::mask

Lens shape Mask.

Definition at line 136 of file psf_xr.h.

mode

operMode XRayPSF::mode

Definition at line 114 of file psf_xr.h.

Ms

double XRayPSF::Ms

Magnification.

Definition at line 190 of file psf_xr.h.

Nix

size_t XRayPSF::Nix

Size of the image in image plane, to be rescaled if needed.

Definition at line 203 of file psf_xr.h.

Niy

size_t XRayPSF::Niy

Definition at line 204 of file psf_xr.h.

Nox

size_t XRayPSF::Nox

X size of the input image (object plane size)

Definition at line 152 of file psf_xr.h.

Noy

size_t XRayPSF::Noy

Y size of the input image (object plane size)

Definition at line 154 of file psf_xr.h.

Noz

size_t XRayPSF::Noz

Z size of the input image (object plane size)

Definition at line 156 of file psf_xr.h.

nThr

int XRayPSF::nThr

Definition at line 215 of file psf_xr.h.

Nzp

double XRayPSF::Nzp

Number of zones in zone plate.

Definition at line 186 of file psf_xr.h.

PSF

MultidimArray<double> XRayPSF::PSF

Working PSF with nonlinear zdim whose slices are the mean PSF for slabs.

Definition at line 125 of file psf_xr.h.

psfGen

ImageGeneric XRayPSF::psfGen

Definition at line 121 of file psf_xr.h.

psfVol

Image<double> XRayPSF::psfVol

3D PSF read from file

Definition at line 123 of file psf_xr.h.

pupileSizeMin

double XRayPSF::pupileSizeMin

Minimum diameter size of the microscope pupile in the lens plane, measured in pixels.

Definition at line 173 of file psf_xr.h.

Rlens

double XRayPSF::Rlens

Lens Aperture Radius.

Definition at line 140 of file psf_xr.h.

slabIndex

std::vector<int> XRayPSF::slabIndex

Z positions in the original PSF Volume to determine de slabs.

Definition at line 131 of file psf_xr.h.

slabThr

double XRayPSF::slabThr

Threshold to separate The volume into slabs to use the same PSF.

Definition at line 129 of file psf_xr.h.

T

Matrix2D<double> XRayPSF::T

Definition at line 133 of file psf_xr.h.

type

PsfType XRayPSF::type

Define the selected PSF generation algorithm.

Definition at line 117 of file psf_xr.h.

verbose

int XRayPSF::verbose

Switch to control verbose mode.

Definition at line 212 of file psf_xr.h.

Zi

double XRayPSF::Zi

Object plane.

Image plane (CCD position)

Definition at line 146 of file psf_xr.h.

Zo

double XRayPSF::Zo

Object plane on Focus (Reference)

Definition at line 142 of file psf_xr.h.

---

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

• xmipp/legacy/libraries/data/psf_xr.h
• xmipp/legacy/libraries/data/psf_xr.cpp