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

Generated on Tue Jun 11 13:46:13 2013 for EMAN2 by  doxygen 1.3.9.1