cmp.cpp

Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "cmp.h"
00037 #include "emdata.h"
00038 #include "ctf.h"
00039 #include "plugins/cmp_template.h"
00040 
00041 using namespace EMAN;
00042 
00043 const string CccCmp::NAME = "ccc";
00044 const string SqEuclideanCmp::NAME = "sqeuclidean";
00045 const string DotCmp::NAME = "dot";
00046 const string TomoDotCmp::NAME = "dot.tomo";
00047 const string QuadMinDotCmp::NAME = "quadmindot";
00048 const string OptVarianceCmp::NAME = "optvariance";
00049 const string PhaseCmp::NAME = "phase";
00050 const string FRCCmp::NAME = "frc";
00051 
00052 template <> Factory < Cmp >::Factory()
00053 {
00054         force_add<CccCmp>();
00055         force_add<SqEuclideanCmp>();
00056         force_add<DotCmp>();
00057         force_add<TomoDotCmp>();
00058         force_add<QuadMinDotCmp>();
00059         force_add<OptVarianceCmp>();
00060         force_add<PhaseCmp>();
00061         force_add<FRCCmp>();
00062 //      force_add<XYZCmp>();
00063 }
00064 
00065 void Cmp::validate_input_args(const EMData * image, const EMData *with) const
00066 {
00067         if (!image) {
00068                 throw NullPointerException("compared image");
00069         }
00070         if (!with) {
00071                 throw NullPointerException("compare-with image");
00072         }
00073 
00074         if (!EMUtil::is_same_size(image, with)) {
00075                 throw ImageFormatException( "images not same size");
00076         }
00077 
00078         float *d1 = image->get_data();
00079         if (!d1) {
00080                 throw NullPointerException("image contains no data");
00081         }
00082 
00083         float *d2 = with->get_data();
00084         if (!d2) {
00085                 throw NullPointerException("compare-with image data");
00086         }
00087 }
00088 
00089 //  It would be good to add code for complex images!  PAP
00090 float CccCmp::cmp(EMData * image, EMData *with) const
00091 {
00092         ENTERFUNC;
00093         if (image->is_complex() || with->is_complex())
00094                 throw ImageFormatException( "Complex images not supported by CMP::CccCmp");
00095         validate_input_args(image, with);
00096 
00097         const float *const d1 = image->get_const_data();
00098         const float *const d2 = with->get_const_data();
00099 
00100         float negative = (float)params.set_default("negative", 1);
00101         if (negative) negative=-1.0; else negative=1.0;
00102 
00103         double avg1 = 0.0, var1 = 0.0, avg2 = 0.0, var2 = 0.0, ccc = 0.0;
00104         long n = 0;
00105         size_t totsize = image->get_xsize()*image->get_ysize()*image->get_zsize();
00106 
00107         bool has_mask = false;
00108         EMData* mask = 0;
00109         if (params.has_key("mask")) {
00110                 mask = params["mask"];
00111                 if(mask!=0) {has_mask=true;}
00112         }
00113 
00114         if (has_mask) {
00115                 const float *const dm = mask->get_const_data();
00116                 for (size_t i = 0; i < totsize; ++i) {
00117                         if (dm[i] > 0.5) {
00118                                 avg1 += double(d1[i]);
00119                                 var1 += d1[i]*double(d1[i]);
00120                                 avg2 += double(d2[i]);
00121                                 var2 += d2[i]*double(d2[i]);
00122                                 ccc += d1[i]*double(d2[i]);
00123                                 n++;
00124                         }
00125                 }
00126         } else {
00127                 for (size_t i = 0; i < totsize; ++i) {
00128                         avg1 += double(d1[i]);
00129                         var1 += d1[i]*double(d1[i]);
00130                         avg2 += double(d2[i]);
00131                         var2 += d2[i]*double(d2[i]);
00132                         ccc += d1[i]*double(d2[i]);
00133                 }
00134                 n = totsize;
00135         }
00136 
00137         avg1 /= double(n);
00138         var1 = var1/double(n) - avg1*avg1;
00139         avg2 /= double(n);
00140         var2 = var2/double(n) - avg2*avg2;
00141         ccc = ccc/double(n) - avg1*avg2;
00142         ccc /= sqrt(var1*var2);
00143         ccc *= negative;
00144         return static_cast<float>(ccc);
00145         EXITFUNC;
00146 }
00147 
00148 
00149 
00150 //float SqEuclideanCmp::cmp(EMData * image, EMData *withorig) const
00151 float SqEuclideanCmp::cmp(EMData *image,EMData * withorig ) const
00152 {
00153         ENTERFUNC;
00154         EMData *with = withorig;
00155         validate_input_args(image, with);
00156 
00157         int zeromask = params.set_default("zeromask",0);
00158         int normto = params.set_default("normto",0);
00159 
00160         if (normto) {
00161                 if (zeromask) with = withorig->process("normalize.toimage",Dict("to",image));
00162                 else with = withorig->process("normalize.toimage",Dict("to",image,"ignore_zero",0));
00163                 with->set_attr("deleteme",1);
00164                 if ((float)(with->get_attr("norm_mult"))<=0) {          // This means the normalization inverted the image, a clear probablity of noise bias, so we undo the normalization
00165                         delete with;
00166                         with=withorig;
00167                 }
00168         }
00169 
00170         const float *const y_data = with->get_const_data();
00171         const float *const x_data = image->get_const_data();
00172         double result = 0.;
00173         float n = 0.0f;
00174         if(image->is_complex() && with->is_complex()) {
00175         // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
00176                 int nx  = with->get_xsize();
00177                 int ny  = with->get_ysize();
00178                 int nz  = with->get_zsize();
00179                 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
00180                 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
00181 
00182                 int ixb = 2*((nx+1)%2);
00183                 int iyb = ny%2;
00184                 //
00185                 if(nz == 1) {
00186                 //  it looks like it could work in 3D, but it is not, really.
00187                 for ( int iz = 0; iz <= nz-1; iz++) {
00188                         double part = 0.;
00189                         for ( int iy = 0; iy <= ny-1; iy++) {
00190                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ix++) {
00191                                                 size_t ii = ix + (iy  + iz * ny)* lsd2;
00192                                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00193                                 }
00194                         }
00195                         for ( int iy = 1; iy <= ny/2-1 + iyb; iy++) {
00196                                 size_t ii = (iy  + iz * ny)* lsd2;
00197                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00198                                 part += (x_data[ii+1] - y_data[ii+1])*double(x_data[ii+1] - y_data[ii+1]);
00199                         }
00200                         if(nx%2 == 0) {
00201                                 for ( int iy = 1; iy <= ny/2-1 + iyb; iy++) {
00202                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00203                                         part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00204                                         part += (x_data[ii+1] - y_data[ii+1])*double(x_data[ii+1] - y_data[ii+1]);
00205                                 }
00206 
00207                         }
00208                         part *= 2;
00209                         part += (x_data[0] - y_data[0])*double(x_data[0] - y_data[0]);
00210                         if(ny%2 == 0) {
00211                                 int ii = (ny/2  + iz * ny)* lsd2;
00212                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00213                         }
00214                         if(nx%2 == 0) {
00215                                 int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00216                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00217                                 if(ny%2 == 0) {
00218                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00219                                         part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00220                                 }
00221                         }
00222                         result += part;
00223                 }
00224                 n = (float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz;
00225 
00226                 }else{ //This 3D code is incorrect, but it is the best I can do now 01/09/06 PAP
00227                 int ky, kz;
00228                 int ny2 = ny/2; int nz2 = nz/2;
00229                 for ( int iz = 0; iz <= nz-1; iz++) {
00230                         if(iz>nz2) kz=iz-nz; else kz=iz;
00231                         for ( int iy = 0; iy <= ny-1; iy++) {
00232                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00233                                 for ( int ix = 0; ix <= lsd2-1; ix++) {
00234                                 // Skip Friedel related values
00235                                 if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00236                                                 size_t ii = ix + (iy  + iz * ny)* lsd2;
00237                                                 result += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00238                                         }
00239                                 }
00240                         }
00241                 }
00242                 n = ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz)/2.0f;
00243                 }
00244         } else {                // real space
00245                 size_t totsize = image->get_xsize()*image->get_ysize()*image->get_zsize();
00246                 if (params.has_key("mask")) {
00247                   EMData* mask;
00248                   mask = params["mask"];
00249                   const float *const dm = mask->get_const_data();
00250                   for (size_t i = 0; i < totsize; i++) {
00251                            if (dm[i] > 0.5) {
00252                                 double temp = x_data[i]- y_data[i];
00253                                 result += temp*temp;
00254                                 n++;
00255                            }
00256                   }
00257                 } 
00258                 else if (zeromask) {
00259                         n=0;
00260                         for (size_t i = 0; i < totsize; i++) {
00261                                 if (x_data[i]==0 || y_data[i]==0) continue;
00262                                 double temp = x_data[i]- y_data[i];
00263                                 result += temp*temp;
00264                                 n++;
00265                         }
00266                         
00267                 }
00268                 else {
00269                   for (size_t i = 0; i < totsize; i++) {
00270                                 double temp = x_data[i]- y_data[i];
00271                                 result += temp*temp;
00272                    }
00273                    n = (float)totsize;
00274                 }
00275         }
00276         result/=n;
00277 
00278         EXITFUNC;
00279         if (with->has_attr("deleteme")) delete with;
00280         float ret = (float)result;
00281         if (!Util::goodf(&ret)) return FLT_MAX;
00282         return ret;
00283 }
00284 
00285 
00286 // Even though this uses doubles, it might be wise to recode it row-wise
00287 // to avoid numerical errors on large images
00288 float DotCmp::cmp(EMData* image, EMData* with) const
00289 {
00290         ENTERFUNC;
00291         validate_input_args(image, with);
00292 
00293         const float *const x_data = image->get_const_data();
00294         const float *const y_data = with->get_const_data();
00295 
00296         int normalize = params.set_default("normalize", 0);
00297         float negative = (float)params.set_default("negative", 1);
00298 
00299         if (negative) negative=-1.0; else negative=1.0;
00300         double result = 0.;
00301         long n = 0;
00302         if(image->is_complex() && with->is_complex()) {
00303         // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
00304                 int nx  = with->get_xsize();
00305                 int ny  = with->get_ysize();
00306                 int nz  = with->get_zsize();
00307                 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
00308                 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
00309 
00310                 int ixb = 2*((nx+1)%2);
00311                 int iyb = ny%2;
00312                 //
00313                 if(nz == 1) {
00314                 //  it looks like it could work in 3D, but does not
00315                 for ( int iz = 0; iz <= nz-1; ++iz) {
00316                         double part = 0.;
00317                         for ( int iy = 0; iy <= ny-1; ++iy) {
00318                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00319                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00320                                         part += x_data[ii] * double(y_data[ii]);
00321                                 }
00322                         }
00323                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00324                                 size_t ii = (iy  + iz * ny)* lsd2;
00325                                 part += x_data[ii] * double(y_data[ii]);
00326                                 part += x_data[ii+1] * double(y_data[ii+1]);
00327                         }
00328                         if(nx%2 == 0) {
00329                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00330                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00331                                         part += x_data[ii] * double(y_data[ii]);
00332                                         part += x_data[ii+1] * double(y_data[ii+1]);
00333                                 }
00334 
00335                         }
00336                         part *= 2;
00337                         part += x_data[0] * double(y_data[0]);
00338                         if(ny%2 == 0) {
00339                                 size_t ii = (ny/2  + iz * ny)* lsd2;
00340                                 part += x_data[ii] * double(y_data[ii]);
00341                         }
00342                         if(nx%2 == 0) {
00343                                 size_t ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00344                                 part += x_data[ii] * double(y_data[ii]);
00345                                 if(ny%2 == 0) {
00346                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00347                                         part += x_data[ii] * double(y_data[ii]);
00348                                 }
00349                         }
00350                         result += part;
00351                 }
00352                 if( normalize ) {
00353                 //  it looks like it could work in 3D, but does not
00354                 double square_sum1 = 0., square_sum2 = 0.;
00355                 for ( int iz = 0; iz <= nz-1; ++iz) {
00356                         for ( int iy = 0; iy <= ny-1; ++iy) {
00357                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00358                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00359                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00360                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00361                                 }
00362                         }
00363                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00364                                 size_t ii = (iy  + iz * ny)* lsd2;
00365                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00366                                 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00367                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00368                                 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00369                         }
00370                         if(nx%2 == 0) {
00371                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00372                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00373                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00374                                         square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00375                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00376                                         square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00377                                 }
00378 
00379                         }
00380                         square_sum1 *= 2;
00381                         square_sum1 += x_data[0] * double(x_data[0]);
00382                         square_sum2 *= 2;
00383                         square_sum2 += y_data[0] * double(y_data[0]);
00384                         if(ny%2 == 0) {
00385                                 int ii = (ny/2  + iz * ny)* lsd2;
00386                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00387                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00388                         }
00389                         if(nx%2 == 0) {
00390                                 int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00391                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00392                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00393                                 if(ny%2 == 0) {
00394                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00395                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00396                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00397                                 }
00398                         }
00399                 }
00400                 result /= sqrt(square_sum1*square_sum2);
00401                 } else  result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
00402 
00403                 } else { //This 3D code is incorrect, but it is the best I can do now 01/09/06 PAP
00404                 int ky, kz;
00405                 int ny2 = ny/2; int nz2 = nz/2;
00406                 for ( int iz = 0; iz <= nz-1; ++iz) {
00407                         if(iz>nz2) kz=iz-nz; else kz=iz;
00408                         for ( int iy = 0; iy <= ny-1; ++iy) {
00409                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00410                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00411                                         // Skip Friedel related values
00412                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00413                                                 size_t ii = ix + (iy  + iz * ny)* lsd2;
00414                                                 result += x_data[ii] * double(y_data[ii]);
00415                                         }
00416                                 }
00417                         }
00418                 }
00419                 if( normalize ) {
00420                 //  still incorrect
00421                 double square_sum1 = 0., square_sum2 = 0.;
00422                 int ky, kz;
00423                 int ny2 = ny/2; int nz2 = nz/2;
00424                 for ( int iz = 0; iz <= nz-1; ++iz) {
00425                         if(iz>nz2) kz=iz-nz; else kz=iz;
00426                         for ( int iy = 0; iy <= ny-1; ++iy) {
00427                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00428                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00429                                         // Skip Friedel related values
00430                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00431                                                 size_t ii = ix + (iy  + iz * ny)* lsd2;
00432                                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00433                                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00434                                         }
00435                                 }
00436                         }
00437                 }
00438                 result /= sqrt(square_sum1*square_sum2);
00439                 } else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz/2);
00440                 }
00441         } else {
00442                 size_t totsize = image->get_xsize() * image->get_ysize() * image->get_zsize();
00443 
00444                 double square_sum1 = 0., square_sum2 = 0.;
00445 
00446                 if (params.has_key("mask")) {
00447                         EMData* mask;
00448                         mask = params["mask"];
00449                         const float *const dm = mask->get_const_data();
00450                         if (normalize) {
00451                                 for (size_t i = 0; i < totsize; i++) {
00452                                         if (dm[i] > 0.5) {
00453                                                 square_sum1 += x_data[i]*double(x_data[i]);
00454                                                 square_sum2 += y_data[i]*double(y_data[i]);
00455                                                 result += x_data[i]*double(y_data[i]);
00456                                         }
00457                                 }
00458                         } else {
00459                                 for (size_t i = 0; i < totsize; i++) {
00460                                         if (dm[i] > 0.5) {
00461                                                 result += x_data[i]*double(y_data[i]);
00462                                                 n++;
00463                                         }
00464                                 }
00465                         }
00466                 } else {
00467                         // this little bit of manual loop unrolling makes the dot product as fast as sqeuclidean with -O2
00468                         for (size_t i=0; i<totsize; i++) result+=x_data[i]*y_data[i];
00469 
00470                         if (normalize) {
00471                                 square_sum1 = image->get_attr("square_sum");
00472                                 square_sum2 = with->get_attr("square_sum");
00473                         } else n = totsize;
00474                 }
00475                 if (normalize) result /= (sqrt(square_sum1*square_sum2)); else result /= n;
00476         }
00477 
00478 
00479         EXITFUNC;
00480         return (float) (negative*result);
00481 }
00482 
00483 
00484 float TomoDotCmp::cmp(EMData * image, EMData *with) const
00485 {
00486         ENTERFUNC;
00487         float threshold = params.set_default("threshold",0.f);
00488         if (threshold < 0.0f) throw InvalidParameterException("The threshold parameter must be greater than or equal to zero");
00489 
00490         if ( threshold > 0) {
00491                 EMData* ccf = params.set_default("ccf",(EMData*) NULL);
00492                 bool ccf_ownership = false;
00493                 if (!ccf) {
00494                         ccf = image->calc_ccf(with);
00495                         ccf_ownership = true;
00496                 }
00497                 bool norm = params.set_default("norm",false);
00498                 if (norm) ccf->process_inplace("normalize");
00499                 int tx = params.set_default("tx",0); int ty = params.set_default("ty",0); int tz = params.set_default("tz",0);
00500                 float best_score = ccf->get_value_at_wrap(tx,ty,tz)/static_cast<float>(image->get_size());
00501                 EMData* ccf_fft = ccf->do_fft();// so cuda works, or else we could do an fft_inplace - though honestly doing an fft inplace is less efficient anyhow
00502                 if (ccf_ownership) delete ccf; ccf = 0;
00503                 ccf_fft->process_inplace("threshold.binary.fourier",Dict("value",threshold));
00504                 float map_sum =  ccf_fft->get_attr("mean");
00505                 if (map_sum == 0.0f) throw UnexpectedBehaviorException("The number of voxels in the Fourier image with an amplitude above your threshold is zero. Please adjust your parameters");
00506                 best_score /= map_sum;
00507                 delete ccf_fft; ccf_fft = 0;
00508                 return -best_score;
00509         } else {
00510                 return -image->dot(with);
00511         }
00512 
00513 
00514 }
00515 
00516 // Even though this uses doubles, it might be wise to recode it row-wise
00517 // to avoid numerical errors on large images
00518 float QuadMinDotCmp::cmp(EMData * image, EMData *with) const
00519 {
00520         ENTERFUNC;
00521         validate_input_args(image, with);
00522 
00523         if (image->get_zsize()!=1) throw InvalidValueException(0, "QuadMinDotCmp supports 2D only");
00524 
00525         int nx=image->get_xsize();
00526         int ny=image->get_ysize();
00527 
00528         int normalize = params.set_default("normalize", 0);
00529         float negative = (float)params.set_default("negative", 1);
00530 
00531         if (negative) negative=-1.0; else negative=1.0;
00532 
00533         double result[4] = { 0,0,0,0 }, sq1[4] = { 0,0,0,0 }, sq2[4] = { 0,0,0,0 } ;
00534 
00535         vector<int> image_saved_offsets = image->get_array_offsets();
00536         vector<int> with_saved_offsets = with->get_array_offsets();
00537         image->set_array_offsets(-nx/2,-ny/2);
00538         with->set_array_offsets(-nx/2,-ny/2);
00539         int i,x,y;
00540         for (y=-ny/2; y<ny/2; y++) {
00541                 for (x=-nx/2; x<nx/2; x++) {
00542                         int quad=(x<0?0:1) + (y<0?0:2);
00543                         result[quad]+=(*image)(x,y)*(*with)(x,y);
00544                         if (normalize) {
00545                                 sq1[quad]+=(*image)(x,y)*(*image)(x,y);
00546                                 sq2[quad]+=(*with)(x,y)*(*with)(x,y);
00547                         }
00548                 }
00549         }
00550         image->set_array_offsets(image_saved_offsets);
00551         with->set_array_offsets(with_saved_offsets);
00552 
00553         if (normalize) {
00554                 for (i=0; i<4; i++) result[i]/=sqrt(sq1[i]*sq2[i]);
00555         } else {
00556                 for (i=0; i<4; i++) result[i]/=nx*ny/4;
00557         }
00558 
00559         float worst=static_cast<float>(result[0]);
00560         for (i=1; i<4; i++) if (static_cast<float>(result[i])<worst) worst=static_cast<float>(result[i]);
00561 
00562         EXITFUNC;
00563         return (float) (negative*worst);
00564 }
00565 
00566 float OptVarianceCmp::cmp(EMData * image, EMData *with) const
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 }
00761 
00762 float PhaseCmp::cmp(EMData * image, EMData *with) const
00763 {
00764         ENTERFUNC;
00765 
00766         int snrweight = params.set_default("snrweight", 0);
00767         int snrfn = params.set_default("snrfn",0);
00768         int ampweight = params.set_default("ampweight",0);
00769         int zeromask = params.set_default("zeromask",0);
00770         float minres = params.set_default("minres",500.0f);
00771         float maxres = params.set_default("maxres",10.0f);
00772 
00773         if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator");
00774 
00775 #ifdef EMAN2_USING_CUDA
00776         if (image->gpu_operation_preferred()) {
00777 //              cout << "Cuda cmp" << endl;
00778                 EXITFUNC;
00779                 return cuda_cmp(image,with);
00780         }
00781 #endif
00782 
00783         EMData *image_fft = NULL;
00784         EMData *with_fft = NULL;
00785 
00786         int ny = image->get_ysize();
00787 //      int np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2;
00788         int np = 0;
00789         vector<float> snr;
00790 
00791         // weighting based on SNR estimate from CTF
00792         if (snrweight) {
00793                 Ctf *ctf = NULL;
00794                 if (!image->has_attr("ctf")) {
00795                         if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
00796                         ctf=with->get_attr("ctf");
00797                 }
00798                 else ctf=image->get_attr("ctf");
00799 
00800                 float ds=1.0f/(ctf->apix*ny);
00801                 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR); // note that this returns ny/2 values
00802                 if(ctf) {delete ctf; ctf=0;}
00803                 np=snr.size();
00804         }
00805         // weighting based on empirical SNR function (not really good)
00806         else if (snrfn==1) {
00807                 np=ny/2;
00808                 snr.resize(np);
00809 
00810                 for (int i = 0; i < np; i++) {
00811                         float x2 = 10.0f*i/np;
00812                         snr[i] = x2 * exp(-x2);
00813                 }
00814         }
00815         // no weighting
00816         else {
00817                 np=ny/2;
00818                 snr.resize(np);
00819 
00820                 for (int i = 0; i < np; i++) snr[i]=1.0;
00821         }
00822 
00823         // Min/max modifications to weighting
00824         float pmin,pmax;
00825         if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;                //cutoff in pixels, assume square
00826         else pmin=0;
00827         if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
00828         else pmax=0;
00829 
00830 //      printf("%f\t%f\n",pmin,pmax);
00831 
00832         // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation
00833         for (int i = 0; i < np; i++) {
00834                 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f;
00835                 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f;
00836 //              printf("%d\t%f\n",i,snr[i]);
00837         }
00838 
00839         if (zeromask) {
00840                 image_fft=image->copy();
00841                 with_fft=with->copy();
00842                 
00843                 if (image_fft->is_complex()) image_fft->do_ift_inplace();
00844                 if (with_fft->is_complex()) with_fft->do_ift_inplace();
00845                 
00846                 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize();
00847                 float *d1=image_fft->get_data();
00848                 float *d2=with_fft->get_data();
00849                 
00850                 for (int i=0; i<sz; i++) {
00851                         if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
00852                 }
00853                 
00854                 image_fft->update();
00855                 with_fft->update();
00856                 image_fft->do_fft_inplace();
00857                 with_fft->do_fft_inplace();
00858                 image_fft->set_attr("free_me",1); 
00859                 with_fft->set_attr("free_me",1); 
00860         }
00861         else {
00862                 if (image->is_complex()) image_fft=image;
00863                 else {
00864                         image_fft=image->do_fft();
00865                         image_fft->set_attr("free_me",1);
00866                 }
00867                 
00868                 if (with->is_complex()) with_fft=with;
00869                 else {
00870                         with_fft=with->do_fft();
00871                         with_fft->set_attr("free_me",1);
00872                 }
00873         }
00874 //      image_fft->ri2ap();
00875 //      with_fft->ri2ap();
00876 
00877         const float *const image_fft_data = image_fft->get_const_data();
00878         const float *const with_fft_data = with_fft->get_const_data();
00879         double sum = 0;
00880         double norm = FLT_MIN;
00881         size_t i = 0;
00882         int nx=image_fft->get_xsize();
00883             ny=image_fft->get_ysize();
00884         int nz=image_fft->get_zsize();
00885         int nx2=image_fft->get_xsize()/2;
00886         int ny2=image_fft->get_ysize()/2;
00887         int nz2=image_fft->get_zsize()/2;
00888 
00889         // This can never happen any more...
00890         if (np==0) {
00891                 for (int z = 0; z < nz; z++){
00892                         for (int y = 0; y < ny; y++) {
00893                                 for (int x = 0; x < nx2; x ++) {
00894                                         float a;
00895                                         if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
00896                                         else a=1.0f;
00897                                         sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
00898                                         norm += a;
00899                                         i += 2;
00900                                 }
00901                         }
00902                 }
00903                 
00904         }
00905         else if (nz==1) {
00906                 for (int y = 0; y < ny; y++) {
00907                         for (int x = 0; x < nx2; x ++) {
00908                                 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y);
00909                                 if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
00910                                 
00911                                 float a;
00912                                 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
00913                                 else a=1.0f;
00914                                 a*=snr[r];
00915                                 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
00916                                 norm += a;
00917                                 i += 2;
00918                         }
00919                 }
00920         }
00921         else {
00922                 for (int z = 0; z < nz; z++){
00923                         for (int y = 0; y < ny; y++) {
00924                                 for (int x = 0; x < nx2; x ++) {
00925                                         int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z);
00926                                         if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
00927                                         
00928                                         float a;
00929                                         if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
00930                                         else a=1.0f;
00931                                         a*=snr[r];
00932                                         sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
00933                                         norm += a;
00934                                         i += 2;
00935                                 } 
00936                         }
00937                 }
00938                 
00939         }
00940 
00941         EXITFUNC;
00942 
00943         if( image_fft->has_attr("free_me") )
00944         {
00945                 delete image_fft;
00946                 image_fft = 0;
00947         }
00948         if( with_fft->has_attr("free_me") )
00949         {
00950                 delete with_fft;
00951                 with_fft = 0;
00952         }
00953 #if 0
00954         return (1.0f - sum / norm);
00955 #endif
00956         return (float)(sum / norm);
00957 }
00958 
00959 #ifdef EMAN2_USING_CUDA
00960 #include "cuda/cuda_cmp.h"
00961 float PhaseCmp::cuda_cmp(EMData * image, EMData *with) const
00962 {
00963         ENTERFUNC;
00964         validate_input_args(image, with);
00965 
00966         typedef vector<EMData*> EMDatas;
00967         static EMDatas hist_pyramid;
00968         static EMDatas norm_pyramid;
00969         static EMData weighting;
00970         static int image_size = 0;
00971 
00972         int size;
00973         EMData::CudaDataLock imagelock(image);
00974         EMData::CudaDataLock withlock(with);
00975 
00976         if (image->is_complex()) {
00977                 size = image->get_xsize();
00978         } else {
00979                 int nx = image->get_xsize()+2;
00980                 nx -= nx%2;
00981                 size = nx*image->get_ysize()*image->get_zsize();
00982         }
00983         if (size != image_size) {
00984                 for(unsigned int i =0; i < hist_pyramid.size(); ++i) {
00985                         delete hist_pyramid[i];
00986                         delete norm_pyramid[i];
00987                 }
00988                 hist_pyramid.clear();
00989                 norm_pyramid.clear();
00990                 int s = size;
00991                 if (s < 1) throw UnexpectedBehaviorException("The image is 0 size");
00992                 int p2 = 1;
00993                 while ( s != 1 ) {
00994                         s /= 2;
00995                         p2 *= 2;
00996                 }
00997                 if ( p2 != size ) {
00998                         p2 *= 2;
00999                         s = p2;
01000                 }
01001                 if (s != 1) s /= 2;
01002                 while (true) {
01003                         EMData* h = new EMData();
01004                         h->set_size_cuda(s); h->to_value(0.0);
01005                         hist_pyramid.push_back(h);
01006                         EMData* n = new EMData();
01007                         n->set_size_cuda(s); n->to_value(0.0);
01008                         norm_pyramid.push_back(n);
01009                         if ( s == 1) break;
01010                         s /= 2;
01011                 }
01012                 int nx = image->get_xsize()+2;
01013                 nx -= nx%2; // for Fourier stuff
01014                 int ny = image->get_ysize();
01015                 int nz = image->get_zsize();
01016                 weighting.set_size_cuda(nx,ny,nz);
01017                 // Size of weighting need only be half this, but does that translate into faster code?
01018                 weighting.set_size_cuda(nx/2,ny,nz);
01019                 float np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2;
01020                 EMDataForCuda tmp = weighting.get_data_struct_for_cuda();
01021                 calc_phase_weights_cuda(&tmp,np);
01022                 //weighting.write_image("phase_wieghts.hdf");
01023                 image_size = size;
01024         }
01025 
01026         EMDataForCuda hist[hist_pyramid.size()];
01027         EMDataForCuda norm[hist_pyramid.size()];
01028 
01029         EMDataForCuda wt = weighting.get_data_struct_for_cuda();
01030         EMData::CudaDataLock lock1(&weighting);
01031         for(unsigned int i = 0; i < hist_pyramid.size(); ++i ) {
01032                 hist[i] = hist_pyramid[i]->get_data_struct_for_cuda();
01033                 hist_pyramid[i]->cuda_lock();
01034                 norm[i] = norm_pyramid[i]->get_data_struct_for_cuda();
01035                 norm_pyramid[i]->cuda_lock();
01036         }
01037 
01038         EMData *image_fft = image->do_fft_cuda();
01039         EMDataForCuda left = image_fft->get_data_struct_for_cuda();
01040         EMData::CudaDataLock lock2(image_fft);
01041         EMData *with_fft = with->do_fft_cuda();
01042         EMDataForCuda right = with_fft->get_data_struct_for_cuda();
01043         EMData::CudaDataLock lock3(image_fft);
01044 
01045         mean_phase_error_cuda(&left,&right,&wt,hist,norm,hist_pyramid.size());
01046         float result;
01047         float* gpu_result = hist_pyramid[hist_pyramid.size()-1]->get_cuda_data();
01048         cudaError_t error = cudaMemcpy(&result,gpu_result,sizeof(float),cudaMemcpyDeviceToHost);
01049         if ( error != cudaSuccess) throw UnexpectedBehaviorException( "CudaMemcpy (host to device) in the phase comparator failed:" + string(cudaGetErrorString(error)));
01050 
01051         delete image_fft; image_fft=0;
01052         delete with_fft; with_fft=0;
01053 
01054         for(unsigned int i = 0; i < hist_pyramid.size(); ++i ) {
01055 //              hist_pyramid[i]->write_image("hist.hdf",-1); // debug
01056 //              norm_pyramid[i]->write_image("norm.hdf",-1); // debug
01057                 hist_pyramid[i]->cuda_unlock();
01058                 norm_pyramid[i]->cuda_unlock();
01059         }
01060 
01061         EXITFUNC;
01062         return result;
01063 
01064 }
01065 
01066 #endif // EMAN2_USING_CUDA
01067 
01068 
01069 float FRCCmp::cmp(EMData * image, EMData * with) const
01070 {
01071         ENTERFUNC;
01072         validate_input_args(image, with);
01073 
01074         int snrweight = params.set_default("snrweight", 0);
01075         int ampweight = params.set_default("ampweight", 0);
01076         int sweight = params.set_default("sweight", 1);
01077         int nweight = params.set_default("nweight", 0);
01078         int zeromask = params.set_default("zeromask",0);
01079         float minres = params.set_default("minres",500.0f);
01080         float maxres = params.set_default("maxres",10.0f);
01081 
01082         if (zeromask) {
01083                 image=image->copy();
01084                 with=with->copy();
01085                 
01086                 int sz=image->get_xsize()*image->get_ysize()*image->get_zsize();
01087                 float *d1=image->get_data();
01088                 float *d2=with->get_data();
01089                 
01090                 for (int i=0; i<sz; i++) {
01091                         if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
01092                 }
01093                 
01094                 image->update();
01095                 with->update();
01096                 image->do_fft_inplace();
01097                 with->do_fft_inplace();
01098                 image->set_attr("free_me",1); 
01099                 with->set_attr("free_me",1); 
01100         }
01101 
01102 
01103         if (!image->is_complex()) {
01104                 image=image->do_fft(); 
01105                 image->set_attr("free_me",1); 
01106         }
01107         if (!with->is_complex()) { 
01108                 with=with->do_fft(); 
01109                 with->set_attr("free_me",1); 
01110         }
01111 
01112         static vector < float >default_snr;
01113 
01114 //      if (image->get_zsize() > 1) {
01115 //              throw ImageDimensionException("2D only");
01116 //      }
01117 
01118 //      int nx = image->get_xsize();
01119         int ny = image->get_ysize();
01120         int ny2=ny/2+1;
01121 
01122         vector < float >fsc;
01123 
01124                 
01125 
01126         fsc = image->calc_fourier_shell_correlation(with,1);
01127 
01128         // The fast hypot here was supposed to speed things up. Little effect
01129 //      if (image->get_zsize()>1) fsc = image->calc_fourier_shell_correlation(with,1);
01130 //      else {
01131 //              double *sxy = (double *)malloc(ny2*sizeof(double)*4);
01132 //              double *sxx = sxy+ny2;
01133 //              double *syy = sxy+2*ny2;
01134 //              double *norm= sxy+3*ny2;
01135 //
01136 //              float *df1=image->get_data();
01137 //              float *df2=with->get_data();
01138 //              int nx2=image->get_xsize();
01139 //
01140 //              for (int y=-ny/2; y<ny/2; y++) {
01141 //                      for (int x=0; x<nx2/2; x++) {
01142 //                              if (x==0 && y<0) continue;      // skip Friedel pair
01143 //                              short r=Util::hypot_fast_int(x,y);
01144 //                              if (r>ny2-1) continue;
01145 //                              int l=x*2+(y<0?ny+y:y)*nx2;
01146 //                              sxy[r]+=df1[l]*df2[l]+df1[l+1]*df2[l+1];
01147 //                              sxx[r]+=df1[l]*df1[l];
01148 //                              syy[r]+=df2[l]*df2[l];
01149 //                              norm[r]+=1.0;
01150 //                      }
01151 //              }
01152 //              fsc.resize(ny2*3);
01153 //              for (int r=0; r<ny2; r++) {
01154 //                      fsc[r]=r*0.5/ny2;
01155 //                      fsc[ny2+r]=sxy[r]/(sqrt(sxx[r])*sqrt(syy[r]));
01156 //                      fsc[ny2*2+r]=norm[r];
01157 //              }
01158 //              free(sxy);
01159 //      }
01160 
01161         vector<float> snr;
01162         if (snrweight) {
01163                 Ctf *ctf = NULL;
01164                 if (!image->has_attr("ctf")) {
01165                         if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
01166                         ctf=with->get_attr("ctf");
01167                 }
01168                 else ctf=image->get_attr("ctf");
01169 
01170                 float ds=1.0f/(ctf->apix*ny);
01171                 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR);
01172                 if(ctf) {delete ctf; ctf=0;}
01173         }
01174 
01175         vector<float> amp;
01176         if (ampweight) amp=image->calc_radial_dist(ny/2,0,1,0);
01177 
01178         // Min/max modifications to weighting
01179         float pmin,pmax;
01180         if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;                //cutoff in pixels, assume square
01181         else pmin=0;
01182         if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
01183         else pmax=0;
01184 
01185         double sum=0.0, norm=0.0;
01186 
01187         for (int i=0; i<ny/2; i++) {
01188                 double weight=1.0;
01189                 if (sweight) weight*=fsc[(ny2)*2+i];
01190                 if (ampweight) weight*=amp[i];
01191                 if (snrweight) weight*=snr[i];
01192                 if (pmin>0) weight*=(tanh(5.0*(i-pmin)/pmin)+1.0)/2.0;
01193                 if (pmax>0) weight*=(1.0-tanh(i-pmax))/2.0;
01194                 
01195                 sum+=weight*fsc[ny2+i];
01196                 norm+=weight;
01197 //              printf("%d\t%f\t%f\n",i,weight,fsc[ny/2+1+i]);
01198         }
01199 
01200         // This performs a weighting that tries to normalize FRC by correcting from the number of particles represented by the average
01201         sum/=norm;
01202         if (nweight && with->get_attr_default("ptcl_repr",0) && sum>=0 && sum<1.0) {
01203                 sum=sum/(1.0-sum);                                                      // convert to SNR
01204                 sum/=(float)with->get_attr_default("ptcl_repr",0);      // divide by ptcl represented
01205                 sum=sum/(1.0+sum);                                                      // convert back to correlation
01206         }
01207 
01208         if (image->has_attr("free_me")) delete image;
01209         if (with->has_attr("free_me")) delete with;
01210 
01211         EXITFUNC;
01212 
01213 
01214         //.Note the negative! This is because EMAN2 follows the convention that
01215         // smaller return values from comparitors indicate higher similarity -
01216         // this enables comparitors to be used in a generic fashion.
01217         return (float)-sum;
01218 }
01219 
01220 void EMAN::dump_cmps()
01221 {
01222         dump_factory < Cmp > ();
01223 }
01224 
01225 map<string, vector<string> > EMAN::dump_cmps_list()
01226 {
01227         return dump_factory_list < Cmp > ();
01228 }
01229 
01230 /* vim: set ts=4 noet: */

Generated on Tue May 25 17:13:31 2010 for EMAN2 by  doxygen 1.4.7