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

Generated on Tue Jul 12 13:49:24 2011 for EMAN2 by  doxygen 1.3.9.1