#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 447 of file cmp.h.
|
Definition at line 450 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 860 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. 00861 { 00862 ENTERFUNC; 00863 validate_input_args(image, with); 00864 00865 int keepzero = params.set_default("keepzero", 1); 00866 int invert = params.set_default("invert",0); 00867 int matchfilt = params.set_default("matchfilt",1); 00868 int matchamp = params.set_default("matchamp",0); 00869 int radweight = params.set_default("radweight",0); 00870 int dbug = params.set_default("debug",0); 00871 00872 size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize(); 00873 00874 00875 EMData *with2=NULL; 00876 if (matchfilt) { 00877 EMData *a = image->do_fft(); 00878 EMData *b = with->do_fft(); 00879 00880 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1); 00881 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1); 00882 00883 float avg=0; 00884 for (size_t i=0; i<a->get_ysize()/2.0f; i++) { 00885 rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i])); 00886 avg+=rfa[i]; 00887 } 00888 00889 avg/=a->get_ysize()/2.0f; 00890 for (size_t i=0; i<a->get_ysize()/2.0f; i++) { 00891 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 00892 } 00893 rfa[0]=0.0; 00894 00895 if (dbug) b->write_image("a.hdf",-1); 00896 00897 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa); 00898 with2=b->do_ift(); 00899 00900 if (dbug) b->write_image("a.hdf",-1); 00901 if (dbug) a->write_image("a.hdf",-1); 00902 00903 /* if (dbug) { 00904 FILE *out=fopen("a.txt","w"); 00905 for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]); 00906 fclose(out); 00907 00908 out=fopen("b.txt","w"); 00909 for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]); 00910 fclose(out); 00911 }*/ 00912 00913 00914 delete a; 00915 delete b; 00916 00917 if (dbug) { 00918 with2->write_image("a.hdf",-1); 00919 image->write_image("a.hdf",-1); 00920 } 00921 00922 // with2->process_inplace("matchfilt",Dict("to",this)); 00923 // x_data = with2->get_data(); 00924 } 00925 00926 // This applies the individual Fourier amplitudes from 'image' and 00927 // applies them to 'with' 00928 if (matchamp) { 00929 EMData *a = image->do_fft(); 00930 EMData *b = with->do_fft(); 00931 size_t size2 = (size_t)a->get_xsize() * a->get_ysize() * a->get_zsize(); 00932 00933 a->ri2ap(); 00934 b->ri2ap(); 00935 00936 const float *const ad=a->get_const_data(); 00937 float * bd=b->get_data(); 00938 00939 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i]; 00940 b->update(); 00941 00942 b->ap2ri(); 00943 with2=b->do_ift(); 00944 //with2->write_image("a.hdf",-1); 00945 delete a; 00946 delete b; 00947 } 00948 00949 const float * x_data; 00950 if (with2) x_data=with2->get_const_data(); 00951 else x_data = with->get_const_data(); 00952 const float *const y_data = image->get_const_data(); 00953 00954 size_t nx = image->get_xsize(); 00955 float m = 0; 00956 float b = 0; 00957 00958 // This will write the x vs y file used to calculate the density 00959 // optimization. This behavior may change in the future 00960 if (dbug) { 00961 FILE *out=fopen("dbug.optvar.txt","w"); 00962 if (out) { 00963 for (size_t i=0; i<size; i++) { 00964 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]); 00965 } 00966 fclose(out); 00967 } 00968 } 00969 00970 00971 Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero); 00972 if (m == 0) { 00973 m = FLT_MIN; 00974 } 00975 b = -b / m; 00976 m = 1.0f / m; 00977 00978 // While negative slopes are really not a valid comparison in most cases, we 00979 // still want to detect these instances, so this if is removed 00980 /* if (m < 0) { 00981 b = 0; 00982 m = 1000.0; 00983 }*/ 00984 00985 double result = 0; 00986 int count = 0; 00987 00988 if (radweight) { 00989 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only"); 00990 if (keepzero) { 00991 for (size_t i = 0,y=0; i < size; y++) { 00992 for (size_t x=0; x<nx; i++,x++) { 00993 if (y_data[i] && x_data[i]) { 00994 #ifdef _WIN32 00995 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0); 00996 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0); 00997 #else 00998 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0); 00999 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0); 01000 #endif 01001 count++; 01002 } 01003 } 01004 } 01005 result/=count; 01006 } 01007 else { 01008 for (size_t i = 0,y=0; i < size; y++) { 01009 for (size_t x=0; x<nx; i++,x++) { 01010 #ifdef _WIN32 01011 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0); 01012 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0); 01013 #else 01014 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0); 01015 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0); 01016 #endif 01017 } 01018 } 01019 result = result / size; 01020 } 01021 } 01022 else { 01023 if (keepzero) { 01024 for (size_t i = 0; i < size; i++) { 01025 if (y_data[i] && x_data[i]) { 01026 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m); 01027 else result += Util::square((x_data[i] * m) + b - y_data[i]); 01028 count++; 01029 } 01030 } 01031 result/=count; 01032 } 01033 else { 01034 for (size_t i = 0; i < size; i++) { 01035 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m); 01036 else result += Util::square((x_data[i] * m) + b - y_data[i]); 01037 } 01038 result = result / size; 01039 } 01040 } 01041 scale = m; 01042 shift = b; 01043 01044 image->set_attr("ovcmp_m",m); 01045 image->set_attr("ovcmp_b",b); 01046 if (with2) delete with2; 01047 EXITFUNC; 01048 01049 #if 0 01050 return (1 - result); 01051 #endif 01052 01053 return static_cast<float>(result); 01054 }
|
|
Implements EMAN::Cmp. Definition at line 459 of file cmp.h. 00460 { 00461 return "Real-space variance after density optimization, self should be noisy and target less noisy. Linear transform applied to density to minimize variance."; 00462 }
|
|
Get the Cmp's name. Each Cmp is identified by a unique name.
Implements EMAN::Cmp. Definition at line 454 of file cmp.h. 00455 {
00456 return NAME;
00457 }
|
|
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 469 of file cmp.h. References EMAN::TypeDict::put(). 00470 { 00471 TypeDict d; 00472 d.put("invert", EMObject::INT, "If set, 'with' is rescaled rather than 'this'. 'this' should still be the noisier image. (default=0)"); 00473 d.put("keepzero", EMObject::INT, "If set, zero pixels will not be adjusted in the linear density optimization. (default=1)"); 00474 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)"); 00475 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)"); 00476 d.put("radweight", EMObject::INT, "Upweight variances closer to the edge of the image. (default=0)"); 00477 d.put("debug", EMObject::INT, "Performs various debugging actions if set."); 00478 return d; 00479 }
|
|
Definition at line 481 of file cmp.h. 00482 {
00483 return scale;
00484 }
|
|
Definition at line 486 of file cmp.h. 00487 {
00488 return shift;
00489 }
|
|
Definition at line 464 of file cmp.h. 00465 { 00466 return new OptVarianceCmp(); 00467 }
|
|
|
|
Definition at line 494 of file cmp.h. Referenced by cmp(). |
|
Definition at line 495 of file cmp.h. Referenced by cmp(). |