VariabilityClass Class Reference

#include <recons_misc.h>

Collaboration diagram for VariabilityClass:

Public Types
enum	t_VAR_status { VAR_none, VAR_measuring, VAR_analyzing }

Public Member Functions
	VariabilityClass (BasicARTParameters *_prm, int _Zoutput_volume_size, int _Youtput_volume_size, int _Xoutput_volume_size)
	Constructor. More...
void	newIteration ()
void	newUpdateVolume (GridVolume *ptr_vol_out, Projection &read_proj)
void	finishAnalysis ()

Public Attributes
t_VAR_status	VAR_state
int	Zoutput_volume_size
int	Youtput_volume_size
int	Xoutput_volume_size
BasicARTParameters *	prm
std::vector< MultidimArray< double > >	VA
	Vector of training vectors. More...
int	N
	Number of updates so far. More...

Detailed Description

Variability structure

Definition at line 77 of file recons_misc.h.

Member Enumeration Documentation

t_VAR_status

enum VariabilityClass::t_VAR_status

Enumerator
VAR_none
VAR_measuring
VAR_analyzing

Definition at line 80 of file recons_misc.h.

{VAR_none, VAR_measuring, VAR_analyzing} t_VAR_status;

Constructor & Destructor Documentation

VariabilityClass()

VariabilityClass::VariabilityClass	(	BasicARTParameters *	_prm,
		int	_Zoutput_volume_size,
		int	_Youtput_volume_size,
		int	_Xoutput_volume_size
	)

Constructor.

Definition at line 396 of file recons_misc.cpp.

{
    prm = _prm;
    Zoutput_volume_size = _Zoutput_volume_size;
    Youtput_volume_size = _Youtput_volume_size;
    Xoutput_volume_size = _Xoutput_volume_size;
    N = 0;
    VAR_state = VAR_none;
}

Member Function Documentation

finishAnalysis()

void VariabilityClass::finishAnalysis

(

)

Finish analysis.

Definition at line 463 of file recons_misc.cpp.

{
    if (VA.size() == 0)
        return;

    // Coocurrence matrix
    int nmax = VA.size();
    MultidimArray<int> coocurrence(nmax, nmax);

    // Study each voxel
    Image<double> SignificantT2, SignificantMaxRatio, SignificantMinRatio;
    int zsize = ZSIZE(VA[0]) / 2;
    int ysize = YSIZE(VA[0]) / 2;
    int xsize = XSIZE(VA[0]) / 2;
    SignificantT2().initZeros(zsize, ysize, xsize);
    SignificantMaxRatio().initZeros(zsize, ysize, xsize);
    SignificantMinRatio().initZeros(zsize, ysize, xsize);
    std::cout << "Classifying voxels ...\n";
    init_progress_bar(MULTIDIM_SIZE(SignificantT2()));
    int counter = 0, nxny = ysize * zsize;
#ifdef MODE7
#define NFEATURES 7
#endif
#ifdef MODE8
#define NFEATURES 8
#endif
#ifdef MODE15
#define NFEATURES 15
#endif

    FOR_ALL_ELEMENTS_IN_ARRAY3D(SignificantT2())
    {
        // Save the data for this voxel
        std::ofstream fh_dat;
        fh_dat.open("PPP.dat");
        if (!fh_dat)
            REPORT_ERROR(ERR_IO_NOTOPEN, "VariabilityClass::finishAnalysis: "
                         "Cannot open PPP.dat for output");
        fh_dat << NFEATURES << " " << nmax << std::endl;
        std::vector< Matrix1D<double> > v;
        v.clear();
        for (int n = 0; n < nmax; n++)
        {
            Matrix1D<double> v_aux(NFEATURES);

#ifdef MODE7

            v_aux(0) = VA[n](k,      i, j + xsize);
            v_aux(1) = VA[n](k, i + ysize,      j);
            v_aux(2) = VA[n](k, i + ysize, j + xsize);
            v_aux(3) = VA[n](k + zsize,      i,      j);
            v_aux(4) = VA[n](k + zsize,      i, j + xsize);
            v_aux(5) = VA[n](k + zsize, i + ysize,      j);
            v_aux(6) = VA[n](k + zsize, i + ysize, j + xsize);
#endif
#ifdef MODE8

            v_aux(0) = VA[n](k,      i,      j);
            v_aux(1) = VA[n](k,      i, j + xsize);
            v_aux(2) = VA[n](k, i + ysize,      j);
            v_aux(3) = VA[n](k, i + ysize, j + xsize);
            v_aux(4) = VA[n](k + zsize,      i,      j);
            v_aux(5) = VA[n](k + zsize,      i, j + xsize);
            v_aux(6) = VA[n](k + zsize, i + ysize,      j);
            v_aux(7) = VA[n](k + zsize, i + ysize, j + xsize);
#endif
#ifdef MODE15

            v_aux(0) = VA[n](k / 2, i / 2,    j / 2);
            v_aux(1) = VA[n](k / 2, i / 2, j / 2 + xsize2);
            v_aux(2) = VA[n](k / 2, i / 2 + ysize2,    j / 2);
            v_aux(3) = VA[n](k / 2, i / 2 + ysize2, j / 2 + xsize2);
            v_aux(4) = VA[n](k / 2 + zsize2,    i / 2,       j / 2);
            v_aux(5) = VA[n](k / 2 + zsize2,    i / 2, j / 2 + xsize2);
            v_aux(6) = VA[n](k / 2 + zsize2, i / 2 + ysize2,       j / 2);
            v_aux(7) = VA[n](k / 2 + zsize2, i / 2 + ysize2, j / 2 + xsize2);
            v_aux(8) = VA[n](k,      i, j + xsize);
            v_aux(9) = VA[n](k, i + ysize,      j);
            v_aux(10) = VA[n](k, i + ysize, j + xsize);
            v_aux(11) = VA[n](k + zsize,      i,      j);
            v_aux(12) = VA[n](k + zsize,      i, j + xsize);
            v_aux(13) = VA[n](k + zsize, i + ysize,      j);
            v_aux(14) = VA[n](k + zsize, i + ysize, j + xsize);
#endif

            v_aux = v_aux * v_aux;
            // COSS: Doesn't work: v_aux=v_aux.sort();

            fh_dat << v_aux.transpose() << std::endl;
            v.push_back(v_aux);
        }
        fh_dat.close();

        // Classify
        if (system("xmipp_fcmeans -din PPP.dat -std::cout PPP -c 2 -saveclusters > /dev/null")==-1)
            REPORT_ERROR(ERR_UNCLASSIFIED,"Cannot open shell");

        // Pick up results
        Matrix2D<double> aux;
        int n_previous;
#define GET_RESULTS(fh,fn,avg,cov,N,idx) \
    fh.open(fn);\
    if (!fh) \
        REPORT_ERROR(ERR_IO_NOTOPEN,(std::string)"VariabilityClass::finishAnalysis: " \
                     "Cannot open "+fn+" for input"); \
    n_previous=-1; \
    while (!fh.eof()) { \
        int n; fh >> n; \
        if (n!=n_previous) { \
            n_previous=n; N++; \
            idx(n)=1; \
            avg+=v[n]; \
            aux.fromVector(v[n]); \
            cov+=aux*aux.transpose(); \
        }\
    } \
    avg/=N; \
    cov/=N; \
    aux.fromVector(avg); \
    cov-=aux*aux.transpose();

        std::ifstream fh_0;
        Matrix1D<double> avg0(NFEATURES);
        Matrix1D<int>    idx0(nmax);
        Matrix2D<double> covariance0(NFEATURES, NFEATURES);
        int N0 = 0;
        GET_RESULTS(fh_0, "PPP.0", avg0, covariance0, N0, idx0);
#ifdef DEBUG

        std::cout << "Class 0 is:\n";
        for (size_t n = 0; n < idx0.size(); n++)
        {
            if (idx0(n))
            {
                int iact_proj = prm->ordered_list(n);
                std::cout << prm->IMG_Inf[iact_proj].fn_proj << std::endl;
            }
        }
#endif

        std::ifstream fh_1;
        Matrix1D<double> avg1(NFEATURES);
        Matrix1D<int>    idx1(nmax);
        Matrix2D<double> covariance1(NFEATURES, NFEATURES);
        int N1 = 0;
        GET_RESULTS(fh_1, "PPP.1", avg1, covariance1, N1, idx1);
#ifdef DEBUG

        std::cout << "Class 1 is:\n";
        for (size_t n = 0; n < idx1.size(); n++)
        {
            if (idx1(n))
            {
                int iact_proj = prm->ordered_list(n);
                std::cout << prm->IMG_Inf[iact_proj].fn_proj << std::endl;
            }
        }
#endif

        Matrix2D<double> T2, covariance;
        if (NFEATURES > 1)
        {
            // Perform T2-Hotelling test
            Matrix1D<double> avg_diff = avg1 - avg0;
            covariance = 1.0 / (N0 + N1 - 2) *
                         ((N0 - 1.0) * covariance0 + (N1 - 1.0) * covariance1);
            covariance *= (1.0 / N0 + 1.0 / N1);
            aux.fromVector(avg_diff);
            T2 = (double)(N0 + N1 - avg_diff.size() - 1) /
                 ((N0 + N1 - 2) * avg_diff.size()) *
                 aux.transpose() * covariance.inv() * aux;
        }
        else
        {
            // Perform t-test
            double variance = ((N0 - 1) * covariance0(0, 0) + (N1 - 1) * covariance1(0, 0)) /
                              (N0 + N1 - 2);
            double t = (avg1(0) - avg0(0)) / sqrt(variance * (1.0 / N0 + 1.0 / N1));
            T2.initZeros(1, 1);
            T2(0, 0) = t;
        }

        // Analysis of the covariance structure
        Matrix1D<double> eigenvalues;
        if (NFEATURES > 1)
        {
            Matrix2D<double> U, V;
            svdcmp(covariance, U, eigenvalues, V);
        }

        // Analysis of the coocurrences
        for (size_t n = 0; n < idx0.size(); n++)
            for (size_t np = n + 1; np < idx0.size(); np++)
                if (idx0(n) && idx0(np))
                    coocurrence(n, np)++;

        for (size_t n = 0; n < idx1.size(); n++)
            for (size_t np = n + 1; np < idx1.size(); np++)
                if (idx1(n) && idx1(np))
                    coocurrence(n, np)++;

        // Keep results
        SignificantT2(k, i, j) = T2(0, 0);
        if (NFEATURES > 1)
        {
            SignificantMinRatio(k, i, j) = eigenvalues(1) / eigenvalues(0);
            SignificantMaxRatio(k, i, j) = eigenvalues(NFEATURES - 1) / eigenvalues(0);
        }
#ifdef DEBUG
        std::cout << "T2 for this classification is " << T2(0, 0) << std::endl;
        std::cout << "Eigenvalues are " << eigenvalues.transpose() << std::endl;
#endif

        if (++counter % nxny == 0)
            progress_bar(counter);
    }
    progress_bar(MULTIDIM_SIZE(SignificantT2()));
    SignificantT2.write(prm->fn_root + "_variability.vol");
    SignificantMinRatio.write(prm->fn_root + "_minratio.vol");
    SignificantMaxRatio.write(prm->fn_root + "_maxratio.vol");
    if (system("rm PPP.dat PPP.cod PPP.vs PPP.err PPP.his PPP.inf PPP.0 PPP.1")==-1)
        REPORT_ERROR(ERR_UNCLASSIFIED,"Cannot open shell");

    for (int n = 0; n < nmax; n++)
        for (int np = n + 1; np < nmax; np++)
            coocurrence(np, n) = coocurrence(n, np);
    Image<double> save;
    typeCast(coocurrence, save());
    save.write(prm->fn_root + "_coocurrence.xmp");
}

newIteration()

void VariabilityClass::newIteration

(

)

Start a new ART iteration.

Definition at line 408 of file recons_misc.cpp.

{}

newUpdateVolume()

void VariabilityClass::newUpdateVolume	(	GridVolume *	ptr_vol_out,
		Projection &	read_proj
	)

Update data with a new volume. The update volume is set to zeros after this function

Definition at line 415 of file recons_misc.cpp.

{
    Image<double> vol_voxels;

    // Convert from basis to voxels
    prm->basis.changeToVoxels(*ptr_vol_out, &(vol_voxels()),
                              Zoutput_volume_size, Youtput_volume_size, Xoutput_volume_size);
    (*ptr_vol_out).initZeros();
    N++;

    // Make the DWT
    MultidimArray<double> DWTV;
#ifdef MODE7

    int DWT_iterations = 2;
    int keep_from_iteration = 0;
#endif
#ifdef MODE8

    int DWT_iterations = 2;
    int keep_from_iteration = 0;
#endif
#ifdef MODE15

    int DWT_iterations = 2;
    int keep_from_iteration = 0;
#endif

    Bilib_DWT(vol_voxels(), DWTV, DWT_iterations);
#ifdef DEBUG

    vol_voxels.write("PPPVariability.vol");
    char c;
    std::cout << "Press any key\n";
    std::cin >> c;
#endif

    // Select the LLL block and keep it
    int x1, x2, y1, y2, z1, z2;
    SelectDWTBlock(keep_from_iteration, DWTV, "000", x1, x2, y1, y2, z1, z2);
    DWTV.selfWindow(z1, y1, x1, z2, y2, x2);
    STARTINGZ(DWTV) = STARTINGY(DWTV) = STARTINGX(DWTV) = 0;
    VA.push_back(DWTV);
}

Member Data Documentation

N

int VariabilityClass::N

Number of updates so far.

Definition at line 91 of file recons_misc.h.

prm

BasicARTParameters* VariabilityClass::prm

Definition at line 85 of file recons_misc.h.

VA

std::vector< MultidimArray<double> > VariabilityClass::VA

Vector of training vectors.

Definition at line 88 of file recons_misc.h.

VAR_state

t_VAR_status VariabilityClass::VAR_state

Definition at line 81 of file recons_misc.h.

Xoutput_volume_size

int VariabilityClass::Xoutput_volume_size

Definition at line 84 of file recons_misc.h.

Youtput_volume_size

int VariabilityClass::Youtput_volume_size

Definition at line 83 of file recons_misc.h.

Zoutput_volume_size

int VariabilityClass::Zoutput_volume_size

Definition at line 82 of file recons_misc.h.

---

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

• xmipp/libraries/reconstruction/recons_misc.h
• xmipp/libraries/reconstruction/recons_misc.cpp