#include <cmp.h>
Inheritance diagram for EMAN::OptVarianceCmp:


Public Member Functions | |
| OptVarianceCmp () | |
| float | cmp (EMData *image, EMData *with) const |
| To compare 'image' with another image passed in through its parameters. | |
| string | get_name () const |
| Get the Cmp's name. | |
| string | get_desc () const |
| TypeDict | get_param_types () const |
| Get Cmp parameter information in a dictionary. | |
| float | get_scale () const |
| float | get_shift () const |
Static Public Member Functions | |
| Cmp * | NEW () |
Static Public Attributes | |
| const string | NAME = "optvariance" |
Private Attributes | |
| float | scale |
| float | shift |
Generally, 'this' should be noisy and 'with' should be less noisy. linear scaling (mx + b) of the densities in 'this' is performed to produce the smallest possible variance between images.
If keepzero is set, then zero pixels are left at zero (b is not added). if matchfilt is set, then 'with' is filtered so its radial power spectrum matches 'this' If invert is set, then (y-b)/m is applied to the second image rather than mx+b to the first.
To modify 'this' to match the operation performed here, scale should be applied first, then b should be added
Definition at line 403 of file cmp.h.
|
|
Definition at line 406 of file cmp.h.
|
|
||||||||||||
|
To compare 'image' with another image passed in through its parameters. An optional transformation may be used to transform the 2 images.
Implements EMAN::Cmp. Definition at line 733 of file cmp.cpp. References EMAN::EMData::ap2ri(), EMAN::EMData::apply_radial_func(), b, EMAN::Util::calc_least_square_fit(), EMAN::EMData::calc_radial_dist(), EMAN::EMData::do_fft(), EMAN::EMData::do_ift(), EMAN::EMData::get_const_data(), EMAN::EMData::get_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ImageDimensionException, nx, EMAN::EMData::ri2ap(), scale, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), shift, square, EMAN::EMData::update(), EMAN::Cmp::validate_input_args(), EMAN::EMData::write_image(), x, and y. 00734 {
00735 ENTERFUNC;
00736 validate_input_args(image, with);
00737
00738 int keepzero = params.set_default("keepzero", 1);
00739 int invert = params.set_default("invert",0);
00740 int matchfilt = params.set_default("matchfilt",1);
00741 int matchamp = params.set_default("matchamp",0);
00742 int radweight = params.set_default("radweight",0);
00743 int dbug = params.set_default("debug",0);
00744
00745 size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00746
00747
00748 EMData *with2=NULL;
00749 if (matchfilt) {
00750 EMData *a = image->do_fft();
00751 EMData *b = with->do_fft();
00752
00753 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1);
00754 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1);
00755
00756 float avg=0;
00757 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00758 rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i]));
00759 avg+=rfa[i];
00760 }
00761
00762 avg/=a->get_ysize()/2.0f;
00763 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00764 if (rfa[i]>avg*10.0) rfa[i]=10.0; // If some particular location has a small but non-zero value, we don't want to overcorrect it
00765 }
00766 rfa[0]=0.0;
00767
00768 if (dbug) b->write_image("a.hdf",-1);
00769
00770 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa);
00771 with2=b->do_ift();
00772
00773 if (dbug) b->write_image("a.hdf",-1);
00774 if (dbug) a->write_image("a.hdf",-1);
00775
00776 /* if (dbug) {
00777 FILE *out=fopen("a.txt","w");
00778 for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]);
00779 fclose(out);
00780
00781 out=fopen("b.txt","w");
00782 for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]);
00783 fclose(out);
00784 }*/
00785
00786
00787 delete a;
00788 delete b;
00789
00790 if (dbug) {
00791 with2->write_image("a.hdf",-1);
00792 image->write_image("a.hdf",-1);
00793 }
00794
00795 // with2->process_inplace("matchfilt",Dict("to",this));
00796 // x_data = with2->get_data();
00797 }
00798
00799 // This applies the individual Fourier amplitudes from 'image' and
00800 // applies them to 'with'
00801 if (matchamp) {
00802 EMData *a = image->do_fft();
00803 EMData *b = with->do_fft();
00804 size_t size2 = (size_t)a->get_xsize() * a->get_ysize() * a->get_zsize();
00805
00806 a->ri2ap();
00807 b->ri2ap();
00808
00809 const float *const ad=a->get_const_data();
00810 float * bd=b->get_data();
00811
00812 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i];
00813 b->update();
00814
00815 b->ap2ri();
00816 with2=b->do_ift();
00817 //with2->write_image("a.hdf",-1);
00818 delete a;
00819 delete b;
00820 }
00821
00822 const float * x_data;
00823 if (with2) x_data=with2->get_const_data();
00824 else x_data = with->get_const_data();
00825 const float *const y_data = image->get_const_data();
00826
00827 size_t nx = image->get_xsize();
00828 float m = 0;
00829 float b = 0;
00830
00831 // This will write the x vs y file used to calculate the density
00832 // optimization. This behavior may change in the future
00833 if (dbug) {
00834 FILE *out=fopen("dbug.optvar.txt","w");
00835 if (out) {
00836 for (size_t i=0; i<size; i++) {
00837 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]);
00838 }
00839 fclose(out);
00840 }
00841 }
00842
00843
00844 Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero);
00845 if (m == 0) {
00846 m = FLT_MIN;
00847 }
00848 b = -b / m;
00849 m = 1.0f / m;
00850
00851 // While negative slopes are really not a valid comparison in most cases, we
00852 // still want to detect these instances, so this if is removed
00853 /* if (m < 0) {
00854 b = 0;
00855 m = 1000.0;
00856 }*/
00857
00858 double result = 0;
00859 int count = 0;
00860
00861 if (radweight) {
00862 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only");
00863 if (keepzero) {
00864 for (size_t i = 0,y=0; i < size; y++) {
00865 for (size_t x=0; x<nx; i++,x++) {
00866 if (y_data[i] && x_data[i]) {
00867 #ifdef _WIN32
00868 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00869 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00870 #else
00871 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00872 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00873 #endif
00874 count++;
00875 }
00876 }
00877 }
00878 result/=count;
00879 }
00880 else {
00881 for (size_t i = 0,y=0; i < size; y++) {
00882 for (size_t x=0; x<nx; i++,x++) {
00883 #ifdef _WIN32
00884 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00885 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00886 #else
00887 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00888 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00889 #endif
00890 }
00891 }
00892 result = result / size;
00893 }
00894 }
00895 else {
00896 if (keepzero) {
00897 for (size_t i = 0; i < size; i++) {
00898 if (y_data[i] && x_data[i]) {
00899 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00900 else result += Util::square((x_data[i] * m) + b - y_data[i]);
00901 count++;
00902 }
00903 }
00904 result/=count;
00905 }
00906 else {
00907 for (size_t i = 0; i < size; i++) {
00908 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00909 else result += Util::square((x_data[i] * m) + b - y_data[i]);
00910 }
00911 result = result / size;
00912 }
00913 }
00914 scale = m;
00915 shift = b;
00916
00917 image->set_attr("ovcmp_m",m);
00918 image->set_attr("ovcmp_b",b);
00919 if (with2) delete with2;
00920 EXITFUNC;
00921
00922 #if 0
00923 return (1 - result);
00924 #endif
00925
00926 return static_cast<float>(result);
00927 }
|
|
|
Implements EMAN::Cmp. Definition at line 415 of file cmp.h. 00416 {
00417 return "Real-space variance after density optimization, self should be noisy and target less noisy. Linear transform applied to density to minimize variance.";
00418 }
|
|
|
Get the Cmp's name. Each Cmp is identified by a unique name.
Implements EMAN::Cmp. Definition at line 410 of file cmp.h. 00411 {
00412 return NAME;
00413 }
|
|
|
Get Cmp parameter information in a dictionary. Each parameter has one record in the dictionary. Each record contains its name, data-type, and description.
Implements EMAN::Cmp. Definition at line 425 of file cmp.h. References EMAN::TypeDict::put(). 00426 {
00427 TypeDict d;
00428 d.put("invert", EMObject::INT, "If set, 'with' is rescaled rather than 'this'. 'this' should still be the noisier image. (default=0)");
00429 d.put("keepzero", EMObject::INT, "If set, zero pixels will not be adjusted in the linear density optimization. (default=1)");
00430 d.put("matchfilt", EMObject::INT, "If set, with will be filtered so its radial power spectrum matches 'this' before density optimization of this. (default=1)");
00431 d.put("matchamp", EMObject::INT, "Takes per-pixel Fourier amplitudes from self and imposes them on the target, but leaves the phases alone. (default=0)");
00432 d.put("radweight", EMObject::INT, "Upweight variances closer to the edge of the image. (default=0)");
00433 d.put("debug", EMObject::INT, "Performs various debugging actions if set.");
00434 return d;
00435 }
|
|
|
Definition at line 437 of file cmp.h. 00438 {
00439 return scale;
00440 }
|
|
|
Definition at line 442 of file cmp.h. 00443 {
00444 return shift;
00445 }
|
|
|
Definition at line 420 of file cmp.h. 00421 {
00422 return new OptVarianceCmp();
00423 }
|
|
|
|
|
|
Definition at line 450 of file cmp.h. Referenced by cmp(). |
|
|
Definition at line 451 of file cmp.h. Referenced by cmp(). |
1.3.9.1