#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 358 of file cmp.h.
|
Definition at line 361 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 563 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. 00564 { 00565 ENTERFUNC; 00566 validate_input_args(image, with); 00567 00568 int keepzero = params.set_default("keepzero", 1); 00569 int invert = params.set_default("invert",0); 00570 int matchfilt = params.set_default("matchfilt",1); 00571 int matchamp = params.set_default("matchamp",0); 00572 int radweight = params.set_default("radweight",0); 00573 int dbug = params.set_default("debug",0); 00574 00575 size_t size = image->get_xsize() * image->get_ysize() * image->get_zsize(); 00576 00577 00578 EMData *with2=NULL; 00579 if (matchfilt) { 00580 EMData *a = image->do_fft(); 00581 EMData *b = with->do_fft(); 00582 00583 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1); 00584 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1); 00585 00586 float avg=0; 00587 for (size_t i=0; i<a->get_ysize()/2.0f; i++) { 00588 rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i])); 00589 avg+=rfa[i]; 00590 } 00591 00592 avg/=a->get_ysize()/2.0f; 00593 for (size_t i=0; i<a->get_ysize()/2.0f; i++) { 00594 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 00595 } 00596 rfa[0]=0.0; 00597 00598 if (dbug) b->write_image("a.hdf",-1); 00599 00600 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa); 00601 with2=b->do_ift(); 00602 00603 if (dbug) b->write_image("a.hdf",-1); 00604 if (dbug) a->write_image("a.hdf",-1); 00605 00606 /* if (dbug) { 00607 FILE *out=fopen("a.txt","w"); 00608 for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]); 00609 fclose(out); 00610 00611 out=fopen("b.txt","w"); 00612 for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]); 00613 fclose(out); 00614 }*/ 00615 00616 00617 delete a; 00618 delete b; 00619 00620 if (dbug) { 00621 with2->write_image("a.hdf",-1); 00622 image->write_image("a.hdf",-1); 00623 } 00624 00625 // with2->process_inplace("matchfilt",Dict("to",this)); 00626 // x_data = with2->get_data(); 00627 } 00628 00629 // This applies the individual Fourier amplitudes from 'image' and 00630 // applies them to 'with' 00631 if (matchamp) { 00632 EMData *a = image->do_fft(); 00633 EMData *b = with->do_fft(); 00634 size_t size2 = a->get_xsize() * a->get_ysize() * a->get_zsize(); 00635 00636 a->ri2ap(); 00637 b->ri2ap(); 00638 00639 const float *const ad=a->get_const_data(); 00640 float * bd=b->get_data(); 00641 00642 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i]; 00643 b->update(); 00644 00645 b->ap2ri(); 00646 with2=b->do_ift(); 00647 //with2->write_image("a.hdf",-1); 00648 delete a; 00649 delete b; 00650 } 00651 00652 const float * x_data; 00653 if (with2) x_data=with2->get_const_data(); 00654 else x_data = with->get_const_data(); 00655 const float *const y_data = image->get_const_data(); 00656 00657 size_t nx = image->get_xsize(); 00658 float m = 0; 00659 float b = 0; 00660 00661 // This will write the x vs y file used to calculate the density 00662 // optimization. This behavior may change in the future 00663 if (dbug) { 00664 FILE *out=fopen("dbug.optvar.txt","w"); 00665 if (out) { 00666 for (size_t i=0; i<size; i++) { 00667 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]); 00668 } 00669 fclose(out); 00670 } 00671 } 00672 00673 00674 Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero); 00675 if (m == 0) { 00676 m = FLT_MIN; 00677 } 00678 b = -b / m; 00679 m = 1.0f / m; 00680 00681 // While negative slopes are really not a valid comparison in most cases, we 00682 // still want to detect these instances, so this if is removed 00683 /* if (m < 0) { 00684 b = 0; 00685 m = 1000.0; 00686 }*/ 00687 00688 double result = 0; 00689 int count = 0; 00690 00691 if (radweight) { 00692 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only"); 00693 if (keepzero) { 00694 for (size_t i = 0,y=0; i < size; y++) { 00695 for (size_t x=0; x<nx; i++,x++) { 00696 if (y_data[i] && x_data[i]) { 00697 #ifdef _WIN32 00698 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0); 00699 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0); 00700 #else 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 #endif 00704 count++; 00705 } 00706 } 00707 } 00708 result/=count; 00709 } 00710 else { 00711 for (size_t i = 0,y=0; i < size; y++) { 00712 for (size_t x=0; x<nx; i++,x++) { 00713 #ifdef _WIN32 00714 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0); 00715 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0); 00716 #else 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 #endif 00720 } 00721 } 00722 result = result / size; 00723 } 00724 } 00725 else { 00726 if (keepzero) { 00727 for (size_t i = 0; i < size; i++) { 00728 if (y_data[i] && x_data[i]) { 00729 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m); 00730 else result += Util::square((x_data[i] * m) + b - y_data[i]); 00731 count++; 00732 } 00733 } 00734 result/=count; 00735 } 00736 else { 00737 for (size_t i = 0; i < size; i++) { 00738 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m); 00739 else result += Util::square((x_data[i] * m) + b - y_data[i]); 00740 } 00741 result = result / size; 00742 } 00743 } 00744 scale = m; 00745 shift = b; 00746 00747 image->set_attr("ovcmp_m",m); 00748 image->set_attr("ovcmp_b",b); 00749 if (with2) delete with2; 00750 EXITFUNC; 00751 00752 #if 0 00753 return (1 - result); 00754 #endif 00755 00756 return static_cast<float>(result); 00757 }
|
|
Implements EMAN::Cmp. Definition at line 370 of file cmp.h. 00371 { 00372 return "Real-space variance after density optimization, self should be noisy and target less noisy. Linear transform applied to density to minimize variance."; 00373 }
|
|
Get the Cmp's name. Each Cmp is identified by a unique name.
Implements EMAN::Cmp. Definition at line 365 of file cmp.h. 00366 {
00367 return NAME;
00368 }
|
|
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::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 }
|
|
Definition at line 392 of file cmp.h. 00393 {
00394 return scale;
00395 }
|
|
Definition at line 397 of file cmp.h. 00398 {
00399 return shift;
00400 }
|
|
Definition at line 375 of file cmp.h. 00376 { 00377 return new OptVarianceCmp(); 00378 }
|
|
|
|
Definition at line 405 of file cmp.h. Referenced by cmp(). |
|
Definition at line 406 of file cmp.h. Referenced by cmp(). |