#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 | |
static Cmp * | NEW () |
Static Public Attributes | |
static 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.
EMAN::OptVarianceCmp::OptVarianceCmp | ( | ) | [inline] |
To compare 'image' with another image passed in through its parameters.
An optional transformation may be used to transform the 2 images.
image | The first image to be compared. | |
with | The second image to be comppared. |
Implements EMAN::Cmp.
Definition at line 733 of file cmp.cpp.
References b, EMAN::Util::calc_least_square_fit(), EMAN::EMData::calc_radial_dist(), EMAN::EMData::do_fft(), ENTERFUNC, EXITFUNC, EMAN::EMData::get_const_data(), EMAN::EMData::get_xsize(), EMAN::EMData::get_ysize(), EMAN::EMData::get_zsize(), ImageDimensionException, nx, EMAN::Cmp::params, EMAN::EMData::ri2ap(), scale, EMAN::EMData::set_attr(), EMAN::Dict::set_default(), shift, square, EMAN::Util::square(), 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 }
string EMAN::OptVarianceCmp::get_desc | ( | ) | const [inline, virtual] |
string EMAN::OptVarianceCmp::get_name | ( | ) | const [inline, virtual] |
TypeDict EMAN::OptVarianceCmp::get_param_types | ( | ) | const [inline, virtual] |
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::EMObject::INT, and 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 }
float EMAN::OptVarianceCmp::get_scale | ( | ) | const [inline] |
float EMAN::OptVarianceCmp::get_shift | ( | ) | const [inline] |
static Cmp* EMAN::OptVarianceCmp::NEW | ( | ) | [inline, static] |
Definition at line 420 of file cmp.h.
References OptVarianceCmp().
00421 { 00422 return new OptVarianceCmp(); 00423 }
const string OptVarianceCmp::NAME = "optvariance" [static] |
float EMAN::OptVarianceCmp::scale [mutable, private] |
float EMAN::OptVarianceCmp::shift [mutable, private] |