Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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 #ifdef EMAN2_USING_CUDA
00042 // Only CCC, DOT  and CCC.TOMO are cuda enabled
00043 #include "cuda/cuda_processor.h"
00044 #include "cuda/cuda_cmp.h"
00045 #endif // EMAN2_USING_CUDA
00046 
00047 using namespace EMAN;
00048 
00049 const string CccCmp::NAME = "ccc";
00050 const string LodCmp::NAME = "lod";
00051 const string SqEuclideanCmp::NAME = "sqeuclidean";
00052 const string DotCmp::NAME = "dot";
00053 const string TomoCccCmp::NAME = "ccc.tomo";
00054 const string QuadMinDotCmp::NAME = "quadmindot";
00055 const string OptVarianceCmp::NAME = "optvariance";
00056 const string PhaseCmp::NAME = "phase";
00057 const string FRCCmp::NAME = "frc";
00058 
00059 template <> Factory < Cmp >::Factory()
00060 {
00061         force_add<CccCmp>();
00062         force_add<LodCmp>();
00063         force_add<SqEuclideanCmp>();
00064         force_add<DotCmp>();
00065         force_add<TomoCccCmp>();
00066         force_add<QuadMinDotCmp>();
00067         force_add<OptVarianceCmp>();
00068         force_add<PhaseCmp>();
00069         force_add<FRCCmp>();
00070 //      force_add<XYZCmp>();
00071 }
00072 
00073 void Cmp::validate_input_args(const EMData * image, const EMData *with) const
00074 {
00075         
00076 #ifdef EMAN2_USING_CUDA  
00077         if (image->getcudarwdata() && with->getcudarwdata()) {
00078                 //no need for futher checking, which will induce an expensive copy from device to host
00079                 return;
00080         }
00081 #endif
00082         if (!image) {
00083                 throw NullPointerException("compared image");
00084         }
00085         if (!with) {
00086                 throw NullPointerException("compare-with image");
00087         }
00088 
00089         if (!EMUtil::is_same_size(image, with)) {
00090                 throw ImageFormatException( "images not same size");
00091         }
00092 
00093         float *d1 = image->get_data();
00094         if (!d1) {
00095                 throw NullPointerException("image contains no data");
00096         }
00097 
00098         float *d2 = with->get_data();
00099         if (!d2) {
00100                 throw NullPointerException("compare-with image data");
00101         }
00102 }
00103 
00104 //  It would be good to add code for complex images!  PAP
00105 float CccCmp::cmp(EMData * image, EMData *with) const
00106 {
00107         ENTERFUNC;
00108         if (image->is_complex() || with->is_complex())
00109                 throw ImageFormatException( "Complex images not supported by CMP::CccCmp");
00110         validate_input_args(image, with);
00111 
00112         const float *const d1 = image->get_const_data();
00113         const float *const d2 = with->get_const_data();
00114 
00115         float negative = (float)params.set_default("negative", 1);
00116         if (negative) negative=-1.0; else negative=1.0;
00117 
00118         double avg1 = 0.0, var1 = 0.0, avg2 = 0.0, var2 = 0.0, ccc = 0.0;
00119         long n = 0;
00120         size_t totsize = image->get_xsize()*image->get_ysize()*image->get_zsize();
00121 
00122         bool has_mask = false;
00123         EMData* mask = 0;
00124         if (params.has_key("mask")) {
00125                 mask = params["mask"];
00126                 if(mask!=0) {has_mask=true;}
00127         }
00128 #ifdef EMAN2_USING_CUDA
00129         if (image->getcudarwdata() && with->getcudarwdata()) {
00130                 //cout << "CUDA ccc cmp" << endl;
00131                 float* maskdata = 0;
00132                 if(has_mask && !mask->getcudarwdata()){
00133                         mask->copy_to_cuda();
00134                         maskdata = mask->getcudarwdata();
00135                 }
00136                 float ccc = ccc_cmp_cuda(image->getcudarwdata(), with->getcudarwdata(), maskdata, image->get_xsize(), image->get_ysize(), image->get_zsize());
00137                 ccc *= negative;
00138                 //cout << "CUDA CCC is: " << ccc << endl;
00139                 return ccc;
00140         }
00141 #endif
00142         if (has_mask) {
00143                 const float *const dm = mask->get_const_data();
00144                 for (size_t i = 0; i < totsize; ++i) {
00145                         if (dm[i] > 0.5) {
00146                                 avg1 += double(d1[i]);
00147                                 var1 += d1[i]*double(d1[i]);
00148                                 avg2 += double(d2[i]);
00149                                 var2 += d2[i]*double(d2[i]);
00150                                 ccc += d1[i]*double(d2[i]);
00151                                 n++;
00152                         }
00153                 }
00154         } else {
00155                 for (size_t i = 0; i < totsize; ++i) {
00156                         avg1 += double(d1[i]);
00157                         var1 += d1[i]*double(d1[i]);
00158                         avg2 += double(d2[i]);
00159                         var2 += d2[i]*double(d2[i]);
00160                         ccc += d1[i]*double(d2[i]);
00161                 }
00162                 n = totsize;
00163         }
00164 
00165         avg1 /= double(n);
00166         var1 = var1/double(n) - avg1*avg1;
00167         avg2 /= double(n);
00168         var2 = var2/double(n) - avg2*avg2;
00169         ccc = ccc/double(n) - avg1*avg2;
00170         ccc /= sqrt(var1*var2);
00171         ccc *= negative;
00172         return static_cast<float>(ccc);
00173         EXITFUNC;
00174 }
00175 
00176 
00177 // Added by JAK 11/12/10
00178 // L^1-norm difference of two maps, after normalization.
00179 float LodCmp::cmp(EMData * image, EMData *with) const
00180 {
00181         ENTERFUNC;
00182         if (image->is_complex() || with->is_complex())
00183                 throw ImageFormatException( "Complex images not yet supported by CMP::LodCmp");
00184         validate_input_args(image, with);
00185 
00186         const float *const d1 = image->get_const_data();
00187         const float *const d2 = with->get_const_data();
00188 
00189         float negative = (float)params.set_default("negative", 1);
00190         if (negative) negative=-1.0; else negative=1.0;
00191 
00192         double avg1 = 0.0, sig1 = 0.0, avg2 = 0.0, sig2 = 0.0, lod = 0.0;
00193         long n = 0;
00194         size_t totsize = image->get_xsize()*image->get_ysize()*image->get_zsize();
00195         size_t i;
00196         
00197         bool has_mask = false;
00198         EMData* mask = 0;
00199         if (params.has_key("mask")) {
00200                 mask = params["mask"];
00201                 if(mask!=0) {has_mask=true;}
00202         }
00203 
00204         int normalize = 0;
00205         if (params.has_key("normalize")) {
00206                 normalize = params["normalize"];
00207         }
00208 
00209         if (normalize) {
00210                 if (has_mask) {
00211                         const float *const dm = mask->get_const_data();
00212                         for (i = 0; i < totsize; ++i) {
00213                                 if (dm[i] > 0.5) {
00214                                         avg1 += double(d1[i]);
00215                                         avg2 += double(d2[i]);
00216                                         n++;
00217                                 }
00218                         }
00219                 } else {
00220                         for (i = 0; i < totsize; ++i) {
00221                                 avg1 += double(d1[i]);
00222                                 avg2 += double(d2[i]);
00223                         }
00224                         n = totsize;
00225                 }
00226 
00227                 avg1 /= double(n);
00228                 avg2 /= double(n);
00229 
00230                 if (has_mask) {
00231                         const float *const dm = mask->get_const_data();
00232                         for (i = 0; i < totsize; ++i) {
00233                                 if (dm[i] > 0.5) {
00234                                         sig1 += fabs(double(d1[i])-avg1);
00235                                         sig2 += fabs(double(d2[i])-avg2);
00236                                 }
00237                         }
00238                 } else {
00239                         for (i = 0; i < totsize; ++i) {
00240                                 sig1 += fabs(double(d1[i])-avg1);
00241                                 sig2 += fabs(double(d2[i])-avg2);
00242                         }
00243                 }
00244         } else {
00245                 avg1 = 0.0; avg2 = 0.0;
00246                 sig1 = 1.0; sig2 = 1.0;
00247         }
00248 
00249         if (has_mask) {
00250                 const float *const dm = mask->get_const_data();
00251                 for (i = 0; i < totsize; ++i) {
00252                         if (dm[i] > 0.5) {
00253                                 lod += fabs((double(d1[i])-avg1)/sig1 - (double(d2[i])-avg2)/sig2);
00254                         }
00255                 }
00256         } else {
00257                 for (i = 0; i < totsize; ++i) {
00258                         lod += fabs((double(d1[i])-avg1)/sig1 - (double(d2[i])-avg2)/sig2);
00259                 }
00260         }
00261         
00262         lod *= (-0.5);
00263         lod *= negative;
00264         return static_cast<float>(lod);
00265         EXITFUNC;
00266 }
00267 
00268 
00269 //float SqEuclideanCmp::cmp(EMData * image, EMData *withorig) const
00270 float SqEuclideanCmp::cmp(EMData *image,EMData * withorig ) const
00271 {
00272         ENTERFUNC;
00273         EMData *with = withorig;
00274         validate_input_args(image, with);
00275 
00276         int zeromask = params.set_default("zeromask",0);
00277         int normto = params.set_default("normto",0);
00278 
00279         if (normto) {
00280                 if (zeromask) with = withorig->process("normalize.toimage",Dict("to",image));
00281                 else with = withorig->process("normalize.toimage",Dict("to",image,"ignore_zero",0));
00282                 with->set_attr("deleteme",1);
00283                 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
00284                         delete with;
00285                         with=withorig;
00286                 }
00287         }
00288 
00289         const float *const y_data = with->get_const_data();
00290         const float *const x_data = image->get_const_data();
00291         double result = 0.;
00292         float n = 0.0f;
00293         if(image->is_complex() && with->is_complex()) {
00294         // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
00295                 int nx  = with->get_xsize();
00296                 int ny  = with->get_ysize();
00297                 int nz  = with->get_zsize();
00298                 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
00299                 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
00300 
00301                 int ixb = 2*((nx+1)%2);
00302                 int iyb = ny%2;
00303                 //
00304                 if(nz == 1) {
00305                 //  it looks like it could work in 3D, but it is not, really.
00306                 for ( int iz = 0; iz <= nz-1; iz++) {
00307                         double part = 0.;
00308                         for ( int iy = 0; iy <= ny-1; iy++) {
00309                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ix++) {
00310                                                 size_t ii = ix + (iy  + iz * ny)* lsd2;
00311                                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00312                                 }
00313                         }
00314                         for ( int iy = 1; iy <= ny/2-1 + iyb; iy++) {
00315                                 size_t ii = (iy  + iz * ny)* lsd2;
00316                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00317                                 part += (x_data[ii+1] - y_data[ii+1])*double(x_data[ii+1] - y_data[ii+1]);
00318                         }
00319                         if(nx%2 == 0) {
00320                                 for ( int iy = 1; iy <= ny/2-1 + iyb; iy++) {
00321                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00322                                         part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00323                                         part += (x_data[ii+1] - y_data[ii+1])*double(x_data[ii+1] - y_data[ii+1]);
00324                                 }
00325 
00326                         }
00327                         part *= 2;
00328                         part += (x_data[0] - y_data[0])*double(x_data[0] - y_data[0]);
00329                         if(ny%2 == 0) {
00330                                 int ii = (ny/2  + iz * ny)* lsd2;
00331                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00332                         }
00333                         if(nx%2 == 0) {
00334                                 int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00335                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00336                                 if(ny%2 == 0) {
00337                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00338                                         part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00339                                 }
00340                         }
00341                         result += part;
00342                 }
00343                 n = (float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz;
00344 
00345                 }else{ //This 3D code is incorrect, but it is the best I can do now 01/09/06 PAP
00346                 int ky, kz;
00347                 int ny2 = ny/2; int nz2 = nz/2;
00348                 for ( int iz = 0; iz <= nz-1; iz++) {
00349                         if(iz>nz2) kz=iz-nz; else kz=iz;
00350                         for ( int iy = 0; iy <= ny-1; iy++) {
00351                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00352                                 for ( int ix = 0; ix <= lsd2-1; ix++) {
00353                                 // Skip Friedel related values
00354                                 if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00355                                                 size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
00356                                                 result += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00357                                         }
00358                                 }
00359                         }
00360                 }
00361                 n = ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz)/2.0f;
00362                 }
00363         } else {                // real space
00364                 size_t totsize = (size_t)image->get_xsize()*image->get_ysize()*image->get_zsize();
00365                 if (params.has_key("mask")) {
00366                   EMData* mask;
00367                   mask = params["mask"];
00368                   const float *const dm = mask->get_const_data();
00369                   for (size_t i = 0; i < totsize; i++) {
00370                            if (dm[i] > 0.5) {
00371                                 double temp = x_data[i]- y_data[i];
00372                                 result += temp*temp;
00373                                 n++;
00374                            }
00375                   }
00376                 } 
00377                 else if (zeromask) {
00378                         n=0;
00379                         for (size_t i = 0; i < totsize; i++) {
00380                                 if (x_data[i]==0 || y_data[i]==0) continue;
00381                                 double temp = x_data[i]- y_data[i];
00382                                 result += temp*temp;
00383                                 n++;
00384                         }
00385                         
00386                 }
00387                 else {
00388                   for (size_t i = 0; i < totsize; i++) {
00389                                 double temp = x_data[i]- y_data[i];
00390                                 result += temp*temp;
00391                    }
00392                    n = (float)totsize;
00393                 }
00394         }
00395         result/=n;
00396 
00397         EXITFUNC;
00398         if (with->has_attr("deleteme")) delete with;
00399         float ret = (float)result;
00400         if (!Util::goodf(&ret)) return FLT_MAX;
00401         return ret;
00402 }
00403 
00404 
00405 // Even though this uses doubles, it might be wise to recode it row-wise
00406 // to avoid numerical errors on large images
00407 float DotCmp::cmp(EMData* image, EMData* with) const
00408 {
00409         ENTERFUNC;
00410         
00411         validate_input_args(image, with);
00412 
00413         int normalize = params.set_default("normalize", 0);
00414         float negative = (float)params.set_default("negative", 1);
00415         if (negative) negative=-1.0; else negative=1.0;
00416 #ifdef EMAN2_USING_CUDA // SO far only works for real images I put CUDA first to avoid running non CUDA overhead (calls to getdata are expensive!!!!)
00417         if(image->is_complex() && with->is_complex()) {
00418         } else {
00419                 if (image->getcudarwdata() && with->getcudarwdata()) {
00420                         //cout << "CUDA dot cmp" << endl;
00421                         float* maskdata = 0;
00422                         bool has_mask = false;
00423                         EMData* mask = 0;
00424                         if (params.has_key("mask")) {
00425                                 mask = params["mask"];
00426                                 if(mask!=0) {has_mask=true;}
00427                         }
00428                         if(has_mask && !mask->getcudarwdata()){
00429                                 mask->copy_to_cuda();
00430                                 maskdata = mask->getcudarwdata();
00431                         }
00432 
00433                         float result = dot_cmp_cuda(image->getcudarwdata(), with->getcudarwdata(), maskdata, image->get_xsize(), image->get_ysize(), image->get_zsize());
00434                         result *= negative;
00435 
00436                         return result;
00437                         
00438                 }
00439         }
00440 #endif
00441         const float *const x_data = image->get_const_data();
00442         const float *const y_data = with->get_const_data();
00443 
00444         double result = 0.;
00445         long n = 0;
00446         if(image->is_complex() && with->is_complex()) {
00447         // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
00448                 int nx  = with->get_xsize();
00449                 int ny  = with->get_ysize();
00450                 int nz  = with->get_zsize();
00451                 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
00452                 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
00453 
00454                 int ixb = 2*((nx+1)%2);
00455                 int iyb = ny%2;
00456                 //
00457                 if(nz == 1) {
00458                 //  it looks like it could work in 3D, but does not
00459                 for ( int iz = 0; iz <= nz-1; ++iz) {
00460                         double part = 0.;
00461                         for ( int iy = 0; iy <= ny-1; ++iy) {
00462                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00463                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00464                                         part += x_data[ii] * double(y_data[ii]);
00465                                 }
00466                         }
00467                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00468                                 size_t ii = (iy  + iz * ny)* lsd2;
00469                                 part += x_data[ii] * double(y_data[ii]);
00470                                 part += x_data[ii+1] * double(y_data[ii+1]);
00471                         }
00472                         if(nx%2 == 0) {
00473                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00474                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00475                                         part += x_data[ii] * double(y_data[ii]);
00476                                         part += x_data[ii+1] * double(y_data[ii+1]);
00477                                 }
00478 
00479                         }
00480                         part *= 2;
00481                         part += x_data[0] * double(y_data[0]);
00482                         if(ny%2 == 0) {
00483                                 size_t ii = (ny/2  + iz * ny)* lsd2;
00484                                 part += x_data[ii] * double(y_data[ii]);
00485                         }
00486                         if(nx%2 == 0) {
00487                                 size_t ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00488                                 part += x_data[ii] * double(y_data[ii]);
00489                                 if(ny%2 == 0) {
00490                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00491                                         part += x_data[ii] * double(y_data[ii]);
00492                                 }
00493                         }
00494                         result += part;
00495                 }
00496                 if( normalize ) {
00497                 //  it looks like it could work in 3D, but does not
00498                 double square_sum1 = 0., square_sum2 = 0.;
00499                 for ( int iz = 0; iz <= nz-1; ++iz) {
00500                         for ( int iy = 0; iy <= ny-1; ++iy) {
00501                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00502                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00503                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00504                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00505                                 }
00506                         }
00507                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00508                                 size_t ii = (iy  + iz * ny)* lsd2;
00509                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00510                                 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00511                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00512                                 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00513                         }
00514                         if(nx%2 == 0) {
00515                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00516                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00517                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00518                                         square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00519                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00520                                         square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00521                                 }
00522 
00523                         }
00524                         square_sum1 *= 2;
00525                         square_sum1 += x_data[0] * double(x_data[0]);
00526                         square_sum2 *= 2;
00527                         square_sum2 += y_data[0] * double(y_data[0]);
00528                         if(ny%2 == 0) {
00529                                 int ii = (ny/2  + iz * ny)* lsd2;
00530                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00531                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00532                         }
00533                         if(nx%2 == 0) {
00534                                 int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00535                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00536                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00537                                 if(ny%2 == 0) {
00538                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00539                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00540                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00541                                 }
00542                         }
00543                 }
00544                 result /= sqrt(square_sum1*square_sum2);
00545                 } else  result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
00546 
00547                 } else { //This 3D code is incorrect, but it is the best I can do now 01/09/06 PAP
00548                 int ky, kz;
00549                 int ny2 = ny/2; int nz2 = nz/2;
00550                 for ( int iz = 0; iz <= nz-1; ++iz) {
00551                         if(iz>nz2) kz=iz-nz; else kz=iz;
00552                         for ( int iy = 0; iy <= ny-1; ++iy) {
00553                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00554                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00555                                         // Skip Friedel related values
00556                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00557                                                 size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
00558                                                 result += x_data[ii] * double(y_data[ii]);
00559                                         }
00560                                 }
00561                         }
00562                 }
00563                 if( normalize ) {
00564                 //  still incorrect
00565                 double square_sum1 = 0., square_sum2 = 0.;
00566                 int ky, kz;
00567                 int ny2 = ny/2; int nz2 = nz/2;
00568                 for ( int iz = 0; iz <= nz-1; ++iz) {
00569                         if(iz>nz2) kz=iz-nz; else kz=iz;
00570                         for ( int iy = 0; iy <= ny-1; ++iy) {
00571                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00572                                 for ( int ix = 0; ix <= lsd2-1; ++ix) {
00573                                         // Skip Friedel related values
00574                                         if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00575                                                 size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
00576                                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00577                                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00578                                         }
00579                                 }
00580                         }
00581                 }
00582                 result /= sqrt(square_sum1*square_sum2);
00583                 } else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz/2);
00584                 }
00585         } else {
00586                 
00587                 size_t totsize = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00588 
00589                 double square_sum1 = 0., square_sum2 = 0.;
00590 
00591                 if (params.has_key("mask")) {
00592                         EMData* mask;
00593                         mask = params["mask"];
00594                         const float *const dm = mask->get_const_data();
00595                         if (normalize) {
00596                                 for (size_t i = 0; i < totsize; i++) {
00597                                         if (dm[i] > 0.5) {
00598                                                 square_sum1 += x_data[i]*double(x_data[i]);
00599                                                 square_sum2 += y_data[i]*double(y_data[i]);
00600                                                 result += x_data[i]*double(y_data[i]);
00601                                         }
00602                                 }
00603                         } else {
00604                                 for (size_t i = 0; i < totsize; i++) {
00605                                         if (dm[i] > 0.5) {
00606                                                 result += x_data[i]*double(y_data[i]);
00607                                                 n++;
00608                                         }
00609                                 }
00610                         }
00611                 } else {
00612                         // this little bit of manual loop unrolling makes the dot product as fast as sqeuclidean with -O2
00613                         for (size_t i=0; i<totsize; i++) result+=x_data[i]*y_data[i];
00614 
00615                         if (normalize) {
00616                                 square_sum1 = image->get_attr("square_sum");
00617                                 square_sum2 = with->get_attr("square_sum");
00618                         } else n = totsize;
00619                 }
00620                 if (normalize) result /= (sqrt(square_sum1*square_sum2)); else result /= n;
00621         }
00622 
00623 
00624         EXITFUNC;
00625         return (float) (negative*result);
00626 }
00627 
00628 // This implements the technique of Mike Schmid where by the cross correlation is normalized
00629 // in an effort to remove the effects of the missing wedge. Somewhat of a heuristic solution, but it seems
00630 // to work. Basiaclly it relies on the observation that 'good' matchs will conentrate the correlation
00631 // signal in the peak, wheras 'bad' correlations will distribute the signal.
00632 // John Flanagan 18-10-2010
00633 
00634 float TomoCccCmp::cmp(EMData * image, EMData *with) const
00635 {
00636         ENTERFUNC;
00637         EMData* ccf = params.set_default("ccf",(EMData*) NULL);
00638         bool ccf_ownership = false;
00639         bool norm = params.set_default("norm",true);
00640         float negative = (float)params.set_default("negative", 1);
00641         if (negative) negative=-1.0; else negative=1.0;
00642         int searchx, searchy, searchz;
00643         
00644 
00645         searchx = params.set_default("searchx",-1); 
00646         searchy = params.set_default("searchy",-1); 
00647         searchz = params.set_default("searchz",-1);
00648         
00649 #ifdef EMAN2_USING_CUDA 
00650         if(image->getcudarwdata() && with->getcudarwdata()){
00651                 if (!ccf) {
00652                         ccf = image->calc_ccf(with);
00653                         ccf_ownership = true;
00654                 }
00655                 //cout << "using CUDA" << endl;
00656                 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
00657                 float best_score = get_value_at_wrap_cuda(ccf->getcudarwdata(), 0, 0, 0, ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
00658                 if(norm) {
00659                         best_score = negative*(best_score - stats.x)/sqrt(stats.y);
00660                 } else {
00661                         best_score = negative*best_score;
00662                 }
00663                 
00664                 if (ccf_ownership) delete ccf; ccf = 0;
00665                 
00666                 return best_score;
00667                 
00668         }
00669 #endif
00670 
00671         if (!ccf) {
00672                 ccf = image->calc_ccf(with);
00673                 ccf_ownership = true;
00674         }
00675         if (norm) ccf->process_inplace("normalize");
00676         
00677         float best_score = ccf->get_value_at_wrap(0,0,0);
00678         if (ccf_ownership) delete ccf; ccf = 0;
00679         
00680         return negative*best_score;
00681 
00682 }
00683 
00684 // Even though this uses doubles, it might be wise to recode it row-wise
00685 // to avoid numerical errors on large images
00686 float QuadMinDotCmp::cmp(EMData * image, EMData *with) const
00687 {
00688         ENTERFUNC;
00689         validate_input_args(image, with);
00690 
00691         if (image->get_zsize()!=1) throw InvalidValueException(0, "QuadMinDotCmp supports 2D only");
00692 
00693         int nx=image->get_xsize();
00694         int ny=image->get_ysize();
00695 
00696         int normalize = params.set_default("normalize", 0);
00697         float negative = (float)params.set_default("negative", 1);
00698 
00699         if (negative) negative=-1.0; else negative=1.0;
00700 
00701         double result[4] = { 0,0,0,0 }, sq1[4] = { 0,0,0,0 }, sq2[4] = { 0,0,0,0 } ;
00702 
00703         vector<int> image_saved_offsets = image->get_array_offsets();
00704         vector<int> with_saved_offsets = with->get_array_offsets();
00705         image->set_array_offsets(-nx/2,-ny/2);
00706         with->set_array_offsets(-nx/2,-ny/2);
00707         int i,x,y;
00708         for (y=-ny/2; y<ny/2; y++) {
00709                 for (x=-nx/2; x<nx/2; x++) {
00710                         int quad=(x<0?0:1) + (y<0?0:2);
00711                         result[quad]+=(*image)(x,y)*(*with)(x,y);
00712                         if (normalize) {
00713                                 sq1[quad]+=(*image)(x,y)*(*image)(x,y);
00714                                 sq2[quad]+=(*with)(x,y)*(*with)(x,y);
00715                         }
00716                 }
00717         }
00718         image->set_array_offsets(image_saved_offsets);
00719         with->set_array_offsets(with_saved_offsets);
00720 
00721         if (normalize) {
00722                 for (i=0; i<4; i++) result[i]/=sqrt(sq1[i]*sq2[i]);
00723         } else {
00724                 for (i=0; i<4; i++) result[i]/=nx*ny/4;
00725         }
00726 
00727         float worst=static_cast<float>(result[0]);
00728         for (i=1; i<4; i++) if (static_cast<float>(result[i])<worst) worst=static_cast<float>(result[i]);
00729 
00730         EXITFUNC;
00731         return (float) (negative*worst);
00732 }
00733 
00734 float OptVarianceCmp::cmp(EMData * image, EMData *with) const
00735 {
00736         ENTERFUNC;
00737         validate_input_args(image, with);
00738 
00739         int keepzero = params.set_default("keepzero", 1);
00740         int invert = params.set_default("invert",0);
00741         int matchfilt = params.set_default("matchfilt",1);
00742         int matchamp = params.set_default("matchamp",0);
00743         int radweight = params.set_default("radweight",0);
00744         int dbug = params.set_default("debug",0);
00745 
00746         size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00747 
00748 
00749         EMData *with2=NULL;
00750         if (matchfilt) {
00751                 EMData *a = image->do_fft();
00752                 EMData *b = with->do_fft();
00753 
00754                 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1);
00755                 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1);
00756 
00757                 float avg=0;
00758                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00759                         rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i]));
00760                         avg+=rfa[i];
00761                 }
00762 
00763                 avg/=a->get_ysize()/2.0f;
00764                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00765                         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
00766                 }
00767                 rfa[0]=0.0;
00768 
00769                 if (dbug) b->write_image("a.hdf",-1);
00770 
00771                 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa);
00772                 with2=b->do_ift();
00773 
00774                 if (dbug) b->write_image("a.hdf",-1);
00775                 if (dbug) a->write_image("a.hdf",-1);
00776 
00777 /*              if (dbug) {
00778                         FILE *out=fopen("a.txt","w");
00779                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]);
00780                         fclose(out);
00781 
00782                         out=fopen("b.txt","w");
00783                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]);
00784                         fclose(out);
00785                 }*/
00786 
00787 
00788                 delete a;
00789                 delete b;
00790 
00791                 if (dbug) {
00792                         with2->write_image("a.hdf",-1);
00793                         image->write_image("a.hdf",-1);
00794                 }
00795 
00796 //              with2->process_inplace("matchfilt",Dict("to",this));
00797 //              x_data = with2->get_data();
00798         }
00799 
00800         // This applies the individual Fourier amplitudes from 'image' and
00801         // applies them to 'with'
00802         if (matchamp) {
00803                 EMData *a = image->do_fft();
00804                 EMData *b = with->do_fft();
00805                 size_t size2 = (size_t)a->get_xsize() * a->get_ysize() * a->get_zsize();
00806 
00807                 a->ri2ap();
00808                 b->ri2ap();
00809 
00810                 const float *const ad=a->get_const_data();
00811                 float * bd=b->get_data();
00812 
00813                 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i];
00814                 b->update();
00815 
00816                 b->ap2ri();
00817                 with2=b->do_ift();
00818 //with2->write_image("a.hdf",-1);
00819                 delete a;
00820                 delete b;
00821         }
00822 
00823         const float * x_data;
00824         if (with2) x_data=with2->get_const_data();
00825         else x_data = with->get_const_data();
00826         const float *const y_data = image->get_const_data();
00827 
00828         size_t nx = image->get_xsize();
00829         float m = 0;
00830         float b = 0;
00831 
00832         // This will write the x vs y file used to calculate the density
00833         // optimization. This behavior may change in the future
00834         if (dbug) {
00835                 FILE *out=fopen("dbug.optvar.txt","w");
00836                 if (out) {
00837                         for (size_t i=0; i<size; i++) {
00838                                 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]);
00839                         }
00840                         fclose(out);
00841                 }
00842         }
00843 
00844 
00845         Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero);
00846         if (m == 0) {
00847                 m = FLT_MIN;
00848         }
00849         b = -b / m;
00850         m = 1.0f / m;
00851 
00852         // While negative slopes are really not a valid comparison in most cases, we
00853         // still want to detect these instances, so this if is removed
00854 /*      if (m < 0) {
00855                 b = 0;
00856                 m = 1000.0;
00857         }*/
00858 
00859         double  result = 0;
00860         int count = 0;
00861 
00862         if (radweight) {
00863                 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only");
00864                 if (keepzero) {
00865                         for (size_t i = 0,y=0; i < size; y++) {
00866                                 for (size_t x=0; x<nx; i++,x++) {
00867                                         if (y_data[i] && x_data[i]) {
00868 #ifdef  _WIN32
00869                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00870                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00871 #else
00872                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00873                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00874 #endif
00875                                                 count++;
00876                                         }
00877                                 }
00878                         }
00879                         result/=count;
00880                 }
00881                 else {
00882                         for (size_t i = 0,y=0; i < size; y++) {
00883                                 for (size_t x=0; x<nx; i++,x++) {
00884 #ifdef  _WIN32
00885                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00886                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00887 #else
00888                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
00889                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
00890 #endif
00891                                 }
00892                         }
00893                         result = result / size;
00894                 }
00895         }
00896         else {
00897                 if (keepzero) {
00898                         for (size_t i = 0; i < size; i++) {
00899                                 if (y_data[i] && x_data[i]) {
00900                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00901                                         else result += Util::square((x_data[i] * m) + b - y_data[i]);
00902                                         count++;
00903                                 }
00904                         }
00905                         result/=count;
00906                 }
00907                 else {
00908                         for (size_t i = 0; i < size; i++) {
00909                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
00910                                 else result += Util::square((x_data[i] * m) + b - y_data[i]);
00911                         }
00912                         result = result / size;
00913                 }
00914         }
00915         scale = m;
00916         shift = b;
00917 
00918         image->set_attr("ovcmp_m",m);
00919         image->set_attr("ovcmp_b",b);
00920         if (with2) delete with2;
00921         EXITFUNC;
00922 
00923 #if 0
00924         return (1 - result);
00925 #endif
00926 
00927         return static_cast<float>(result);
00928 }
00929 
00930 float PhaseCmp::cmp(EMData * image, EMData *with) const
00931 {
00932         ENTERFUNC;
00933 
00934         int snrweight = params.set_default("snrweight", 0);
00935         int snrfn = params.set_default("snrfn",0);
00936         int ampweight = params.set_default("ampweight",0);
00937         int zeromask = params.set_default("zeromask",0);
00938         float minres = params.set_default("minres",500.0f);
00939         float maxres = params.set_default("maxres",2.0f);
00940 
00941         if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator");
00942 
00943         EMData *image_fft = NULL;
00944         EMData *with_fft = NULL;
00945 
00946         int ny = image->get_ysize();
00947 //      int np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2;
00948         int np = 0;
00949         vector<float> snr;
00950 
00951         // weighting based on SNR estimate from CTF
00952         if (snrweight) {
00953                 Ctf *ctf = NULL;
00954                 if (!image->has_attr("ctf")) {
00955                         if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
00956                         ctf=with->get_attr("ctf");
00957                 }
00958                 else ctf=image->get_attr("ctf");
00959 
00960                 float ds=1.0f/(ctf->apix*ny);
00961                 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR); // note that this returns ny/2 values
00962                 if(ctf) {delete ctf; ctf=0;}
00963                 np=snr.size();
00964         }
00965         // weighting based on empirical SNR function (not really good)
00966         else if (snrfn==1) {
00967                 np=ny/2;
00968                 snr.resize(np);
00969 
00970                 for (int i = 0; i < np; i++) {
00971                         float x2 = 10.0f*i/np;
00972                         snr[i] = x2 * exp(-x2);
00973                 }
00974         }
00975         // no weighting
00976         else {
00977                 np=ny/2;
00978                 snr.resize(np);
00979 
00980                 for (int i = 0; i < np; i++) snr[i]=1.0;
00981         }
00982 
00983         // Min/max modifications to weighting
00984         float pmin,pmax;
00985         if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;                //cutoff in pixels, assume square
00986         else pmin=0;
00987         if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
00988         else pmax=0;
00989 
00990 //      printf("%f\t%f\n",pmin,pmax);
00991 
00992         // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation
00993         for (int i = 0; i < np; i++) {
00994                 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f;
00995                 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f;
00996 //              printf("%d\t%f\n",i,snr[i]);
00997         }
00998 
00999         if (zeromask) {
01000                 image_fft=image->copy();
01001                 with_fft=with->copy();
01002                 
01003                 if (image_fft->is_complex()) image_fft->do_ift_inplace();
01004                 if (with_fft->is_complex()) with_fft->do_ift_inplace();
01005                 
01006                 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize();
01007                 float *d1=image_fft->get_data();
01008                 float *d2=with_fft->get_data();
01009                 
01010                 for (int i=0; i<sz; i++) {
01011                         if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
01012                 }
01013                 
01014                 image_fft->update();
01015                 with_fft->update();
01016                 image_fft->do_fft_inplace();
01017                 with_fft->do_fft_inplace();
01018                 image_fft->set_attr("free_me",1); 
01019                 with_fft->set_attr("free_me",1); 
01020         }
01021         else {
01022                 if (image->is_complex()) image_fft=image;
01023                 else {
01024                         image_fft=image->do_fft();
01025                         image_fft->set_attr("free_me",1);
01026                 }
01027                 
01028                 if (with->is_complex()) with_fft=with;
01029                 else {
01030                         with_fft=with->do_fft();
01031                         with_fft->set_attr("free_me",1);
01032                 }
01033         }
01034 //      image_fft->ri2ap();
01035 //      with_fft->ri2ap();
01036 
01037         const float *const image_fft_data = image_fft->get_const_data();
01038         const float *const with_fft_data = with_fft->get_const_data();
01039         double sum = 0;
01040         double norm = FLT_MIN;
01041         size_t i = 0;
01042 //      int nx=image_fft->get_xsize();
01043             ny=image_fft->get_ysize();
01044         int nz=image_fft->get_zsize();
01045         int nx2=image_fft->get_xsize()/2;
01046         int ny2=image_fft->get_ysize()/2;
01047 //      int nz2=image_fft->get_zsize()/2;
01048 
01049         // This can never happen any more...
01050         if (np==0) {
01051                 for (int z = 0; z < nz; z++){
01052                         for (int y = 0; y < ny; y++) {
01053                                 for (int x = 0; x < nx2; x ++) {
01054                                         float a;
01055                                         if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01056                                         else a=1.0f;
01057                                         sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01058                                         norm += a;
01059                                         i += 2;
01060                                 }
01061                         }
01062                 }
01063                 
01064         }
01065         else if (nz==1) {
01066                 for (int y = 0; y < ny; y++) {
01067                         for (int x = 0; x < nx2; x ++) {
01068                                 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y);
01069                                 if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
01070                                 
01071                                 float a;
01072                                 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01073                                 else a=1.0f;
01074                                 a*=snr[r];
01075                                 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01076                                 norm += a;
01077                                 i += 2;
01078                         }
01079                 }
01080         }
01081         else {
01082                 for (int z = 0; z < nz; z++){
01083                         for (int y = 0; y < ny; y++) {
01084                                 for (int x = 0; x < nx2; x ++) {
01085                                         int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z);
01086                                         if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
01087                                         
01088                                         float a;
01089                                         if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01090                                         else a=1.0f;
01091                                         a*=snr[r];
01092                                         sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01093                                         norm += a;
01094                                         i += 2;
01095                                 } 
01096                         }
01097                 }
01098                 
01099         }
01100 
01101         EXITFUNC;
01102 
01103         if( image_fft->has_attr("free_me") )
01104         {
01105                 delete image_fft;
01106                 image_fft = 0;
01107         }
01108         if( with_fft->has_attr("free_me") )
01109         {
01110                 delete with_fft;
01111                 with_fft = 0;
01112         }
01113 #if 0
01114         return (1.0f - sum / norm);
01115 #endif
01116         return (float)(sum / norm);
01117 }
01118 
01119 float FRCCmp::cmp(EMData * image, EMData * with) const
01120 {
01121         ENTERFUNC;
01122         validate_input_args(image, with);
01123 
01124         int snrweight = params.set_default("snrweight", 0);
01125         int ampweight = params.set_default("ampweight", 0);
01126         int sweight = params.set_default("sweight", 1);
01127         int nweight = params.set_default("nweight", 0);
01128         int zeromask = params.set_default("zeromask",0);
01129         float minres = params.set_default("minres",500.0f);
01130         float maxres = params.set_default("maxres",2.0f);
01131 
01132         vector < float >fsc;
01133         bool use_cpu = true;
01134 
01135         if (use_cpu) {
01136                 if (zeromask) {
01137                         image=image->copy();
01138                         with=with->copy();
01139                 
01140                         int sz=image->get_xsize()*image->get_ysize()*image->get_zsize();
01141                         float *d1=image->get_data();
01142                         float *d2=with->get_data();
01143                 
01144                         for (int i=0; i<sz; i++) {
01145                                 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
01146                         }
01147                 
01148                         image->update();
01149                         with->update();
01150                         image->do_fft_inplace();
01151                         with->do_fft_inplace();
01152                         image->set_attr("free_me",1); 
01153                         with->set_attr("free_me",1); 
01154                 }
01155 
01156 
01157                 if (!image->is_complex()) {
01158                         image=image->do_fft(); 
01159                         image->set_attr("free_me",1); 
01160                 }
01161                 if (!with->is_complex()) { 
01162                         with=with->do_fft(); 
01163                         with->set_attr("free_me",1); 
01164                 }
01165 
01166                 fsc = image->calc_fourier_shell_correlation(with,1);
01167         }
01168         
01169         int ny = image->get_ysize();
01170         int ny2=ny/2+1;
01171 
01172         // The fast hypot here was supposed to speed things up. Little effect
01173 //      if (image->get_zsize()>1) fsc = image->calc_fourier_shell_correlation(with,1);
01174 //      else {
01175 //              double *sxy = (double *)malloc(ny2*sizeof(double)*4);
01176 //              double *sxx = sxy+ny2;
01177 //              double *syy = sxy+2*ny2;
01178 //              double *norm= sxy+3*ny2;
01179 //
01180 //              float *df1=image->get_data();
01181 //              float *df2=with->get_data();
01182 //              int nx2=image->get_xsize();
01183 //
01184 //              for (int y=-ny/2; y<ny/2; y++) {
01185 //                      for (int x=0; x<nx2/2; x++) {
01186 //                              if (x==0 && y<0) continue;      // skip Friedel pair
01187 //                              short r=Util::hypot_fast_int(x,y);
01188 //                              if (r>ny2-1) continue;
01189 //                              int l=x*2+(y<0?ny+y:y)*nx2;
01190 //                              sxy[r]+=df1[l]*df2[l]+df1[l+1]*df2[l+1];
01191 //                              sxx[r]+=df1[l]*df1[l];
01192 //                              syy[r]+=df2[l]*df2[l];
01193 //                              norm[r]+=1.0;
01194 //                      }
01195 //              }
01196 //              fsc.resize(ny2*3);
01197 //              for (int r=0; r<ny2; r++) {
01198 //                      fsc[r]=r*0.5/ny2;
01199 //                      fsc[ny2+r]=sxy[r]/(sqrt(sxx[r])*sqrt(syy[r]));
01200 //                      fsc[ny2*2+r]=norm[r];
01201 //              }
01202 //              free(sxy);
01203 //      }
01204 
01205         vector<float> snr;
01206         if (snrweight) {
01207                 Ctf *ctf = NULL;
01208                 if (!image->has_attr("ctf")) {
01209                         if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
01210                         ctf=with->get_attr("ctf");
01211                 }
01212                 else ctf=image->get_attr("ctf");
01213 
01214                 float ds=1.0f/(ctf->apix*ny);
01215                 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR);
01216                 if(ctf) {delete ctf; ctf=0;}
01217         }
01218 
01219         vector<float> amp;
01220         if (ampweight) amp=image->calc_radial_dist(ny/2,0,1,0);
01221 
01222         // Min/max modifications to weighting
01223         float pmin,pmax;
01224         if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;                //cutoff in pixels, assume square
01225         else pmin=0;
01226         if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
01227         else pmax=0;
01228 
01229         double sum=0.0, norm=0.0;
01230 
01231         for (int i=0; i<ny/2; i++) {
01232                 double weight=1.0;
01233                 if (sweight) weight*=fsc[(ny2)*2+i];
01234                 if (ampweight) weight*=amp[i];
01235                 if (snrweight) weight*=snr[i];
01236 //              if (snrweight)  {
01237 //                      if (snr[i]>0) weight*=sqrt(snr[i]);
01238 //                      else weight=0;
01239 //              }
01240 //if(snr[i]<0) printf("snr[%d] = %1.5g\n",i,snr[i]);
01241                 if (pmin>0) weight*=(tanh(5.0*(i-pmin)/pmin)+1.0)/2.0;
01242                 if (pmax>0) weight*=(1.0-tanh(i-pmax))/2.0;
01243                 
01244                 sum+=weight*fsc[ny2+i];
01245                 norm+=weight;
01246 //              printf("%d\t%f\t%f\n",i,weight,fsc[ny/2+1+i]);
01247         }
01248 
01249         // This performs a weighting that tries to normalize FRC by correcting from the number of particles represented by the average
01250         sum/=norm;
01251         if (nweight && with->get_attr_default("ptcl_repr",0) && sum>=0 && sum<1.0) {
01252                 sum=sum/(1.0-sum);                                                      // convert to SNR
01253                 sum/=(float)with->get_attr_default("ptcl_repr",0);      // divide by ptcl represented
01254                 sum=sum/(1.0+sum);                                                      // convert back to correlation
01255         }
01256 
01257         if (image->has_attr("free_me")) delete image;
01258         if (with->has_attr("free_me")) delete with;
01259 
01260         EXITFUNC;
01261 
01262 
01263         if (!Util::goodf(&sum)) sum=-3.0e38;
01264         //.Note the negative! This is because EMAN2 follows the convention that
01265         // smaller return values from comparitors indicate higher similarity -
01266         // this enables comparitors to be used in a generic fashion.
01267         return (float)-sum;
01268 }
01269 
01270 void EMAN::dump_cmps()
01271 {
01272         dump_factory < Cmp > ();
01273 }
01274 
01275 map<string, vector<string> > EMAN::dump_cmps_list()
01276 {
01277         return dump_factory_list < Cmp > ();
01278 }
01279 
01280 /* vim: set ts=4 noet: */

Generated on Thu Nov 17 12:43:44 2011 for EMAN2 by  doxygen 1.3.9.1