EMAN::FourierReconstructor Class Reference

Fourier space 3D reconstruction The Fourier reconstructor is designed to work in an iterative fashion, where similarity ("quality") metrics are used to determine if a slice should be inserted into the 3D in each subsequent iteration. More...

#include <reconstructor.h>

Inheritance diagram for EMAN::FourierReconstructor:

[legend]

Collaboration diagram for EMAN::FourierReconstructor:

[legend]

List of all members.

Public Member Functions

FourierReconstructor ()

Default constructor calls load_default_settings().

virtual ~FourierReconstructor ()

Deconstructor calls free_memory().

virtual void setup ()

Setup the Fourier reconstructor relies on the correct parameters throws a variety of exceptions if the parameters are unworkable
Exceptions:

InvalidValueException when the x_in parameter is odd or less than zero

InvalidValueException when the y_in parameter is odd or less than zero

InvalidValueException when the x_out parameter is odd or less than zero

InvalidValueException when the y_out parameter is odd or less than zero

InvalidValueException when the z_out parameter is odd or less than zero

InvalidValueException when the pad parameter is less than the greatest dimension of the input images (the greater of x_in and y_in) Note the restraint that the Fourier volumes should be even will be lifted once debugging of the the xform.fourierorigin processor is performed, but however, when this happens a rigorous testing should be performed because no testing has been done in this regard.

virtual int insert_slice (const EMData *const slice, const Transform &t3d=Transform())

Insert a slice into a 3D volume, in a given orientation.

virtual int determine_slice_agreement (const EMData *const input_slice, const Transform &t3d, const unsigned int num_particles_in_slice=1)

Determine slice agreement with the current reconstruction.

virtual int insert_slice (const EMData *const slice, const Transform3D &t3d)

Insert an image slice to the reconstructor.

virtual int determine_slice_agreement (const EMData *const input_slice, const Transform3D &t3d, const unsigned int num_particles_in_slice=1)

virtual EMData * finish ()

Get the reconstructed volume Peforms Fourier inversion on a potentially large volume and may take several minutes.

virtual string get_name () const

Get the unique name of the reconstructor.

virtual string get_desc () const

Get the one line description of the reconstructor.

virtual TypeDict get_param_types () const

Get the parameter types of this object.

virtual float get_score (const unsigned int idx)

Get the quality score that has been determined for this slice this is used in e2make3d.py to print information to std out.

virtual float get_norm (const unsigned int idx)

Get the normalization value that has been determined for this slice this is used in e2make3d.py to print information to std out.

virtual void zero_memory ()

Sends the pixels in tmp_data and image to zero Convenience only.

Static Public Member Functions

static Reconstructor * NEW ()

Factory incorporation uses the pointer of this function.

Protected Member Functions

EMData * preprocess_slice (const EMData *const slice, const Transform &t=Transform())

Preprocess the slice prior to insertion into the 3D volume this Fourier tranforms the slice and make sure all the pixels are in the right position.

void load_default_settings ()

Load default settings.

void free_memory ()

Frees the memory owned by this object (but not parent objects) Deletes the FourierPixelInserter3D pointer.

void load_inserter ()

Load the pixel inserter based on the information in params.

void load_interp_FRC_calculator ()

Load the pixel inserter based on the information in params.

void do_insert_slice_work (const EMData *const input_slice, const Transform &euler)

A function to perform the nuts and bolts of inserting an image slice.

void print_stats (const vector< InterpolatedFRC::QualityScores > &scores)

print stats is called internally at various points in the reconstruction routine and is for the benefit of people using e2make3d.py

Protected Attributes

vector
< InterpolatedFRC::QualityScores > quality_scores

Quality scores vectors for storing normalization constants, and SNR normalized Fourier ring correlation scores.

vector
< InterpolatedFRC::QualityScores > prev_quality_scores

unsigned int image_idx

Index used internally to index into the quality scores vectors at different points ( in insert_slice and determine_slice_agreement).

FourierPixelInserter3D * inserter

A pixel inserter pointer which inserts pixels into the 3D volume using one of a variety of insertion methods.

InterpolatedFRC * interp_FRC_calculator

A pixel inserter pointer which inserts pixels into the 3D volume using one of a variety of insertion methods.

bool slice_insertion_flag

Internal flags used to perform memory zeroing and normalization transparently.

bool slice_agreement_flag

float x_scale_factor

Used for scaling frequency axes when any of the x_out, y_out or z_out parameters are specified.

float y_scale_factor

float z_scale_factor

int max_input_dim

Keeps track of the maximum dimension size.

int output_x

Keeps track of the eventual size of the output real space volume.

int output_y

int output_z

Private Member Functions

FourierReconstructor (const FourierReconstructor &that)

Disallow copy construction.

FourierReconstructor & operator= (const FourierReconstructor &)

Disallow assignment.

Detailed Description

Fourier space 3D reconstruction The Fourier reconstructor is designed to work in an iterative fashion, where similarity ("quality") metrics are used to determine if a slice should be inserted into the 3D in each subsequent iteration.

The client creates a Fourier reconstructor to insert real images into a 3D volume. The return image is a real space image

This reconstructor is based on EMAN1's Fourier reconstructor with a handful of modifications including 1. - Fourier ring correlation (FRC) as opposed to the mean phase residual is used to estimate slice quality. The FRC of the slice in the 3D volume is determined - but the slice is removed from the 3D volume before doing this so the score reflects the extent to which the slice agrees with the contribution of the other slices in the 3D volume. The FRC is converted to SNR using the relationship described by Penczek ( Three-dimensional spectral signal to noise ratio for a class of reconstruction algorithms, JSB 2002 138 (24-46) ) FRC = S/(sqrt(S+N1)sqrt(S+N2)) Where N1 is the noise in the slice of the 3D volume and N2 is the noise in the image slice being inserted. We make the assumption that the noise in the 3D volume is 0 (N1=0) to get FRC^2 = SNR/(1+SNR) which gives a spectral SNR plot - we then divide each SNR value by the number of particles in the class average (seeing as SNR should scale linearly with the number of particles) to get the estimated SNR per contributing particle in this class average. If the particles that have been averaged are not homogenous this score should be low etc. The scaled SNR curve is then converted back to a FRC curve and integrated. This integral is the similarity metric, and depends on how far information extends to in Fourier space - typical values range from 0.05 to 0.2, but can vary substantially depending on the data.

2 - Uses half of the memory used by EMAN1's equivalent reconstruction algorithm

Fourier reconstructor usage

        Reconstructor* r = Factory<Reconstructor>::get("fourier", params);
        r->setup();
        for k in 0:num_iterations-1
                // First do a round of slice quality (metric) determination - only possible if a 3D volume has
                // already been generated (k>0)
                if ( k  > 0 )
                        // Determine the agreement of the slices with the previous reconstructed volume (in memory)
                        for i in 0:num_slices-1
                                r->determine_slice_agreement(image[i], image[i].euler_orientation);

                // Insert the slices into the 3D volume
                // Will decide not to insert the slice if the its "quality" is not good enough
                for i in 0:num_slices-1
                        int failure = r->insert_slice(image[i], image[i].euler_orientation);
                        if ( failure ) cout << "Slice was not inserted due to poor quality" << endl;

        // Get the resulting volume
        EMData* result = r->finish();
        result->write_image("threed.mrc");

Definition at line 348 of file reconstructor.h.

Constructor & Destructor Documentation

EMAN::FourierReconstructor::FourierReconstructor ( ) [inline]

Default constructor calls load_default_settings().

Definition at line 354 of file reconstructor.h.

References load_default_settings().

Referenced by NEW().

00354                                        : image_idx(0), inserter(0), interp_FRC_calculator(0), slice_insertion_flag(true),
00355                 slice_agreement_flag(false), x_scale_factor(0.0), y_scale_factor(0.0), z_scale_factor(0.0),
00356                 max_input_dim(0), output_x(0), output_y(0), output_z(0) { load_default_settings(); }

virtual EMAN::FourierReconstructor::~FourierReconstructor ( ) [inline, virtual]

Deconstructor calls free_memory().

Definition at line 361 of file reconstructor.h.

References free_memory().

00361 { free_memory(); }

EMAN::FourierReconstructor::FourierReconstructor ( const FourierReconstructor & that ) [private]

Disallow copy construction.

Member Function Documentation

void FourierReconstructor::setup ( ) [virtual]

Setup the Fourier reconstructor relies on the correct parameters throws a variety of exceptions if the parameters are unworkable

Exceptions:

	InvalidValueException	when the x_in parameter is odd or less than zero
	InvalidValueException	when the y_in parameter is odd or less than zero
	InvalidValueException	when the x_out parameter is odd or less than zero
	InvalidValueException	when the y_out parameter is odd or less than zero
	InvalidValueException	when the z_out parameter is odd or less than zero
	InvalidValueException	when the pad parameter is less than the greatest dimension of the input images (the greater of x_in and y_in) Note the restraint that the Fourier volumes should be even will be lifted once debugging of the the xform.fourierorigin processor is performed, but however, when this happens a rigorous testing should be performed because no testing has been done in this regard.

Implements EMAN::Reconstructor.

Reimplemented in EMAN::BaldwinWoolfordReconstructor.

Definition at line 328 of file reconstructor.cpp.

References EMAN::ReconstructorVolumeData::image, InvalidValueException, is_fftodd(), load_inserter(), load_interp_FRC_calculator(), max_input_dim, EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, output_x, output_y, output_z, EMAN::FactoryBase::params, EMAN::EMData::set_complex(), EMAN::Dict::set_default(), EMAN::EMData::set_fftodd(), EMAN::EMData::set_ri(), EMAN::EMData::set_size(), EMAN::ReconstructorVolumeData::tmp_data, EMAN::EMData::to_zero(), EMAN::EMData::update(), x_scale_factor, y_scale_factor, and z_scale_factor.

Referenced by EMAN::BaldwinWoolfordReconstructor::setup().

00329 {
00330         // default setting behavior - does not override if the parameter is already set
00331         params.set_default("mode",2);
00332 
00333         int x_size = params["x_in"];
00334         if ( x_size < 0 ) throw InvalidValueException(x_size, "x size of images must be greater than 0");
00335         int y_size = params["y_in"];
00336         if ( y_size < 0 ) throw InvalidValueException(y_size, "y size of images must be greater than 0");
00337 
00338         if ( x_size > y_size ) max_input_dim = x_size;
00339         else max_input_dim = y_size;
00340 
00341         // This is a helical adaptation - FIXME explain
00342         bool helical_special_behavior = false;
00343         if ( helical_special_behavior )
00344         {
00345                 if ( x_size > y_size )
00346                 {
00347                         if ( (int) params["xsample"] == 0 ) params["xsample"] = y_size;
00348                         if ( (int) params["zsample"] == 0 ) params["zsample"] = x_size;
00349                 }
00350                 if ( y_size > x_size )
00351                 {
00352                         if ( (int) params["ysample"] == 0 ) params["ysample"] = x_size;
00353                         if ( (int) params["zsample"] == 0 ) params["zsample"] = y_size;
00354                 }
00355         }
00356 
00357         if ( (int) params["xsample"] != 0  ) output_x = (int) params["xsample"];
00358         else output_x = x_size;
00359 
00360         if ( (int) params["ysample"] != 0  ) output_y = (int) params["ysample"];
00361         else output_y = y_size;
00362 
00363         if ( (int) params["zsample"] != 0  ) output_z = (int) params["zsample"];
00364         else
00365         {
00366                 if ( x_size == y_size ) output_z = x_size;
00367                 else if ( x_size > y_size ) output_z = x_size;
00368                 else output_z = y_size;
00369         }
00370 
00371         nx = output_x;
00372         ny = output_y;
00373         nz = output_z;
00374 
00375         // Adjust nx if for Fourier transform even odd issues
00376         bool is_fftodd = (nx % 2 == 1);
00377         // The Fourier transform requires one extra pixel in the x direction,
00378         // which is two spaces in memory, one each for the complex and the
00379         // real components
00380         nx += 2-is_fftodd;
00381 
00382         if ( (int) params["zsample"] == 0  ) nz = max_input_dim;
00383 
00384         if ( nz < max_input_dim ) z_scale_factor = (float) nz / (float) max_input_dim;
00385         if ( ny < max_input_dim ) y_scale_factor = (float) ny / (float) max_input_dim;
00386         if ( nx < max_input_dim ) x_scale_factor = (float) nx / (float) max_input_dim;
00387 
00388         // Odd dimension support is here atm, but not above.
00389         image = new EMData();
00390         image->set_size(nx, ny, nz);
00391         image->set_complex(true);
00392         image->set_fftodd(is_fftodd);
00393         image->set_ri(true);
00394 
00395         tmp_data = new EMData();
00396         tmp_data->set_size(nx/2, ny, nz);
00397         tmp_data->to_zero();
00398         tmp_data->update();
00399 
00400         load_inserter();
00401         load_interp_FRC_calculator();
00402 
00403         if ( (bool) params["quiet"] == false )
00404         {
00405                 cout << "3D Fourier dimensions are " << nx << " " << ny << " " << nz << endl;
00406                 cout << "You will require approximately " << setprecision(3) << (nx*ny*nz*sizeof(float)*1.5)/1000000000.0 << "GB of memory to reconstruct this volume" << endl;
00407                 cout << "Scale factors are " << x_scale_factor << " " << y_scale_factor << " " << z_scale_factor << endl;
00408                 cout << "Max input dim is " << max_input_dim << endl;
00409         }
00410 }

int FourierReconstructor::insert_slice	(	const EMData *const	slice,
		const Transform &	t3d = `Transform()`
	)			`[virtual]`

Insert a slice into a 3D volume, in a given orientation.

Returns:: 0 if successful, 1 otherwise

Parameters:

	slice	the image slice to be inserted into the 3D volume
	t3d	the Transform that stores the image Euler angles

Exceptions:

	NullPointerException	if the input EMData pointer is null
	ImageFormatException	if the image is complex as opposed to real

Reimplemented from EMAN::Reconstructor.

Reimplemented in EMAN::BaldwinWoolfordReconstructor.

Definition at line 481 of file reconstructor.cpp.

References do_insert_slice_work(), EMAN::EMData::get_attr(), get_norm(), EMAN::EMData::has_attr(), image_idx, ImageFormatException, EMAN::EMData::is_complex(), EMAN::EMData::mult(), NullPointerException, EMAN::FactoryBase::params, preprocess_slice(), quality_scores, EMAN::Transform::set_mirror(), EMAN::Transform::set_scale(), EMAN::Transform::set_trans(), slice_agreement_flag, slice_insertion_flag, and zero_memory().

00482 {
00483         // Are these exceptions really necessary? (d.woolford)
00484         if (!input_slice) throw NullPointerException("EMData pointer (input image) is NULL");
00485 
00486         // I could also test to make sure the image is the right dimensions...
00487         if (input_slice->is_complex()) throw ImageFormatException("The image is complex, expecting real");
00488 
00489         // Get the proprecessed slice - there are some things that always happen to a slice,
00490         // such as as Fourier conversion and optional padding etc.
00491         //
00492         // We must use only the rotational component of the transform, scaling, translation and mirroring
00493         // are not implemented in Fourier space
00494         Transform * rotation;
00495         if ( input_slice->has_attr("xform.projection") ) {
00496                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
00497         } else {
00498                 rotation = new Transform(arg); // assignment operator
00499         }
00500 
00501         EMData* slice = preprocess_slice( input_slice, *rotation);
00502 
00503         // Catch the first time a slice is inserted to do some things....
00504         if ( slice_insertion_flag == true )
00505         {
00506                 // Zero the memory in the associated EMData objects
00507                 zero_memory();
00508                 // Reset the image_idx to zero
00509                 image_idx = 0;
00510                 // Set flags appropriately
00511                 slice_agreement_flag = true;
00512                 slice_insertion_flag = false;
00513 
00514                 // should be a if verbose here
00515                 //if ( prev_quality_scores.size() != 0 ) print_stats( prev_quality_scores );
00516         }
00517 
00518         // quality_scores.size() is zero on the first run, so this enforcement of slice quality does not take
00519         // place until after the first round of slice insertion
00520         if ( quality_scores.size() != 0 )
00521         {
00522                 if ( quality_scores[image_idx].get_snr_normed_frc_integral() < (float) params["hard"] )
00523                 {
00524                         image_idx++;
00525                         return 1;
00526                 }
00527                 slice->mult(1.f/quality_scores[image_idx].get_norm());
00528 //              cout << "Norm multiplied by " << 1.f/quality_scores[image_idx].get_norm() << endl;
00529                 image_idx++;
00530         }
00531 
00532         // Finally to the pixel wise slice insertion
00533         rotation->set_scale(1.0);
00534         rotation->set_mirror(false);
00535         rotation->set_trans(0,0,0);
00536         do_insert_slice_work(slice, *rotation);
00537 
00538         if(rotation) {delete rotation; rotation=0;}
00539         delete slice;
00540 
00541 //      image->update();
00542         return 0;
00543 }

int FourierReconstructor::determine_slice_agreement	(	const EMData *const	input_slice,
		const Transform &	t3d,
		const unsigned int	num_particles_in_slice = `1`
	)			`[virtual]`

Determine slice agreement with the current reconstruction.

Returns:: 0 if successful, 1 otherwise

Parameters:

	input_slice	the image slice to compared against the 3D volume
	t3d	the Transform that the image Euler angles
	num_particles_in_slice	the number of particles in the slice - used to determine the SNR normalized FSC, defaults to 1.

Exceptions:

	NullPointerException	if the input EMData pointer is null
	ImageFormatException	if the image is complex as opposed to real

Reimplemented from EMAN::Reconstructor.

Definition at line 622 of file reconstructor.cpp.

References EMAN::InterpolatedFRC::continue_frc_calc(), dt, EMAN::InterpolatedFRC::finish(), EMAN::EMData::get_attr(), EMAN::EMData::get_data(), EMAN::InterpolatedFRC::QualityScores::get_norm(), get_norm(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::has_attr(), image_idx, ImageFormatException, interp_FRC_calculator, EMAN::EMData::is_complex(), EMAN::EMData::is_fftodd(), LOGERR, max_input_dim, EMAN::EMData::mult(), EMAN::ReconstructorVolumeData::normalize_threed(), NullPointerException, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, preprocess_slice(), prev_quality_scores, quality_scores, EMAN::InterpolatedFRC::reset(), EMAN::Dict::set_default(), EMAN::Transform::set_mirror(), EMAN::InterpolatedFRC::QualityScores::set_norm(), EMAN::Transform::set_scale(), EMAN::Transform::set_trans(), slice_agreement_flag, slice_insertion_flag, EMAN::Util::square(), square, weight, x, x_scale_factor, y, y_scale_factor, and z_scale_factor.

00623 {
00624         // Are these exceptions really necessary? (d.woolford)
00625         if (!input_slice) {
00626                 LOGERR("Insertion of NULL slice in FourierReconstructor::insert_slice");
00627                 throw NullPointerException("EMData pointer (input image) is NULL");
00628         }
00629 
00630         // I could also test to make sure the image is the right dimensions...
00631         if (input_slice->is_complex()) {
00632                 LOGERR("Do not Fourier transform the image before it is passed to determine_slice_agreement in FourierReconstructor, this is performed internally");
00633                 throw ImageFormatException("The image is complex, expecting real");
00634         }
00635 
00636         // The first time determine_slice_agreement is called (in each iteration) some things need to happen...
00637         if ( slice_agreement_flag == true )
00638         {
00639                 // Normalize the real data using the normalization values in the normalization volume
00640                 // This is important because the InterpolatedFRC objects assumes this behaviour
00641                 normalize_threed();
00642                 // Reset the image_idx to index into prev_quality_scores correctly
00643                 image_idx = 0;
00644                 // Make a copy of the previous quality scores - the first time around this will be like copying nothing
00645                 prev_quality_scores = quality_scores;
00646                 // Now clear the old scores so they can be redetermined
00647                 quality_scores.clear();
00648                 // Reset the flags
00649                 slice_insertion_flag = true;
00650                 slice_agreement_flag = false;
00651         }
00652 
00653         // Get the proprecessed slice - there are some things that always happen to a slice,
00654         // such as as Fourier conversion and optional padding etc.
00655         Transform * rotation;
00656         if ( input_slice->has_attr("xform.projection") ) {
00657                 rotation = (Transform*) (input_slice->get_attr("xform.projection")); // assignment operator
00658 
00659         } else {
00660                 rotation = new Transform(t3d); // assignment operator
00661         }
00662         EMData* slice = preprocess_slice( input_slice, *rotation );
00663 
00664         // quality_scores.size() is zero on the first run, so this enforcement of slice quality does not take
00665         // place until after the first round of slice insertion
00666         if (  prev_quality_scores.size() != 0 )
00667         {
00668                 // The slice must be multiplied by the normalization value used in insert or else the calculations will
00669                 // be inconsistent
00670                 slice->mult(1.f/prev_quality_scores[image_idx].get_norm());
00671         }
00672 
00673         // Reset zeros the associated memory in the ifrc
00674         interp_FRC_calculator->reset();
00675 
00676         int rl = Util::square( max_input_dim / 2);
00677 
00678         float dt[2];
00679 
00680         float y_scale = 0, x_scale = 0;
00681 
00682         int y_in = slice->get_ysize();
00683         int x_in = slice->get_xsize();
00684         if (input_slice->is_fftodd()) x_in -= 1;
00685         else x_in -= 2;
00686 
00687         if ( y_in != x_in )
00688         {
00689                 if ( x_in > y_in ) y_scale = (float) x_in / (float) y_in;
00690                 else x_scale = (float) y_in / (float) x_in;
00691         }
00692         float *dat = slice->get_data();
00693 
00694         float weight = params.set_default("weight",1.0f);
00695 
00696 
00697         rotation->set_scale(1.0);
00698         rotation->set_mirror(false);
00699         rotation->set_trans(0,0,0);
00700 
00701         for (int y = 0; y < slice->get_ysize(); y++) {
00702                 for (int x = 0; x < slice->get_xsize() / 2; x++) {
00703 
00704                         float rx = (float) x;
00705                         float ry = (float) y;
00706 
00707                         if ( y_in != x_in )
00708                         {
00709                                 if ( x_in > y_in ) ry *= y_scale;
00710                                 else rx *= x_scale;
00711                         }
00712 
00713                         if ((rx * rx + Util::square(ry - max_input_dim / 2)) > rl)
00714                                 continue;
00715 
00716 //                      float xx = (float) (rx * t3d[0][0] + (ry - max_input_dim / 2) * t3d[1][0]);
00717 //                      float yy = (float) (rx * t3d[0][1] + (ry - max_input_dim / 2) * t3d[1][1]);
00718 //                      float zz = (float) (rx * t3d[0][2] + (ry - max_input_dim / 2) * t3d[1][2]);
00719 
00720                         Vec3f coord(rx,(ry - max_input_dim / 2),0);
00721                         coord = coord*(*rotation); // transpose multiplication
00722                         float xx = coord[0];
00723                         float yy = coord[1];
00724                         float zz = coord[2];
00725 
00726                         float cc = 1;
00727 
00728                         if (xx < 0) {
00729                                 xx = -xx;
00730                                 yy = -yy;
00731                                 zz = -zz;
00732                                 cc = -1.0;
00733                         }
00734 
00735                         if ( z_scale_factor != 0 ) zz *= z_scale_factor;
00736                         if ( y_scale_factor != 0 ) yy *= y_scale_factor;
00737                         if ( x_scale_factor != 0 ) xx *= x_scale_factor;
00738 
00739                         yy += ny / 2;
00740                         zz += nz / 2;
00741 
00742                         int idx = x * 2 + y * (slice->get_xsize());
00743                         dt[0] = dat[idx];
00744                         dt[1] = cc * dat[idx+1];
00745 
00746 
00747                         if ( prev_quality_scores.size() != 0 )
00748                         {
00749                                 // If the slice was not inserted into the 3D volume in the previous round of slice insertion
00750                                 // then the weight must be set to zero, or else the result produced by the InterpolatedFRC will be wrong.
00751                                 // This is because the InterpolatedFRC subtracts the incoming pixels from the 3D volume in order
00752                                 // to estimate quality using only data from the other slices in the volume. If the slice was not previously inserted
00753                                 // then it should not be subtracted prior to quality estimation....
00754                                 if ( prev_quality_scores[image_idx].get_snr_normed_frc_integral() < (float) params["hard"] )
00755                                 {
00756                                         weight = 0;
00757                                 }
00758                         }
00759                         // FIXME: this could be replaced in favor of a class implementation with no switch statement, similar to the
00760                         // inserter in the Fourier reconstructor do_insert_slice_work method
00761 
00762                         interp_FRC_calculator->continue_frc_calc(xx, yy, zz, dt, weight);
00763                 }
00764         }
00765 
00766         InterpolatedFRC::QualityScores q_scores = interp_FRC_calculator->finish( num_particles_in_slice );
00767         // Print the quality scores here for debug information
00768         //q_scores.debug_print();
00769 
00770         if (prev_quality_scores.size() != 0 )
00771         {
00772                 // If a previous normalization value has been calculated then average it with the one that
00773                 // was just calculated, this will cause the normalization values to converge more rapidly.
00774                 q_scores.set_norm((q_scores.get_norm() + prev_quality_scores[image_idx].get_norm())/2.0f);
00775                 image_idx++;
00776         }
00777 
00778         quality_scores.push_back(q_scores);
00779 
00780         if(rotation) {delete rotation; rotation=0;}
00781         delete slice;
00782         return 0;
00783 }

virtual int EMAN::FourierReconstructor::insert_slice	(	const EMData *const	slice,
		const Transform3D &	euler
	)			`[inline, virtual]`

Insert an image slice to the reconstructor.

To insert multiple image slices, call this function multiple times.

Parameters:

	slice	Image slice.
	euler	Euler angle of this image slice.

Returns:: 0 if OK. 1 if error.

Reimplemented from EMAN::Reconstructor.

Reimplemented in EMAN::BaldwinWoolfordReconstructor.

Definition at line 398 of file reconstructor.h.

References UnexpectedBehaviorException.

00398                                                                                              {
00399                         throw UnexpectedBehaviorException("No support for Transform3D in insert_slice, use Transform instead");
00400                 };

virtual int EMAN::FourierReconstructor::determine_slice_agreement	(	const EMData *const	const,
		const Transform3D &	,
		const unsigned int	= `1`
	)			`[inline, virtual]`

Returns:

Parameters:

	input_slice
	arg
	num_particles_in_slice

Exceptions:

Reimplemented from EMAN::Reconstructor.

Definition at line 402 of file reconstructor.h.

References UnexpectedBehaviorException.

00402                                                                                                                                                                 {
00403                         throw UnexpectedBehaviorException("No support for Transform3D in determine_slice_agreement, use Transform instead");
00404                 };

EMData * FourierReconstructor::finish ( ) [virtual]

Get the reconstructed volume Peforms Fourier inversion on a potentially large volume and may take several minutes.

Returns:: The real space reconstructed volume

Implements EMAN::Reconstructor.

Reimplemented in EMAN::BaldwinWoolfordReconstructor.

Definition at line 849 of file reconstructor.cpp.

References EMAN::EMData::clip_inplace(), EMAN::EMData::depad(), EMAN::EMData::do_ift_inplace(), EMAN::EMData::get_data(), EMAN::EMData::get_fft_amplitude(), EMAN::ReconstructorVolumeData::image, in, is_fftodd(), norm(), EMAN::ReconstructorVolumeData::normalize_threed(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, output_x, output_y, output_z, EMAN::FactoryBase::params, EMAN::EMData::process_inplace(), rdata, sqrt(), EMAN::Util::square_sum(), EMAN::ReconstructorVolumeData::tmp_data, EMAN::EMData::update(), and EMAN::EMData::write_image().

00850 {
00851         float *norm = tmp_data->get_data();
00852         float *rdata = image->get_data();
00853 
00854         if (params["dlog"]) {
00855                 size_t size = nx*ny*nz;
00856                 for (size_t i = 0; i < size; i += 2) {
00857                         float d = norm[i];
00858                         if (d == 0) {
00859                                 rdata[i] = 0;
00860                                 rdata[i + 1] = 0;
00861                         }
00862                         else {
00863                                 float in = norm[i + 1] / norm[i];
00864                                 float cin = Util::square_sum(rdata[i], rdata[i + 1]);
00865                                 rdata[i] *= sqrt(in / cin);
00866                                 rdata[i + 1] *= sqrt(in / cin);
00867                         }
00868                 }
00869         }
00870         else {
00871                 normalize_threed();
00872         }
00873 
00874 //      tmp_data->write_image("density.mrc");
00875 
00876         // we may as well delete the tmp data now... it saves memory and the calling program might
00877         // need memory after it gets the return volume.
00878         // If this delete didn't happen now, it would happen when the deconstructor was called,
00879         if ( tmp_data != 0 )
00880         {
00881                 delete tmp_data;
00882                 tmp_data = 0;
00883         }
00884 
00885 
00886         if ( params["3damp"]) {
00887                 EMData* fftimage = image->get_fft_amplitude ();
00888                 fftimage->write_image("threed_fft_amp.mrc");
00889                 delete fftimage;
00890         }
00891 
00892         //
00893         image->process_inplace("xform.fourierorigin.tocorner");
00894         image->do_ift_inplace();
00895         image->depad();
00896         image->process_inplace("xform.phaseorigin.tocenter");
00897 
00898         // If the image was padded it should be the original size, as the client would expect
00899         //  I blocked the rest, it is almost certainly incorrect  PAP 07/31/08
00900         // No, it's not incorrect. You are wrong. You have the meaning of nx mixed up. DSAW 09/23/cd
00901         bool is_fftodd = (nx % 2 == 1);
00902         if ( (nx-2*(!is_fftodd)) != output_x || ny != output_y || nz != output_z )
00903         {
00904                 FloatPoint origin( (nx-output_x)/2, (ny-output_y)/2, (nz-output_z)/2 );
00905                 FloatSize region_size( output_x, output_y, output_z);
00906                 Region clip_region( origin, region_size );
00907                 image->clip_inplace( clip_region );
00908         }
00909 
00910         // Should be an "if (verbose)" here or something
00911         //print_stats(quality_scores);
00912 
00913         image->update();
00914 
00915         return image;
00916 }

virtual string EMAN::FourierReconstructor::get_name ( ) const [inline, virtual]

Get the unique name of the reconstructor.

Implements EMAN::FactoryBase.

Reimplemented in EMAN::BaldwinWoolfordReconstructor.

Definition at line 415 of file reconstructor.h.

00416                 {
00417                         return "fourier";
00418                 }

virtual string EMAN::FourierReconstructor::get_desc ( ) const [inline, virtual]

Get the one line description of the reconstructor.

Implements EMAN::FactoryBase.

Reimplemented in EMAN::BaldwinWoolfordReconstructor.

Definition at line 422 of file reconstructor.h.

00423                 {
00424                         return "Reconstruction via direct Fourier methods using a combination of nearest neighbour and Gaussian kernels";
00425                 }

static Reconstructor* EMAN::FourierReconstructor::NEW ( ) [inline, static]

Factory incorporation uses the pointer of this function.

Returns:: a Reconstructor pointer to a newly allocated FourierReconstructor

Reimplemented in EMAN::BaldwinWoolfordReconstructor.

Definition at line 430 of file reconstructor.h.

References FourierReconstructor().

00431                 {
00432                         return new FourierReconstructor();
00433                 }

virtual TypeDict EMAN::FourierReconstructor::get_param_types ( ) const [inline, virtual]

Get the parameter types of this object.

Returns:: a TypeDict detailing all of the acceptable (and necessary) parameters

Implements EMAN::FactoryBase.

Reimplemented in EMAN::BaldwinWoolfordReconstructor.

Definition at line 438 of file reconstructor.h.

References EMAN::EMObject::BOOL, EMAN::EMObject::EMDATA, EMAN::EMObject::FLOAT, EMAN::EMObject::INT, EMAN::TypeDict::put(), and EMAN::EMObject::STRING.

00439                 {
00440                         TypeDict d;
00441                         d.put("x_in", EMObject::INT, "Necessary. The x dimension of the input images.");
00442                         d.put("y_in", EMObject::INT, "Necessary. The y dimension of the input images.");
00443                         d.put("zsample", EMObject::INT, "Optional. The z dimension (Fourier sampling) of the reconstructed volume, very useful for tomographic reconstruction. Works for general volumes.");
00444                         d.put("ysample", EMObject::INT, "Optional. The y dimension (Fourier sampling) of the reconstructed volume, works for general volumes. Not commonly specified.");
00445                         d.put("xsample", EMObject::INT, "Optional. The x dimension (Fourier sampling) of the reconstructed volume, works for general volumes. Not commonly specified.");
00446                         d.put("mode", EMObject::INT, "Optional. Fourier pixel insertion mode [1-8] - mode 2 is default.");
00447                         d.put("hard", EMObject::FLOAT, "Optional. The quality metric threshold. Default is 0 (off).");
00448                         d.put("sym", EMObject::STRING, "Optional. The symmetry of the reconstructed volume, c?, d?, oct, tet, icos, h?. Default is c1");
00449                         d.put("quiet", EMObject::BOOL, "Optional. Toggles writing useful information to standard out. Default is false.");
00450                         d.put("3damp", EMObject::BOOL, "Optional. Toggles writing the 3D FFT amplitude image. Default is false.");
00451                         d.put("weight", EMObject::FLOAT, "Weight of the slice that is being inserted. Default is 1.0.");
00452                         d.put("start_model", EMObject::EMDATA, "Start model");
00453                         d.put("start_model_weight", EMObject::FLOAT, "start_model_weight");
00454                         return d;
00455                 }

virtual float EMAN::FourierReconstructor::get_score ( const unsigned int idx ) [inline, virtual]

Get the quality score that has been determined for this slice this is used in e2make3d.py to print information to std out.

Returns:: the snr normalized Fourier ring correlation score

Parameters:

idx

the index of the slice in the quality_scores vector

Exceptions:

GenericException(just

throw) when the idx is beyond the range of the quality_scores vector

Reimplemented from EMAN::Reconstructor.

Definition at line 463 of file reconstructor.h.

References quality_scores, and UnexpectedBehaviorException.

00463 { if ( quality_scores.size() > idx ) return quality_scores[idx].get_snr_normed_frc_integral(); else throw UnexpectedBehaviorException("The requested index was beyond the length of the quality scores vector."); }

virtual float EMAN::FourierReconstructor::get_norm ( const unsigned int idx ) [inline, virtual]

Get the normalization value that has been determined for this slice this is used in e2make3d.py to print information to std out.

Returns:: the normalization value

Parameters:

idx

the index of the slice in the quality_scores vector

Exceptions:

GenericException(throw)

when the idx is beyond the range of the quality_scores vector

Reimplemented from EMAN::Reconstructor.

Definition at line 471 of file reconstructor.h.

References quality_scores, and UnexpectedBehaviorException.

Referenced by determine_slice_agreement(), and insert_slice().

00471 { if ( quality_scores.size() > idx ) return quality_scores[idx].get_norm();  else throw UnexpectedBehaviorException("The requested index was beyond the length of the quality scores vector."); }

void FourierReconstructor::zero_memory ( ) [virtual]

Sends the pixels in tmp_data and image to zero Convenience only.

Reimplemented from EMAN::ReconstructorVolumeData.

Definition at line 448 of file reconstructor.cpp.

References EMAN::EMData::add(), EMAN::EMData::do_fft_inplace(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), EMAN::ReconstructorVolumeData::image, ImageDimensionException, EMAN::EMData::is_complex(), EMAN::EMData::is_shuffled(), EMAN::EMData::mult(), EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, EMAN::EMData::process_inplace(), EMAN::Dict::set_default(), EMAN::ReconstructorVolumeData::tmp_data, and EMAN::EMData::to_zero().

Referenced by insert_slice().

00449 {
00450         if (tmp_data != 0 ) tmp_data->to_zero();
00451         if (image != 0 ) image->to_zero();
00452 
00453         EMData* start_model = params.set_default("start_model",(EMData*)0);
00454 
00455         if (start_model != 0) {
00456                 if (!start_model->is_complex()) {
00457                         start_model->do_fft_inplace();
00458                         start_model->process_inplace("xform.phaseorigin.tocenter");
00459                 }
00460 
00461                 if (start_model->get_xsize() != nx || start_model->get_ysize() != ny || start_model->get_zsize() != nz ) {
00462                         throw ImageDimensionException("The dimensions of the start_model are incorrect");
00463                 }
00464 
00465                 if (!start_model->is_shuffled()) {
00466                         start_model->process_inplace("xform.fourierorigin.tocenter");
00467                 }
00468 
00469                 memcpy(image->get_data(),start_model->get_data(),nx*ny*nz*sizeof(float));
00470 
00471                 float start_model_weight =  params.set_default("start_model_weight",1.0f);
00472 
00473                 if (start_model_weight != 1.0) {
00474                         image->mult(start_model_weight);
00475                 }
00476 
00477                 tmp_data->add(start_model_weight);
00478         }
00479 }

EMData * FourierReconstructor::preprocess_slice	(	const EMData *const	slice,
		const Transform &	t = `Transform()`
	)			`[protected]`

Preprocess the slice prior to insertion into the 3D volume this Fourier tranforms the slice and make sure all the pixels are in the right position.

Returns:: the processed slice

Parameters:

	slice	the slice to be prepocessed
	t	transform

Exceptions:

InvalidValueException

when the specified padding value is less than the size of the images

Definition at line 412 of file reconstructor.cpp.

References EMAN::EMData::do_fft_inplace(), EMAN::Transform::get_mirror(), EMAN::Transform::get_scale(), EMAN::Transform::get_trans_2d(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::mult(), EMAN::EMData::process(), EMAN::EMData::process_inplace(), EMAN::Transform::set_rotation(), and sqrt().

Referenced by determine_slice_agreement(), and insert_slice().

00413 {
00414         // Shift the image pixels so the real space origin is now located at the phase origin (at the bottom left of the image)
00415         EMData* return_slice = 0;
00416         Transform tmp(t);
00417         tmp.set_rotation(Dict("type","eman")); // resets the rotation to 0 implicitly
00418 
00419         Vec2f trans = tmp.get_trans_2d();
00420         float scale = tmp.get_scale();
00421         bool mirror = tmp.get_mirror();
00422         if (trans[0] != 0 || trans[1] != 0 || scale != 1.0 ) {
00423                 return_slice = slice->process("math.transform",Dict("transform",&tmp));
00424         } else if ( mirror == true ) {
00425                 return_slice = slice->process("xform.flip",Dict("axis","x"));
00426         }
00427 
00428         if (return_slice == 0) {
00429                 return_slice = slice->process("xform.phaseorigin.tocorner");
00430         } else {
00431                 return_slice->process_inplace("xform.phaseorigin.tocorner");
00432         }
00433 
00434 //      EMData* return_slice = slice->process("normalize.edgemean");
00435 //      return_slice->process_inplace("xform.phaseorigin.tocorner");
00436 
00437         // Fourier transform the slice
00438         return_slice->do_fft_inplace();
00439 
00440         return_slice->mult((float)sqrt(1.0f/(return_slice->get_ysize())*return_slice->get_xsize()));
00441 
00442         // Shift the Fourier transform so that it's origin is in the center (bottom) of the image.
00443         return_slice->process_inplace("xform.fourierorigin.tocenter");
00444 
00445         return return_slice;
00446 }

void EMAN::FourierReconstructor::load_default_settings ( ) [inline, protected]

Load default settings.

Reimplemented in EMAN::BaldwinWoolfordReconstructor.

Definition at line 487 of file reconstructor.h.

Referenced by FourierReconstructor().

00488                 {
00489 //                      params.put("x_in", = 0;
00490 //                      params.put("y_in", = 0;
00491 //                      params["zsample"] = 0;
00492 //                      params["ysample"] = 0;
00493 //                      params["xsample"] = 0;
00494 //                      params["mode"] = 2;
00495 //                      params["hard"] = 0.05;
00496 //                      params["sym"] = "c1";
00497 //                      params["quiet"] = true;
00498 //                      params["3damp"] = false;
00499                 }

void FourierReconstructor::free_memory ( ) [protected]

Frees the memory owned by this object (but not parent objects) Deletes the FourierPixelInserter3D pointer.

Reimplemented from EMAN::ReconstructorVolumeData.

Definition at line 258 of file reconstructor.cpp.

References inserter, and interp_FRC_calculator.

Referenced by ~FourierReconstructor().

00259 {
00260         if ( inserter != 0 )
00261         {
00262                 delete inserter;
00263                 inserter = 0;
00264         }
00265         if ( interp_FRC_calculator != 0 )
00266         {
00267                 delete interp_FRC_calculator;
00268                 interp_FRC_calculator = 0;
00269         }
00270 }

void FourierReconstructor::load_inserter ( ) [protected]

Load the pixel inserter based on the information in params.

Definition at line 274 of file reconstructor.cpp.

References EMAN::EMData::get_data(), EMAN::ReconstructorVolumeData::image, EMAN::FourierPixelInserter3D::init(), inserter, EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, and EMAN::ReconstructorVolumeData::tmp_data.

Referenced by do_insert_slice_work(), and setup().

00275 {
00276         // ints get converted to strings byt the Dict object here
00277 
00278         stringstream ss;
00279         ss << (int)params["mode"];
00280         string mode;
00281         ss >> mode;
00282 //      ss
00283 //      string mode = (string)params["mode"];
00284         Dict parms;
00285         parms["rdata"] = image->get_data();
00286         parms["norm"] = tmp_data->get_data();
00287         parms["nx"] = nx;
00288         parms["ny"] = ny;
00289         parms["nz"] = nz;
00290 
00291         if ( inserter != 0 )
00292         {
00293                 delete inserter;
00294         }
00295 
00296         inserter = Factory<FourierPixelInserter3D>::get(mode, parms);
00297         inserter->init();
00298 }

void FourierReconstructor::load_interp_FRC_calculator ( ) [protected]

Load the pixel inserter based on the information in params.

Definition at line 302 of file reconstructor.cpp.

References EMAN::EMData::get_data(), EMAN::ReconstructorVolumeData::image, interp_FRC_calculator, EMAN::ReconstructorVolumeData::nx, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, EMAN::ReconstructorVolumeData::tmp_data, x_scale_factor, y_scale_factor, and z_scale_factor.

Referenced by setup().

00303 {
00304         Dict init_parms;
00305         init_parms["rdata"] = image->get_data();
00306         init_parms["norm"] = tmp_data->get_data();
00307         init_parms["nx"] = nx;
00308         init_parms["ny"] = ny;
00309         init_parms["nz"] = nz;
00310 
00311         if ( x_scale_factor != 0 ) init_parms["x_scale"] = 1.0/x_scale_factor;
00312         if ( y_scale_factor != 0 ) init_parms["y_scale"] = 1.0/y_scale_factor;
00313         if ( z_scale_factor != 0 ) init_parms["z_scale"] = 1.0/z_scale_factor;
00314 
00315         if ( interp_FRC_calculator != 0 )
00316         {
00317                 delete interp_FRC_calculator;
00318         }
00319 
00320         stringstream ss;
00321         ss << (int)params["mode"];
00322         string mode;
00323         ss >> mode;
00324 
00325         interp_FRC_calculator = Factory<InterpolatedFRC>::get(mode, init_parms);
00326 }

void FourierReconstructor::do_insert_slice_work	(	const EMData *const	input_slice,
		const Transform &	euler
	)			`[protected]`

A function to perform the nuts and bolts of inserting an image slice.

Parameters:

	input_slice	the slice to insert into the 3D volume
	euler	a transform3D storing the slice euler angle

Definition at line 545 of file reconstructor.cpp.

References dt, EMAN::EMData::get_data(), EMAN::FactoryBase::get_name(), EMAN::Symmetry3D::get_symmetries(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::FourierPixelInserter3D::insert_pixel(), inserter, EMAN::EMData::is_fftodd(), load_inserter(), max_input_dim, EMAN::ReconstructorVolumeData::ny, EMAN::ReconstructorVolumeData::nz, EMAN::FactoryBase::params, EMAN::Util::square(), square, weight, x, x_scale_factor, y, y_scale_factor, and z_scale_factor.

Referenced by insert_slice().

00546 {
00547         // Reload the inserter if the mode has changed
00548         string mode = (string) params["mode"];
00549         if ( mode != inserter->get_name() )     load_inserter();
00550 
00551         int rl = Util::square( max_input_dim / 2);
00552 
00553         float dt[2];
00554 
00555         float y_scale = 0, x_scale = 0;
00556 
00557         int y_in = input_slice->get_ysize();
00558         int x_in = input_slice->get_xsize();
00559         // Adjust the dimensions to account for odd and even ffts
00560         if (input_slice->is_fftodd()) x_in -= 1;
00561         else x_in -= 2;
00562 
00563         if ( y_in != x_in )
00564         {
00565                 if ( x_in > y_in ) y_scale = (float) x_in / (float) y_in;
00566                 else x_scale = (float) y_in / (float) x_in;
00567         }
00568 
00569         float *dat = input_slice->get_data();
00570         vector<Transform> syms = Symmetry3D::get_symmetries((string)params["sym"]);
00571         float weight = params.set_default("weight",1.0f);
00572 
00573         for ( vector<Transform>::const_iterator it = syms.begin(); it != syms.end(); ++it ) {
00574                 Transform t3d = arg*(*it);
00575                 for (int y = 0; y < input_slice->get_ysize(); y++) {
00576                         for (int x = 0; x < input_slice->get_xsize() / 2; x++) {
00577 
00578                                 float rx = (float) x;
00579                                 float ry = (float) y;
00580 
00581                                 if ( y_in != x_in )
00582                                 {
00583                                         if ( x_in > y_in ) ry *= y_scale;
00584                                         else rx *= x_scale;
00585                                 }
00586 
00587                                 if ((rx * rx + Util::square(ry - max_input_dim / 2)) > rl)
00588                                         continue;
00589 
00590                                 Vec3f coord(rx,(ry - max_input_dim / 2),0);
00591                                 coord = coord*t3d; // transpose multiplication
00592                                 float xx = coord[0];
00593                                 float yy = coord[1];
00594                                 float zz = coord[2];
00595 
00596                                 float cc = 1;
00597 
00598                                 if (xx < 0) {
00599                                         xx = -xx;
00600                                         yy = -yy;
00601                                         zz = -zz;
00602                                         cc = -1.0;
00603                                 }
00604 
00605                                 if ( z_scale_factor != 0 ) zz *= z_scale_factor;
00606                                 if ( y_scale_factor != 0 ) yy *= y_scale_factor;
00607                                 if ( x_scale_factor != 0 ) xx *= x_scale_factor;
00608 
00609                                 yy += ny / 2;
00610                                 zz += nz / 2;
00611 
00612                                 int idx = x * 2 + y * (input_slice->get_xsize());
00613                                 dt[0] = dat[idx];
00614                                 dt[1] = cc * dat[idx+1];
00615 
00616                                 inserter->insert_pixel(xx,yy,zz,dt,weight);
00617                         }
00618                 }
00619         }
00620 }

void FourierReconstructor::print_stats ( const vector< InterpolatedFRC::QualityScores > & scores ) [protected]

print stats is called internally at various points in the reconstruction routine and is for the benefit of people using e2make3d.py

Definition at line 785 of file reconstructor.cpp.

References EMAN::FactoryBase::params, prev_quality_scores, and sqrt().

00786 {
00787         if ( prev_quality_scores.size() == 0 )
00788         {
00789                 //cout << "No quality scores present in FourierReconstructor::print_stats, nothing to print" << endl;
00790                 return;
00791         }
00792 
00793         unsigned int size = scores.size();
00794 
00795         unsigned int contributing_images = 0;
00796         for( unsigned int i = 0; i < size; ++ i )
00797         {
00798                 if (scores[i].get_snr_normed_frc_integral() < (float) params["hard"])
00799                         continue;
00800                 else contributing_images++;
00801         }
00802 
00803         double* norm_frc = new double[contributing_images];
00804         double* frc = new double[contributing_images];
00805         double* norm_snr = new double[contributing_images];
00806 
00807         unsigned int idx = 0;
00808         for( unsigned int i = 0; i < size; ++ i )
00809         {
00810                 if (scores[i].get_snr_normed_frc_integral() < (float) params["hard"])
00811                         continue;
00812 
00813                 norm_frc[idx] = scores[i].get_snr_normed_frc_integral();
00814                 frc[idx] = scores[i].get_frc_integral();
00815                 norm_snr[idx] = scores[i].get_normed_snr_integral();
00816 
00817                 ++idx;
00818         }
00819 
00820         double mean = gsl_stats_mean(norm_frc, 1, contributing_images);
00821         double variance = gsl_stats_variance_m(norm_frc, 1, contributing_images, mean);
00822 
00823         cout << "Normalized FRC mean " << mean << " std dev " << sqrt(variance) << endl;
00824 
00825         mean = gsl_stats_mean(frc, 1, contributing_images);
00826         variance = gsl_stats_variance_m(frc, 1, contributing_images, mean);
00827 
00828         cout << "FRC mean " << mean << " std dev " << sqrt(variance) << endl;
00829 
00830         mean = gsl_stats_mean(norm_snr, 1, contributing_images);
00831         variance = gsl_stats_variance_m(norm_snr, 1, contributing_images, mean);
00832         cout << "SNR mean " << mean << " std dev " << sqrt(variance) << endl;
00833 
00834         double c0, c1, cov00, cov01, cov11, sumsq;
00835         gsl_fit_linear (norm_frc, 1, frc, 1, contributing_images, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
00836         cout << "The correlation between frc and norm_frc is " << c0 << " + " << c1 << "x" << endl;
00837 
00838         gsl_fit_linear (norm_frc, 1, norm_snr, 1, contributing_images, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
00839         cout << "The correlation between norm_snr and norm_frc is " << c0 << " + " << c1 << "x" << endl;
00840 
00841         gsl_fit_linear (norm_snr, 1, frc, 1, contributing_images, &c0, &c1, &cov00, &cov01, &cov11, &sumsq);
00842         cout << "The correlation between frc and norm_snr is " << c0 << " + " << c1 << "x" << endl;
00843 
00844         delete [] norm_frc;
00845         delete [] frc;
00846         delete [] norm_snr;
00847 }