EMAN2
cmp.cpp
Go to the documentation of this file.
00001 /*
00002  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00003  * Copyright (c) 2000-2006 Baylor College of Medicine
00004  *
00005  * This software is issued under a joint BSD/GNU license. You may use the
00006  * source code in this file under either license. However, note that the
00007  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00008  * so you are responsible for compliance with the licenses of these packages
00009  * if you opt to use BSD licensing. The warranty disclaimer below holds
00010  * in either instance.
00011  *
00012  * This complete copyright notice must be included in any revised version of the
00013  * source code. Additional authorship citations may be added, but existing
00014  * author citations must be preserved.
00015  *
00016  * This program is free software; you can redistribute it and/or modify
00017  * it under the terms of the GNU General Public License as published by
00018  * the Free Software Foundation; either version 2 of the License, or
00019  * (at your option) any later version.
00020  *
00021  * This program is distributed in the hope that it will be useful,
00022  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00023  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00024  * GNU General Public License for more details.
00025  *
00026  * You should have received a copy of the GNU General Public License
00027  * along with this program; if not, write to the Free Software
00028  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00029  *
00030  * */
00031 
00032 #include "cmp.h"
00033 #include "emdata.h"
00034 #include "ctf.h"
00035 #include "plugins/cmp_template.h"
00036 #undef max
00037 #include <climits>
00038 
00039 #ifdef EMAN2_USING_CUDA
00040 // Only CCC, DOT  and CCC.TOMO are cuda enabled
00041 #include "cuda/cuda_processor.h"
00042 #include "cuda/cuda_cmp.h"
00043 #endif // EMAN2_USING_CUDA
00044 
00045 using namespace EMAN;
00046 
00047 const string CccCmp::NAME = "ccc";
00048 const string LodCmp::NAME = "lod";
00049 const string SqEuclideanCmp::NAME = "sqeuclidean";
00050 const string DotCmp::NAME = "dot";
00051 const string TomoCccCmp::NAME = "ccc.tomo";
00052 const string TomoFscCmp::NAME = "fsc.tomo";
00053 const string QuadMinDotCmp::NAME = "quadmindot";
00054 const string OptVarianceCmp::NAME = "optvariance";
00055 const string OptSubCmp::NAME = "optsub";
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<TomoFscCmp>();
00067         force_add<QuadMinDotCmp>();
00068         force_add<OptVarianceCmp>();
00069         force_add<OptSubCmp>();
00070         force_add<PhaseCmp>();
00071         force_add<FRCCmp>();
00072 //      force_add<XYZCmp>();
00073 }
00074 
00075 void Cmp::validate_input_args(const EMData * image, const EMData *with) const
00076 {
00077         
00078 #ifdef EMAN2_USING_CUDA  
00079         if (image->getcudarwdata() && with->getcudarwdata()) {
00080                 //no need for futher checking, which will induce an expensive copy from device to host
00081                 return;
00082         }
00083 #endif
00084         if (!image) {
00085                 throw NullPointerException("compared image");
00086         }
00087         if (!with) {
00088                 throw NullPointerException("compare-with image");
00089         }
00090 
00091         if (!EMUtil::is_same_size(image, with)) {
00092                 throw ImageFormatException( "images not same size");
00093         }
00094 
00095         float *d1 = image->get_data();
00096         if (!d1) {
00097                 throw NullPointerException("image contains no data");
00098         }
00099 
00100         float *d2 = with->get_data();
00101         if (!d2) {
00102                 throw NullPointerException("compare-with image data");
00103         }
00104 }
00105 
00106 //  It would be good to add code for complex images!  PAP
00107 float CccCmp::cmp(EMData * image, EMData *with) const
00108 {
00109         ENTERFUNC;
00110         if (image->is_complex() || with->is_complex())
00111                 throw ImageFormatException( "Complex images not supported by CMP::CccCmp");
00112         validate_input_args(image, with);
00113 
00114         const float *const d1 = image->get_const_data();
00115         const float *const d2 = with->get_const_data();
00116 
00117         float negative = (float)params.set_default("negative", 1);
00118         if (negative) negative=-1.0; else negative=1.0;
00119 
00120         double avg1 = 0.0, var1 = 0.0, avg2 = 0.0, var2 = 0.0, ccc = 0.0;
00121         long n = 0;
00122         size_t totsize = image->get_xsize()*image->get_ysize()*image->get_zsize();
00123 
00124         bool has_mask = false;
00125         EMData* mask = 0;
00126         if (params.has_key("mask")) {
00127                 mask = params["mask"];
00128                 if(mask!=0) {has_mask=true;}
00129         }
00130 #ifdef EMAN2_USING_CUDA
00131         if (image->getcudarwdata() && with->getcudarwdata()) {
00132                 //cout << "CUDA ccc cmp" << endl;
00133                 float* maskdata = 0;
00134                 if(has_mask && !mask->getcudarwdata()){
00135                         mask->copy_to_cuda();
00136                         maskdata = mask->getcudarwdata();
00137                 }
00138                 float ccc = ccc_cmp_cuda(image->getcudarwdata(), with->getcudarwdata(), maskdata, image->get_xsize(), image->get_ysize(), image->get_zsize());
00139                 ccc *= negative;
00140                 //cout << "CUDA CCC is: " << ccc << endl;
00141                 return ccc;
00142         }
00143 #endif
00144         if (has_mask) {
00145                 const float *const dm = mask->get_const_data();
00146                 for (size_t i = 0; i < totsize; ++i) {
00147                         if (dm[i] > 0.5) {
00148                                 avg1 += double(d1[i]);
00149                                 var1 += d1[i]*double(d1[i]);
00150                                 avg2 += double(d2[i]);
00151                                 var2 += d2[i]*double(d2[i]);
00152                                 ccc += d1[i]*double(d2[i]);
00153                                 n++;
00154                         }
00155                 }
00156         } else {
00157                 for (size_t i = 0; i < totsize; ++i) {
00158                         avg1 += double(d1[i]);
00159                         var1 += d1[i]*double(d1[i]);
00160                         avg2 += double(d2[i]);
00161                         var2 += d2[i]*double(d2[i]);
00162                         ccc += d1[i]*double(d2[i]);
00163                 }
00164                 n = totsize;
00165         }
00166 
00167         avg1 /= double(n);
00168         var1 = var1/double(n) - avg1*avg1;
00169         avg2 /= double(n);
00170         var2 = var2/double(n) - avg2*avg2;
00171         ccc = ccc/double(n) - avg1*avg2;
00172         ccc /= sqrt(var1*var2);
00173         if (!Util::goodf(&ccc)) ccc=-2.0;               // Steve - if one image was 0, this returned nan, which messes certain things up. -2.0 is out of range, and should serve as a proxy
00174         ccc *= negative;
00175         return static_cast<float>(ccc);
00176         EXITFUNC;
00177 }
00178 
00179 
00180 // Added by JAK 11/12/10
00181 // L^1-norm difference of two maps, after normalization.
00182 float LodCmp::cmp(EMData * image, EMData *with) const
00183 {
00184         ENTERFUNC;
00185         if (image->is_complex() || with->is_complex())
00186                 throw ImageFormatException( "Complex images not yet supported by CMP::LodCmp");
00187         validate_input_args(image, with);
00188 
00189         const float *const d1 = image->get_const_data();
00190         const float *const d2 = with->get_const_data();
00191 
00192         float negative = (float)params.set_default("negative", 1);
00193         if (negative) negative=-1.0; else negative=1.0;
00194 
00195         double avg1 = 0.0, sig1 = 0.0, avg2 = 0.0, sig2 = 0.0, lod = 0.0;
00196         long n = 0;
00197         size_t totsize = image->get_xsize()*image->get_ysize()*image->get_zsize();
00198         size_t i;
00199         
00200         bool has_mask = false;
00201         EMData* mask = 0;
00202         if (params.has_key("mask")) {
00203                 mask = params["mask"];
00204                 if(mask!=0) {has_mask=true;}
00205         }
00206 
00207         int normalize = 0;
00208         if (params.has_key("normalize")) {
00209                 normalize = params["normalize"];
00210         }
00211 
00212         if (normalize) {
00213                 if (has_mask) {
00214                         const float *const dm = mask->get_const_data();
00215                         for (i = 0; i < totsize; ++i) {
00216                                 if (dm[i] > 0.5) {
00217                                         avg1 += double(d1[i]);
00218                                         avg2 += double(d2[i]);
00219                                         n++;
00220                                 }
00221                         }
00222                 } else {
00223                         for (i = 0; i < totsize; ++i) {
00224                                 avg1 += double(d1[i]);
00225                                 avg2 += double(d2[i]);
00226                         }
00227                         n = totsize;
00228                 }
00229 
00230                 avg1 /= double(n);
00231                 avg2 /= double(n);
00232 
00233                 if (has_mask) {
00234                         const float *const dm = mask->get_const_data();
00235                         for (i = 0; i < totsize; ++i) {
00236                                 if (dm[i] > 0.5) {
00237                                         sig1 += fabs(double(d1[i])-avg1);
00238                                         sig2 += fabs(double(d2[i])-avg2);
00239                                 }
00240                         }
00241                 } else {
00242                         for (i = 0; i < totsize; ++i) {
00243                                 sig1 += fabs(double(d1[i])-avg1);
00244                                 sig2 += fabs(double(d2[i])-avg2);
00245                         }
00246                 }
00247         } else {
00248                 avg1 = 0.0; avg2 = 0.0;
00249                 sig1 = 1.0; sig2 = 1.0;
00250         }
00251 
00252         if (has_mask) {
00253                 const float *const dm = mask->get_const_data();
00254                 for (i = 0; i < totsize; ++i) {
00255                         if (dm[i] > 0.5) {
00256                                 lod += fabs((double(d1[i])-avg1)/sig1 - (double(d2[i])-avg2)/sig2);
00257                         }
00258                 }
00259         } else {
00260                 for (i = 0; i < totsize; ++i) {
00261                         lod += fabs((double(d1[i])-avg1)/sig1 - (double(d2[i])-avg2)/sig2);
00262                 }
00263         }
00264         
00265         lod *= (-0.5);
00266         lod *= negative;
00267         return static_cast<float>(lod);
00268         EXITFUNC;
00269 }
00270 
00271 
00272 //float SqEuclideanCmp::cmp(EMData * image, EMData *withorig) const
00273 float SqEuclideanCmp::cmp(EMData *image,EMData * withorig ) const
00274 {
00275         ENTERFUNC;
00276         EMData *with = withorig;
00277         validate_input_args(image, with);
00278 
00279         int zeromask = params.set_default("zeromask",0);
00280         int normto = params.set_default("normto",0);
00281 
00282         if (normto) {
00283                 if (zeromask) with = withorig->process("normalize.toimage",Dict("to",image));
00284                 else with = withorig->process("normalize.toimage",Dict("to",image,"ignore_zero",0));
00285                 with->set_attr("deleteme",1);
00286                 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
00287                         delete with;
00288                         with=withorig;
00289                 }
00290         }
00291 
00292         const float *const y_data = with->get_const_data();
00293         const float *const x_data = image->get_const_data();
00294         double result = 0.;
00295         float n = 0.0f;
00296         if(image->is_complex() && with->is_complex()) {
00297         // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
00298                 int nx  = with->get_xsize();
00299                 int ny  = with->get_ysize();
00300                 int nz  = with->get_zsize();
00301                 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
00302                 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
00303 
00304                 int ixb = 2*((nx+1)%2);
00305                 int iyb = ny%2;
00306                 //
00307                 if(nz == 1) {
00308                 //  it looks like it could work in 3D, but it does not, really.
00309                 for ( int iz = 0; iz <= nz-1; iz++) {
00310                         double part = 0.;
00311                         for ( int iy = 0; iy <= ny-1; iy++) {
00312                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ix++) {
00313                                                 size_t ii = ix + (iy  + iz * ny)* lsd2;
00314                                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00315                                 }
00316                         }
00317                         for ( int iy = 1; iy <= ny/2-1 + iyb; iy++) {
00318                                 size_t ii = (iy  + iz * ny)* lsd2;
00319                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00320                                 part += (x_data[ii+1] - y_data[ii+1])*double(x_data[ii+1] - y_data[ii+1]);
00321                         }
00322                         if(nx%2 == 0) {
00323                                 for ( int iy = 1; iy <= ny/2-1 + iyb; iy++) {
00324                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00325                                         part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00326                                         part += (x_data[ii+1] - y_data[ii+1])*double(x_data[ii+1] - y_data[ii+1]);
00327                                 }
00328 
00329                         }
00330                         part *= 2;
00331                         part += (x_data[0] - y_data[0])*double(x_data[0] - y_data[0]);
00332                         if(ny%2 == 0) {
00333                                 int ii = (ny/2  + iz * ny)* lsd2;
00334                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00335                         }
00336                         if(nx%2 == 0) {
00337                                 int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00338                                 part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00339                                 if(ny%2 == 0) {
00340                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00341                                         part += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00342                                 }
00343                         }
00344                         result += part;
00345                 }
00346                 n = (float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz;
00347 
00348                 } else { //This 3D code is incorrect, but it is the best I can do now 01/09/06 PAP
00349                 int ky, kz;
00350                 int ny2 = ny/2; int nz2 = nz/2;
00351                 for ( int iz = 0; iz <= nz-1; iz++) {
00352                         if(iz>nz2) kz=iz-nz; else kz=iz;
00353                         for ( int iy = 0; iy <= ny-1; iy++) {
00354                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00355                                 for ( int ix = 0; ix <= lsd2-1; ix++) {
00356                                 // Skip Friedel related values
00357                                 if(ix>0 || (kz>=0 && (ky>=0 || kz!=0))) {
00358                                                 size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
00359                                                 result += (x_data[ii] - y_data[ii])*double(x_data[ii] - y_data[ii]);
00360                                         }
00361                                 }
00362                         }
00363                 }
00364                 n = ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz)/2.0f;
00365                 }
00366         } else {                // real space
00367                 size_t totsize = (size_t)image->get_xsize()*image->get_ysize()*image->get_zsize();
00368                 if (params.has_key("mask")) {
00369                   EMData* mask;
00370                   mask = params["mask"];
00371                   const float *const dm = mask->get_const_data();
00372                   for (size_t i = 0; i < totsize; i++) {
00373                            if (dm[i] > 0.5) {
00374                                 double temp = x_data[i]- y_data[i];
00375                                 result += temp*temp;
00376                                 n++;
00377                            }
00378                   }
00379                 } 
00380                 else if (zeromask) {
00381                         n=0;
00382                         for (size_t i = 0; i < totsize; i++) {
00383                                 if (x_data[i]==0 || y_data[i]==0) continue;
00384                                 double temp = x_data[i]- y_data[i];
00385                                 result += temp*temp;
00386                                 n++;
00387                         }
00388                         
00389                 }
00390                 else {
00391                   for (size_t i = 0; i < totsize; i++) {
00392                                 double temp = x_data[i]- y_data[i];
00393                                 result += temp*temp;
00394                    }
00395                    n = (float)totsize;
00396                 }
00397         }
00398         result/=n;
00399 
00400         EXITFUNC;
00401         if (with->has_attr("deleteme")) delete with;
00402         float ret = (float)result;
00403         if (!Util::goodf(&ret)) return FLT_MAX;
00404         return ret;
00405 }
00406 
00407 
00408 // Even though this uses doubles, it might be wise to recode it row-wise
00409 // to avoid numerical errors on large images
00410 float DotCmp::cmp(EMData* image, EMData* with) const
00411 {
00412         ENTERFUNC;
00413         
00414         validate_input_args(image, with);
00415 
00416         int normalize = params.set_default("normalize", 0);
00417         float negative = (float)params.set_default("negative", 1);
00418         if (negative) negative=-1.0; else negative=1.0;
00419 #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!!!!)
00420         if(image->is_complex() && with->is_complex()) {
00421         } else {
00422                 if (image->getcudarwdata() && with->getcudarwdata()) {
00423                         //cout << "CUDA dot cmp" << endl;
00424                         float* maskdata = 0;
00425                         bool has_mask = false;
00426                         EMData* mask = 0;
00427                         if (params.has_key("mask")) {
00428                                 mask = params["mask"];
00429                                 if(mask!=0) {has_mask=true;}
00430                         }
00431                         if(has_mask && !mask->getcudarwdata()){
00432                                 mask->copy_to_cuda();
00433                                 maskdata = mask->getcudarwdata();
00434                         }
00435 
00436                         float result = dot_cmp_cuda(image->getcudarwdata(), with->getcudarwdata(), maskdata, image->get_xsize(), image->get_ysize(), image->get_zsize());
00437                         result *= negative;
00438 
00439                         return result;
00440                         
00441                 }
00442         }
00443 #endif
00444         const float *const x_data = image->get_const_data();
00445         const float *const y_data = with->get_const_data();
00446 
00447         double result = 0.;
00448         long n = 0;
00449         if(image->is_complex() && with->is_complex()) {
00450         // Implemented by PAP  01/09/06 - please do not change.  If in doubts, write/call me.
00451                 int nx  = with->get_xsize();
00452                 int ny  = with->get_ysize();
00453                 int nz  = with->get_zsize();
00454                 nx = (nx - 2 + with->is_fftodd()); // nx is the real-space size of the input image
00455                 int lsd2 = (nx + 2 - nx%2) ; // Extended x-dimension of the complex image
00456 
00457                 int ixb = 2*((nx+1)%2);
00458                 int iyb = ny%2;
00459                 //
00460                 if(nz == 1) {
00461                 //  it looks like it could work in 3D, but does not
00462                 for ( int iz = 0; iz <= nz-1; ++iz) {
00463                         double part = 0.;
00464                         for ( int iy = 0; iy <= ny-1; ++iy) {
00465                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00466                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00467                                         part += x_data[ii] * double(y_data[ii]);
00468                                 }
00469                         }
00470                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00471                                 size_t ii = (iy  + iz * ny)* lsd2;
00472                                 part += x_data[ii] * double(y_data[ii]);
00473                                 part += x_data[ii+1] * double(y_data[ii+1]);
00474                         }
00475                         if(nx%2 == 0) {
00476                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00477                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00478                                         part += x_data[ii] * double(y_data[ii]);
00479                                         part += x_data[ii+1] * double(y_data[ii+1]);
00480                                 }
00481 
00482                         }
00483                         part *= 2;
00484                         part += x_data[0] * double(y_data[0]);
00485                         if(ny%2 == 0) {
00486                                 size_t ii = (ny/2  + iz * ny)* lsd2;
00487                                 part += x_data[ii] * double(y_data[ii]);
00488                         }
00489                         if(nx%2 == 0) {
00490                                 size_t ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00491                                 part += x_data[ii] * double(y_data[ii]);
00492                                 if(ny%2 == 0) {
00493                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00494                                         part += x_data[ii] * double(y_data[ii]);
00495                                 }
00496                         }
00497                         result += part;
00498                 }
00499                 if( normalize ) {
00500                 //  it looks like it could work in 3D, but does not
00501                 double square_sum1 = 0., square_sum2 = 0.;
00502                 for ( int iz = 0; iz <= nz-1; ++iz) {
00503                         for ( int iy = 0; iy <= ny-1; ++iy) {
00504                                 for ( int ix = 2; ix <= lsd2 - 1 - ixb; ++ix) {
00505                                         size_t ii = ix + (iy  + iz * ny)* lsd2;
00506                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00507                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00508                                 }
00509                         }
00510                         for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00511                                 size_t ii = (iy  + iz * ny)* lsd2;
00512                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00513                                 square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00514                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00515                                 square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00516                         }
00517                         if(nx%2 == 0) {
00518                                 for ( int iy = 1; iy <= ny/2-1 + iyb; ++iy) {
00519                                         size_t ii = lsd2 - 2 + (iy  + iz * ny)* lsd2;
00520                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00521                                         square_sum1 += x_data[ii+1] * double(x_data[ii+1]);
00522                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00523                                         square_sum2 += y_data[ii+1] * double(y_data[ii+1]);
00524                                 }
00525 
00526                         }
00527                         square_sum1 *= 2;
00528                         square_sum1 += x_data[0] * double(x_data[0]);
00529                         square_sum2 *= 2;
00530                         square_sum2 += y_data[0] * double(y_data[0]);
00531                         if(ny%2 == 0) {
00532                                 int ii = (ny/2  + iz * ny)* lsd2;
00533                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00534                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00535                         }
00536                         if(nx%2 == 0) {
00537                                 int ii = lsd2 - 2 + (0  + iz * ny)* lsd2;
00538                                 square_sum1 += x_data[ii] * double(x_data[ii]);
00539                                 square_sum2 += y_data[ii] * double(y_data[ii]);
00540                                 if(ny%2 == 0) {
00541                                         int ii = lsd2 - 2 +(ny/2  + iz * ny)* lsd2;
00542                                         square_sum1 += x_data[ii] * double(x_data[ii]);
00543                                         square_sum2 += y_data[ii] * double(y_data[ii]);
00544                                 }
00545                         }
00546                 }
00547                 result /= sqrt(square_sum1*square_sum2);
00548                 } 
00549                         else  result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
00550                         // This concludes the 2D part.
00551                 } else {
00552                 // The 3D code
00553                 int ky, kz;
00554                 int ny2 = ny/2; int nz2 = nz/2;
00555                 for ( int iz = 0; iz <= nz-1; ++iz) {
00556                         if(iz>nz2) kz=iz-nz; else kz=iz;
00557                         for ( int iy = 0; iy <= ny-1; ++iy) {
00558                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00559                                 for (int ix = 0; ix <= lsd2-1; ++ix) {
00560                                                 int p = 2;
00561                                                 if  (ix <= 1 || ix >= lsd2-2 && nx%2 == 0)   p = 1;
00562                                                 size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
00563                                                 result += p * x_data[ii] * double(y_data[ii]);
00564                                 }
00565                         }
00566                 }
00567                 if( normalize ) {
00568                 double square_sum1 = 0., square_sum2 = 0.;
00569                 int ky, kz;
00570                 int ny2 = ny/2; int nz2 = nz/2;
00571                 for ( int iz = 0; iz <= nz-1; ++iz) {
00572                         if(iz>nz2) kz=iz-nz; else kz=iz;
00573                         for ( int iy = 0; iy <= ny-1; ++iy) {
00574                                 if(iy>ny2) ky=iy-ny; else ky=iy;
00575                                 for (int ix = 0; ix <= lsd2-1; ++ix) {
00576                                                 int p = 2;
00577                                                 if  (ix <= 1 || ix >= lsd2-2 && nx%2 == 0)   p = 1;
00578                                                 size_t ii = ix + (iy  + iz * ny)* (size_t)lsd2;
00579                                                 square_sum1 += p * x_data[ii] * double(x_data[ii]);
00580                                                 square_sum2 += p * y_data[ii] * double(y_data[ii]);
00581                                 }
00582                         }
00583                 }
00584                 result /= sqrt(square_sum1*square_sum2);
00585                 } 
00586                         else result /= ((float)nx*(float)ny*(float)nz*(float)nx*(float)ny*(float)nz);
00587                 }
00588         } else {
00589                 // This part is for when two images are real, which is much easier because 2-D or 3-D
00590                 // is the same code.
00591                 size_t totsize = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00592 
00593                 double square_sum1 = 0., square_sum2 = 0.;
00594 
00595                 if (params.has_key("mask")) {
00596                         EMData* mask;
00597                         mask = params["mask"];
00598                         const float *const dm = mask->get_const_data();
00599                         if (normalize) {
00600                                 for (size_t i = 0; i < totsize; i++) {
00601                                         if (dm[i] > 0.5) {
00602                                                 square_sum1 += x_data[i]*double(x_data[i]);
00603                                                 square_sum2 += y_data[i]*double(y_data[i]);
00604                                                 result += x_data[i]*double(y_data[i]);
00605                                         }
00606                                 }
00607                         } else {
00608                                 for (size_t i = 0; i < totsize; i++) {
00609                                         if (dm[i] > 0.5) {
00610                                                 result += x_data[i]*double(y_data[i]);
00611                                                 n++;
00612                                         }
00613                                 }
00614                         }
00615                 } else {
00616                         // this little bit of manual loop unrolling makes the dot product as fast as sqeuclidean with -O2
00617                         for (size_t i=0; i<totsize; i++) result+=x_data[i]*y_data[i];
00618 
00619                         if (normalize) {
00620                                 square_sum1 = image->get_attr("square_sum");
00621                                 square_sum2 = with->get_attr("square_sum");
00622                         } else n = totsize;
00623                 }
00624                 if (normalize) result /= (sqrt(square_sum1*square_sum2)); else result /= n;
00625         }
00626 
00627 
00628         EXITFUNC;
00629         return (float) (negative*result);
00630 }
00631 
00632 // This implements the technique of Mike Schmid where by the cross correlation is normalized
00633 // in an effort to remove the effects of the missing wedge. Somewhat of a heuristic solution, but it seems
00634 // to work. Basically it relies on the observation that 'good' matchs will conentrate the correlation
00635 // signal in the peak, wheras 'bad' correlations will distribute the signal.
00636 // John Flanagan 18-10-2010
00637 
00638 float TomoCccCmp::cmp(EMData * image, EMData *with) const
00639 {
00640         ENTERFUNC;
00641         EMData* ccf = params.set_default("ccf",(EMData*) NULL);
00642         bool ccf_ownership = false;
00643         bool norm = params.set_default("norm",true);
00644         float negative = (float)params.set_default("negative", 1);
00645         if (negative) negative=-1.0; else negative=1.0;
00646         
00647 #ifdef EMAN2_USING_CUDA 
00648         if(image->getcudarwdata() && with->getcudarwdata()){
00649                 if (!ccf) {
00650                         ccf = image->calc_ccf(with);
00651                         ccf_ownership = true;
00652                 }
00653                 //cout << "using CUDA" << endl;
00654                 float2 stats = get_stats_cuda(ccf->getcudarwdata(), ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
00655                 float best_score = get_value_at_wrap_cuda(ccf->getcudarwdata(), 0, 0, 0, ccf->get_xsize(), ccf->get_ysize(), ccf->get_zsize());
00656                 if(norm) {
00657                         best_score = negative*(best_score - stats.x)/sqrt(stats.y);
00658                 } else {
00659                         best_score = negative*best_score;
00660                 }
00661                 
00662                 if (ccf_ownership) delete ccf; ccf = 0;
00663                 
00664                 EXITFUNC;
00665                 return best_score;
00666                 
00667         }
00668 #endif
00669 
00670         if (!ccf) {
00671                 ccf = image->calc_ccf(with);
00672                 ccf_ownership = true;
00673         }
00674         if (norm) ccf->process_inplace("normalize");
00675         
00676         float best_score = ccf->get_value_at_wrap(0,0,0);
00677         if (ccf_ownership) delete ccf; ccf = 0;
00678         
00679         return negative*best_score;
00680 
00681 }
00682 
00683 float TomoFscCmp::cmp(EMData * image, EMData *with) const
00684 {
00685         ENTERFUNC;
00686         bool usecpu = 1;
00687         bool del_imagefft = 0;
00688         bool del_withfft = 0;
00689         float score = 1.0f;
00690         
00691         //get parameters
00692         if (!image->has_attr("spt_wedge_mean") || !image->has_attr("spt_wedge_sigma"))  throw InvalidCallException("Rubbish!!! Image Subtomogram does not have mena and/or sigma amps metadata");
00693         // BAD DESIGN!!!! The fact that I have to load attrs into a variable before I can do operations on them is silly
00694         
00695         //get threshold information
00696         float image_meanwedgeamp = image->get_attr("spt_wedge_mean");
00697         float image_sigmawedgeamp = image->get_attr("spt_wedge_sigma");
00698         
00699         // Use with amp stats if avaliable otherwise effectivly ignore
00700         float with_meanwedgeamp = image_meanwedgeamp;
00701         float with_sigmawedgeamp = image_sigmawedgeamp;
00702         if (with->has_attr("spt_wedge_mean") && with->has_attr("spt_wedge_sigma"))
00703         {
00704                 with_meanwedgeamp = with->get_attr("spt_wedge_mean");
00705                 with_sigmawedgeamp = with->get_attr("spt_wedge_sigma");
00706         }
00707         
00708         // Find threshold
00709         float sigmas = params.set_default("sigmas",5.0f);
00710         float img_amp_thres = pow(image_meanwedgeamp + sigmas*image_sigmawedgeamp, 2.0f);
00711         float with_amp_thres = pow(with_meanwedgeamp + sigmas*with_sigmawedgeamp, 2.0f);
00712         
00713         // take negative of score
00714         float negative = (float)params.set_default("negative", 1.0f);
00715         if (negative) negative=-1.0; else negative=1.0;
00716         //get apix, use param apix, if not specified use apix_x, if this is not specified then apix=1.0
00717         float apix = params.set_default("apix",image->get_attr_default("apix_x", 1.0f));
00718         //get min and max res
00719         float minres = params.set_default("minres",std::numeric_limits<float>::max());
00720         float maxres = params.set_default("maxres", 0.0f);
00721         
00722         //Check to ensure that images are complex
00723         EMData* image_fft = image;
00724         EMData* with_fft = with;
00725         
00726 #ifdef EMAN2_USING_CUDA 
00727         //do CUDA FFT, does not allow minres, maxres yet
00728         if(EMData::usecuda == 1 && image->getcudarwdata() && with->getcudarwdata()) {
00729                 if(!image->is_complex()){
00730                         del_imagefft = 1;
00731                         image_fft = image->do_fft_cuda();
00732                 }
00733                 if(!with->is_complex()){
00734                         del_withfft = 1;
00735                         with_fft = with->do_fft_cuda();
00736                 }
00737                 score = fsc_tomo_cmp_cuda(image_fft->getcudarwdata(), with_fft->getcudarwdata(), img_amp_thres, with_amp_thres, 0.0, 0.0, image_fft->get_xsize(), image_fft->get_ysize(), image_fft->get_zsize());
00738                 usecpu = 0;
00739         }
00740 #endif  
00741         if(usecpu){
00742                 if(!image->is_complex()){
00743                         del_imagefft = 1;
00744                         image_fft = image->do_fft();
00745                 }
00746                 if(!with->is_complex()){
00747                         del_withfft = 1;
00748                         with_fft = with->do_fft();
00749                 }
00750                 
00751                 //loop over all voxels
00752                 int count = 0;
00753                 double sum_imgamp_sq = 0.0;
00754                 double sum_withamp_sq = 0.0;
00755                 double cong = 0.0;
00756                 float* img_data = image_fft->get_data();
00757                 float* with_data = with_fft->get_data();
00758         
00759                 int nx  = image_fft->get_xsize();
00760                 int ny  = image_fft->get_ysize();
00761                 int nz  = image_fft->get_zsize();
00762                 int ny2 = ny/2;
00763                 int nz2 = nz/2;
00764         
00765                 //compute FSC
00766                 int ii, kz, ky;
00767                 for (int iz = 0; iz < nz; iz++) {
00768                         if(iz > nz2) kz = nz-iz; else kz=iz;
00769                         for (int iy = 0; iy < ny; iy++) {
00770                                 if(iy > ny2) ky = ny-iy; else ky=iy;
00771                                 for (int ix = 0; ix < nx; ix+=2) {
00772                                         //compute spatial frequency and convert to resolution stat
00773                                         float freq = std::sqrt(kz*kz + ky*ky + ix*ix/4.0f)/float(nz);
00774                                         float reso = apix/freq;
00775                                         
00776                                         //only look within a resolution domain
00777                                         if(reso < minres && reso > maxres){
00778                                                 ii = ix + (iy  + iz * ny)* nx;
00779                                                 float img_r = img_data[ii];
00780                                                 float img_i = img_data[ii+1];
00781                                                 float with_r = with_data[ii];
00782                                                 float with_i = with_data[ii+1];
00783                                                 double img_amp_sq = img_r*img_r + img_i*img_i;
00784                                                 double with_amp_sq = with_r*with_r + with_i*with_i;
00785 
00786                                                 if((img_amp_sq >  img_amp_thres) && (with_amp_sq >  with_amp_thres)){
00787                                                         count ++;
00788                                                         sum_imgamp_sq += img_amp_sq;
00789                                                         sum_withamp_sq += with_amp_sq;
00790                                                         cong += img_r*with_r + img_i*with_i;
00791                                                 }
00792                                         }
00793                                 }
00794                         }
00795                 }
00796                 
00797                 if(count > 0){ 
00798                         score = (float)(cong/sqrt(sum_imgamp_sq*sum_withamp_sq));
00799                 }else{
00800                         score = 1.0f;
00801                 }
00802         }
00803         
00804         //avoid mem leaks
00805         if(del_imagefft) delete image_fft;
00806         if(del_withfft) delete with_fft;
00807         
00808         return negative*score;
00809         EXITFUNC;
00810 }
00811 
00812 // Even though this uses doubles, it might be wise to recode it row-wise
00813 // to avoid numerical errors on large images
00814 float QuadMinDotCmp::cmp(EMData * image, EMData *with) const
00815 {
00816         ENTERFUNC;
00817         validate_input_args(image, with);
00818 
00819         if (image->get_zsize()!=1) throw InvalidValueException(0, "QuadMinDotCmp supports 2D only");
00820 
00821         int nx=image->get_xsize();
00822         int ny=image->get_ysize();
00823 
00824         int normalize = params.set_default("normalize", 0);
00825         float negative = (float)params.set_default("negative", 1);
00826 
00827         if (negative) negative=-1.0; else negative=1.0;
00828 
00829         double result[4] = { 0,0,0,0 }, sq1[4] = { 0,0,0,0 }, sq2[4] = { 0,0,0,0 } ;
00830 
00831         vector<int> image_saved_offsets = image->get_array_offsets();
00832         vector<int> with_saved_offsets = with->get_array_offsets();
00833         image->set_array_offsets(-nx/2,-ny/2);
00834         with->set_array_offsets(-nx/2,-ny/2);
00835         int i,x,y;
00836         for (y=-ny/2; y<ny/2; y++) {
00837                 for (x=-nx/2; x<nx/2; x++) {
00838                         int quad=(x<0?0:1) + (y<0?0:2);
00839                         result[quad]+=(*image)(x,y)*(*with)(x,y);
00840                         if (normalize) {
00841                                 sq1[quad]+=(*image)(x,y)*(*image)(x,y);
00842                                 sq2[quad]+=(*with)(x,y)*(*with)(x,y);
00843                         }
00844                 }
00845         }
00846         image->set_array_offsets(image_saved_offsets);
00847         with->set_array_offsets(with_saved_offsets);
00848 
00849         if (normalize) {
00850                 for (i=0; i<4; i++) result[i]/=sqrt(sq1[i]*sq2[i]);
00851         } else {
00852                 for (i=0; i<4; i++) result[i]/=nx*ny/4;
00853         }
00854 
00855         float worst=static_cast<float>(result[0]);
00856         for (i=1; i<4; i++) if (static_cast<float>(result[i])<worst) worst=static_cast<float>(result[i]);
00857 
00858         EXITFUNC;
00859         return (float) (negative*worst);
00860 }
00861 
00862 float OptVarianceCmp::cmp(EMData * image, EMData *with) const
00863 {
00864         ENTERFUNC;
00865         validate_input_args(image, with);
00866 
00867         int keepzero = params.set_default("keepzero", 1);
00868         int invert = params.set_default("invert",0);
00869         int matchfilt = params.set_default("matchfilt",1);
00870         int matchamp = params.set_default("matchamp",0);
00871         int radweight = params.set_default("radweight",0);
00872         int dbug = params.set_default("debug",0);
00873 
00874         size_t size = (size_t)image->get_xsize() * image->get_ysize() * image->get_zsize();
00875 
00876 
00877         EMData *with2=NULL;
00878         if (matchfilt) {
00879                 EMData *a = image->do_fft();
00880                 EMData *b = with->do_fft();
00881 
00882                 vector <float> rfa=a->calc_radial_dist(a->get_ysize()/2,0.0f,1.0f,1);
00883                 vector <float> rfb=b->calc_radial_dist(b->get_ysize()/2,0.0f,1.0f,1);
00884 
00885                 float avg=0;
00886                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00887                         rfa[i]=(rfb[i]==0?0.0f:(rfa[i]/rfb[i]));
00888                         avg+=rfa[i];
00889                 }
00890 
00891                 avg/=a->get_ysize()/2.0f;
00892                 for (size_t i=0; i<a->get_ysize()/2.0f; i++) {
00893                         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
00894                 }
00895                 rfa[0]=0.0;
00896 
00897                 if (dbug) b->write_image("a.hdf",-1);
00898 
00899                 b->apply_radial_func(0.0f,1.0f/a->get_ysize(),rfa);
00900                 with2=b->do_ift();
00901 
00902                 if (dbug) b->write_image("a.hdf",-1);
00903                 if (dbug) a->write_image("a.hdf",-1);
00904 
00905 /*              if (dbug) {
00906                         FILE *out=fopen("a.txt","w");
00907                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfa[i]);
00908                         fclose(out);
00909 
00910                         out=fopen("b.txt","w");
00911                         for (int i=0; i<a->get_ysize()/2.0; i++) fprintf(out,"%d\t%f\n",i,rfb[i]);
00912                         fclose(out);
00913                 }*/
00914 
00915 
00916                 delete a;
00917                 delete b;
00918 
00919                 if (dbug) {
00920                         with2->write_image("a.hdf",-1);
00921                         image->write_image("a.hdf",-1);
00922                 }
00923 
00924 //              with2->process_inplace("matchfilt",Dict("to",this));
00925 //              x_data = with2->get_data();
00926         }
00927 
00928         // This applies the individual Fourier amplitudes from 'image' and
00929         // applies them to 'with'
00930         if (matchamp) {
00931                 EMData *a = image->do_fft();
00932                 EMData *b = with->do_fft();
00933                 size_t size2 = (size_t)a->get_xsize() * a->get_ysize() * a->get_zsize();
00934 
00935                 a->ri2ap();
00936                 b->ri2ap();
00937 
00938                 const float *const ad=a->get_const_data();
00939                 float * bd=b->get_data();
00940 
00941                 for (size_t i=0; i<size2; i+=2) bd[i]=ad[i];
00942                 b->update();
00943 
00944                 b->ap2ri();
00945                 with2=b->do_ift();
00946 //with2->write_image("a.hdf",-1);
00947                 delete a;
00948                 delete b;
00949         }
00950 
00951         const float * x_data;
00952         if (with2) x_data=with2->get_const_data();
00953         else x_data = with->get_const_data();
00954         const float *const y_data = image->get_const_data();
00955 
00956         size_t nx = image->get_xsize();
00957         float m = 0;
00958         float b = 0;
00959 
00960         // This will write the x vs y file used to calculate the density
00961         // optimization. This behavior may change in the future
00962         if (dbug) {
00963                 FILE *out=fopen("dbug.optvar.txt","w");
00964                 if (out) {
00965                         for (size_t i=0; i<size; i++) {
00966                                 if ( !keepzero || (x_data[i] && y_data[i])) fprintf(out,"%g\t%g\n",x_data[i],y_data[i]);
00967                         }
00968                         fclose(out);
00969                 }
00970         }
00971 
00972 
00973         Util::calc_least_square_fit(size, x_data, y_data, &m, &b, keepzero);
00974         if (m == 0) {
00975                 m = FLT_MIN;
00976         }
00977         b = -b / m;
00978         m = 1.0f / m;
00979 
00980         // While negative slopes are really not a valid comparison in most cases, we
00981         // still want to detect these instances, so this if is removed
00982 /*      if (m < 0) {
00983                 b = 0;
00984                 m = 1000.0;
00985         }*/
00986 
00987         double  result = 0;
00988         int count = 0;
00989 
00990         if (radweight) {
00991                 if (image->get_zsize()!=1) throw ImageDimensionException("radweight option is 2D only");
00992                 if (keepzero) {
00993                         for (size_t i = 0,y=0; i < size; y++) {
00994                                 for (size_t x=0; x<nx; i++,x++) {
00995                                         if (y_data[i] && x_data[i]) {
00996 #ifdef  _WIN32
00997                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
00998                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
00999 #else
01000                                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
01001                                                 else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
01002 #endif
01003                                                 count++;
01004                                         }
01005                                 }
01006                         }
01007                         result/=count;
01008                 }
01009                 else {
01010                         for (size_t i = 0,y=0; i < size; y++) {
01011                                 for (size_t x=0; x<nx; i++,x++) {
01012 #ifdef  _WIN32
01013                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(_hypot((float)x,(float)y)+nx/4.0);
01014                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(_hypot((float)x,(float)y)+nx/4.0);
01015 #else
01016                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m)*(hypot((float)x,(float)y)+nx/4.0);
01017                                         else result += Util::square((x_data[i] * m) + b - y_data[i])*(hypot((float)x,(float)y)+nx/4.0);
01018 #endif
01019                                 }
01020                         }
01021                         result = result / size;
01022                 }
01023         }
01024         else {
01025                 if (keepzero) {
01026                         for (size_t i = 0; i < size; i++) {
01027                                 if (y_data[i] && x_data[i]) {
01028                                         if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
01029                                         else result += Util::square((x_data[i] * m) + b - y_data[i]);
01030                                         count++;
01031                                 }
01032                         }
01033                         result/=count;
01034                 }
01035                 else {
01036                         for (size_t i = 0; i < size; i++) {
01037                                 if (invert) result += Util::square(x_data[i] - (y_data[i]-b)/m);
01038                                 else result += Util::square((x_data[i] * m) + b - y_data[i]);
01039                         }
01040                         result = result / size;
01041                 }
01042         }
01043         scale = m;
01044         shift = b;
01045 
01046         image->set_attr("ovcmp_m",m);
01047         image->set_attr("ovcmp_b",b);
01048         if (with2) delete with2;
01049         EXITFUNC;
01050 
01051 #if 0
01052         return (1 - result);
01053 #endif
01054 
01055         return static_cast<float>(result);
01056 }
01057 
01058 float PhaseCmp::cmp(EMData * image, EMData *with) const
01059 {
01060         ENTERFUNC;
01061 
01062         int snrweight = params.set_default("snrweight", 0);
01063         int snrfn = params.set_default("snrfn",0);
01064         int ampweight = params.set_default("ampweight",0);
01065         int zeromask = params.set_default("zeromask",0);
01066         float minres = params.set_default("minres",500.0f);
01067         float maxres = params.set_default("maxres",10.0f);
01068 
01069         if (snrweight && snrfn) throw InvalidCallException("SNR weight and SNRfn cannot both be set in the phase comparator");
01070 
01071         EMData *image_fft = NULL;
01072         EMData *with_fft = NULL;
01073 
01074         int ny = image->get_ysize();
01075 //      int np = (int) ceil(Ctf::CTFOS * sqrt(2.0f) * ny / 2) + 2;
01076         int np = 0;
01077         vector<float> snr;
01078 
01079         // weighting based on SNR estimate from CTF
01080         if (snrweight) {
01081                 Ctf *ctf = NULL;
01082                 if (!image->has_attr("ctf")) {
01083                         if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
01084                         ctf=with->get_attr("ctf");
01085                 }
01086                 else ctf=image->get_attr("ctf");
01087 
01088                 float ds=1.0f/(ctf->apix*ny);
01089                 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR); // note that this returns ny/2 values
01090                 for (int i=0; i<snr.size(); i++) {
01091                         if (snr[i]<=0) snr[i]=0.001;            // make sure that points don't get completely excluded due to SNR estimation issues, or worse, contribute with a negative weight
01092                 }
01093                 if(ctf) {delete ctf; ctf=0;}
01094                 np=snr.size();
01095         }
01096         // weighting based on empirical SNR function (not really good)
01097         else if (snrfn==1) {
01098                 np=ny/2;
01099                 snr.resize(np);
01100 
01101                 for (int i = 0; i < np; i++) {
01102                         float x2 = 10.0f*i/np;
01103                         snr[i] = x2 * exp(-x2);
01104                 }
01105         }
01106         // no weighting
01107         else {
01108                 np=ny/2;
01109                 snr.resize(np);
01110 
01111                 for (int i = 0; i < np; i++) snr[i]=1.0;
01112         }
01113 
01114         // Min/max modifications to weighting
01115         float pmin,pmax;
01116         if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;                //cutoff in pixels, assume square
01117         else pmin=0;
01118         if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
01119         else pmax=0;
01120 
01121 //      printf("%f\t%f\n",pmin,pmax);
01122 
01123         // We use 'soft' edges for the Fourier cutoffs to minimize issues with pixel interpolation
01124         for (int i = 0; i < np; i++) {
01125                 if (pmin>0) snr[i]*=(tanh(5.0f*(i-pmin)/pmin)+1.0f)/2.0f;
01126                 if (pmax>0) snr[i]*=(1.0f-tanh(i-pmax))/2.0f;
01127 //              printf("%d\t%f\n",i,snr[i]);
01128         }
01129 
01130         if (zeromask) {
01131                 image_fft=image->copy();
01132                 with_fft=with->copy();
01133                 
01134                 if (image_fft->is_complex()) image_fft->do_ift_inplace();
01135                 if (with_fft->is_complex()) with_fft->do_ift_inplace();
01136                 
01137                 int sz=image_fft->get_xsize()*image_fft->get_ysize()*image_fft->get_zsize();
01138                 float *d1=image_fft->get_data();
01139                 float *d2=with_fft->get_data();
01140                 
01141                 for (int i=0; i<sz; i++) {
01142                         if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
01143                 }
01144                 
01145                 image_fft->update();
01146                 with_fft->update();
01147                 image_fft->do_fft_inplace();
01148                 with_fft->do_fft_inplace();
01149                 image_fft->set_attr("free_me",1); 
01150                 with_fft->set_attr("free_me",1); 
01151         }
01152         else {
01153                 if (image->is_complex()) image_fft=image;
01154                 else {
01155                         image_fft=image->do_fft();
01156                         image_fft->set_attr("free_me",1);
01157                 }
01158                 
01159                 if (with->is_complex()) with_fft=with;
01160                 else {
01161                         with_fft=with->do_fft();
01162                         with_fft->set_attr("free_me",1);
01163                 }
01164         }
01165 //      image_fft->ri2ap();
01166 //      with_fft->ri2ap();
01167 
01168         const float *const image_fft_data = image_fft->get_const_data();
01169         const float *const with_fft_data = with_fft->get_const_data();
01170         double sum = 0;
01171         double norm = FLT_MIN;
01172         size_t i = 0;
01173 //      int nx=image_fft->get_xsize();
01174             ny=image_fft->get_ysize();
01175         int nz=image_fft->get_zsize();
01176         int nx2=image_fft->get_xsize()/2;
01177         int ny2=image_fft->get_ysize()/2;
01178 //      int nz2=image_fft->get_zsize()/2;
01179 
01180         // This can never happen any more...
01181         if (np==0) {
01182                 for (int z = 0; z < nz; z++){
01183                         for (int y = 0; y < ny; y++) {
01184                                 for (int x = 0; x < nx2; x ++) {
01185                                         float a;
01186                                         if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01187                                         else a=1.0f;
01188                                         sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01189                                         norm += a;
01190                                         i += 2;
01191                                 }
01192                         }
01193                 }
01194                 
01195         }
01196         else if (nz==1) {
01197                 for (int y = 0; y < ny; y++) {
01198                         for (int x = 0; x < nx2; x ++) {
01199                                 int r=Util::hypot_fast_int(x,y>ny/2?ny-y:y);
01200                                 if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
01201                                 
01202                                 float a;
01203                                 if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01204                                 else a=1.0f;
01205                                 a*=snr[r];
01206                                 sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01207                                 norm += a;
01208                                 i += 2;
01209                         }
01210                 }
01211         }
01212         else {
01213                 for (int z = 0; z < nz; z++){
01214                         for (int y = 0; y < ny; y++) {
01215                                 for (int x = 0; x < nx2; x ++) {
01216                                         int r=(int)Util::hypot3(x,y>ny/2?ny-y:y,z>nz/2?nz-z:z);
01217                                         if (r>=ny2) { i+=2; continue; }         // we only have snr values to the box radius
01218                                         
01219                                         float a;
01220                                         if (ampweight) a= (float)hypot(with_fft_data[i],with_fft_data[i+1]);
01221                                         else a=1.0f;
01222                                         a*=snr[r];
01223                                         sum += Util::angle_err_ri(image_fft_data[i],image_fft_data[i+1],with_fft_data[i],with_fft_data[i+1]) * a;
01224                                         norm += a;
01225                                         i += 2;
01226                                 } 
01227                         }
01228                 }
01229                 
01230         }
01231 
01232         EXITFUNC;
01233 
01234         if( image_fft->has_attr("free_me") )
01235         {
01236                 delete image_fft;
01237                 image_fft = 0;
01238         }
01239         if( with_fft->has_attr("free_me") )
01240         {
01241                 delete with_fft;
01242                 with_fft = 0;
01243         }
01244 #if 0
01245         return (1.0f - sum / norm);
01246 #endif
01247         return (float)(sum / norm);
01248 }
01249 
01250 float FRCCmp::cmp(EMData * image, EMData * with) const
01251 {
01252         ENTERFUNC;
01253         validate_input_args(image, with);
01254 
01255         int snrweight = params.set_default("snrweight", 0);
01256         int ampweight = params.set_default("ampweight", 0);
01257         int sweight = params.set_default("sweight", 1);
01258         int nweight = params.set_default("nweight", 0);
01259         int zeromask = params.set_default("zeromask",0);
01260         float minres = params.set_default("minres",500.0f);
01261         float maxres = params.set_default("maxres",10.0f);
01262 
01263         vector < float >fsc;
01264         bool use_cpu = true;
01265 
01266         if (use_cpu) {
01267                 if (zeromask) {
01268                         image=image->copy();
01269                         with=with->copy();
01270                 
01271                         int sz=image->get_xsize()*image->get_ysize()*image->get_zsize();
01272                         float *d1=image->get_data();
01273                         float *d2=with->get_data();
01274                 
01275                         for (int i=0; i<sz; i++) {
01276                                 if (d1[i]==0.0 || d2[i]==0.0) { d1[i]=0.0; d2[i]=0.0; }
01277                         }
01278                 
01279                         image->update();
01280                         with->update();
01281                         image->do_fft_inplace();
01282                         with->do_fft_inplace();
01283                         image->set_attr("free_me",1); 
01284                         with->set_attr("free_me",1); 
01285                 }
01286 
01287 
01288                 if (!image->is_complex()) {
01289                         image=image->do_fft(); 
01290                         image->set_attr("free_me",1); 
01291                 }
01292                 if (!with->is_complex()) { 
01293                         with=with->do_fft(); 
01294                         with->set_attr("free_me",1); 
01295                 }
01296 
01297                 fsc = image->calc_fourier_shell_correlation(with,1);
01298         }
01299         
01300         int ny = image->get_ysize();
01301         int ny2=ny/2+1;
01302 
01303         // The fast hypot here was supposed to speed things up. Little effect
01304 //      if (image->get_zsize()>1) fsc = image->calc_fourier_shell_correlation(with,1);
01305 //      else {
01306 //              double *sxy = (double *)malloc(ny2*sizeof(double)*4);
01307 //              double *sxx = sxy+ny2;
01308 //              double *syy = sxy+2*ny2;
01309 //              double *norm= sxy+3*ny2;
01310 //
01311 //              float *df1=image->get_data();
01312 //              float *df2=with->get_data();
01313 //              int nx2=image->get_xsize();
01314 //
01315 //              for (int y=-ny/2; y<ny/2; y++) {
01316 //                      for (int x=0; x<nx2/2; x++) {
01317 //                              if (x==0 && y<0) continue;      // skip Friedel pair
01318 //                              short r=Util::hypot_fast_int(x,y);
01319 //                              if (r>ny2-1) continue;
01320 //                              int l=x*2+(y<0?ny+y:y)*nx2;
01321 //                              sxy[r]+=df1[l]*df2[l]+df1[l+1]*df2[l+1];
01322 //                              sxx[r]+=df1[l]*df1[l];
01323 //                              syy[r]+=df2[l]*df2[l];
01324 //                              norm[r]+=1.0;
01325 //                      }
01326 //              }
01327 //              fsc.resize(ny2*3);
01328 //              for (int r=0; r<ny2; r++) {
01329 //                      fsc[r]=r*0.5/ny2;
01330 //                      fsc[ny2+r]=sxy[r]/(sqrt(sxx[r])*sqrt(syy[r]));
01331 //                      fsc[ny2*2+r]=norm[r];
01332 //              }
01333 //              free(sxy);
01334 //      }
01335 
01336         vector<float> snr;
01337         if (snrweight) {
01338                 Ctf *ctf = NULL;
01339                 if (!image->has_attr("ctf")) {
01340                         if (!with->has_attr("ctf")) throw InvalidCallException("SNR weight with no CTF parameters");
01341                         ctf=with->get_attr("ctf");
01342                 }
01343                 else ctf=image->get_attr("ctf");
01344 
01345                 float ds=1.0f/(ctf->apix*ny);
01346                 snr=ctf->compute_1d(ny,ds,Ctf::CTF_SNR);
01347                 for (int i=0; i<snr.size(); i++) {
01348                         if (snr[i]<=0) snr[i]=0.001;            // make sure that points don't get completely excluded due to SNR estimation issues, or worse, contribute with a negative weight
01349                 }
01350                 if(ctf) {delete ctf; ctf=0;}
01351         }
01352 
01353         vector<float> amp;
01354         if (ampweight) amp=image->calc_radial_dist(ny/2,0,1,0);
01355 
01356         // Min/max modifications to weighting
01357         float pmin,pmax;
01358         if (minres>0) pmin=((float)image->get_attr("apix_x")*image->get_ysize())/minres;                //cutoff in pixels, assume square
01359         else pmin=0;
01360         if (maxres>0) pmax=((float)image->get_attr("apix_x")*image->get_ysize())/maxres;
01361         else pmax=0;
01362 
01363         double sum=0.0, norm=0.0;
01364 
01365         for (int i=0; i<ny/2; i++) {
01366                 double weight=1.0;
01367                 if (sweight) weight*=fsc[(ny2)*2+i];
01368                 if (ampweight) weight*=amp[i];
01369                 if (snrweight) weight*=snr[i];
01370 //              if (snrweight)  {
01371 //                      if (snr[i]>0) weight*=sqrt(snr[i]);
01372 //                      else weight=0;
01373 //              }
01374 //if(snr[i]<0) printf("snr[%d] = %1.5g\n",i,snr[i]);
01375                 if (pmin>0) weight*=(tanh(5.0*(i-pmin)/pmin)+1.0)/2.0;
01376                 if (pmax>0) weight*=(1.0-tanh(i-pmax))/2.0;
01377                 
01378                 sum+=weight*fsc[ny2+i];
01379                 norm+=weight;
01380 //              printf("%d\t%f\t%f\n",i,weight,fsc[ny/2+1+i]);
01381         }
01382 
01383         // This performs a weighting that tries to normalize FRC by correcting from the number of particles represented by the average
01384         sum/=norm;
01385         if (nweight && with->get_attr_default("ptcl_repr",0) && sum>=0 && sum<1.0) {
01386                 sum=sum/(1.0-sum);                                                      // convert to SNR
01387                 sum/=(float)with->get_attr_default("ptcl_repr",0);      // divide by ptcl represented
01388                 sum=sum/(1.0+sum);                                                      // convert back to correlation
01389         }
01390 
01391         if (image->has_attr("free_me")) delete image;
01392         if (with->has_attr("free_me")) delete with;
01393 
01394         EXITFUNC;
01395 
01396         if (!Util::goodf(&sum)) sum=-2.0;       // normally should be >-1.0
01397 
01398         //.Note the negative! This is because EMAN2 follows the convention that
01399         // smaller return values from comparitors indicate higher similarity -
01400         // this enables comparitors to be used in a generic fashion.
01401         return (float)-sum;
01402 }
01403 
01404 float OptSubCmp::cmp(EMData * image, EMData * with) const
01405 {
01406         ENTERFUNC;
01407         validate_input_args(image, with);
01408 
01409 //      int snrweight = params.set_default("snrweight", 0);
01410 //      int ampweight = params.set_default("ampweight", 0);
01411 //      int sweight = params.set_default("sweight", 1);
01412 //      int nweight = params.set_default("nweight", 0);
01413         int ctfweight = params.set_default("ctfweight",0);
01414         int zeromask = params.set_default("zeromask",0);
01415 //      int swap = params.set_default("swap",0);
01416         float minres = params.set_default("minres",200.0f);
01417         float maxres = params.set_default("maxres",10.0f);
01418         EMData *mask = params.set_default("mask",(EMData *)NULL);
01419         
01420 //      float ds=1.0f/((float)image->get_attr("apix_x")*(int)image->get_ysize());
01421         float apix=(float)image->get_attr("apix_x");
01422 
01423         // Sometimes we will get the "raw" image with CTF as image and sometimes as with, we always want to subtract the less noisy
01424         // reference, so if one has CTF parameters (even if ctfweight isn't used) we always pass it in as the primary
01425         EMData *diff;
01426         if (image->has_attr("ctf")) diff=image->process("math.sub.optimal",Dict("ref",with,"return_presigma",1,"low_cutoff_frequency",apix/minres ,"high_cutoff_frequency",apix/maxres,"ctfweight",ctfweight));
01427         else diff=with->process("math.sub.optimal",Dict("ref",image,"return_presigma",1,"low_cutoff_frequency",apix/minres ,"high_cutoff_frequency",apix/maxres,"ctfweight",ctfweight));
01428 
01429         if (mask!=NULL) diff->mult(*mask);
01430         if (zeromask) {
01431                 EMData *tmp=with->process("threshold.notzero");
01432                 diff->mult(*tmp);
01433                 delete tmp;
01434         }
01435         
01436 //      diff->process_inplace("filter.highpass.tophat",Dict("cutoff_freq",(float)1.0/minres));
01437 //      diff->process_inplace("filter.lowpass.tophat",Dict("cutoff_freq",(float)1.0/maxres));
01438         
01439         float ret=(float)diff->get_attr("sigma")/(float)diff->get_attr("sigma_presub");
01440         delete diff;
01441         return ret;
01442         
01443         
01444 //      EMData *diff=image->process("math.sub.optimal",Dict("ref",with,"return_fft",1));
01445 //      
01446 //      // Very expensive to use a mask, since it requires another ift/fft pair
01447 //      if (mask!=NULL) {
01448 //              EMData *tmp = diff->do_ift();
01449 //              tmp->mult(*mask);
01450 //              delete diff;
01451 //              diff=tmp->do_fft();
01452 //              delete tmp;
01453 //      }
01454 //      
01455 //      // This gives us basically the 1-D power spectrum of what's left after subtraction (and masking)
01456 //      vector<float> dist=diff->calc_radial_dist(diff->get_ysize()/2,0.0f,1.0f,1);
01457 //      int s0=int(floor(1.0/(minres*ds)));
01458 //      int s1=int(ceil(1.0/(maxres*ds)));
01459 //      if (s0<2) s0=2;
01460 //      
01461 //      if (s1<=s0) throw InvalidCallException("OptSubCmp error. minres must be greater than maxres.");
01462 // //   printf("%d %d\n",s0,s1);
01463 //      
01464 //      double sum=0.0f,sum2=0.0f;
01465 //      for (int i=s0; i<s1; i++) sum+=dist[i]*i;
01466 //      sum/=image->get_size()*s1;
01467 //      
01468 //      delete diff;
01469 //      return sum;
01470 }
01471 
01472 
01473 void EMAN::dump_cmps()
01474 {
01475         dump_factory < Cmp > ();
01476 }
01477 
01478 map<string, vector<string> > EMAN::dump_cmps_list()
01479 {
01480         return dump_factory_list < Cmp > ();
01481 }
01482 
01483 /* vim: set ts=4 noet: */