Blobs

Collaboration diagram for Blobs:

Classes
struct	blobtype

Macros
#define	blob_val(r, blob)   kaiser_value(r, blob.radius, blob.alpha, blob.order)
#define	blob_proj(r, blob)   kaiser_proj(r, blob.radius, blob.alpha, blob.order)
#define	blob_Fourier_val(w, blob)   kaiser_Fourier_value(w, blob.radius, blob.alpha, blob.order)
#define	blob_mass(blob)   basvolume(blob.radius, blob.alpha, blob.order,3)

Functions
double	kaiser_value (double r, double a, double alpha, int m)
double	kaiser_proj (double r, double a, double alpha, int m)
double	kaiser_Fourier_value (double w, double a, double alpha, int m)
double	basvolume (double a, double alpha, int m, int n)
double	in_zeroarg (int n)
double	inph_zeroarg (int n)
double	i_nph (int n, double x)
double	i_n (int n, double x)
double	blob_freq_zero (struct blobtype b)
double	blob_att (double w, struct blobtype b)
double	blob_ops (double w, struct blobtype b)
double	optimal_CC_grid_relative_size (struct blobtype b)
double	optimal_BCC_grid_relative_size (struct blobtype b)
double	optimal_FCC_grid_relative_size (struct blobtype b)
blobtype	best_blob (double alpha_0, double alpha_F, double inc_alpha, double a_0, double a_F, double inc_a, double w, double *target_att, int target_length=1)
void	footprint_blob (ImageOver &blobprint, const struct blobtype &blob, int istep=50, int normalise=0)
double	sum_blob_Grid (const struct blobtype &blob, const Grid &grid, const Matrix2D< double > *D=nullptr)
void	voxel_volume_shape (const GridVolume &vol_blobs, const struct blobtype &blob, const Matrix2D< double > *D, Matrix1D< int > &corner1, Matrix1D< int > &size)
void	blobs2voxels (const GridVolume &vol_blobs, const struct blobtype &blob, MultidimArray< double > *vol_voxels, const Matrix2D< double > *D=nullptr, int threads=1, int Zdim=0, int Ydim=0, int Xdim=0)
void	blobs2space_coefficients (const GridVolume &vol_blobs, const struct blobtype &blob, MultidimArray< double > *vol_coefs)
void	voxels2blobs (const MultidimArray< double > *vol_voxels, const struct blobtype &blob, GridVolume &vol_blobs, int grid_type, double grid_relative_size, double lambda, const MultidimArray< double > *vol_mask=nullptr, const Matrix2D< double > *D=nullptr, double final_error_change=0.01, int tell=0, double R=-1, int threads=1)
void	ART_voxels2blobs_single_step (GridVolume &vol_in, GridVolume *vol_out, const struct blobtype &blob, const Matrix2D< double > *D, double lambda, MultidimArray< double > *theo_vol, const MultidimArray< double > *read_vol, MultidimArray< double > *corr_vol, const MultidimArray< double > *mask_vol, double &mean_error, double &max_error, int eq_mode=VARTK, int threads=1)

Variables
constexpr int	SHOW_CONVERSION = 1
constexpr int	VARTK = 1
constexpr int	VMAXARTK = 2

Detailed Description

Macro Definition Documentation

blob_Fourier_val

#define blob_Fourier_val	(		w,
			blob
	)		kaiser_Fourier_value(w, blob.radius, blob.alpha, blob.order)

Fourier transform of a blob. This function returns the value of the Fourier transform of the blob at a given frequency (w). This frequency must be normalized by the sampling rate. For instance, for computing the Fourier Transform of a blob at 1/Ts (Ts in Amstrongs) you must provide the frequency Tm/Ts, where Tm is the sampling rate.

The Fourier Transform can be computed only for blobs with m=2 or m=0.

Definition at line 174 of file blobs.h.

blob_mass

#define blob_mass

(

blob

)

basvolume(blob.radius, blob.alpha, blob.order,3)

Formula for a volume integral of a blob (n is the blob dimension)

Definition at line 181 of file blobs.h.

blob_proj

#define blob_proj	(		r,
			blob
	)		kaiser_proj(r, blob.radius, blob.alpha, blob.order)

Blob projection. This function returns the value of the blob line integral through a straight line which passes at a distance 'r' (in Universal System units) from the center of the blob. Remember that a blob is spherically symmetrycal so the only parameter to know this blob line integral is its distance to the center of the blob. It doesn't matter if this distance is larger than the real blob spatial extension, in this case the function returns 0. \ Ex:

struct blobtype blob; blob.radius = 2; blob.order = 2; blob.alpha = 3.6;
Matrix1D<double> v=vectorR3(1,1,1);
std::cout << "Blob line integral through (1,1,1) = " << blob_proj(v.mod(),blob)
     << std::endl;

Definition at line 161 of file blobs.h.

blob_val

#define blob_val	(		r,
			blob
	)		kaiser_value(r, blob.radius, blob.alpha, blob.order)

Blob value. This function returns the value of a blob at a given distance from its center (in Universal System units). The distance must be always positive. Remember that a blob is spherically symmetrycal so the only parameter to know the blob value at a point is its distance to the center of the blob. It doesn't matter if this distance is larger than the real blob spatial extension, in this case the function returns 0 as blob value. \ Ex:

struct blobtype blob; blob.radius = 2; blob.order = 2; blob.alpha = 3.6;
Matrix1D<double> v=vectorR3(1,1,1);
std::cout << "Blob value at (1,1,1) = " << blob_val(v.mod(),blob) << std::endl;

Definition at line 139 of file blobs.h.

Function Documentation

ART_voxels2blobs_single_step()

void ART_voxels2blobs_single_step	(	GridVolume &	vol_in,
		GridVolume *	vol_out,
		const struct blobtype &	blob,
		const Matrix2D< double > *	D,
		double	lambda,
		MultidimArray< double > *	theo_vol,
		const MultidimArray< double > *	read_vol,
		MultidimArray< double > *	corr_vol,
		const MultidimArray< double > *	mask_vol,
		double &	mean_error,
		double &	max_error,
		int	eq_mode = VARTK,
		int	threads = 1
	)

Voxels —> Blobs, single step. This is the basic step of the voxels to blobs conversion. It applies an ART step to the blob volume and produces its output in a different (or the same) volume. You must provide the blob structure, lambda and deformation matrix (typical for crystals). The theoretical volume and correction volume are resized inside and are output volumes meaning the conversion from blobs to voxels before this step and the correction applied to the blob volume. The read volume is the one to reproduce with blobs and the mask volume is used to select only one region to adjust.

If the mask is NULL then it is assumed that all voxels in the input voxel volume must be adjusted by the blobs. If the input voxel volume is NULL the it is assumed that it is equal to 0 and only those corresponding to the 1 positions within the mask must be adjusted. An exception is thrown if the mask and voxels volume are both NULL.

Then mean and maximum error committed for each blob are returned.

Valid eq_modes are \VARTK: ART by blocks for volumes. \VMAXARTK: update with the maximum update.

Definition at line 1094 of file blobs.cpp.

{

    // Resize output volumes ................................................
    if (read_vol != nullptr)
    {
        (*theo_vol).initZeros(*read_vol);
    }
    else if (mask_vol != nullptr)
    {
        (*theo_vol).initZeros(*mask_vol);
    }
    else
    {
        REPORT_ERROR(ERR_MATRIX_EMPTY,
                     "ART_voxels2blobs_single_step: Mask and voxel volumes are empty");
    }
    (*corr_vol).initZeros(*theo_vol);

    auto * th_ids = (pthread_t *)malloc( threads * sizeof( pthread_t));
    auto * threads_d = (ThreadBlobsToVoxels *) malloc ( threads * sizeof( ThreadBlobsToVoxels ) );

    // Translate actual blob volume to voxels ...............................
    for (size_t i = 0; i < vol_in.VolumesNo(); i++)
    {
        int min_distance = (int)ceil((2*(vol_in.grid(i)).relative_size ) / blob.radius ) + 1;

        slices_status = new int[(int)(ZZ(vol_in.grid(i).highest)-ZZ(vol_in.grid(i).lowest)+1)]();
        slices_processed = 0;

        for( int c = 0 ; c < threads ; c++ )
        {
            threads_d[c].vol_blobs = &(vol_in(i)());
            threads_d[c].grid = &(vol_in.grid(i));
            threads_d[c].blob = &blob;
            threads_d[c].vol_voxels = theo_vol;
            threads_d[c].D = D;
            threads_d[c].istep = 50;
            threads_d[c].vol_corr = corr_vol;
            threads_d[c].vol_mask = mask_vol;
            threads_d[c].FORW = true;
            threads_d[c].eq_mode = eq_mode;
            threads_d[c].thread_id = c;
            threads_d[c].threads_num = threads;
            threads_d[c].min_separation = min_distance;

            pthread_create( (th_ids+c), nullptr, blobs2voxels_SimpleGrid, (void *)(threads_d+c) );
        }

        // Wait for threads to finish
        for( int c = 0 ; c < threads ; c++ )
        {
            pthread_join(*(th_ids+c),nullptr);
        }

        delete[] slices_status;
        //        blobs2voxels_SimpleGrid(vol_in(i)(), vol_in.grid(i), blob, theo_vol, D,
        //                                50, corr_vol, mask_vol, FORWARD, eq_mode);
#ifdef DEBUG

        std::cout << "Blob grid no " << i << " stats: ";
        vol_in(i)().printStats();
        std::cout << std::endl;
        std::cout << "So far vol stats: ";
        (*theo_vol)().printStats();
        std::cout << std::endl;
#endif

    }

    // Now normalise the resulting volume ..................................
    double norm = sum_blob_Grid(blob, vol_in.grid(), D); // Aqui tambien hay que multiplicar ****!!!!
    FOR_ALL_ELEMENTS_IN_ARRAY3D(*theo_vol)
    A3D_ELEM(*theo_vol, k, i, j) /= norm;

#ifdef DEBUG

    Image<double> save, save2;
    save() = *theo_vol;
    save.write("PPPtheovol.vol");
    std::cout << "Theo stats:";
    save().printStats();
    std::cout << std::endl;
    save() = *corr_vol;
    save.write("PPPcorr2vol.vol");
    save2().resize(save());
#endif

    // Compute differences ..................................................
    mean_error = 0;
    double read_val;
    if (read_vol != nullptr)
        read_val = A3D_ELEM(*read_vol, 0, 0, 0);
    else
        read_val = 0;
    max_error = ABS(read_val - A3D_ELEM(*theo_vol, 0, 0, 0));

    double diff;
    int N = 0;
    FOR_ALL_ELEMENTS_IN_ARRAY3D(*theo_vol)
    {
        // Compute difference volume and error
        if (read_vol != nullptr)
            read_val = A3D_ELEM(*read_vol, k, i, j);
        else
            read_val = 0;

        if (mask_vol == nullptr)
        {
            diff = read_val - A3D_ELEM(*theo_vol, k, i, j);
            N++;
        }
        else
            if (A3D_ELEM(*mask_vol, k, i, j) == 1)
            {
                diff = read_val - A3D_ELEM(*theo_vol, k, i, j);
                N++;
            }
            else
                diff = 0;

        max_error = XMIPP_MAX(max_error, ABS(diff));
        mean_error += diff * diff;
#ifdef DEBUG

        save(k, i, j) = diff;
        save2(k, i, j) = read_val;
#endif


        // Compute the correction volume
        if (ABS(A3D_ELEM(*corr_vol, k, i, j)) < 1e-2)
            A3D_ELEM(*corr_vol, k, i, j) = SGN(A3D_ELEM(*corr_vol, k, i, j));
        A3D_ELEM(*corr_vol, k, i, j) =
            lambda * diff / A3D_ELEM(*corr_vol, k, i, j);
    }
#ifdef DEBUG
    save.write("PPPdiffvol.vol");
    std::cout << "Diff stats:";
    save().printStats();
    std::cout << std::endl;
    save2.write("PPPreadvol.vol");
    std::cout << "Read stats:";
    save2().printStats();
    std::cout << std::endl;
    save() = *corr_vol;
    save.write("PPPcorrvol.vol");
    std::cout << "Corr stats:";
    save().printStats();
    std::cout << std::endl;
#endif

    mean_error /= XMIPP_MAX(N, 1); // At worst, divided by 1

    // Backprojection of correction volume ..................................
    for (size_t i = 0; i < vol_in.VolumesNo(); i++)
    {
        slices_status = new int[(int)(ZZ(vol_out->grid(i).highest)-ZZ(vol_out->grid(i).lowest)+1)]();
        slices_processed = 0;

        for( int c = 0 ; c < threads ; c++ )
        {
            threads_d[c].vol_blobs = &((*vol_out)(i)());
            threads_d[c].grid = &(vol_out->grid(i));
            threads_d[c].blob = &blob;
            threads_d[c].vol_voxels = theo_vol;
            threads_d[c].D = D;
            threads_d[c].istep = 50;
            threads_d[c].vol_corr = corr_vol;
            threads_d[c].vol_mask = mask_vol;
            threads_d[c].FORW = false;
            threads_d[c].eq_mode = eq_mode;
            threads_d[c].thread_id = c;
            threads_d[c].threads_num = threads;
            threads_d[c].min_separation = 1;

            pthread_create( (th_ids+c), nullptr, blobs2voxels_SimpleGrid, (void *)(threads_d+c) );
        }

        // Wait for threads to finish
        for( int c = 0 ; c < threads ; c++ )
        {
            pthread_join(*(th_ids+c), nullptr);
        }
        delete[] slices_status;
        //        blobs2voxels_SimpleGrid((*vol_out)(i)(), (*vol_out).grid(i), blob,
        //                                theo_vol, D, 50, corr_vol, mask_vol, BACKWARD, eq_mode);
#ifdef DEBUG

        std::cout << "Blob grid no " << i << " stats: ";
        vol_in(i)().printStats();
        std::cout << std::endl;
#endif

    }
    free(threads_d);
    free(th_ids);
#ifdef DEBUG
    char c;
    std::cout << "Press any key to continue\n";
    std::cin >> c;
#endif
}

basvolume()

double basvolume	(	double	a,
		double	alpha,
		int	m,
		int	n
	)

Function actually computing the blob integral

Definition at line 215 of file blobs.cpp.

{
    double  hn;
    double tpi;
    double v;
    hn = 0.5 * n;
    tpi = 2.0 * PI;

    if (alpha == 0.0)
    {
        if ((n / 2)*2 == n)           /* n even                               */
            v = pow(tpi, hn) * in_zeroarg(n / 2 + m) / in_zeroarg(m);
        else                        /* n odd                                */
            v = pow(tpi, hn) * inph_zeroarg(n / 2 + m) / in_zeroarg(m);

    }
    else
    {                        /* alpha > 0.0                          */
        if ((n / 2)*2 == n)           /* n even                               */
            v = pow(tpi / alpha, hn) * i_n(n / 2 + m, alpha) / i_n(m, alpha);
        else                        /* n odd                                */
            v = pow(tpi / alpha, hn) * i_nph(n / 2 + m, alpha) / i_n(m, alpha);
    }

    return v * pow(a, (double)n);
}

best_blob()

blobtype best_blob	(	double	alpha_0,
		double	alpha_F,
		double	inc_alpha,
		double	a_0,
		double	a_F,
		double	inc_a,
		double	w,
		double *	target_att,
		int	target_length = 1
	)

Choose the best blob within a region. You must select a resolution limit, the maximum attenuation allowed at that frequency and then the blob with less computational cost is returned. m=2 always. The resolution criterion is given by w (remind that this w=Tm*w(cont)).

If there is no blob meeting the specification then alpha=a=-1.

Definition at line 363 of file blobs.cpp.

{
    blobtype retval;
    retval.order = 2;
    int alpha_size = FLOOR((alpha_F - alpha_0) / inc_alpha + 1);
    int a_size = FLOOR((a_F - a_0) / inc_a + 1);
    Image<double> att;
    Image<double> ops;
    att().resize(alpha_size, a_size);
    ops().resize(alpha_size, a_size);
    int i;
    int j;
    double a;
    double alpha;
    double best_a = -1;
    double best_alpha = -1;
    double best_ops = 1e10;
    double best_att = 0;
    for (i = 0, alpha = alpha_0; i < alpha_size; alpha += inc_alpha, i++)
        for (j = 0, a = a_0; j < a_size; a += inc_a, j++)
        {
            retval.radius = a;
            retval.alpha = alpha;
            A2D_ELEM(att(), i, j) = blob_att(w, retval);
            A2D_ELEM(ops(), i, j) = blob_ops(w, retval);
            if (j > 0)
                for (int n = target_length - 1; n >= 0; n--)
                    if (A2D_ELEM(att(), i, j - 1) > target_att[n] &&
                        A2D_ELEM(att(), i, j) < target_att[n])
                    {
                        A2D_ELEM(att(), i, j) = 0;
                        if (A2D_ELEM(ops(), i, j - 1) < best_ops &&
                            A2D_ELEM(att(), i, j - 1) >= best_att)
                        {
                            best_alpha = alpha;
                            best_a = a - inc_a;
                            best_ops = A2D_ELEM(ops(), i, j - 1);
                            best_att = target_att[n];
                        }
                    }
        }
#ifdef DEBUG
    att.write((std::string)"att" + floatToString(w) + ".xmp");
    ops.write((std::string)"ops" + floatToString(w) + ".xmp");
#endif

    retval.radius = best_a;
    retval.alpha = best_alpha;
    return retval;
}

blob_att()

double blob_att	(	double	w,
		struct blobtype	b
	)

Attenuation of a blob. The Fourier transform of the blob at w is the Fourier transform at w=0 multiplied by the attenuation. This is the value returned. Remind that the frequency must be normalized by the sampling rate. Ie, Tm*w(cont)

Definition at line 336 of file blobs.cpp.

{
    return blob_Fourier_val(w, b) / blob_Fourier_val(0, b);
}

blob_freq_zero()

double blob_freq_zero

(

struct blobtype

b

)

Blob pole. This is the normalized frequency at which the blob goes to 0.

Definition at line 330 of file blobs.cpp.

{
    return sqrt(b.alpha*b.alpha + 6.9879*6.9879) / (2*PI*b.radius);
}

blob_ops()

double blob_ops	(	double	w,
		struct blobtype	b
	)

Number of operations for a blob. This is a number proportional to the number of operations that ART would need to make a reconstruction with this blob.

Definition at line 342 of file blobs.cpp.

{
    return pow(b.alpha*b.alpha + 6.9879*6.9879, 1.5) / b.radius;
}

blobs2space_coefficients()

void blobs2space_coefficients	(	const GridVolume &	vol_blobs,
		const struct blobtype &	blob,
		MultidimArray< double > *	vol_coefs
	)

Blob coefficients as a voxel volume. This function returns a volume with the blob coefficients in their right position in space. The voxel size is g/2

Definition at line 1029 of file blobs.cpp.

{

    // Compute vol_coefs shape
    Matrix1D<int> corner1;
    Matrix1D<int> size;
    voxel_volume_shape(vol_blobs, blob, nullptr, corner1, size);
    double g_2 = vol_blobs.grid(0).relative_size / 2;
    XX(corner1) = (int)(FLOOR(XX(corner1)) / g_2);
    XX(size) = (int)(CEIL(XX(size)) / g_2);
    YY(corner1) = (int)(FLOOR(YY(corner1)) / g_2);
    YY(size) = (int)(CEIL(YY(size)) / g_2);
    ZZ(corner1) = (int)(FLOOR(ZZ(corner1)) / g_2);
    ZZ(size) = (int)(CEIL(ZZ(size)) / g_2);
    (*vol_coefs).initZeros(ZZ(size), YY(size), XX(size));
    STARTINGX(*vol_coefs) = XX(corner1);
    STARTINGY(*vol_coefs) = YY(corner1);
    STARTINGZ(*vol_coefs) = ZZ(corner1);

    // Set all blob coefficients at the right position
    for (size_t n = 0; n < vol_blobs.VolumesNo(); n++)
    {
        auto ZZ_lowest = (int)ZZ(vol_blobs.grid(n).lowest);
        auto YY_lowest = (int)YY(vol_blobs.grid(n).lowest);
        auto XX_lowest = (int)XX(vol_blobs.grid(n).lowest);
        auto ZZ_highest = (int)ZZ(vol_blobs.grid(n).highest);
        auto YY_highest = (int)YY(vol_blobs.grid(n).highest);
        auto XX_highest = (int)XX(vol_blobs.grid(n).highest);
        for (int k = ZZ_lowest; k <= ZZ_highest; k++)
            for (int i = YY_lowest; i <= YY_highest; i++)
                for (int j = XX_lowest; j <= XX_highest; j++)
                {
                    Matrix1D<double> grid_index(3);
                    Matrix1D<double> univ_position(3),
                    coef_position(3);
                    VECTOR_R3(grid_index, j, i, k);
                    vol_blobs.grid(n).grid2universe(grid_index, univ_position);
                    V3_BY_CT(coef_position, univ_position, 1.0 / g_2);
                    A3D_ELEM((*vol_coefs),
                             (int)ZZ(coef_position),
                             (int)YY(coef_position),
                             (int)XX(coef_position) ) = A3D_ELEM(vol_blobs(n)(), k, i, j);
#ifdef DEBUG

                    std::cout << "Blob value at (" << j << "," << i << ","
                    << k << ") (" << XX(univ_position)
                    << "," << YY(univ_position) << ","
                    << ZZ(univ_position) << ") ("
                    << XX(coef_position) << "," << YY(coef_position) << ","
                    << ZZ(coef_position) << ") --> "
                    << A3D_ELEM(vol_blobs(n)(), (int)k, (int)i, (int)j) << std::endl;
#endif

                }
    }
}

blobs2voxels()

void blobs2voxels	(	const GridVolume &	vol_blobs,
		const struct blobtype &	blob,
		MultidimArray< double > *	vol_voxels,
		const Matrix2D< double > *	D = nullptr,
		int	threads = 1,
		int	Zdim = 0,
		int	Ydim = 0,
		int	Xdim = 0
	)

Blobs —> Voxels. The voxel size is defined between two coordinates (Gcorner1 and Gcorner2) which accomplish

XX(Gcorner1)=MIN(XX(SGcorner1(i)));
YY(Gcorner1)=MIN(YY(SGcorner1(i)));
ZZ(Gcorner1)=MIN(ZZ(SGcorner1(i)));

XX(Gcorner2)=MAX(XX(SGcorner2(i)));
YY(Gcorner2)=MAX(YY(SGcorner2(i)));
ZZ(Gcorner2)=MAX(ZZ(SGcorner2(i)));

where SGcorner1(i) and SGcorner2(i) are the lowest and highest corners of each subgrid. These corners are expressed in Universal coordinates.

This quantity is then enlarged to be an integer voxel position and to take into account that the former computed Gcorner1 and 2 are the furthest centers of blobs, ie, there will be voxels even further affected by these blobs.

Gcorner1 = CEILnD (Gcorner1 - blob.radius);
Gcorner2 = FLOORnD(Gcorner2 + blob.radius);

However, you might give a size, usually set to 0, i.e., no external size. If no size is provided a size is produced such that all blob centers fit into the output volume.

D is a 3x3 matrix specifying a volume deformation. It means that the sample at logical position (i,j,k) is really placed at space position D*(i,j,k)'.

Definition at line 926 of file blobs.cpp.

{

    // Resize and set starting corner .......................................
    if (Zdim == 0 || Ydim == 0 || Xdim == 0)
    {
        Matrix1D<int> size;
        Matrix1D<int> corner;
        voxel_volume_shape(vol_blobs, blob, D, corner, size);
        (*vol_voxels).initZeros(ZZ(size), YY(size), XX(size));
        STARTINGX(*vol_voxels) = XX(corner);
        STARTINGY(*vol_voxels) = YY(corner);
        STARTINGZ(*vol_voxels) = ZZ(corner);
    }
    else
    {
        (*vol_voxels).initZeros(Zdim, Ydim, Xdim);
        (*vol_voxels).setXmippOrigin();
    }

    auto * th_ids = new pthread_t [threads];
    auto * threads_d = new ThreadBlobsToVoxels [threads];

    // Convert each subvolume ...............................................
    for (size_t i = 0; i < vol_blobs.VolumesNo(); i++)
    {
        int min_distance = (int)ceil((2*(vol_blobs.grid(i)).relative_size ) / blob.radius ) + 1;
        slices_status = new int[(int)(ZZ(vol_blobs.grid(i).highest)-ZZ(vol_blobs.grid(i).lowest)+1)]();
        slices_processed = 0;

        for( int c = 0 ; c < threads ; c++ )
        {
            threads_d[c].vol_blobs = &(vol_blobs(i)());
            threads_d[c].grid = &(vol_blobs.grid(i));
            threads_d[c].blob = &blob;
            threads_d[c].vol_voxels = vol_voxels;
            threads_d[c].D = D;
            threads_d[c].istep = 50;
            threads_d[c].vol_corr = nullptr;
            threads_d[c].vol_mask = nullptr;
            threads_d[c].FORW = true;
            threads_d[c].eq_mode = VARTK;
            threads_d[c].thread_id = c;
            threads_d[c].threads_num = threads;
            threads_d[c].min_separation = min_distance;

            pthread_create( (th_ids+c), nullptr, blobs2voxels_SimpleGrid, (void *)(threads_d+c) );
        }

        // Wait for threads to finish
        for( int c = 0 ; c < threads ; c++ )
        {
            pthread_join(*(th_ids+c),nullptr);
        }

#ifdef DEBUG
        std::cout << "Blob grid no " << i << " stats: ";
        vol_blobs(i)().printStats();
        std::cout << std::endl;
        std::cout << "So far vol stats: ";
        (*vol_voxels).printStats();
        std::cout << std::endl;
        Image<double> save;
        save() = *vol_voxels;
        save.write((std::string)"PPPvoxels" + integerToString(i));
#endif

        delete[] slices_status;
    }

    // Now normalise the resulting volume ..................................
    double inorm = 1.0 / sum_blob_Grid(blob, vol_blobs.grid(), D); // Aqui tambien hay que multiplicar ****!!!!
    FOR_ALL_ELEMENTS_IN_ARRAY3D(*vol_voxels)
    A3D_ELEM(*vol_voxels, k, i, j) *= inorm;

    // Set voxels outside interest region to minimum value .................
    double R = vol_blobs.grid(0).get_interest_radius();
    if (R != -1)
    {
        double R2 = (R - 6) * (R - 6);

        // Compute minimum value within sphere
        double min_val = A3D_ELEM(*vol_voxels, 0, 0, 0);
        FOR_ALL_ELEMENTS_IN_ARRAY3D(*vol_voxels)
        if (j*j + i*i + k*k <= R2 - 4)
            min_val = XMIPP_MIN(min_val, A3D_ELEM(*vol_voxels, k, i, j));

        // Substitute minimum value
        R2 = (R - 2) * (R - 2);
        FOR_ALL_ELEMENTS_IN_ARRAY3D(*vol_voxels)
        if (j*j + i*i + k*k >= R2)
            A3D_ELEM(*vol_voxels, k, i, j) = min_val;
    }

    delete[] th_ids;
    delete[] threads_d;
}

footprint_blob()

void footprint_blob	(	ImageOver &	blobprint,
		const struct blobtype &	blob,
		int	istep = 50,
		int	normalise = 0
	)

Footprint of a blob. In the actual implementation of the 3D reconstruction one of the most important things to save time is to have a precomputed blob projection. As it is spherically symmetrical all projections from different directions of the same blob happen to be the same. This function makes this precalculation storing it in an oversampled image. This is done so because the blob footprint might not be centered with respect to the universal grid needing so, the footprint at non-integer positions. As parameters to this function you may give the sampling rate in each direction, and an optional normalisation (usually disabled) at the end of the computation which divides the whole image by the sum of the whole image. \ Ex:

ImageOver             blobprint;
struct blobtype       blob;
blob.radius = 1.7;
blob.order  = 2;
blob.alpha  = 3.6;
footprint_blob(blobprint, blob);

As the blob radius is 1.7, the function will create a normalised footprint of logical corners (-2,-2) to (2,2) (in the universal coord. system) (this size is the minimum integer number of pixels which contain entirely the blob), with 50 samples each pixel. The real size of the footprint is (550)x(550).

Definition at line 417 of file blobs.cpp.

{
    // Resize output image and redefine the origin of it
    int footmax = CEIL(blob.radius);
    blobprint.init(-footmax, footmax, istep, -footmax, footmax, istep);

    // Run for Imge class indexes
    for (int i = STARTINGY(blobprint()); i <= FINISHINGY(blobprint()); i++)
        for (int j = STARTINGX(blobprint()); j <= FINISHINGX(blobprint()); j++)
        {
            // Compute oversampled index and blob value
            double vi;
            double ui;
            IMG2OVER(blobprint, i, j, vi, ui);
            double r = sqrt(vi * vi + ui * ui);
            IMGPIXEL(blobprint, i, j) = blob_proj(r, blob);
        }

    // Adjust the footprint structure
    if (normalise)
        blobprint() /= blobprint().sum();
}

i_n()

double i_n	(	int	n,
		double	x
	)

Bessel function I_n (x), n = 0, 1, 2, ... Use ONLY for small values of n

Definition at line 244 of file blobs.cpp.

{
    int i;
    double i_ns1;
    double i_n;
    double i_np1;

    if (n == 0)
        return bessi0(x);
    if (n == 1)
        return bessi1(x);
    if (x == 0.0)
        return 0.0;
    i_ns1 = bessi0(x);
    i_n   = bessi1(x);

    for (i = 1; i < n; i++)
    {
        i_np1 = i_ns1 - (2 * i) / x * i_n;
        i_ns1 = i_n;
        i_n   = i_np1;
    }
    return i_n;
}

i_nph()

double i_nph	(	int	n,
		double	x
	)

Bessel function I_(n+1/2) (x), n = 0, 1, 2, ...

Definition at line 270 of file blobs.cpp.

{
    int i;
    double r2dpix;
    double i_ns1;
    double i_n;
    double i_np1;

    if (x == 0.0)
        return 0.0;
    r2dpix = sqrt(2.0 / (PI * x));
    i_ns1 = r2dpix * cosh(x);
    i_n   = r2dpix * sinh(x);

    for (i = 1; i <= n; i++)
    {
        i_np1 = i_ns1 - (2 * i - 1) / x * i_n;
        i_ns1 = i_n;
        i_n   = i_np1;
    }
    return i_n;
}

in_zeroarg()

double in_zeroarg

(

int

n

)

Limit (z->0) of (1/z)^n I_n(z) (needed by basvolume)

Definition at line 294 of file blobs.cpp.

{
    int i;
    double fact;
    fact = 1.0;
    for (i = 1; i <= n; i++)
    {
        fact *= 0.5 / i;
    }
    return fact;
}

inph_zeroarg()

double inph_zeroarg

(

int

n

)

Limit (z->0) of (1/z)^(n+1/2) I_(n+1/2) (z) (needed by basvolume)

Definition at line 307 of file blobs.cpp.

{
    int i;
    double fact;
    fact = 1.0;
    for (i = 1; i <= n; i++)
    {
        fact *= 1.0 / (2 * i + 1.0);
    }
    return fact*sqrt(2.0 / PI);
}

kaiser_Fourier_value()

double kaiser_Fourier_value	(	double	w,
		double	a,
		double	alpha,
		int	m
	)

Function actually computing the blob Fourier transform.

Definition at line 144 of file blobs.cpp.

{
    if (m != 2 && m !=0)
        REPORT_ERROR(ERR_VALUE_INCORRECT, "m out of range in kaiser_Fourier_value()");
    double sigma = sqrt(ABS(alpha * alpha - (2. * PI * a * w) * (2. * PI * a * w)));
    if (m == 2)
    {
        if (2.*PI*a*w > alpha)
            return  pow(2.*PI, 3. / 2.)*pow(a, 3.)*pow(alpha, 2.)*bessj3_5(sigma)
                    / (bessi0(alpha)*pow(sigma, 3.5));
        else
            return  pow(2.*PI, 3. / 2.)*pow(a, 3.)*pow(alpha, 2.)*bessi3_5(sigma)
                    / (bessi0(alpha)*pow(sigma, 3.5));
    }
    else if (m == 0)
    {
        if (2*PI*a*w > alpha)
            return  pow(2.*PI, 3. / 2.)*pow(a, 3)*bessj1_5(sigma)
                    / (bessi0(alpha)*pow(sigma, 1.5));
        else
            return  pow(2.*PI, 3. / 2.)*pow(a, 3)*bessi1_5(sigma)
                    / (bessi0(alpha)*pow(sigma, 1.5));
    }
    else
        REPORT_ERROR(ERR_ARG_INCORRECT,"Invalid blob order");
}

kaiser_proj()

double kaiser_proj	(	double	r,
		double	a,
		double	alpha,
		int	m
	)

Function actually computing the blob projection.

Definition at line 94 of file blobs.cpp.

{
    double sda;
    double sdas;
    double w;
    double arg;
    double p;

    sda = s / a;
    sdas = sda * sda;
    w = 1.0 - sdas;
    if (w > 1.0e-10)
    {
        arg = alpha * sqrt(w);
        if (m == 0)
        {
            if (alpha == 0.0)
                p = 2.0 * a * sqrt(w);
            else
                p = (2.0 * a / alpha) * sinh(arg) / bessi0(alpha);

        }
        else if (m == 1)
        {
            if (alpha == 0.0)
                p = 2.0 * a * w * sqrt(w) * (2.0 / 3.0);
            else
                p = (2.0 * a / alpha) * sqrt(w) * (cosh(arg) - sinh(arg) / arg)
                    / bessi1(alpha);

        }
        else if (m == 2)
        {
            if (alpha == 0.0)
                p = 2.0 * a * w * w * sqrt(w) * (8.0 / 15.0);
            else
                p = (2.0 * a / alpha) * w *
                    ((3.0 / (arg * arg) + 1.0) * sinh(arg) - (3.0 / arg) * cosh(arg)) / bessi2(alpha);
        }
        else
            REPORT_ERROR(ERR_VALUE_INCORRECT, "m out of range in kaiser_proj()");

    }
    else
        p = 0.0;

    return p;
}

kaiser_value()

double kaiser_value	(	double	r,
		double	a,
		double	alpha,
		int	m
	)

Function actually computing the blob value.

Definition at line 37 of file blobs.cpp.

{
    double rda;
    double rdas;
    double arg;
    double w;

    rda = r / a;
    if (rda <= 1.0)
    {
        rdas = rda * rda;
        arg = alpha * sqrt(1.0 - rdas);
        if (m == 0)
        {
            w = bessi0(arg) / bessi0(alpha);
        }
        else if (m == 1)
        {
            w = sqrt (1.0 - rdas);
            if (alpha != 0.0)
                w *= bessi1(arg) / bessi1(alpha);
        }
        else if (m == 2)
        {
            w = sqrt (1.0 - rdas);
            w = w * w;
            if (alpha != 0.0)
                w *= bessi2(arg) / bessi2(alpha);
        }
        else if (m == 3)
        {
            w = sqrt (1.0 - rdas);
            w = w * w * w;
            if (alpha != 0.0)
                w *= bessi3(arg) / bessi3(alpha);
        }
        else if (m == 4)
        {
            w = sqrt (1.0 - rdas);
            w = w * w * w *w;
            if (alpha != 0.0)
                w *= bessi4(arg) / bessi4(alpha);
        }
        else
            REPORT_ERROR(ERR_VALUE_INCORRECT, "m out of range in kaiser_value()");

    }
    else
        w = 0.0;

    return w;
}

optimal_BCC_grid_relative_size()

double optimal_BCC_grid_relative_size

(

struct blobtype

b

)

Optimal BCC grid spastd::cing. This function returns the optimal grid relative size for the blob selected.

Definition at line 352 of file blobs.cpp.

{
    return sqrt(8.0) / (2*blob_freq_zero(b));
}

optimal_CC_grid_relative_size()

double optimal_CC_grid_relative_size

(

struct blobtype

b

)

Optimal CC grid spastd::cing. This function returns the optimal grid relative size for the blob selected.

Definition at line 348 of file blobs.cpp.

{
    return 1 / blob_freq_zero(b);
}

optimal_FCC_grid_relative_size()

double optimal_FCC_grid_relative_size

(

struct blobtype

b

)

Optimal FCC grid spastd::cing. This function returns the optimal grid relative size for the blob selected.

Definition at line 356 of file blobs.cpp.

{
    return sqrt(27.0) / (4*blob_freq_zero(b));
}

sum_blob_Grid()

double sum_blob_Grid	(	const struct blobtype &	blob,
		const Grid &	grid,
		const Matrix2D< double > *	D = nullptr
	)

Sum of a single blob over a grid. As a normalisation factor, the sum of the blob values over a given grid is needed. What this function does is put a blob at coordinate (0,0,0) and sums all the blob values at points of the grid which are inside the blob. It doesn't matter if the grid is compound or not. \ Ex:

// Blob definition
struct blobtype       blob;
blob.radius = 2;
blob.order  = 2;
blob.alpha  = 3.6;

// Grid definition
Grid BCCgrid;
BCCgrid=BCC_grid(1.41,vectorR3(-5,-5,-5),vector_R3( 5, 5, 5));

std::cout << "The sum of a single blob over the grid is " <<
     << sum_blob_grid(blob, BCCgrid);

D is a 3x3 matrix specifying a volume deformation. It means that the sample at logical position (i,j,k) is really placed at space position D*(i,j,k)'.

Definition at line 320 of file blobs.cpp.

{
    double sum = 0;
    for (size_t i = 0; i < grid.GridsNo(); i++)
        sum += sum_blob_SimpleGrid(blob, grid(i), D);
    return sum;
}

voxel_volume_shape()

void voxel_volume_shape	(	const GridVolume &	vol_blobs,
		const struct blobtype &	blob,
		const Matrix2D< double > *	D,
		Matrix1D< int > &	corner1,
		Matrix1D< int > &	size
	)

Voxel shape for a blob volume. Given a blob volume this function returns the logical origin and size for the minimum voxel volume which holds it. See blobs2voxels for an explanation of the limit and V parameters.

Definition at line 851 of file blobs.cpp.

{
    Matrix1D<double>  Gcorner1(3);
    Matrix1D<double>  Gcorner2(3);     // lowest and highest coord.

    corner1.resize(3);
    size.resize(3);

    // Look for the lowest and highest volume coordinate
    vol_blobs.grid().voxel_corners(Gcorner1, Gcorner2, D);

    // Add blob radius in each direction, and find the furthest integer
    // samples => compute size
    Gcorner1 -= blob.radius;
    Gcorner1.selfCEIL();
    Gcorner2 += blob.radius;
    Gcorner2.selfFLOOR();

    XX(size) = (int)XX(Gcorner2) - (int)XX(Gcorner1) + 1;
    YY(size) = (int)YY(Gcorner2) - (int)YY(Gcorner1) + 1;
    ZZ(size) = (int)ZZ(Gcorner2) - (int)ZZ(Gcorner1) + 1;
#ifdef DEBUG

    std::cout << "Gcorner1  " << Gcorner1.transpose() << std::endl;
    std::cout << "Gcorner2  " << Gcorner2.transpose() << std::endl;
    std::cout << "Size of voxel volume " << (int)ZZ(size) << " x "
    << (int)YY(size) << " x " << (int)XX(size) << std::endl;
#endif

#ifdef NEVER_DEFINED
    // In principle this limitation has been substituted by a direct
    // specification of the output volume size, and should no longer
    // be valid. However it is a beautiful piece of code to be removed
    // already
    if (limit != 0 && XX(size) != limit)
    {
        double diff = ((double)XX(size) - (double)limit) / 2;
        if (diff == (int)diff)
        {
            Gcorner1 += diff;
            Gcorner2 -= diff;
        }
        else
        {
            Gcorner1 += (diff - 0.5);
            Gcorner2 -= (diff + 0.5);
        }
        XX(size) = (int)XX(Gcorner2) - (int)XX(Gcorner1) + 1;
        YY(size) = (int)YY(Gcorner2) - (int)YY(Gcorner1) + 1;
        ZZ(size) = (int)ZZ(Gcorner2) - (int)ZZ(Gcorner1) + 1;
#ifdef DEBUG

        std::cout << "Limiting to " << limit << " diff = " << diff << std::endl;
        std::cout << "New Gcorner1  " << Gcorner1.transpose() << std::endl;
        std::cout << "New Gcorner2  " << Gcorner2.transpose() << std::endl;
#endif

    }
#endif

    typeCast(Gcorner1, corner1);

#ifdef DEBUG

    std::cout << "Final size of voxel volume " << (int)ZZ(size) << " x "
    << (int)YY(size) << " x " << (int)XX(size) << std::endl;
    std::cout << "Corner1= " << corner1.transpose() << std::endl;
#endif
}

voxels2blobs()

void voxels2blobs	(	const MultidimArray< double > *	vol_voxels,
		const struct blobtype &	blob,
		GridVolume &	vol_blobs,
		int	grid_type,
		double	grid_relative_size,
		double	lambda,
		const MultidimArray< double > *	vol_mask = nullptr,
		const Matrix2D< double > *	D = nullptr,
		double	final_error_change = 0.01,
		int	tell = 0,
		double	R = -1,
		int	threads = 1
	)

Voxels —> Blobs. The voxels to blobs procedure is an iterative process resolved by ART as an iterative technique to solve an equation system. The blobs to voxels process is the BASIC way to solve this problem. This applies the whole ART process. You can set a debugging level by setting the flag SHOW_CONVERSION in the tell argument.

If the mask is NULL then it is assumed that all voxels in the input voxel volume must be adjusted by the blobs. If the input voxel volume is NULL the it is assumed that it is equal to 0 and only those corresponding to the 1 positions within the mask must be adjusted.

R is the interest radius. If it is -1 two superimposed grids are created, otherwise a single grid with tilted axis is.

Valid grid types are defined in Some useful grids

Definition at line 1316 of file blobs.cpp.

{
    Image<double> theo_vol;
    Image<double> corr_vol;
    double mean_error;
    double mean_error_1;
    double max_error;
    int it = 1;

    tell = SHOW_CONVERSION;

    // Resize output volume .................................................
    Grid grid_blobs;
    if (R == -1)
    {
        Matrix1D<double> corner1(3);
        Matrix1D<double> corner2(3);
        XX(corner1) = STARTINGX(*vol_voxels);
        YY(corner1) = STARTINGY(*vol_voxels);
        ZZ(corner1) = STARTINGZ(*vol_voxels);
        XX(corner2) = FINISHINGX(*vol_voxels);
        YY(corner2) = FINISHINGY(*vol_voxels);
        ZZ(corner2) = FINISHINGZ(*vol_voxels);

        switch (grid_type)
        {
        case CC:
                        grid_blobs = Create_CC_grid(grid_relative_size, corner1, corner2);
            break;
        case FCC:
                        grid_blobs = Create_FCC_grid(grid_relative_size, corner1, corner2);
            break;
        case BCC:
                        grid_blobs = Create_BCC_grid(grid_relative_size, corner1, corner2);
            break;
        }
    }
    else
{
        switch (grid_type)
        {
        case CC:
                        grid_blobs = Create_CC_grid(grid_relative_size, R);
            break;
        case FCC:
                        grid_blobs = Create_FCC_grid(grid_relative_size, R);
            break;
        case BCC:
                        grid_blobs = Create_BCC_grid(grid_relative_size, R);
            break;
        }
    }
    vol_blobs.adapt_to_grid(grid_blobs);

    // Convert ..............................................................
    std::cout << "Converting voxel volume to blobs\n";
    ART_voxels2blobs_single_step(vol_blobs, &vol_blobs, blob, D, lambda,
                                 &(theo_vol()), vol_voxels,
                                 &(corr_vol()), vol_mask, mean_error_1, max_error, VARTK, threads);
    if (tell && SHOW_CONVERSION)
        std::cout << "   Finished iteration: " << it++
        << " Mean Error= " << mean_error_1
        << " Max_Error= " << max_error
        << std::endl;
    else
        std::cout << "0%";
    std::cout.flush();
#ifdef DEBUG

    theo_vol.write("PPPtheo.vol");
    corr_vol.write("PPPcorr.vol");
    std::cout << "Press any key\n";
    char c;
    std::cin >> c;
#endif

    double change;
    bool end_condition;
    do
{
        ART_voxels2blobs_single_step(vol_blobs, &vol_blobs, blob, D, lambda,
                                     &(theo_vol()), vol_voxels,
                                     &(corr_vol()), vol_mask, mean_error, max_error, VARTK);
#ifdef DEBUG

        theo_vol.write("PPPtheo.vol");
        corr_vol.write("PPPcorr.vol");
        std::cout << "Press any key\n";
        char c;
        std::cin >> c;
#endif

        change = ABS(mean_error - mean_error_1) / mean_error_1;
        mean_error_1 = mean_error;
        if (tell && SHOW_CONVERSION)
            std::cout << "   Finished iteration: " << it++
            << " Mean Error= " << mean_error
            << " Max_Error= " << max_error
            << std::endl;
        else
        {
            printf("\r");
            std::cout << 100 - change << "%";
        }
        std::cout.flush();
        end_condition = change <= final_error_change;
    }
    while (!end_condition);
    std::cout << std::endl;
}

Variable Documentation

SHOW_CONVERSION

constexpr int SHOW_CONVERSION = 1

Definition at line 353 of file blobs.h.

VARTK

constexpr int VARTK = 1

Definition at line 378 of file blobs.h.

VMAXARTK

constexpr int VMAXARTK = 2

Definition at line 379 of file blobs.h.