#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 358 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 566 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.
00567 { 00568 ENTERFUNC; 00569 validate_input_args(image, with); 00570 00571 int keepzero = params.set_default("keepzero", 1); 00572 int invert = params.set_default("invert",0); 00573 int matchfilt = params.set_default("matchfilt",1); 00574 int matchamp = params.set_default("matchamp",0); 00575 int radweight = params.set_default("radweight",0); 00576 int dbug = params.set_default("debug",0); 00577 00578 size_t size = image->get_xsize() * image->get_ysize() * image->get_zsize(); 00579 00580 00581 EMData *with2=NULL; 00582 if (matchfilt) { 00583 EMData *a = image->do_fft(); 00584 EMData *b = with->do_fft(); 00585 00586 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1); 00587 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1); 00588 00589 float avg=0; 00590 for (size_t i=0; i<a->get_ysize()/2.0f; i++) { 00591 rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i])); 00592 avg+=rfa[i]; 00593 } 00594 00595 avg/=a->get_ysize()/2.0f; 00596 for (size_t i=0; i<a->get_ysize()/2.0f; i++) { 00597 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 00598 } 00599 rfa[0]=0.0; 00600 00601 if (dbug) b->write_image("a.hdf",-1); 00602 00603 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa); 00604 with2=b->do_ift(); 00605 00606 if (dbug) b->write_image("a.hdf",-1); 00607 if (dbug) a->write_image("a.hdf",-1); 00608 00609 /* if (dbug) { 00610 FILE *out=fopen("a.txt","w"); 00611 for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]); 00612 fclose(out); 00613 00614 out=fopen("b.txt","w"); 00615 for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]); 00616 fclose(out); 00617 }*/ 00618 00619 00620 delete a; 00621 delete b; 00622 00623 if (dbug) { 00624 with2->write_image("a.hdf",-1); 00625 image->write_image("a.hdf",-1); 00626 } 00627 00628 // with2->process_inplace("matchfilt",Dict("to",this)); 00629 // x_data = with2->get_data(); 00630 } 00631 00632 // This applies the individual Fourier amplitudes from 'image' and 00633 // applies them to 'with' 00634 if (matchamp) { 00635 EMData *a = image->do_fft(); 00636 EMData *b = with->do_fft(); 00637 size_t size2 = a->get_xsize() * a->get_ysize() * a->get_zsize(); 00638 00639 a->ri2ap(); 00640 b->ri2ap(); 00641 00642 const float *const ad=a->get_const_data(); 00643 float * bd=b->get_data(); 00644 00645 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i]; 00646 b->update(); 00647 00648 b->ap2ri(); 00649 with2=b->do_ift(); 00650 //with2->write_image("a.hdf",-1); 00651 delete a; 00652 delete b; 00653 } 00654 00655 const float * x_data; 00656 if (with2) x_data=with2->get_const_data(); 00657 else x_data = with->get_const_data(); 00658 const float *const y_data = image->get_const_data(); 00659 00660 size_t nx = image->get_xsize(); 00661 float m = 0; 00662 float b = 0; 00663 00664 // This will write the x vs y file used to calculate the density 00665 // optimization. This behavior may change in the future 00666 if (dbug) { 00667 FILE *out=fopen("dbug.optvar.txt","w"); 00668 if (out) { 00669 for (size_t i=0; i<size; i++) { 00670 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]); 00671 } 00672 fclose(out); 00673 } 00674 } 00675 00676 00677 Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero); 00678 if (m == 0) { 00679 m = FLT_MIN; 00680 } 00681 b = -b / m; 00682 m = 1.0f / m; 00683 00684 // While negative slopes are really not a valid comparison in most cases, we 00685 // still want to detect these instances, so this if is removed 00686 /* if (m < 0) { 00687 b = 0; 00688 m = 1000.0; 00689 }*/ 00690 00691 double result = 0; 00692 int count = 0; 00693 00694 if (radweight) { 00695 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only"); 00696 if (keepzero) { 00697 for (size_t i = 0,y=0; i < size; y++) { 00698 for (size_t x=0; x<nx; i++,x++) { 00699 if (y_data[i] && x_data[i]) { 00700 #ifdef _WIN32 00701 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0); 00702 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0); 00703 #else 00704 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0); 00705 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0); 00706 #endif 00707 count++; 00708 } 00709 } 00710 } 00711 result/=count; 00712 } 00713 else { 00714 for (size_t i = 0,y=0; i < size; y++) { 00715 for (size_t x=0; x<nx; i++,x++) { 00716 #ifdef _WIN32 00717 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0); 00718 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0); 00719 #else 00720 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0); 00721 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0); 00722 #endif 00723 } 00724 } 00725 result = result / size; 00726 } 00727 } 00728 else { 00729 if (keepzero) { 00730 for (size_t i = 0; i < size; i++) { 00731 if (y_data[i] && x_data[i]) { 00732 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m); 00733 else result += Util::square((x_data[i] * m) + b - y_data[i]); 00734 count++; 00735 } 00736 } 00737 result/=count; 00738 } 00739 else { 00740 for (size_t i = 0; i < size; i++) { 00741 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m); 00742 else result += Util::square((x_data[i] * m) + b - y_data[i]); 00743 } 00744 result = result / size; 00745 } 00746 } 00747 scale = m; 00748 shift = b; 00749 00750 image->set_attr("ovcmp_m",m); 00751 image->set_attr("ovcmp_b",b); 00752 if (with2) delete with2; 00753 EXITFUNC; 00754 00755 #if 0 00756 return (1 - result); 00757 #endif 00758 00759 return static_cast<float>(result); 00760 }
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 380 of file cmp.h.
References EMAN::EMObject::INT, and EMAN::TypeDict::put().
00381 { 00382 TypeDict d; 00383 d.put("invert", EMObject::INT, "If set, 'with' is rescaled rather than 'this'. 'this' should still be the noisier image. (default=0)"); 00384 d.put("keepzero", EMObject::INT, "If set, zero pixels will not be adjusted in the linear density optimization. (default=1)"); 00385 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)"); 00386 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)"); 00387 d.put("radweight", EMObject::INT, "Upweight variances closer to the edge of the image. (default=0)"); 00388 d.put("debug", EMObject::INT, "Performs various debugging actions if set."); 00389 return d; 00390 }
float EMAN::OptVarianceCmp::get_scale | ( | ) | const [inline] |
float EMAN::OptVarianceCmp::get_shift | ( | ) | const [inline] |
static Cmp* EMAN::OptVarianceCmp::NEW | ( | ) | [inline, static] |
Definition at line 375 of file cmp.h.
References OptVarianceCmp().
00376 { 00377 return new OptVarianceCmp(); 00378 }
const string OptVarianceCmp::NAME = "optvariance" [static] |
float EMAN::OptVarianceCmp::scale [mutable, private] |
float EMAN::OptVarianceCmp::shift [mutable, private] |